三維各向異性功能梯度材料的有限體積法

宣領(lǐng)寬，明平劍，龔京風(fēng)，張文平，韓春旭

(1.哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院，150001哈爾濱;2.中國艦船研究設(shè)計(jì)中心，430064武漢)

功能梯度材料(FGM)一般是由兩種或兩種以上材料復(fù)合而成，各組分材料的體積分?jǐn)?shù)在空間位置上是連續(xù)變化的，其宏觀材料特性表現(xiàn)出梯度(逐漸變化)的性質(zhì)［1］.FGM擁有高強(qiáng)度、耐熱性、耐磨性等優(yōu)勢(shì)，文獻(xiàn)［2-3］對(duì)FGM的發(fā)展進(jìn)行了詳細(xì)介紹.由于實(shí)際工程的需要，F(xiàn)GM可能是各向同性的也有可能是各向異性的.2000年以后，國內(nèi)外學(xué)者圍繞各向異性FGM的力學(xué)問題開展了相關(guān)研究.常用計(jì)算方法有解析法與數(shù)值法.李翔宇［4］給出了軸對(duì)稱的橫觀各向同性FGM圓板和環(huán)板的靜力學(xué)問題解析解，黃德進(jìn)等［5-6］推導(dǎo)出各向異性FGM平面梁的彈性靜力學(xué)解析解;雖然解析法能夠給出精確解，并且能夠靈活分析規(guī)律性的內(nèi)在聯(lián)系，但難以處理工程上的復(fù)雜結(jié)構(gòu)問題，因此在工程實(shí)際中經(jīng)常采用數(shù)值法.

目前，有限元法(FEM)是計(jì)算力學(xué)領(lǐng)域發(fā)展最成熟、應(yīng)用最廣泛的數(shù)值方法，近年來發(fā)展較為緩慢，與此同時(shí)，學(xué)者們?nèi)匀徊煌５匕l(fā)展新的數(shù)值方法.邊界元法(BEM)曾在各向同性線彈性問題中獲得了成功，因此有學(xué)者試圖將BEM應(yīng)用于各向異性、非均勻材料的線彈性問題.BEM需要求解格林函數(shù)，而材料的各向異性會(huì)使彈性矩陣中的系數(shù)增加，導(dǎo)致構(gòu)造格林函數(shù)變得更加復(fù)雜，而且材料的非均性更會(huì)加劇這一問題.Albuquerque等［7］將BEM應(yīng)用于二維各向異性材料的動(dòng)力學(xué)問題，Criado等［8-9］導(dǎo)出指數(shù)形式的FGM線彈性體的格林函數(shù)，并基于BEM求解指數(shù)形式FGM的各向同性線彈性問題.然而到目前為止，尚未看到有學(xué)者將BEM應(yīng)用于各向異性FGM的線彈性問題.此外，無網(wǎng)格法由于良好的適應(yīng)性及避免劃分網(wǎng)格等優(yōu)勢(shì)得到了廣泛關(guān)注，Sladek等［10-11］將無網(wǎng)格伽遼金法(MLPG)應(yīng)用于求解二維、三維各向異性FGM的線彈性問題，然而，MLPG法對(duì)于復(fù)雜工程問題也存在計(jì)算量大、效率低等缺點(diǎn).

另一方面，在流體動(dòng)力學(xué)(CFD)領(lǐng)域內(nèi)十分流行的有限體積法(FVM)，也逐漸被應(yīng)用于結(jié)構(gòu)問題的求解，這么做的最終目標(biāo)是采用統(tǒng)一方法求解多物理場(chǎng)(如流固耦合)問題，避免混合方法帶來的一系列問題(如收斂困難)［12-13］.一般認(rèn)為FEM對(duì)傳統(tǒng)的結(jié)構(gòu)自伴問題具有更高的精度，然而對(duì)于偏微分方程二者精度區(qū)別并不大［13］，而且在許多應(yīng)用中二者是等價(jià)的［14］.事實(shí)上結(jié)構(gòu)與流體控制方程相同，不同的只是本構(gòu)方程［15］.FVM已被成功應(yīng)用于結(jié)構(gòu)問題，Wheel［16-17］在分析應(yīng)力集中問題時(shí)，基于相同網(wǎng)格類型FVM獲得結(jié)果的精度比FEM高，并且基于FVM解決了FEM難以解決的彈性力學(xué)中的“閉鎖”現(xiàn)象;Fallah［18］、Demird?ic＇等［19］也將 FVM 應(yīng)用于求解結(jié)構(gòu)問題，但到目前，尚未看到將FVM應(yīng)用于各向異性FGM線彈性問題的報(bào)道.本文采用非結(jié)構(gòu)有限體積法研究正交各向異性立方體結(jié)構(gòu)的靜態(tài)特性，將計(jì)算結(jié)果與其他數(shù)值方法結(jié)果對(duì)比，驗(yàn)證本文方法的正確性;模擬橫觀各向同性FGM轉(zhuǎn)盤結(jié)構(gòu)從啟動(dòng)到穩(wěn)定過程中的動(dòng)態(tài)特性，研究密度、彈性模量沿徑向變化對(duì)其力學(xué)性能的影響.

1 數(shù)學(xué)模型

1.1 控制方程

各向異性線彈性體的平衡方程［20］可表示為

各向異性線彈性體的本構(gòu)方程［21］的張量形式可表示為

式中:Cijkl為彈性張量，共81個(gè)分量;εkl為應(yīng)變張量．式(2)的矩陣形式為

式中:D為彈性矩陣;ε為應(yīng)變向量．

正交各向異性材料的彈性矩陣為

彈性軸與z軸一致的橫觀各向同性材料彈性矩陣D可表示為

邊界條件可表示為

式中:Γd為給定位移的邊界;Γt為給定牽引力的邊界．T可表示為

式中n=(nx，ny，nz)為邊界的單位外法線矢量．

1.2 數(shù)值離散

三維計(jì)算域采用四面體網(wǎng)格劃分，如圖1所示為任意的四面體網(wǎng)格單元C1234，O為C1234的重心，E12、E13、E14分別為邊 12、13、14 的中點(diǎn)，O2、O3、O4分別為 134、124、123 的中點(diǎn)，空間多邊形E12O4OO3、E13O2OO4、E14O3OO2為圍繞節(jié)點(diǎn) 1 的控制體積在C1234中的邊界，在圍繞節(jié)點(diǎn)1的控制體上對(duì)式(1)進(jìn)行體積分可得

圖1 三維非結(jié)構(gòu)網(wǎng)格

應(yīng)力與材料屬性均定義在單元中心上，并且認(rèn)為其在單元C1234內(nèi)部是均勻的，對(duì)式(4)的右端第1項(xiàng)應(yīng)用高斯定理可得

式中:四面體單元1234對(duì)節(jié)點(diǎn)1的貢獻(xiàn)可以表示為σ1·A1，其中A1為空間多邊形E12O4E13O2E14O3E12(由3個(gè)四邊形組成)的面積矢量．因此有

對(duì)四面體單元C1234的其他節(jié)點(diǎn)2，3，4，同理可得

每個(gè)單元內(nèi)的應(yīng)力可由應(yīng)變通過式(3)求得，在任意單元C1234中，應(yīng)變?yōu)椋?2］

式中:φ為ux、uy、uz;a1=;b1=

;(x2，y2，z2)、(x3，y3，z3)、(x4，y4，z4)分別為節(jié)點(diǎn)2、3、4 的空間坐標(biāo);Vc為四面體C1234的體積．式(5)～(7)不僅可以用來計(jì)算Ai，而且其相當(dāng)于得到應(yīng)變向量ε，可以利用式(3)計(jì)算單元內(nèi)的應(yīng)力，因此可以將ai、bi、ci一次計(jì)算并儲(chǔ)存起來，這樣不僅使用起來方便，而且能夠減少計(jì)算量、提高計(jì)算速度．

對(duì)于位移邊界Γd可將邊界節(jié)點(diǎn)的位移直接等于給定位移，對(duì)于固支邊界則有=0．

若考慮靜力學(xué)問題，則式(8)的左端為0.將式(8)整理為以位移為求解變量的線性方程組，可以直接求解或迭代求解，這里采用INTEL IMSL的LSARG求解器直接求解.

由于本文計(jì)算得到的應(yīng)力均位于單元上，若想知道某一節(jié)點(diǎn)(如圖1的節(jié)點(diǎn)1)上的應(yīng)力，可由下式獲得

2 數(shù)值應(yīng)用

2.1 立方體結(jié)構(gòu)

靜態(tài)分析如圖2所示，正交各向異性FGM、邊長(zhǎng)為10 m的立方體結(jié)構(gòu)，來驗(yàn)證本文方法的正確性.邊界條件為:x=0 m、y=0 m與z=0 m平面均為法向簡(jiǎn)支，z=10 m平面受到法向均布拉力σ=1 N/m2，其他各面自由.立方體為正交各向異性FGM材料，材料對(duì)應(yīng)的彈性矩陣中元素為:c11=236.4×102MPa，c12=63.64 MPa，c13=18.18 MPa;c22=119.7 MPa，c23=15.15 MPa;c33=42.42×(1+z/10)MPa;c44=80 MPa;c55=c66=16 MPa.計(jì)算域采用四面體網(wǎng)格劃分，表1給出本文采用的4種網(wǎng)格的相關(guān)信息.圖3為AB邊上的位移ux與uy，F(xiàn)VM、FEM結(jié)果由本文方法、ANSYS 12.0采用相同網(wǎng)格Grid 4得到，可以看出本文FVM結(jié)果與FEM及文獻(xiàn)［11］中的MLPG結(jié)果吻合良好.以Grid 4的計(jì)算結(jié)果作為參考，將基于4種網(wǎng)格所得的邊AB上的節(jié)點(diǎn)位移uy列于表2，可以看出隨著網(wǎng)格數(shù)的增加，Grid3、Grid4的計(jì)算結(jié)果基本吻合;采用相同網(wǎng)格時(shí)，從表2可以看出FVM與FEM的計(jì)算結(jié)果幾乎完全一致.

圖2 邊長(zhǎng)為10 m的立方體

表1 計(jì)算網(wǎng)格信息

圖3 邊AB上的位移

表2 不同網(wǎng)格時(shí)位移uy計(jì)算結(jié)果對(duì)比(10-8m)

2.2 轉(zhuǎn)盤結(jié)構(gòu)

研究如圖4所示圍繞z軸以角速度ω旋轉(zhuǎn)的轉(zhuǎn)盤，幾何尺寸為:內(nèi)徑rin=0.1 m，外徑rout=0.24 m，厚度h=0.06 m，模擬轉(zhuǎn)盤從啟動(dòng)到穩(wěn)定過程中的動(dòng)態(tài)特性，角速度按照式(9)規(guī)律變化.轉(zhuǎn)盤內(nèi)側(cè)固支，其他自由，采用16 461個(gè)四面體網(wǎng)格劃分，節(jié)點(diǎn)數(shù)為3 872，最短邊長(zhǎng) ΔLmin=8.711 mm，時(shí)間步長(zhǎng)可依據(jù) CFL 條件［21］Δt≤ΔLmin/(cp)max選?。?/p>

圖4 轉(zhuǎn)盤結(jié)構(gòu)

考慮以下3種均勻橫觀各向同性材料的轉(zhuǎn)盤，其材料屬性分別為

以下計(jì)算中如無特殊說明，時(shí)間步長(zhǎng)均采用Δt=1.5 μs.3種材料的其他屬性均相同，泊松比vxy=vyz=vzx=0.3，密度 ρ=8 000 kg/m3，對(duì)應(yīng)的剪切模量可計(jì)算出:

由旋轉(zhuǎn)產(chǎn)生的單位質(zhì)量體力可表示為

驗(yàn)證本文方法計(jì)算動(dòng)態(tài)問題時(shí)的正確性.圖5為轉(zhuǎn)盤材料為M1時(shí)、t=0.015 s時(shí)，F(xiàn)VM計(jì)算結(jié)果與ANSYS計(jì)算結(jié)果的對(duì)比，可看出本文方法與ANSYS計(jì)算結(jié)果吻合良好，轉(zhuǎn)盤的徑向位移ur、軸向位移uz與徑向應(yīng)力σr均關(guān)于其中位面z=0.03 m對(duì)稱，徑向最大位移位于轉(zhuǎn)盤的外徑(r=0.24 m)，軸向的最大位移大概位于r=0.131 m處的上、下表面上，徑向最大應(yīng)力位于轉(zhuǎn)盤的內(nèi)徑r=0.1 m.圖6為中位面上各監(jiān)測(cè)點(diǎn)徑向應(yīng)力σr隨時(shí)間的變化，計(jì)算結(jié)果均與FEM吻合良好，隨著速度的均勻增加，應(yīng)力增加的越來越快，當(dāng)速度保持不變后，各點(diǎn)的應(yīng)力也保持不變.

圖5 t＝0.015 s時(shí)FVM計(jì)算結(jié)果與 FEM(ANSYS)計(jì)算結(jié)果對(duì)比

圖7為3種均勻橫觀各向同性材料不同點(diǎn)的時(shí)間響應(yīng)，本文FVM計(jì)算結(jié)果與ANSYS計(jì)算結(jié)果吻合良好.從圖7可以看出彈性模量的減少均會(huì)造成對(duì)應(yīng)方向上的位移增大.

圖6 中位面上各監(jiān)測(cè)點(diǎn)徑向應(yīng)力σr隨時(shí)間的變化

圖7 3種材料不同點(diǎn)的時(shí)間響應(yīng)

對(duì)比FVM與ANSYS計(jì)算消耗，以轉(zhuǎn)盤材料M3為例，由ANSYS模態(tài)計(jì)算可得轉(zhuǎn)盤的最小固有頻率為1 665 Hz，則其ANSYS所采用時(shí)間步長(zhǎng)應(yīng)滿足 Δt≤1/(20f)≈30 μs.表3為 FVM 與ANSYS計(jì)算消耗對(duì)比.在網(wǎng)格與時(shí)間步長(zhǎng)均相同的前提下，F(xiàn)VM消耗較小的內(nèi)存并獲得較快的計(jì)算速度，其主要原因是FVM采用顯格式推進(jìn)求解且不需要求解大型線性方程組，而ANSYS在每一個(gè)時(shí)間內(nèi)都需要建立存儲(chǔ)剛度矩陣并求解大型線性方程組，這將會(huì)占用大量的內(nèi)存及計(jì)算時(shí)間.在該算例中，雖然顯格式推進(jìn)求解限制 FVM采用較小的時(shí)間步長(zhǎng)，但與ANSYS采用較大的時(shí)間步長(zhǎng)相比，F(xiàn)VM仍然消耗較少的內(nèi)存及計(jì)算時(shí)間.

表3 計(jì)算消耗對(duì)比

最后，考慮材料屬性沿徑向變化對(duì)轉(zhuǎn)盤的影響，材料屬性的變化規(guī)律可表示為

式中:β1、β2分別為梯度常數(shù);ρ0=8 000 kg/m3;E0=100 GPa;Ez=200 GPa，其他屬性保持不變.考慮4個(gè)不同梯度常數(shù)即0.0，0.1，0.3，0.5，當(dāng)β1=β2=0.0對(duì)應(yīng)均勻橫觀各向同性材料M3.如圖8為中位面上r=0.24 m處的徑向位移，本文FVM計(jì)算結(jié)果與ANSYS計(jì)算結(jié)果吻合良好.當(dāng)β2保持不變，可以看出隨著徑向彈性模量或梯度常數(shù)β1的增大，轉(zhuǎn)盤的徑向位移減小(圖8(a));當(dāng)β1保持不變，隨著徑向密度或梯度常數(shù)β2的增大，轉(zhuǎn)盤的徑向位移增大(圖8(b));當(dāng)徑向密度與彈性模量按相同趨勢(shì)增大時(shí)，轉(zhuǎn)盤的徑向位移增大(圖8(c));與彈性模量變化相比，密度沿徑向變化對(duì)轉(zhuǎn)盤徑向位移變化有更大的影響.

圖8 中位面上r＝0.24 m處的徑向位移

3 結(jié)論

1)計(jì)算結(jié)果與其他數(shù)值計(jì)算結(jié)果吻合良好，表明本文FVM的正確性.

2)與ANSYS相比，本文方法需要較小的時(shí)間步長(zhǎng)，但占用較小內(nèi)存且具有較快的計(jì)算速度.

3)研究材料屬性沿徑向變化對(duì)轉(zhuǎn)盤的影響.結(jié)果表明:隨著徑向彈性模量的增大，轉(zhuǎn)盤的徑向位移減小;隨著徑向密度的增大，轉(zhuǎn)盤的徑向位移增大;當(dāng)徑向密度與彈性模量按相同趨勢(shì)增大時(shí)，轉(zhuǎn)盤的徑向位移增大;與彈性模量相比，密度沿徑向變化對(duì)轉(zhuǎn)盤徑向位移變化有更大的影響.

［1］王保林，韓杰才，張幸紅.非均勻材料力學(xué)［M］.北京:科學(xué)出版社，2003.

［2］MARKWORTH A J，RAMESH K S，PARKS W P.Review:modeling studies applied to functionally graded materials［J］.Journal of Materials Science，1995，30(9):2183-2193.

［3］BIRMAN V，BYRD L W.Modeling and analysis of functionally graded materialsand structures ［J］.Applied Mechanics Reviews，2007，60(5):83-96.

［4］李翔宇.橫觀各向同性功能梯度圓板和環(huán)板的軸對(duì)稱問題［D］.杭州:浙江大學(xué)，2007.

［5］黃德進(jìn)，丁皓江，陳偉球.線性分布載荷作用下功能梯度各向異性懸臂梁的解析解［J］.應(yīng)用數(shù)學(xué)和力學(xué)，2007，28:763-768.

［6］黃德進(jìn).各向異性功能梯度平面梁的彈性力學(xué)解［D］.杭州:浙江大學(xué)，2009.

［7］ALBUQUERQUE E L，SOLLERO P，ALIABADI M H.The boundary element method applied to time dependent problems in anisotropic materials［J］.International Journal of Solids and Structures，2002，39(5):1405-1422.

［8］CRIADO R，ORTIZ J E，MANTIC V，et al.Green’s function evaluation for three-dimensional exponentially graded elasticity［J］.International Journal for Numerical Methods in Engineering，2008，74(10):1560-1591.

［9］CRIADO R，ORTIZ J E，MANTIC V，et al.Boundary elementanalysis of three-dimensional exponentially graded isotropic elastic solids［J］.CMES:Computer Modeling in Engineering＆Sciences，2007，22(2):151-164.

［10］SLADEK J，SLADEK V，ZHANG C H.Stress analysis in anisotropic functionally graded materials by the MLPG method［J］. Engineering Analysis with Boundary Elements，2005，29(6):597-609.

［11］SLADEK J，SLADEK V，SOLEK P.Elastic analysis in 3D anisotropic functionally graded solids by the MLPG［J］.CMES:Computer Modeling in Engineering＆Sciences，2009，12(1):35-36.

［12］MANEERATANA K.Development of the finite volume method for non-linear structural applications［D］.London:The University of London，2000.

［13］SLONE A K，BAILEY C，CROSS M.Dynamic solid mechanics using finite volume methods［J］.Applied Mathematical Modelling，2003，27(2):69-87.

［14］IDELSOHN S R，O?ATE E.Finite volumes and finite elements:two‘good friends’［J］.International Journal for Numerical Methods in Engineering，1994，37(19):3323-3341.

［16］WHEEL M A.A geometrically versatile finite volume formulation for plane elastostatic stress analysis［J］.Journal of Strain Analysis for Engineering Design，1996，31(2):111-116.

［17］WHEEL M A.A mixed finite volume formulation for determining the small strain deformation of incompressible materials［J］.International Journal for Numerical Methods in Engineering，1999，44(12):1843-1861.

［18］FALLAH N.A cell vertex and cell centered finite volume method for plate bending analysis［J］.Computational Methods Applied in Mechanical Engineering，2004，193(33/35):3457-3470.

［20］程祖依.彈性動(dòng)力學(xué)基礎(chǔ)［M］.武漢:中國地質(zhì)大學(xué)出版社，1989:98-106.

［21］陳明祥.彈塑性力學(xué)［M］.北京:科學(xué)出版社，2007:61-68.

［22］甘舜仙.有限元技術(shù)與程序［M］.北京:北京理工大學(xué)出版社，1988:152-154.

［23］林曉.固體應(yīng)力波的數(shù)值解法［M］.北京:科學(xué)出版社，2008.