基于細觀有限元方法的復合材料橫向力學性能分析

劉萬雷，常新龍，張曉軍，張 磊

(第二炮兵工程大學，西安 710025)

?

基于細觀有限元方法的復合材料橫向力學性能分析

劉萬雷，常新龍，張曉軍，張 磊

(第二炮兵工程大學，西安 710025)

橫向斷裂是制約復合材料結構設計的關鍵點，傳統細觀模型因為不能充分考慮組分性能、體積分數和纖維形狀及分布情況而不能有效預測材料橫向力學性能。采用改進的隨機序列吸收算法建立具有隨機纖維分布的復合材料代表性體積單胞模型，考慮基體破壞和界面脫粘兩種失效模式和固化過程中產生的殘余應力，對模型在橫向拉、壓、剪3種載荷下的力學行為進行仿真計算。分析了不同界面強度對復合材料力學性能的影響規律。仿真結果與實驗數據對比表明：橫向模量預測誤差在7%以內，壓縮和剪切的強度誤差在8%以內，結果一致性較好，表明該模型能夠有效預測復合材料橫向力學性能。

復合材料；細觀有限元；殘余應力；界面強度；橫向力學性能

復合材料在受到多軸載荷作用時，往往首先出現垂直于纖維方向的橫向斷裂，這也成為制約復合材料結構設計的重要因素。然而，復合材料橫向破壞與其基體和界面性能密切相關，傳統的解析方法因不能充分考慮組分性能、體積分數和纖維形狀及分布情況[1]，在應用時受到很大限制。

隨著有限元技術和細觀力學的發展使得建立復雜模型分析復合材料失效機理成為可能。Gonzalez等[2]建立了具有隨機纖維分布的單向復合材料RVE單胞，采用細觀有限元方法預測了復合材料在橫向壓縮載荷作用下的應力-應變響應曲線，并分析了界面強度對橫向壓縮性能的影響規律。Vaughan等[3,4]采用實驗與數值計算相結合的方法，提出了一種生成高體積分數復合材料RVE單胞的新算法，RVE單胞中纖維之間距離的分布規律與實際細觀結構一致性較好。Lu等[5]通過細觀有限元方法研究了界面性能對三維編織復合材料單軸拉伸強度的影響，仿真結果與實驗數據吻合較好。Yang等[6]利用細觀計算力學方法研究了單向復合材料在橫向拉伸和壓縮時的力學特性，生成了纖維隨機分布的RVE單胞，并主要考慮了基體塑性變形和界面脫粘的影響，研究了界面性能對于復合材料損傷特性的影響規律。雷友鋒等[7]利用細觀有限元方法對三維RVE單胞進行了有效模量的計算，與部分實驗結果具有較好的一致性，但其采用的RVE單胞中纖維是均勻分布的，與實際情況出入較大；王兆清等[8]采用重心有限元法對單向纖維復合材料橫截面有效彈性模量進行了數值分析，與傳統有限元方法相比，該方法計算效率較高；徐強等[9]通過Python語言對ABAQUS進行了二次開發，發展了一種含界面的單向復合材料漸進損傷模型，研究了纖維排列方式及體積含量對材料橫向拉伸力學性能的影響。

然而，以往的研究中較少考慮固化過程中殘余熱應力對復合材料性能的影響，Hojo等[10], Maligno等[11]研究表明殘余熱應力的存在使得復合材料加載之前就會在界面處出現損傷，因此研究復合材料橫向力學性能時有必要將這一因素考慮在內。本工作采用細觀有限元方法研究了纖維增強復合材料橫向力學性能。建立了具有隨機纖維分布的復合材料細觀單胞模型，模型中考慮橫向加載時基體塑性破壞和界面脫粘兩種典型的損傷模式。采用摩爾-庫侖塑性模型和Ductile失效準則相結合模擬基體塑性破壞，用Cohesive單元模擬界面脫粘。模型中同時考慮了固化過程中產生的殘余應力對復合材料橫向力學性能的影響。最后用ABAQUS/Explicit程序對細觀單胞模型在橫向拉伸、壓縮、剪切3種載荷作用下的力學響應進行了仿真計算，采用細觀均勻化方法預測了復合材料橫向力學性能，并分析了不同界面強度對復合材料橫向力學性能的影響規律。

1 復合材料細觀力學模型

復合材料橫向力學性能主要受基體塑性變形和界面脫粘兩方面的影響，本工作將纖維視為橫觀各向同性的彈性材料，將基體看成是各向同性的彈塑性材料,組分性能如表1所示。用摩爾-庫侖塑性模型[12]描述基體的非線性行為，該模型中認為材料屈服發生在剪應力作用的平面，當其達到一個臨界值時，發生塑性變形，這個臨界值依賴于平面的法向應力，用下式表示：

τ=c-σtanφ

(1)

式中c和φ分別代表黏聚力和摩擦角，這兩個常數控制著材料的塑性行為。黏聚力c代表純剪切情況下材料的屈服應力，摩擦角φ代表靜水壓力對材料塑性的影響。φ=0時，摩爾-庫侖模型退化為Tresca模型，當φ=90°時，演化為Rankine模型。這兩個參數都可以通過拉伸和壓縮強度獲得：

(2)

固體斷裂面遵循摩爾-庫侖準則，并且單軸壓縮時的斷裂角度α與摩擦角φ之間的關系如下：

α=45°+φ/2

(3)

典型的環氧基體材料中有50°<α<60°，這樣10<φ<30°。一旦φ固定，則對應的c值也就確定了，由文獻[13]知基體拉伸和壓縮強度分別為80MPa和120MPa，故可以確定摩擦角φ=12°，c=50MPa。假設c和φ為常數，獨立于累積塑性應變，摩爾庫侖模型的屈服面以最大和最小主應力表示如下：

F(σΙ,σIII)=(σΙ-σIII)+(σΙ+σIII)sinφ-2ccosφ=0

(4)

(5)

表1 玻璃纖維/環氧樹脂復合材料組分性能[13]

Table 1 Constituent properties of glass fiber/epoxy resin matrix composite[13]

MaterialE/GPaG/GPaνXt/MPaXc/MPaS/MPaFiber80330.2021501450—Matrix3.351.480.358012075

復合材料橫向斷裂常常是由于纖維/基體界面脫粘和基體微裂紋等細觀損傷積累引起的，這些細觀損傷與許多因素有關，其中影響最大的因素是不同組分之間的黏合力。本工作采用一種基于物理機制的零厚度Cohesive界面單元及Traction-Separation本構模擬纖維/基體界面損傷，并認為界面單元的拉伸和剪切強度相等[14]。二維界面單元中存在兩個應力分量(tn，ts)和兩個位移分量(δn，δs)，初始時，界面應力和位移為線性關系，當界面應力達到損傷初始準則后，損傷萌生，界面的應力和位移關系不再保持線性變化，采用二次應力準則作為初始失效準則：

(6)

復合材料橫向加載時，法線方向存在受拉和受壓兩種狀態，當受壓時，無論是否存在損傷，其壓縮剛度保持為初始剛度；當受拉時，其剛度根據損傷程度進行折減(如圖1所示，其中K代表界面剛度)：

(7)

式中：δn代表法向位移，δ0代表初始損傷位移，δf代表失效位移，δmax代表積分點處計算得到的最大有效位移，d為損傷參數。

采用斷裂能準則判斷界面最終失效：

(8)

圖1 Cohesive單元牽引拉伸失效規律

2 細觀有限元模型建立

2.1 RVE單胞生成

本工作采用二維方形RVE單胞計算復合材料橫向力學性能。RVE單胞是由隨機排列的圓形纖維嵌入在基體中組成的，仿真中一個重要的問題就是確定RVE單胞的尺寸，一方面需要盡量縮小RVE單胞的尺寸減小計算量，同時又需要保證RVE單胞中包含足夠的信息能獨立代表復合材料性能，由文獻[15]可知包含30根纖維的RVE單胞可以滿足此條件。

采用改進的隨機順序吸收算法[16](算法如圖2所示)生成63μm×63μm的方形RVE單胞，單胞中包含30根隨機排列的纖維，纖維直徑為5μm，體積分數為60%，誤差控制在1%以內。包含30根纖維的RVE單胞如圖3所示。

圖2 改進隨機順序吸收算法流程圖

圖3 RVE單胞周期性邊界條件施加示意圖

2.2 周期性邊界條件的施加

含周期性單胞的材料在外載荷作用下，其應力和應變場呈現出連續性和周期性，因此需要對RVE單胞施加周期性邊界條件，確保得到合理的細觀應力、應變分布。對圖3所示的二維單胞施加周期性邊界條件[17]的方程如下：

(9)

(10)

對于四個角節點OABC施加如下約束：O點全約束，A，C節點可以根據施加載荷的需要確定施加的約束；B點約束方程為：

uB=uA+uC

(11)

將上述方程用腳本語言python編寫成子程序供ABAQUS軟件調用，完成對RVE單胞周期性邊界條件的施加。

2.3 模型的網格劃分及界面單元的建立

采用四邊形單元CPE4R單元和三角形單元CPE3對RVE單胞進行網格劃分，在基體和纖維邊界添加二維Cohesive單元COH2D4，添加Cohesive單元的程序流程如圖4所示。

圖4 Cohesive單元添加算法

3 橫向載荷作用下的力學性能分析

復合材料固化降溫過程中由于纖維和基體的熱膨脹系數不同，會在纖維/基體界面處產生殘余應力，為了讓仿真結果更接近實際情況，在計算過程中需要考慮殘余應力的影響。纖維和基體的熱膨脹系數分別取為15×10-6/K，45×10-6/K，固化溫度從175℃降到室溫25℃[17]。

通過對RVE單胞施加不同的初始邊界條件，可以實現橫向拉伸、壓縮和剪切加載。計算結束后提取單胞應力、應變，然后采用均勻化方法[1]得到宏觀尺度的應力應變響應，計算公式如下：

(12)

通過式(12)可以確定不同載荷下復合材料宏觀的應力應變關系，根據應力應變關系可以求出復合材料的橫向模量、強度等力學性能參數。

采用ABAQUS/Explicit顯示有限元分析軟件(可以有效避免數值計算的收斂性問題)對所建立的單胞模型進行拉、壓、剪計算(界面強度取值為50MPa，界面剛度為108MPa/m)，并將計算結果與文獻中實驗結果[13]對比，如表2所示。可以看出，仿真計算結果與實驗結果較為一致，表明本工作所建的模型適合于復合材料橫向力學性能預測。

表2 復合材料橫向力學性能預測結果與實驗結果[13]對比

Table 2 Comparison between simulation results and experimental results[13]

ResultEt/GPaEc/GPaG/GPaXt/MPaXc/MPaS/MPaSimulation16.5317.175.4638.92122.5076.16Experiment17.7017.705.8335.00114.0072Error6.61%3.00%6.78%11.20%7.46%5.78%

在復合材料加工制作過程中，纖維/基體界面強度會受到多種因素的影響，且可以通過不同實驗手段提高界面強度，因此有必要研究不同界面強度對復合材料橫向力學性能的影響。本工作對25，50，80MPa(約等于基體拉伸強度)和100MPa四種界面強度下的復合材料橫向力學性能變化規律進行了對比分析。

3.1 拉伸載荷

在圖3中點A處施加X方向位移載荷，固定點O，限制點C的轉動自由度和X方向位移自由度，即可完成拉伸位移的施加。圖5中給出了3種不同界面強度下，橫向拉伸應力與應變之間的關系，“w/”代表考慮殘余應力的影響，“w/o”代表不考慮殘余應力影響的情況。從圖5中可以看出，相比不考慮殘余應力的情況而言，存在殘余應力時復合材料拉伸強度略有提高，這是由于基體熱膨脹系數相對纖維較大，固化過程中會在界面處產生壓應力，抵消了一部分拉伸載荷的影響；隨著界面強度的提高，復合材料拉伸強度隨之提高，但當界面強度接近基體強度后，界面強度的提高對復合材料拉伸強度的影響隨之減弱。

圖5 不同界面強度下橫向拉伸應力-應變關系

圖6中給出了界面強度為50MPa時復合材料的損傷演化過程。由圖中可以看出當界面強度低于基體強度時，纖維/基體之間的界面先于基體發生破壞，在基體達到拉伸強度之前便發生了界面脫粘，原有界面變成了自由表面，內部應變能得到釋放。隨著界面損傷的積累，引起基體破壞，并最終導致材料整體失效，最終斷裂形態與圖6(c)的實驗結果相似。當復合材料界面強度高于基體強度(即界面強度為100MPa)時損傷演化過程為：裂紋首先在基體中產生，隨著損傷積累逐漸擴展到界面處，并最終引起斷裂破壞。

圖6 橫向拉伸斷裂仿真與實驗結果[18]對比 (a)初始損傷；(b)最終斷裂；(c)實驗結果

3.2 壓縮載荷

圖7中給出了不同界面強度下復合材料橫向壓縮的應力-應變曲線，可以看出，界面強度相同時，由于界面處殘余壓應力的影響，使得復合材料橫向壓縮強度略有降低。

圖7 不同界面強度下橫向壓縮斷裂計算結果

從圖7中可以看出復合材料橫向壓縮失效之前都表現出明顯的非線性行為，這是由于在壓縮載荷作用下，首先發生界面脫粘，隨后出現基體裂紋，最終導致材料發生整體破壞。橫向拉伸載荷作用下引起界面脫粘的是基體與纖維間的法向應力,而壓縮載荷作用下界面脫粘是切向應力引起的。隨著界面強度的提高，復合材料壓縮斷裂強度也隨之增高，當界面強度從25MPa上升到50MPa時，復合材料壓縮斷裂強度提高了約22%，而隨著界面強度的進一步提高，斷裂強度提高則不明顯，這是由于當界面強度增大到一定程度時，壓縮時首先出現基體破壞，隨著損傷積累才出現界面破壞，材料的性能主要由基體力學性能決定。標準界面下復合材料斷裂角度如圖8所示，與實驗觀察到的角度較為一致[19]。

3.3 剪切載荷

不同界面強度情況下，復合材料剪切性能的變化規律與壓縮變化規律相似。由圖9中可以看出殘余應力對復合材料剪切性能影響不大。復合材料在剪切載荷作用下,法向應力與切向應力的綜合作用導致部分界面早于基體發生失效。隨著載荷的繼續增加,界面損傷不再擴展,基體開始進入損傷階段。當界面強度明顯低于基體強度時,較弱的界面強度脫粘導致纖維與基體形成過多的既有表面而造成缺少足夠的相互約束以維持整個單胞的剪切變形。因此,在線性段過后,隨著加載的繼續,平均應力有段較長的平緩上升過程。

圖8 橫向壓縮斷裂仿真與實驗結果[19]對比 (a)初始損傷；(b)最終斷裂；(c)實驗結果

圖9 不同界面強度下橫向壓縮斷裂計算結果

4 結論

(1)建立了具有隨機纖維分布的周期性代表單胞模型，通過摩爾-庫侖塑性模型與Ductile damage失效準則結合模擬基體損傷演化過程，并將固化過程中的殘余應力考慮在內，實現了復合材料橫向模量、強度等力學參數的預測。

(2)將仿真結果與實驗數據對比表明：復合材料橫向模量誤差在7%以內，橫向壓縮和剪切強度誤差在8%以內，拉伸強度為11.2%，精度滿足工程應用需求。

(3)固化過程中產生的殘余應力導致復合材料橫向拉伸強度略有提高，壓縮強度略有降低，對剪切強度的影響不大。

(4)當界面強度低于基體強度時，復合材料橫向強度由界面強度決定，隨著界面強度的提高，復合材料橫向強度有明顯提高；當界面強度高于基體強度時，復合材料橫向強度主要由基體性能決定，界面強度的變化影響不大。

[1] 杜善義,王彪. 復合材料細觀力學[M]. 北京: 科學出版社, 1998.

DU S Y，WANG B. Micromechanics of Composite [M].Beijing：Science Press China，1988.

[2] GONZALEZ C, LLORCA J. Mechanical behavior of unidirectional fiber-reinforced polymers under transverse compression: microscopic mechanisms and modeling[J]. Composites Science and Technology, 2007,67(13):2795-2806.

[3] VAUGHAN T J, McCARTHY C T. A combined experimental-numerical approach for generating statistically equivalent fibre distributions for high strength laminated composite materials [J]. Composites Science and Technology,2010,70(2): 291-297.

[4] VAUGHAN T J, McCARTHY C T. Micromechanical modeling of the transverse damage behaviour in fibre reinforced composites [J]. Composites Science and Technology，2011,71(3):388-396.

[5] LU Z X, WANG C Y, XIA B, et al. Effect of interfacial properties on the uniaxial tensile behavior of three-dimensional braided composites [J]. Computational Materials Science，2013, 79: 547-557.

[6] YANG L, YAN Y, LIU Y J, et al. Microscopic failure mechanisms of fiber-reinforced polymer composites under transverse tension and compression[J]. Composites Science and Technology，2012,72(15):1818-1825.

[7] 雷友峰, 魏德明, 高德平. 細觀力學有限元法預測復合材料宏觀有效性能模量[J]. 燃氣渦輪實驗與研究, 2003, 16(3):11-15.

LEI Y F, WEI D M, GAO D P. Predicting macroscopic effective elastic moduli of composites by micromechanics FEM [J] . Gas Turbine Experiment and Research, 2003, 16(3): 11-15.

[8] 王兆清, 張景濤, 李淑萍. 計算復合材料有效彈性模量的重心有限元方法[J].復合材料學報，2007, 24(6): 173-179.

WANG Z Q, ZHANG J T, LI S P. Barycentric finite element method for predicting the effective elastic moduli of composite materials [J]. Acta Materiae Compositae Sinica, 2007, 24(6): 173-179.

[9] 徐強，盧子興. 纖維排列方式對單向復合材料橫向拉伸強度的影響[C]//中國力學大會2011暨錢學森誕辰100周年紀念大會論文集，哈爾濱: 中國力學學會, 2011: 1-10.

XU Q, LU Z X. Effects of fiber distributions on the transverse properties of unidirectional composites[C]//CCTAM2011， Harbin: CSTAM, 2011: 1-10.

[10] HOJO M, MIZUNO M, HOBBIEBRUNKEN T, et al. Effect of fiber array irregularities on microscopic interfacial normal stress states of transversely loaded UD-CFRP from viewpoint of failure initiation[J]. Compos Sci Technol, 2009, 69(11-12):1726-1734.

[11] MALIGNO A R, WARRIOR N A, LONG A C. Finite element investigations on the microstructure of fibre-reinforced composites[J]. Express Polym Lett, 2008, 2(9):665-676.

[12] MENETREY P H, WILLAM K J. Triaxial Failure criterion for concrete and its generalization[J]. ACI Structural Journal, 1995,92(6):311-318.

[13] SODEN P D, HINTON M J, KADDOUR A S. Lamina properties, lay-up configurations and loading conditions for a range of fiber-reinforced composite laminates [J]. Composites Science and Technology，1998,58(7):1011-1022.

[14] PITKETHLY M J，FAVRE J P，GAUR U，et al. A round robin programme on interfacial test methods[J]. Composites Science and Technology，1993,48(1-4):205-214.

[15] SEGURADO J, LLORCA J. A numerical approximation to the elastic properties of sphere-reinforced composites [J]. J Mech Phys Solids,2002,50(10):2107-2121.

[16] XIA Z H, ZHANG Y F, ELLYIN F. A unified periodical boundary conditions for representative volume elements of composites and applications [J]. International Journal of Solids and Structures, 2003, 40(8): 1907-1921.

[17] O’HIGGINS R M, McCARTHY C T, McCARTHY M A. Identification of damage and plasticity parameters for continuum damage mechanics modelling of carbon and glass fibre-reinforced composite materials[J]. Strain，2011,47(1):105-115.

[18] WANG X Q, ZHANG J F, WANG Z Q. Effects of interphase properties in unidirectional fiber reinforced composite materials[J]. Materials & Design， 2011, 32(6): 3486-3492.

[19] BUDIANSKY B, FLECK N A. Compressive failure of fiber composites[J]. Mech Phys Solids,1993,41(1):183-211.

--------------------●

Analysis of Composite Transverse Mechanical Properties Based on Micromechanical Finite Element Method

LIU Wan-lei,CHANG Xin-long,ZHANG Xiao-jun,ZHANG Lei

(The Second Artillery Engineering University, Xi’an 710025, China)

Transverse fracture often occurs early in the loading history and is one of the key issues limiting the composite structural design. However, the mechanical behavior under transverse loading can not be represented by traditional micromechanical model which does not consider the influence of constituent properties, fiber volume fraction and fiber distribution. A new computational micromechanics finite element method, which the microstructure was idealized as a random dispersion of parallel fibers embedded in the polymeric matrix using improved random sequence absorption algorithm, was presented. The plasticity of matrix and interface decohesion of the composite were included in this model and the residual stress caused by the cooling of composite after the curing process was also taken into consideration. The transverse tension, compression and shear of composite were analyzed by the micromechanical finite element method. Compared with the experiment results, the prediction errors of transverse module were less than 7%, and the transverse compression strength and shear strength were less than 8%. The results demonstrate that the method proposed here can be used to predict the composite transverse behavior.

composite;micromechanical finite element;residual stress;interface strength;transverse mechanical property

10.11868/j.issn.1001-4381.2016.11.018

TB33

A

1001-4381(2016)11-0107-07

國家自然科學基金項目(11302249)

2015-05-25；

2016-07-27

常新龍(1965-)，男，教授，博士生導師，博士，主要從事失效物理與可靠性方面的研究工作，聯系地址：陜西省西安灞橋區洪慶鎮同心路2號3602分隊(710025)，E-mail:xinlongch@sina.com