基于S-R和分解定理的無網格法在功能梯度板中的應用

宋彥琦,石博康,李向上

(中國礦業大學(北京)力學與建筑工程學院,北京 100083)

隨著科技的不斷發展,工程環境愈加復雜嚴苛,人們對材料的性能要求也隨之提升.由于傳統的材料已經不能滿足科研和工程要求,研究者將目光轉向了復合材料.材料屬性的不連續性導致傳統復合材料有明顯的缺陷,易于失效.而功能梯度復合材料(functionally graded material,FGM)在微觀結構上連續變化,能夠消除不同組分材料的界面效應,具有高強度、高韌性,可減少應力集中等力學性能[1],被廣泛應用于機械、化工、電子、核工程和航空航天等領域.在某些嚴苛的工程環境中,功能梯度板(functionally graded plate,FGP)常會產生非線性變形,因此研究FGP板的非線性變形特征是有意義的.

國內外學者對功能梯度復合材料板的非線性彎曲進行了大量研究.賈金正等[2]采用物理中面的概念,研究了橫縱荷載作用下梯度指數和長高比對功能梯度梁變形情況的影響.王雪等[3]采用逐步加載技術,研究了熱-機械荷載聯合作用下功能梯度梁的荷載-撓度曲線.Khabbaz等[4]基于高階剪切變形理論研究了功能梯度板的非線性彎曲問題.上述研究均引入了Green應變張量,Green應變張量雖然不會在大變形時產生虛假應變,但沒有與之對應的轉動張量,故該理論是有缺陷的[5].陳至達[5]提出了S-R和分解定理,克服了Green應變張量的缺點,能夠確定唯一的應變.

另外,部分學者基于有限元法對功能梯度板的非線性變形進行了研究.Hassan等[6]研究了SiC-Al功能梯度板的非線性彎曲問題,并計算了塑性區的大小.Nam等[7]基于高階剪切變形理論,利用多邊形有限元法求解了Al-Al2O3板的非線性彎曲問題.利用有限單元法求解大變形問題時,網格會產生嚴重的變形,嚴重影響了計算結果的精度.而無網格法可以部分或者徹底地消除網格,保證結果的精度[8-9].1994年,Belytschko等[10]提出無網格伽遼金(element-free Galerkin,EFG)法后,無網格法得到了迅速的發展,并被應用到眾多領域.葉志強等[11]和劉世宇等[12]利用基于全方向導數的無網格法,對二維翼型繞流問題進行了求解.Li等[13]利用EFG法對金屬成型過程中的接觸沖擊大變形問題進行了分析.張弛等[14]利用無網格法對變厚度功能梯度材料圓板的自由振動問題進行了分析.陳海昌等[15]基于一階剪切變形理論,應用無網格法對不銹鋼-陶瓷功能梯度板彎曲問題進行了計算.陳祖芳等[16]、羅丹[17]和宋彥琦等[18]推導了2D-S-R-EFG離散方程,求解了二維梁結構的非線性變形問題,宋彥琦等[19]又提出了三維情況下的S-R無網格法.

綜上所述,目前研究普遍應用的Green應變張量沒有相應的轉動張量,而有限元法在求解大變形問題時嚴重影響結果的精度.另外,S-R和分解定理在功能梯度板非線性變形的問題中的應用較少.如果將S-R和分解定理與無網格法結合起來,就可以建立更加完善的大變形問題求解方法.因此,本工作基于S-R和分解定理,推導出3D-S-R-EFG離散方程,并利用MATLAB編寫了S-R無網格法程序,對功能梯度材料板的非線性變形彎曲進行了研究.

1 功能梯度材料模型

功能梯度板由陶瓷和金屬混合制成,長度為a,寬度為b,厚度為h,板的上表面受均布荷載q,示意圖如圖1所示.板的材料參數沿板的厚度方向變化,符合冪律型分布,

圖1 功能梯度板示意圖Fig.1 Schematic diagram of functionally graded plate

式中:P為板內任一點的材料參數(彈性模量E、泊松比μ、剪切模量G等);Pc和Pm分別是陶瓷和金屬的材料參數;Vc為陶瓷的體積分數;n為陶瓷的體積分數指數.

圖2給出了體積分數指數n變化時,陶瓷的體積分數隨厚度變化的情況.

圖2 體積分數隨板的厚度變化情況Fig.2 Variation of the volume fraction through the thickness of the plate

2 S-R和分解定理

采用自然拖帶系法[5]描述物體的運動,選取兩個參考系:固定參考系{Xi}(i=1,2,3)和嵌含在變形體中的自然拖帶坐標系{xi}(i=1,2,3).自然拖帶坐標系如圖3所示.質點變形前后的位置矢量分別為r和R,位移矢量為u,

圖3 自然拖帶坐標系Fig.3 The co-moving coordinate

定義變形前后的基標矢量:

對式(2)求導,可得

式中:

其中為一階逆變分量uj對xi的協變導數為初始拖帶系中的第二類Christoffel符號.

由上分析可以得到變形前后基矢的轉換關系為

式中:F是局部坐標系基矢的變換函數;為Kronecker符號.

S-R和分解定理如下:給定一個物理可能的位移函數,此函數在形變體點集域內是單值連續的,處處有一階導數,則此運動變換可以唯一分解為正交與對稱兩個子變換的直和.正交變換體現點集的轉動,而對稱變換體現點集的形變[5],即

式中:S表示應變張量;R表示轉動張量.

3 基于更新拖帶坐標系法的增量變分方程

3.1 增量變分方程

利用虛功率原理并結合增量方法,可得

式中:t+ΔtV j‖i表示在t+Δt時刻的速度t+ΔtV i關于t+Δt時刻拖帶系協變基矢t+Δtgi的協變導數;t+ΔtR為外力虛功,

考慮準靜力學問題,根據速度梯度和應力的線性近似有

將式(12)、(13)代入式(10),并結合Δuj|k=ΔttV j‖k,可得到

3.2 更新拖帶坐標系法

圖4所示為更新拖帶坐標系法.由圖4可知:在初始時刻t0,質點位于點P,其初始拖帶坐標系的基矢為在t時刻,質點移動到點P′處,其初始拖帶系基矢變為定義為t時刻的初始拖帶系的基矢,即在t～t+Δt時刻的增量步上,對初始拖帶系進行了一次更新,則t時刻瞬時拖帶系tgi的tV i‖j可表示為

圖4 更新拖帶坐標系法Fig.4 The updated co-moving coordinate method

將式(15)代入式(14),可得

式中:

引入速率型物性方程

式中:

利用式(17),可得

忽略體力矩的影響,并令初始系為直線直角坐標系,可得[20]

則式(23)可簡化為

式(25)即可作為全局弱勢的無網格法的增量變分方程.

4 基于全局弱勢的無網格Galerkin法

根據無網格Galerkin法,利用移動最小二乘法(moving least squares,MLS)對t時刻初始拖帶系下節點的速度tV進行插值,則有

式中:tφk為t時刻節點k處的MLS形函數;tVk是t時刻k節點的速度參數.采用罰函數法施加邊界條件,對式(25)進行修正,可得

對于tV i‖j,有

將式(26)代入式(27),可得

故可將式(17)寫為

式中:

對于[Bs]和[Br],有

注意到矩陣[Bs]和[Br]中含有,令初始拖帶系和直線直角坐標系同胚,則有這大大降低了計算的復雜度,也體現了更新拖帶坐標系法的優勢.

將式(30)代入式(27),并利用Δt[V]=[Δu],可得

式中:

5 數值算例

5.1 各向同性板的非線性彎曲

表1為本工作中無網格法數值計算的參數設定.為了檢驗本工作的方法和程序的正確性,對受均布荷載的四邊簡支方形薄板進行計算.令FGP兩種材料的彈性模量均為5.38×1010Pa,泊松比均為0.3.根據式(1),該板為各向同性板.板的長度和寬度均為25.4 cm,厚度為2.54 cm,節點分布為13×13×3.該算例被廣泛應用于驗證薄板非線性彎曲計算方法的準確性.

表1 無網格法參數設定Table 1 Essential parameters for S-R-EFG

圖5給出了板中面中點的無量綱撓度與荷載的關系.從圖中可以看出:在初始的小變形階段,3種方法的計算結果吻合度高,不同理論對結果的影響較小,這在一定程度上說明了本方法的準確性;在非線性(大變形)階段,本方法的計算結果與Zhu等[21]和Zhao等[22]的結果相比較小.因為Zhu等[21]采用了一階剪切變形理論和基于滑動Kriging插值的無網格局部Petrov-Galerkin法,而Zhao等[22]采用了基于一階剪切變形理論的無網格里茲法,但一階剪切變形理論會引入一些假設,如假定沿厚度方向的正應變為0,降低了結果的準確性.此外,文獻[21-22]均引入了von K′arm′an應變來體現非線性變形,而von K′arm′an應變從本質上來講,仍然屬于Green應變.Green應變張量沒有推導出相應的轉動張量,是存在理論缺陷的.本工作中的S-R無網格法推導了轉動速率和速度的關系,并將其引入到了系統方程中,考慮了轉動對變形的影響,故能得出更準確的結果.根據以上討論,該算例能夠說明本工作中S-R無網格法的準確性,并且也可以解釋與其他結果的差異.

圖5 各向同性板的無量綱撓度和荷載的關系Fig.5 Non-dimensional central deflection versus load parameter of isotropic plate

5.2 功能梯度板的非線性彎曲

表2為算例中FGP材料的力學參數[21].

表2 功能梯度板的力學參數Table 2 Mechanical parameters of FGP

圖6給出了FGP在四邊簡支約束下無量綱撓度和無量綱荷載之間的關系.與文獻[21-23]中的結果進行對比,可以看出,在采用13×13×3節點進行計算時,本工作的方法能夠得到收斂的結果.與文獻[21-22]相同,文獻[23]中也引入了von K′arm′an應變來體現非線性,故本工作中的結算結果與其他結果的差異也是可以解釋的.

圖6 功能梯度板的無量綱撓度和荷載的關系Fig.6 Non-dimensional central deflection versus load parameter of FGP

下面研究體積分數指數對FGP非線性彎曲的影響.圖7所示為FGP的撓度和體積分數指數n的關系.由圖7可以看出,隨著n的增加,FGP的剛度下降,板的無量綱撓度增加,這與Zhu等[21]、Zhao等[22]和Do等[23]得出的規律是相同的.圖8所示為在不同體積分數指數下板中面中點的無量綱撓度和寬厚比的關系.從圖中可以看出:當FGP的寬厚比從10增加到15時,板的中面撓度迅速減小;但是當板的寬厚比增加到20時,板的中面撓度只產生了很小的變化.因此可以預測,隨著板的寬厚比增加,板的中面無量綱撓度會趨近于一個定值.

圖7 不同體積分數指數下功能梯度板的無量綱撓度和荷載的關系Fig.7 Non-dimensional central deflection versus load parameter of FGP for the various volume fraction exponents

圖8 不同寬厚比情況下功能梯度板的無量綱撓度和荷載的關系Fig.8 Non-dimensional central deflection versus load parameter of FGP for different lengthto-thickness ratios

6 結論

本工作以S-R和分解定理為基礎,基于全局弱勢的無網格伽遼金法得到了三維情況下的S-R-EFG離散方程,利用編制的3D-S-R-EFG程序求解了Ti-6Al-4V-Al2O3板的大撓度彎曲問題,證明了該方法在求解FGP非線性彎曲問題時的收斂性,并研究了體積分數指數n和寬厚比對FGP的無量綱撓度的影響.研究發現,板的無量綱撓度會隨著n的增加而增加,而隨著板的寬厚比的增加,板的無量綱撓度會趨近一個定值.這表明了S-R無網格法在求解FGP大撓度彎曲問題時的合理性.