基于S-R和分解定理的無網(wǎng)格法在功能梯度板中的應(yīng)用

宋彥琦,石博康,李向上

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

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

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

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

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

1 功能梯度材料模型

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

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

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

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

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

2 S-R和分解定理

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

圖3 自然拖帶坐標(biāo)系Fig.3 The co-moving coordinate

定義變形前后的基標(biāo)矢量:

對式(2)求導(dǎo),可得

式中:

其中為一階逆變分量uj對xi的協(xié)變導(dǎo)數(shù)為初始拖帶系中的第二類Christoffel符號.

由上分析可以得到變形前后基矢的轉(zhuǎn)換關(guān)系為

式中:F是局部坐標(biāo)系基矢的變換函數(shù);為Kronecker符號.

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

式中:S表示應(yīng)變張量;R表示轉(zhuǎn)動張量.

3 基于更新拖帶坐標(biāo)系法的增量變分方程

3.1 增量變分方程

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

式中:t+ΔtV j‖i表示在t+Δt時刻的速度t+ΔtV i關(guān)于t+Δt時刻拖帶系協(xié)變基矢t+Δtgi的協(xié)變導(dǎo)數(shù);t+ΔtR為外力虛功,

考慮準靜力學(xué)問題,根據(jù)速度梯度和應(yīng)力的線性近似有

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

3.2 更新拖帶坐標(biāo)系法

圖4所示為更新拖帶坐標(biāo)系法.由圖4可知:在初始時刻t0,質(zhì)點位于點P,其初始拖帶坐標(biāo)系的基矢為在t時刻,質(zhì)點移動到點P′處,其初始拖帶系基矢變?yōu)槎x為t時刻的初始拖帶系的基矢,即在t～t+Δt時刻的增量步上,對初始拖帶系進行了一次更新,則t時刻瞬時拖帶系tgi的tV i‖j可表示為

圖4 更新拖帶坐標(biāo)系法Fig.4 The updated co-moving coordinate method

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

式中:

引入速率型物性方程

式中:

利用式(17),可得

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

則式(23)可簡化為

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

4 基于全局弱勢的無網(wǎng)格Galerkin法

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

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

對于tV i‖j,有

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

故可將式(17)寫為

式中:

對于[Bs]和[Br],有

注意到矩陣[Bs]和[Br]中含有,令初始拖帶系和直線直角坐標(biāo)系同胚,則有這大大降低了計算的復(fù)雜度,也體現(xiàn)了更新拖帶坐標(biāo)系法的優(yōu)勢.

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

式中:

5 數(shù)值算例

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

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

表1 無網(wǎng)格法參數(shù)設(shè)定Table 1 Essential parameters for S-R-EFG

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

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

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

表2為算例中FGP材料的力學(xué)參數(shù)[21].

表2 功能梯度板的力學(xué)參數(shù)Table 2 Mechanical parameters of FGP

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

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

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

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

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

6 結(jié)論

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