馮云青,侯磊,2

(1.上海大學數學系,上海200444;2.上海高校計算科學E-研究院,上海200030)

構造高階插值函數求解表面沖擊碰撞問題

馮云青1,侯磊1,2

(1.上海大學數學系,上海200444;2.上海高校計算科學E-研究院,上海200030)

本文主要采用Lagrange雙三次形函數來構造插值函數,從而利用有限元方法求解表面沖擊碰撞問題中的耦合方程組.為了避免龍格現象,采用Lobatto點構造插值節(jié)點.本文不僅采用高階的形函數,也使用兩種不同的數值積分方法,結合這兩點,從而提高數值解的精度.根據以上的理論分析結果,本中使用Matlab編程模擬發(fā)生碰撞時材料的形變和應力變化.

有限元;Lagrange;正交積分

0 引言

汽車電子化的飛速發(fā)展,對汽車產業(yè)的相關研究帶來了前所未有的現代化沖擊.由于汽車結構發(fā)展的多樣性和復雜性,維修設備的現代化發(fā)展,汽車修復技術需要變革.為了能更好地變革修復技術,進而對研究汽車碰撞的模擬實驗也提出了更多樣性的模擬和更好的仿真效果要求.因此,汽車碰撞問題以及修復工作的研究越來越需要高技術的引入.目前對于汽車碰撞試驗的數據主要來源于實驗室實物模擬與模型數值計算結合.采用數學模型與高性能計算來虛擬化試驗,受到了廣泛的關注.本文主要研究汽車碰撞安全相關問題,以高速碰撞時接觸表面的形變與應力分布為主要研究對象,模擬實物高速碰撞后,復雜接觸表面的形態(tài).在汽車制造的材料選擇中,蜂窩多孔材料由于其質量輕、抗壓抗彎性能高而受到廣泛關注,但這種材料微觀結構復雜,材料性質各異,導致對其的數值模擬困難重重.本文試圖采用非線性流固耦合方程組,對蜂窩多孔材料所構建的復雜接觸表面進行數值模擬.該流固耦合方程組[1-2]為：

其中ρ是流體密度,ξ是切變因子,η是粘度,λ是松弛因子.該耦合方程是建立在高速碰撞的條件下,固體材料在一段時間內發(fā)生了巨大形變,因此認為它同時具有了流體性質.耦合方程組的第一式是Cauchy守恒方程,用于計算受到巨大沖擊時應力場τ分布產生的彈塑性材料大變形,該變形由速度場u描述.第二式是P-T/T方程,用它來描述粘彈塑性材料應力變化時復雜接觸表面的外延和簡單剪切速率的阻力.P-T/T是建立在Maxwell方程[2-4]的基礎上的具有瞬時拉伸特征的應力計算方程.

本文采用有限元和有限差分相結合的方式求解該流固耦合方程組,對于空間變量采用有限元方法進行離散,再用三種差分方法對時間變量進行離散,最終得到大型線性方程組[1,5-6].在對空間變量采用有限元方法離散時,引入高階的形函數——Lagrange雙三次形函數來構造插值函數,避免單元出現過于剛化的現象,從而減小了空間有限元離散所帶來的誤差.同時,當插值函數的階數提高后,為了避免可能出現的龍格現象,本文用Lobatto多項式構造不等距的插值節(jié)點,且在文章[1]的基礎上加入更多的插值節(jié)點使網格剖分更細,從而進一步提高有限元解的精度.在具體計算單元剛度矩陣時,由于被積函數的復雜性,需要引入積分點進行數值積分,選用較少的積分點數是另一種克服單元剛化的手段.本文是采用Gauss和Gauss-Lobatto兩種積分點進行數值積分.對于積分點數的選擇是有要求的,若采用過多的Gauss或Gauss-Lobatto點求積,單元剛度矩陣會因為過度剛性而變得更加糟糕,因為附加的Gauss或Gauss-Lobatto點捕捉的是剛度矩陣更高次的項.這些高次項能抵抗低次項所不能抵抗的變形,因而使單元變得剛性化.同時,更多點的積分也增加了計算量.除此以外,剛度矩陣的高精度經常導致有限元解的低精度[7].然而,若用太少的Gauss或Gauss-Lobatto點可能導致更壞結果,產生不穩(wěn)定性等等.因而在有限元問題中,一個單元的最佳積分法則往往也是值得思考和探索的.文本由于對每個單元用矩形剖分后有4×4個插值節(jié)點,為此分別采用3×3和4×4個Gauss點,4×4個Gauss-Lobatto點在每個單元上進行數值積分[8-9].

1 形函數的構造

為了提高有限元解的收斂速度或減小有限元解的誤差,通常采用兩種方式,一種是適當提高單元形函數中多項式階數,另一種是采用更精細的單元網格來提高有限元結果的精度.構造恰當的有限元插值函數的第一步是確定插值節(jié)點.本文采用Lobatto多項式的零點作為插值節(jié)點,構造Lagrange16點雙三次多項式形函數.由于插值形函數階數的提高,本文為避免插值函數出現震蕩即龍格現象,采用不等距的插值節(jié)點.

首先考慮在區(qū)間[-h,h]上的Lobatto多項式序列[10]為：

由此得到該多項式的零點為：

本文下面考慮的問題都是在標準元素[-1,1]×[-1,1]的矩形區(qū)域上的.將該矩形區(qū)域四等分為2×2個單元,記為單元①、單元②、單元③、單元④.再利用以上求得的Lobatto多項式的零點作為插值節(jié)點做矩形剖分,可得到每個單元上的不等距插值節(jié)點為4×4個.再對節(jié)點進行編號,如圖1.1,首先對每個單元上的16個節(jié)點進行局部編號,它們的編號稱為節(jié)點的局部編碼.如圖1.2,其次對區(qū)域內所有的插值節(jié)點按照節(jié)約內存的方式進行總體編號,這樣的編號稱為節(jié)點的總體編碼[1,10].

圖1 .1 單元網格插值節(jié)點編號Fig.1.1 The number of interpolation nodes of element mesh

圖1 .2 總體網格插值節(jié)點編號Fig.1.2 The number of interpolation nodes of global mesh

編碼完成后,可列出構造方程離散格式所需要的Lagrange16點雙三次形函數.可知在區(qū)間[-1,1]上的Lobatto多項式有4個零點,則一維的Lagrange形函數構造如下：

由此可得二維上的Lagrange16點雙三次形函數,表示為：

其中Li(η)也是一維Lagrange形函數.

本文在每個單元上利用16個Lobatto插值節(jié)點形成單元剛度矩陣,然后再將單元剛度矩陣合成總剛矩陣.由上面的剖分可以得到總剛矩陣為49×49矩陣.剛度矩陣的計算跟選取的有限元形函數有關.基于位移的有限元法可能存在單元出現過于剛化的現象,為了克服這一現象,可采用高階的形函數.本文構造的是Lagrange16點雙三次形函數.單元內任一點的位移是由假定的形函數及單元結點位移插值得到.通常高階單元比低階單元更加光滑,因為高階單元采用更高階的形函數,其變形更光滑,相當于外加給單元的約束變少了.這里具體的剛度矩陣的計算是用Matlab編寫程序實現的,可參考后面的數值結果.

2 耦合方程組的離散

本節(jié)將對流固耦合方程組(0.1)進行有限元——有限差分離散.設方程組(0.1)中, u=iu1(x,y;t)+ju2(x,y;t)為速度場(矢量),其中的u1(x,y;t),u2(x,y;t)分別為速度u的x分量和y分量；為應力場(張量),其中的τ11(x,y;t),τ12(x,y;t), τ21(x,y;t),τ22(x,y;t)分別為應力τ對應的四個分量并且τ12(x,y;t)=τ21(x,y;t).故該耦合方程組實際上包含五個分量方程.將耦合方程組(0.1)按上述五個分量展開得到五個分量方程為：

對于要計算出的u,τ且u∈H,τ∈H,其求解是有困難的.通常做法是將原問題(0.1)轉化為積分形式,降低解對光滑性的要求,此時積分形式對應的解稱為廣義解并稱廣義解所在的空間為容許空間.再選取容許空間的有限維子空間,在該子空間內對廣義解進一步逼近.故本文選取容許空間H的有限維子空間Vh∈H為：

稱這樣的Vh為檢驗空間,其中?m(x,y)(m=1,2,···,16)為檢驗函數空間的基,各?m(x,y)在標準網格上的表達形式同?m(ξ,η).此時H被檢驗空間Vh所替代,問題轉化為求近似解使得

其中,L為對應耦合方程組(0.1)的微分算子.進一步,使用有限元方法,將這一過程模塊化,且確保解的存在唯一性.

下文對空間變量(x,y)以及時間變量t分開處理.對于以上五個方程的空間變量,采用有限元方法進行離散.以下均是在標準網格e=[-1,1]×[-1,1]上,采用Lagrange16點雙三次形函數,實現的空間有限元半離散過程.這里相當于推導了在一個單元上的空間變量的離散過程,寫出的剛度矩陣均為標準網格下的單元剛度矩陣,本文中共有四個這樣的單元,如圖1.2,一共有49個總體編碼.作者在用Matlab編寫時,將標準網格上的離散結果通過坐標變換推導到實際的四個單元上,需要將節(jié)點的局部編碼和總體節(jié)點編碼進行轉換,再將每個單元的剛度矩陣合成總體剛度矩陣.

令u1(xm,ym;tn)表示當t=tn時,速度場的x分量在該節(jié)點上的數值.同理可推得其他分量,這里不一一列舉.令為以?m為插值基的函數,有限元解的插值基函數表達式為：

其他同理可得.

插值基函數的偏導數可表示為：

其他也同理可得.

其中?mx、?my分別表示對?m求x、y方向的偏導.式子(2.10)—(2.14)分別對應(2.1)—(2.5)五個方程的半離散,是對空間變量的有限元離散.需要注意的是,這里是定義在標準網格上的一個單元的離散過程.

對于時間變量的離散,本文采用三種差分方法[1,5].這里僅以方程(2.1)半離散后的有限元格式(2.10)為例,給出三種時間差分方法的全離散過程.先將半離散格式(2.10)改寫為：

推導出對應的Euler格式為：

對應的半隱式Crank-Nicolson格式為：

對應的Adams二步顯格式為：

由此,方程(2.1)的空間變量和時間變量均被離散,得到了三種不同的全離散格式.其余四個方程,均可按照這三種差分格式全離散,不再贅述.三種全離散格式的程序設計見第4節(jié).

3 Gauss-Lobatto數值積分

由第2節(jié)可以看出,在對剛度矩陣的計算中被積函數是較為復雜的,因此求其積分需要采用數值積分的方法.本文采用Gauss和Gauss-Lobatto型兩種求積方法,Gauss-Lobatto求積公式與Gauss求積公式的區(qū)別在于其要求被積區(qū)間的兩個端點均固定,因此所有積分點中的端點都被預先固定.眾所周知,數值積分計算并不總能得到精確解,通常只能得到近似解,如何才能得到較高精度的近似解是一個值得研究的問題[8,12],這與積分點的選取密不可分.在實際計算時會發(fā)現當積分點數增加時,數值積分結果能收斂到精確解,然而其收斂過程并不是單調的.除此之外,一般情況下,為了克服剛度矩陣的剛化現象會采用較少的單元積分點數計算單元矩陣,避免過多的積分點捕捉剛度矩陣中更多的高次項使其剛性變得更加糟糕.為此,選取恰當的單元積分點數是需要經過數值試驗的.本文選取的形函數是Lagrange16點雙三次形函數,采用3×3和4×4個Gauss點,4×4個Gauss-Lobatto點在每個單元上進行數值積分試驗,確保解的準確性.

4 數值結果

綜合以上對耦合方程組(0.1)有限元——有限差分的全離散結果以及數值積分點的選取,運用Matlab編程,模擬在[0,2]×[0,2]區(qū)域中,受不同外力后,網格的速度場和應力場的變化.本文參考有限元軟件Dyna的后處理模式,用網格的形變表示速度場的變化,給各個小單元標識不同顏色表示應力場的變化.這里為了提高解的精度,更好地體現形變和應力變化,本文設步長h=0.5.圖4.1是利用9點Gauss積分得到的三種有限元——有限差分形式的形變圖.其中?t為時間步長,u10、u20為初始速度.9點Gauss積分對應的應力變化圖4.2如下.其中色塊的顏色越深,表示受到的應力越大;每幅圖中的每個小圖分別對應各個時間的應力變化.三幅大圖分別表示用Lagrange-Euler、Lagrange-Crank-Nicolson和Lagrange-Adams這三種網格剖分方式計算得到的應力圖.為了顯示每種網格剖分對應的應力變化,下圖分別選取兩個時間節(jié)點t=1,t=2的應力圖.

圖4.3是利用16點Gauss積分得到的三種有限元——有限差分形式的形變圖.如圖4.4表示16點Gauss積分對應三種剖分格式的應力變化.為展示每種網格剖分的應力變化情況,本文對各個剖分格式下的計算結果分別選取三個時間點的應力圖,從左向右依次對應t=1,t=2,t=3.

圖4 .1 9點Gauss積分,?t=1,u10=0.1,u20=0Fig.4.1 Gauss integration with 9 points,?t=1,u10=0.1,u20=0

圖4 .2 9點Gauss積分對應三種剖分格式的應力變化圖Fig.4.2 The stress distribution calculated by Gauss integration with 9 points of three kinds of mesh

圖4.3 16點Gauss積分,?t=1,u10=0.1,u20=0Fig.4.3 Gauss integration with 16 points,?t=1,u10=0.1,u20=0

圖4 .5是利用16點Gauss-Lobatto積分得到的三種有限元——有限差分形式的形變圖.圖4.6中三幅圖分別表示在Lagrange-Euler、Lagrange-Crank-Nicolson以及Lagrange-Adams剖分格式下,利用16點Gauss-Lobatto積分求得的應力變化圖.其中Lagrange-Euler和Lagrange-Adams剖分格式中截取了兩個時間點的應力變化圖,從左往右分別對應t=1,t=2.Lagrange-Crank-Nicolson格式截取了三個時間點的應力變化圖，即t=1,t=2,t=3.

圖4 .4 16點Gauss積分對應三種剖分格式的應力變化圖Fig.4.4 The stress distribution calculated by Gauss integration with 16 points of three kinds of mesh

圖4 .5 16點Gauss-Lobatto積分,?t=1,u10=0.1,u20=0Fig.4.5 Gauss-Lobatto integration with 16 points,?t=1,u10=0.1,u20=0

圖4 .6 16點Gauss-Lobatto積分對應三種剖分格式的應力變化圖Fig.4.6 The stress distribution calculated by Gauss-Lobatto integration with 16 points of three kinds of mesh

5 結論

本文采用有限元方法對空間變量進行離散,主要運用Lagrange16點雙三次形函數構造插值函數.再結合三種不同的有限差分對時間變量進行離散,最終求解耦合方程組(0.1)的數值解.文中對剛度矩陣中數值積分的求解,采用不同類型的求積公式進行數值模擬.最終發(fā)現,9點的Gauss積分和16點的Gauss-Lobatto積分效果相差不大;16點的Gauss積分更好地體現了形變程度,精確度更高.總體上這三種數值積分方式都可以說明采用Lagrange16點雙三次形函數提高了耦合方程組數值解的精度.

6 致謝

作者在寫這篇文章時受到各位專家學者的鼓勵和指導,特在此對他們的熱情幫助和指導表示感謝.

[1]李涵靈.復雜表面接觸問題的高性能有限元求解及其數據后處理[D].上海：上海大學,2014.

[2]PHAN-TH IEN N.A nonlinear network viscoelastic model[J].Journal of Rheology,1978,22(3)：259-283.

[3]TANNER R I.Engineering Rheology[M].London：Oxford University Press,2000.

[4]THIEN N P,TANNER R I.A new constitutive equation derived from network theory[J].Journal of Non Newtonian Fluid Mechanics,1977(2)：353-365.

[5]李立康.數值計算方法[M].上海：復旦大學出版社,1999.

[6]LIN Q,ZHOU J.Superconvergence in high order galerkin finite element methods[J].Computer Methods in Applied Mechanics and Engineering,1999,196(3)：3779-3784.

[7]王N成,邵敏.有限單元法基本原理和數值方法[M].北京：清華大學出版社,1997.

[8]HELENBROOK B T.On the existence of explicit hp-finite element methods using Gauss-Lobatto integration on the triangle[J].SIAM Journal on Numerical Analysis,2009,47(2)：1304-1318.

[9]XU Y.On Gauss-Lobatto integration on the triangle[J].SIAM Journal on Numerical Analysis,2011,49(2)：541-548.

[10]李開泰,黃艾香,黃慶懷.有限元方法及其應用[M].北京：科學出版社,2006.

[11]HOU L,NASSEHI V.Evaluation of stress-effective flow in rubber mixing[J].Non linear Analysis Theory Methods and Applications,2001,47(3)：1809-1820.

[12]BOYD J P.A numerical comparison of seven grids for polynomial interpolation on the interval[J].Computers and Mathematics with Applications,1999,38(3)：35-50.

(責任編輯：林磊)

High order interpolation function for surface contact problem

FENG Yun-qing1,HOU Lei1,2
(1.Department of Mathematics,Shanghai University,Shanghai 200444,China; 2.Computational Sciences E-institute of Shanghai Universities,Shanghai 200030,China)

This paper mainly adopts Lagrange bicubic shape function to construct interpolation function and uses finite element method to solve the coup ling equations of surface contact.The Lobatto points are used to construct the interpolation nodes to avoid the Runge phenomenon.Higher shape functions and two different numerical integration methods are adopted to imp rove the accuracy of the numerical solution.According to the above analysis,this article uses Matlab program to simulate the deformation and stress changes in surface contact problem.

finite element;Lagrange;orthogonal integration

O 241

A

10.3969/j.issn.1000-5641.2016.03.002

1000-5641(2016)03-0009-12

2015-05

國家自然科學基金項目(11271247)

馮云青,女,碩士研究生,研究方向為有限元軟件應用.E-mail：yqyzsh@126.com.

侯磊,男,教授,博士,研究方向為非線性有限元、自適應方法及應用. E-mail：houlei@shu.edu.cn.