宣領寬，明平劍，龔京風，張文平，韓春旭

(1.哈爾濱工程大學動力與能源工程學院，150001哈爾濱;2.中國艦船研究設計中心，430064武漢)

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

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

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

1 數(shù)學模型

1.1 控制方程

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

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

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

式中:D為彈性矩陣;ε為應變向量．

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

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

邊界條件可表示為

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

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

1.2 數(shù)值離散

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

圖1 三維非結構網(wǎng)格

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

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

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

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

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

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

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

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

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

2 數(shù)值應用

2.1 立方體結構

靜態(tài)分析如圖2所示，正交各向異性FGM、邊長為10 m的立方體結構，來驗證本文方法的正確性.邊界條件為:x=0 m、y=0 m與z=0 m平面均為法向簡支，z=10 m平面受到法向均布拉力σ=1 N/m2，其他各面自由.立方體為正交各向異性FGM材料，材料對應的彈性矩陣中元素為: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.計算域采用四面體網(wǎng)格劃分，表1給出本文采用的4種網(wǎng)格的相關信息.圖3為AB邊上的位移ux與uy，F(xiàn)VM、FEM結果由本文方法、ANSYS 12.0采用相同網(wǎng)格Grid 4得到，可以看出本文FVM結果與FEM及文獻［11］中的MLPG結果吻合良好.以Grid 4的計算結果作為參考，將基于4種網(wǎng)格所得的邊AB上的節(jié)點位移uy列于表2，可以看出隨著網(wǎng)格數(shù)的增加，Grid3、Grid4的計算結果基本吻合;采用相同網(wǎng)格時，從表2可以看出FVM與FEM的計算結果幾乎完全一致.

圖2 邊長為10 m的立方體

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

圖3 邊AB上的位移

表2 不同網(wǎng)格時位移uy計算結果對比(10-8m)

2.2 轉盤結構

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

圖4 轉盤結構

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

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

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

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

圖5 t＝0.015 s時FVM計算結果與 FEM(ANSYS)計算結果對比

圖7為3種均勻橫觀各向同性材料不同點的時間響應，本文FVM計算結果與ANSYS計算結果吻合良好.從圖7可以看出彈性模量的減少均會造成對應方向上的位移增大.

圖6 中位面上各監(jiān)測點徑向應力σr隨時間的變化

圖7 3種材料不同點的時間響應

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

表3 計算消耗對比

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

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

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

3 結論

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

2)與ANSYS相比，本文方法需要較小的時間步長，但占用較小內存且具有較快的計算速度.

3)研究材料屬性沿徑向變化對轉盤的影響.結果表明:隨著徑向彈性模量的增大，轉盤的徑向位移減小;隨著徑向密度的增大，轉盤的徑向位移增大;當徑向密度與彈性模量按相同趨勢增大時，轉盤的徑向位移增大;與彈性模量相比，密度沿徑向變化對轉盤徑向位移變化有更大的影響.

［1］王保林，韓杰才，張幸紅.非均勻材料力學［M］.北京:科學出版社，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)板的軸對稱問題［D］.杭州:浙江大學，2007.

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

［6］黃德進.各向異性功能梯度平面梁的彈性力學解［D］.杭州:浙江大學，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］程祖依.彈性動力學基礎［M］.武漢:中國地質大學出版社，1989:98-106.

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

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

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