用于非均質復合材料應力分析的交錯網格有限體積法

宣領寬,龔京風,張文平,明平劍

(哈爾濱工程大學動力與能源工程學院, 150001, 哈爾濱)

用于非均質復合材料應力分析的交錯網格有限體積法

宣領寬,龔京風,張文平,明平劍

(哈爾濱工程大學動力與能源工程學院, 150001, 哈爾濱)

針對非均質復合材料的應力問題,發展了一種交錯網格有限體積法(SCV-FVM)。該方法基于非結構網格離散線彈性平衡方程,采用交錯網格技術將材料的空間變化引入離散過程,從而不需要顯式處理復合材料交界面。用SCV-FVM對宏觀非均質復合材料應力場進行了數值模擬,結果與理論解吻合良好。與其他數值結果的對比表明,SCV-FVM能夠避免材料物性突變引起的牽引力方向的應力數值波動及不連續現象,但是難以捕捉垂直于牽引力方向的應力跳躍現象,可以通過加密交界面網格來改善計算結果。用SCV-FVM對微觀非均質復合材料應力場進行了數值模擬,結果不存在由物性參數空間變化引起的數值不連續現象及應力集中現象,表明SCV-FVM適合對微觀非均質材料進行應力分析。

非均質復合材料;層合材料;功能梯度材料;交錯網格技術;有限體積法

層合材料、功能梯度材料、涂層材料等復合材料已被廣泛應用于工程實際中。采用不同的制備工藝得到的復合材料可能存在界面形貌的波動[1]、微尺度顆粒的隨機分布[2-3]等。

關于均質(即不存在物性突變)復合材料熱力性能的研究已有較多報道[4-14]。非均質復合材料包括微觀非均質結構,如具有微結構的功能梯度材料;宏觀非均質結構,如層合材料、包含問題等。文獻[15]指出,基于連續介質理論的有限元法不能有效地直接求解非均質問題,如無限體中含微觀結構的問題。為此,文獻[15]建立了基于均勻化理論的確定復合材料結構應力場的方法,用均質的宏觀結構和非均質的具有周期性分布的微觀結構來描述原結構。文獻[3]基于高階功能梯度理論(HOTFGM),對具有微結構的熱障涂層的熱力性能進行了數值分析,并指出:基于有限元法得到的應力場存在不合理的應力集中現象。另一方面,文獻[4-6]基于HOTFGM發展了有限體積理論(FVT),并用于研究經典包含問題,結果表明:用FVT計算得到的沿牽引力方向的應力與理論解吻合良好,但計算得到的垂直于牽引力方向的應力存在數值不連續現象,并且難以得到收斂解。文獻[16]根據分界面法向應力連續和位移切向梯度連續的情況,提出了3種顯式處理材料分界面的途徑,改進了格心型有限體積法(CC-FVM),計算結果表明,改進后的CC-FVM能夠避免不合理的應力波動。

作者曾提出了一種新的數值方法——交錯網格有限體積法(SCV-FVM)。該方法采用網格有限體積法(CV-FVM)離散控制方程,利用交錯網格技術將物性參數的空間變化引入離散過程。SCV-FVM現已被成功應用于求解均質功能梯度材料及層合材料的熱傳導問題[7]。本文進一步將SCV-FVM推廣應用于宏觀和微觀非均質復合材料的應力問題研究。

1 數學模型

1.1 控制方程

考慮二維線彈性材料的靜態應力問題。根據控制體Ω內的力平衡建立控制方程

(1)

式中:σ為應力張量;g為單位質量體積力矢量;V為Ω的體積;S為Ω的邊界;n為垂直于邊界S的單位外法線矢量。將方程(1)寫為張量形式

(2)

式中:σij為σ在垂直于i方向的微元面沿j方向的分量;gi為g沿i方向的分量;nj為n沿j方向的分量。

線彈性體的本構方程為

σij=2Gεij+λεkkδij

(3)

式中:δij為克羅尼克爾符號,當i=j時δij=1,當i≠j時δij=0;拉姆系數G與λ可由彈性模量E和泊松比μ根據式(4)和式(5)計算:

對于平面應變

(4)

對于平面應力

(5)

εij為應變張量ε的分量,其表達式為

(6)

將式(3)和式(6)代入式(2),得到待解控制方程

(7)

式中:ui為位移矢量u沿i方向的分量。

對于固支邊界SD,位移為0;對于載荷邊界SN,給定邊界載荷矢量f0;自由邊界SF是載荷邊界的特殊情況,即f0=0;簡支邊界可以由固支邊界和自由邊界組合得到。

1.2 數值離散

采用三節點三角形(T3)和四節點四邊形(Q4)網格單元劃分計算域。SCV-FVM與CC-FVM最基本的區別在于控制體的建立。如圖1所示,CC-FVM以網格單元作為控制體,變量定義在單元中心(實心圓點),如網格單元15,而SCV-FVM則圍繞單元節點(空心圓點)依次連接相鄰單元中心和邊長中點(實心圓點)建立控制體,網格單元由實線圍成,控制體由虛線圍成,如圍繞內部節點a建立的控制體1-2-3-4-5-6-7-8-9-10,圍繞邊界節點b建立的控制體11-12-13-b-14-4-3-2。

基于交錯網格的思想,將待解變量定義在單元節點上,物性參數定義在單元中心。假設待解變量在控制體內均勻分布,物性參數在單元內均勻分布,則物性參數在控制體內是變化的(如圖1所示),從而可將物性參數的空間變化自然地引入離散過程。

圖1 控制體示意圖

基于CV-FVM離散方程(7),可得

(8)

式中:nc為當前節點周圍的單元總數;nα為第α個單元內的節點數;Nα β為型函數,型函數的導數在不同單元中的積分可參考文獻[7];下角標α表示第α個單元中心的變量值,αβ表示第α個單元內第β個節點上的變量值。

對于固支邊界,保持節點位移為0。對于載荷邊界,將邊界力代入方程(8),得

(9)

耦合求解不同方向的位移,從而可以一次計算得到整個位移場,避免迭代。

2 方法驗證及應用

2.1 雙層材料圓盤問題

考慮如圖2所示的雙層材料圓盤[16],采用平面應力假設。內表面半徑ri=0.05m,受到的均勻壓力pi=1 MPa;外表面半徑ro=2ri,為自由邊界。兩層材料分界面的半徑rm=1.4ri;內層材料的泊松比μi=0.35,外層材料的泊松比μo=0.3;外層材料的彈性模量Eo是內層材料彈性模量Ei的10倍。由于對稱性,取1/4圓盤作為計算域(如圖2a所示),徑向劃分25個均勻網格,周向劃分60個均勻網格。在文獻[16]中,計算模型徑向有50個均勻網格,周向有120個均勻網格。

(a)幾何模型

(b)非均勻網格

圖3為采用傳統CC-FVM、文獻[16]的改進型CC-FVM和本文的SCV-FVM獲得的雙層材料圓盤應力計算結果比較。從圖3a可以看到,傳統CC-FVM得到的徑向應力σr在分界面附近存在不正確的數值波動,由改進型CC-FVM得到的結果與理論解吻合良好,而本文的SCV-FVM不需要特殊處理分界面,得到的σr不存在數值波動。

圖3b為采用不同方法得到的周向應力σθ曲線。由于在分界面存在材料突變,σθ曲線應有跳躍現象。由改進型SCV-FVM得到的σθ在遠離分界面的區域與理論解一致,但在分界面上計算值存在一定誤差。CC-FVM以位移為求解變量,根據本構關系計算應力場,空間一個點僅計算一個應力值,因此得到的分界面上的應力是一個平均值。當細化交界面附近網格,即采用如圖2b所示的非均勻網格時,計算得到的σθ曲線有明顯改善。非均勻網格與均勻網格采用相同的網格數,非均勻網格的最小網格尺寸為1 mm,在交界面處,交界面兩側沿徑向的網格數分別為10和15,所以交界面附近網格更細密(見圖2b)。

(a)徑向應力σr

(b)周向應力σθ

2.2 包含問題

考慮如圖4a所示的經典包含問題[1,6]。方形結構(2L=30 m)基體為環氧樹脂,物性參數為Em=4.9 GPa,μm=0.34。包含區域(r=1 m)為玻璃纖維,物性參數為Ef=69.0 GPa,μf=0.2。方形結構的左面和右面受到沿x方向的均勻拉力p=100 MPa,上、下表面自由。由于結構的對稱性,取其1/4作為計算域(見圖4b),x=0和y=0為簡支邊界,采用平面應變假設。

(a)幾何結構及邊界條件

(b)計算域劃分

(a)σx(MPa)

(b)σy(MPa)

(a)σx曲線

(b)σy曲線

文獻[6]計算了完整的方形結構,采用4×2 600個四邊形網格劃分計算域,而本文僅采用1 200個四邊形網格劃分計算域,如圖4b所示,網格分布與文獻[6]的類似。圖5為基于本文SCV-FVM計算得到的應力云圖,圖6為基于不同方法計算得到的應力曲線。由圖6可見:采用有限元法得到的應力曲線不連續;采用FVT得到的應力曲線沿牽引力方向是連續的(見圖6a),但垂直于牽引力方向則是不連續的(見圖6b);采用本文SCV-FVM得到的應力曲線與理論解吻合良好,不存在數值不連續現象。需要指出的是,在環氧樹脂和玻璃纖維的分界面存在物性參數的突變,垂直于牽引力方向的應力存在跳躍。與2.1節中的問題類似,因為SCV-FVM計算得到的應力在分界面上是一個平均值,因此采用SCV-FVM計算應力跳躍問題時需要加密分界面兩側的網格。

2.3 微結構涂層問題

利用SCV-FVM求解如圖7所示的4層微結構涂層的應力場。用均勻四邊形網格劃分計算域,網格為邊長20 μm的正方形。為了避免結構的剛性平移,底面中間3個節點固支,其余節點簡支。涂層頂端受到沿x方向變化的壓力,兩側自由,采用平面應變條件。

圖7 4層微結構涂層示意圖

涂層長L=2 mm,寬W=1 mm,共有4層。底層L1為純Fe,厚度Δy1=0.2 mm;第2層L2為黏結劑CoCrAlY,厚度Δy2=0.1 mm;第4層為純ZrO2,厚度Δy4=0.2 mm。各涂層材料的物性參數見表1。第3層為梯度層,厚度Δy3=0.5mm,材料由CoCrAlY逐漸變為ZrO2。ZrO2的質量分數沿y軸變化

(10)

則梯度層的物性參數

P=PZφZ+PC(1-φZ)

(11)

式中:P代表表1中的參數k、E、μ;下角標Z代表ZrO2,C代表CoCrAlY;指數m控制梯度層中材料質量分數的變化規律。

表1 涂層材料的物性參數

(a)m=2 (b)m=4

(a)位移uy

(b)應力σx

梯度層L3中ZrO2顆粒隨機分布,但其統計平均質量分數滿足式(10)。假設顆粒為邊長20μm的正方形,與網格尺寸相同,考慮如圖8所示的2種梯度層結構。圖9為計算得到的位移uy和應力σx的云圖,圖10為計算得到的應力沿y=0.5mm的分布曲線。由圖9和圖10可知,由于微結構顆粒的隨機分布,位移和應力在梯度層存在不均勻分布,位移場的不均勻度小于應力場的不均勻度。另外,隨著m的增大,涂層的變形和應力幅值減小,同時不均勻度降低。

由于微尺度顆粒的存在,使梯度層的物性參數變化劇烈。文獻[2]中采用有限元法計算類似的微尺度問題時,存在不合理的應力集中現象,而本文方法的計算結果不存在數值應力集中問題,這進一步驗證了本文方法對非均勻材料的適用性。

(a)σx曲線

(b)σy曲線

3 結 論

本文采用交錯網格技術提出了SCV-FVM。該方法圍繞節點建立控制體,將物性參數定義在單元中心,待解變量定義在單元節點上,從而能夠考慮非均質材料的空間變化。

采用SCV-FVM對宏觀非均質復合材料的應力場進行了數值模擬,通過與解析解的對比,驗證了結果的正確性。與CC-FVM相比,SCV-FVM不需要顯式處理交界面,就能避免物性突變引起的交界面應力數值波動。與HOTFGM及有限元法相比,SCV-FVM能夠避免材料空間變化引起的不合理的應力不連續現象。通過分析發現,SCV-FVM在空間任一點僅計算一個應力值,計算得到的材料交界面處的應力是一個平均值,因此難以捕捉由材料物性突變引起的垂直于牽引力方向的應力跳躍現象,但可以通過加密網格來改善計算結果。

由于SCV-FVM不需要顯式地處理復合材料交界面,因而能夠用于分析具有隨機分布微觀結構的非均質材料問題。對微觀非均質涂層的數值模擬結果表明,用SCV-FVM計算的應力場不存在數值不連續現象及應力集中問題。

[1] PINDERA M J, ABOUDI J, ARNOLD S M.Analysis of the spallation mechanism suppression in plasma-sprayed TBCs through the use of heterogeneous bond coat architectures [J].International Journal of Plasticity, 2005, 21(6): 1061-1096.

[2] ZHONG Y, BANSAL Y, PINDERA M J.Efficient reformulation of the thermal higher-order theory for FGMs with locally variable thermal conductivity [J].International Journal of Computational Engineering Science, 2004, 5(4): 795-831.

[3] BANSAL Y, PINDERA M J.Efficient reformulation of the thermoelastic higher-order theory for functionally graded materials [J].Journal of Thermal Stresses, 2003, 26(11/12): 1055-1092.

[4] CAVALCANTE M A A, MARQUES S P C, PIN- DERA M J.Parametric formulation of the finite- volume theory for functionally graded materials: part IAnalysis [J].ASME Journal of Applied Mechanics, 2007, 74(5): 935-945.

[5] CAVALCANTE M A A, MARQUES S P C, PINDERA M J.Parametric formulation of the finite- volume theory for functionally graded materials: part II Nnumerical results [J].ASME Journal of Applied Mechanics, 2007, 74(5): 946-957.

[6] CAVALCANTE M A A, MARQUES S P C, PINDERA M J.Computational aspects of the parametric finite-volume theory for functionally graded materials [J].Computational Materials Science, 2008, 44(2): 422-438.

[7] GONG Jingfeng, XUAN Lingkuan, MING Pingjian, et al.An unstructured finite-volume method for transient heat conduction analysis of multilayer functionally graded materials with mixed grids [J].Numerical Heat Transfer: Part B, 2013, 63(3): 222-247.

[8] GONG Jingfeng, MING Pingjian, XUAN Lingkuan, et al.Thermoelastic analysis of three-dimensional functionally graded rotating disks based on finite volume method [J/OL].Proceedings of the Institution of Mechanical Engineers: Part C Journal of Mechanical Engineering Science, 2013, 227(12).DOI:10.1177/0954406213489933.

[9] PENG Xulong, LI Xianfang.Thermal stress in rotating functionally graded hollow circular disks [J].Composite Structures, 2010, 92(8): 1896-1904.

[10]PENG Xulong, LI Xianfang.Elastic analysis of rotating functionally graded polar orthotropic disks [J].International Journal of Mechanical Sciences, 2012, 60(1): 84-91.

[11]CARRERA E, BRISCHETTO S, ROBALDO A.Variable kinematic model for the analysis of functionally graded material plates [J].AIAA Journal, 2008, 46(1): 194-203.

[12]GIUNTA G, BELOUETTAR S, CARRERA E.Analysis of FGM beams by means of classical and advanced theories [J].Mechanics of Advanced Materials and Structures, 2010, 17(8): 622-635.

[13]仲政, 于濤.功能梯度懸臂梁彎曲問題的解析解 [J].同濟大學學報: 自然科學版, 2006, 34(4): 443-447.ZHONG Zheng, YU Tao.Analytical bending solution of functionally graded cantilever-beam [J].Journal of Tongji University: Natural Science, 2006, 34(4): 443-447.

[14]于濤, 仲政.均布載荷作用下功能梯度懸臂梁彎曲問題的解析解 [J].固體力學學報, 2006, 27(1): 15-20.YU Tao, ZHONG Zheng.A general solution of a clamped functionally graded cantilever-beam under uniform loading [J].Acta Mechanica Solida Sinica, 2006, 27(1): 15-20.

[15]劉書田, 程耿東.復合材料應力分析的均勻化方法 [J].力學學報, 1997, 29(3): 306-313.LIU Shutian, CHENG Gengdong.Homogenization method of stress analysis of composite structures [J].Acta Mechanica Sinica, 1997, 29(3): 306-313.

(編輯 葛趙青)

StaggeredCell-VertexFiniteVolumeMethodforAnalyzingStressinHeterogeneousCompositeMaterials

XUAN Lingkuan,GONG Jingfeng,ZHANG Wenping,MING Pingjian

(College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China)

A staggered cell-vertex finite volume method (SCV-FVM) is developed for stress analysis in heterogeneous composite materials.The linear elastic equilibrium equation is discretized based on unstructured grids.The staggered grid technique is adopted to introduce the material variation into the discretization, so it is unnecessary to treat the material interfaces explicitly.SCV-FVM is taken to simulate the elastic fields of macrostructures and the results agree well with the analytical ones.Comparisons between different numerical results show that SCV-FVM is able to avoid numerical oscillation of the stress along traction direction, however it is hard for SCV-FVM to capture the stress jump vertical to the traction direction, which can be improved by refining the mesh around the material interface.The elastic performance of multi-layer composites with microstructures is discussed via SCV-FVM, and the numerical discontinuity and stress concentration due to variation of physical parameters do not appear in the predicted results, which demonstrates the feasibility of SCV-FVM for stress analysis of heterogeneous composite materials.

heterogeneous composite material; laminated material; functionally graded material; staggered grid technique; finite volume method

10.7652/xjtuxb201403022

2013-07-09。

宣領寬(1987—),男,博士生;張文平(通信作者),男,教授。

中央高?；究蒲袠I務費專項資金資助項目(HEUCF130302)。

時間: 2013-12-25

O343.7

:A

:0253-987X(2014)03-0121-07

網絡出版地址: http:∥www.cnki.net/kcms/detail/61.1069.T.20131225.1702.006.html