隨機多尺度短切碳纖維復合結構模型中RVE 尺寸效應和方向模量的均一化響應

邱伊健, 鄭萍, 程香平, 郭玉松

(1. 江西省科學院 應用物理研究所，江西 南昌 330001; 2. 南昌師范學院，江西 南昌 330032)

0 引言

手工鋪設、長絲纏繞、蒸壓和真空袋處理技 術[1]，通常與連續(xù)纖維增強熱固性復合材料相關，它通常應用于短期或一次性使用的產品，并且對產品質量和性能的要求較高。然而，當需要大量的此類產品時，這些技術是有限制的。對于當前許多應用來說，強度或剛度峰值的提高并不是這種復合材料的主要需求。而是在滿足特定強度和剛度需求情況下，人們往往希望批量生產這種產品，使得用于連續(xù)纖維的加工技術在機器時間和人工時間方面顯得極其冗長和昂貴。通常聚合物基復合材料中的基體不具有承受很大載荷的能力，其主要作用之一是將負載從一根纖維轉移到另一根纖維。它影響復合材料的抗壓、抗剪強度，也有助于提高復合材料的抗斷裂性能和能量吸收能力。

短纖維增強樹脂基復合材料(SFERC)[2-4]由于其良好的經濟性和相對可靠的機械性能已經被廣泛關注。近年來，由于復合材料制造工藝的進步，短切纖維增強彈性復合材料在中低端工程結構中的力學性能能夠滿足大多數使用需求，在某種程度上能夠替代連續(xù)纖維增強復合材料(CFRC)[5-6]。由于短切纖維在基體中復合狀態(tài)的無序性，并呈現卷曲 狀[7]，使得這種材料在細觀層面存在極其復雜的結構。如何將SFERC 中復雜的細觀結構建立準確的均一化方法，是一項具有挑戰(zhàn)性的工作。

已有大量研究工作提出SFERC 的均質化方 法[8-13]，傳統(tǒng)的方法是依靠實驗技術來表征材料的有效力學性能。這些性能用適當的現象學本構模型來表達時，可以作為傳統(tǒng)工程材料力學分析的基礎。然而，現代方法論傾向于用均一化方案和微觀力學模型結合來取代這種現象學的方法，它是通過在宏觀和微觀兩個尺度下同時求解，以降低求解非線性和時變響應的兩個尺度相關計算代價的解耦方法。復合材料均一化理論最早由Hill[9]和Ogden[10]提出。但是由于復合材料均一化場是復雜結構的非線性場，在受到外部載荷時場量的非線性變化限制了其實際應用。隨后Castaneda[11-12]提出了2 階均一化方案并由Kouznetsova 等[13]成熟應用。由于平均場理論組分均一化響應的限制，Chen 等[14]采用全場法得到了復合結構的局部應力和應變場響應。全場法比平均場法具有更高的微觀場分辨率。最近Caylak 等[15]研究了結合平均場和全場均一化方法的框架下建立多態(tài)不確定模型，用模糊隨機變量發(fā)展多態(tài)不確定平均場和全場方法，在纖維體積分數不確定幾何形狀和材料參數不確定的情況下，確定纖維增強復合材料的隨機有效性能。

除此之外，人們還將各種其他領域的研究方法引入復合材料的細觀力學研究[16-20]，例如機器學習和人工智能記憶法[18,21-24]，即便它們目前與傳統(tǒng)的直接數值計算法相比并沒有帶來切實的計算效益，但是可以作為未來本研究領域的一個新方向。

與直立纖維相比，卷曲纖維的隨機離散化計算很難采用直接的有限元分析(FEA)方法來處理。直立纖維可以直接離散到代表體積元(RVE)結構中，而卷曲纖維由于其模量在各部位是不一致的，需要先將卷曲纖維的模量進行局部均一化，然后對均一化完成后的局部樣本進行離散。傳統(tǒng)均一化方法通常采用一個尺度，即RVE 尺度，直接把短切纖維離散到RVE 結構中，然后對RVE 進行均一化。這種離散算法最典型的是2 階或4 階隨機纖維取向張量的統(tǒng)計離散法。而本文要解決的問題在于卷曲纖維模量是隨著局部坐標位置不斷變化的，它不能像直立纖維那樣直接通過3 個方向的正向模量和3 個方向的剪切模量來表征。

本文對于卷曲纖維的方向性模量需要在纖維局部進行模量路徑上微分轉置處理，將偏軸模量轉置到正軸模量的分量上，然后對卷曲纖維進行離散。與傳統(tǒng)算法不同的是，本文的創(chuàng)新點在于在有限元模型的構建中采用雙尺度，即單胞(UC)和RVE 尺度同時進行分析。首先將卷曲纖維的模量均一化到UC 尺度，獲取UC 結構的均一化模量，然后采用本文提出的算法將UC 尺度模型離散到RVE 結構中，進行二次均一化。為了使得雙尺度均一化的實現，本文采用一種新型的隨機場，稱為空間方向敏感(HE)算法[25]，將空間隨機場分為兩部分，其中一部分為確定場量，另一部分為隨機場量。

1 多尺度均一化實驗

1.1 卷曲短切纖維彈性場微分的坐標變換

為描述單根卷曲短切纖維結構模量，把短切纖維假設為 Hsiao-Daniel 正弦波形式的卷曲狀 態(tài)[8]。其卷曲周期為ω，振幅為A，如圖1 所示。圖1 中，θ為軸與Ox軸的旋轉角度，軸表示局部纖維長度方向，軸、軸分別表示與長度方向垂直的纖維截面方向。

圖1 卷曲短纖維長度微分的坐標變換Fig. 1 Coordinate transformation of length differential of CCF

纖維的波動性方程在Oxyz坐標系(后文簡稱主坐標系)下可以定義為

式中：u為纖維波動的縱坐標值。

在卷曲纖維的任意一端截取一段長度為dx的波動段，并且在平行于這段纖維主方向的切線上建立轉置纖維局部坐標系(后文簡稱轉置坐標系)。其中坐標軸平行于切線方向，坐標軸垂直于切線方向，而坐標軸與Oyz軸保持一致。因此主坐標系與轉置坐標系下的應力轉換可以表示為

式中：σxx、σyy、σzz分別為主坐標系3 個方向的正應力；τyz、τxz、τxy分別為主坐標系3 個方向的剪應力；Tij為轉置矩陣(i、j為矩陣分量的下角標)；分別為轉置標系的3 個正應力；分別為轉置坐標系的3 個方向剪應力。

這里的Tij-1為矩陣Tij的逆矩陣，并且

同理，相應的應變轉換方程可以表示為

式中：εx、εy、εz和分別為主坐標系和轉置坐標系3 個方向的正應變；γx、γy、γz和分別為主坐標系和轉置坐標系下 3 個方向的剪應變。

通過在式(4)兩邊分別左乘對角矩陣Rij，

將左邊正應變和剪應變化簡成一致形式，右邊左乘一個單位矩陣E=Rij-1Rij，將右邊也簡化成一致形式，得到

因此主坐標系下的應力應變關系可以表示為

式中：S和Sij分別為主坐標系和轉置坐標系下的柔量矩陣，二者存在如下轉換關系：

在轉置坐標系下，纖維可以看成主軸方向平行于Ofxf軸的橫觀各向異性材料，對于橫觀各向異性的柔量矩陣采用Hill 矩陣可以寫成如下形式：

1.2 短切纖維的隨機場

1.2.1 Digimat 原生算法

Digimat 軟件采用2 階張量取向填充算法，它是一種以概率隨機取向形式統(tǒng)計的一種纖維離散方法，這種算法一般用于直立纖維，并且纖維結構中不附帶基體相的情況。對于短纖維復合材料，定義微觀結構的最重要的隨機變量是纖維體積分數vf、纖維取向和纖維長度lf。由圖2(a)可知，可以根據兩個纖維角φf和θf來定義單根纖維取向。

圖2 纖維取向Fig. 2 Fiber orientation

圖2 中，Ω為復合材料的結構域，P(q)為 第q根纖維的方向矢量，x1′、x2′、x3′分別為纖維方向矢量分解的3 個正交向量，e1′、e2′、e3′分別為纖維方向矢量分解的3 個正交基向量。

式中：P1f、P2f、P3f分別為纖維取向的3 個正交方向矢量分量；α1f、α2f、α3f分別為纖維與總體坐標系所呈的3 個Euler 角。

需要說明的是，使用歐拉角或方向單位向量時，3 個變量中只有2 個是獨立的。所有4 個獨立的隨機變量Pf、φf、θf和lf(或者分別是歐拉角或方向矢量分量中的2 個)需要提供適當的概率分布，來表示纖維體積分數、取向和長度的不確定性。

對于2 階取向張量填充算法(2 階纖維張量取向描述法)來填充聚集無序短纖維的取向問題，因為它們可以與纖維取向的變換矩陣和概率分布直接相關。這里簡要介紹一下張量取向aij的定義的推導，首先考慮包含根據圖2(b)中有限數量的纖維的任意參考體積Ω，其中包含了有限數目nf根短纖維，如圖2(b)所示。任意纖維號為q的取向是根據其單獨的纖維取向向量P(q)定義的，其向量組成P(q)根據式(9)來定義。因此，對于參考體積為Ω，纖維總數為nf的2 階纖維張量取向由式(10)定義：

式(10)用于域Ω內的所有纖維。如果Ω內的纖維數目接近無窮大，即nf→∞，則用連續(xù)定義 式(11)來代替2 階纖維張量取向的離散定義(見 式(10))：

式中：f(Pf)為纖維取向Pf出現的概率密度分布函數?？紤]到概率分布中期望值的定義：

根據式(12)，2 階纖維張量取向的分量aij定義了纖維(歸一化)方向矢量分量乘積的期望值，因此參考式(9)，即乘積的期望值，分別為離散和連續(xù)纖維取向分布的歐拉角。由于式(10)和式(11)中乘積的交換性，2 階張量取向是對稱的，即aij=aji。因此，對于完全排列的纖維，全部朝向x1方向，2 階纖維張量取向的形式為

而

則表示纖維取向完全無序。

1.2.2 HE 算法

在獲取短切纖維卷曲彈性應力應變場變換之后，還需要考慮卷曲纖維在基體相中的離散分布。因此需要在體積元中建立合適的纖維隨機場來描述離散纖維強化基體復合結構的彈性本構。本文通過RVE 結構來研究離散短纖維在基體中的離散行為，其坐標系為Orxryrzr，而單根短切纖維通過UC結構來表征，其坐標系用Ouxuyuzu表示，它內部包含了單根短切纖維，與附著于纖維壁呈一定厚度的界面層，以及部分UC 結構中包裹纖維及其界面層的基體材料(見圖3)。

圖3 單胞纖維在RVE 中的離散角度Fig. 3 Discrete angle of single cell fiber in RVE

圖3 中，α為UC 坐標系旋轉矢量與主坐標系Ox軸的夾角，β為旋轉矢量在主坐標系Oyz平面上的投影與Oy軸的夾角。

本文中介紹一種新型隨機場模型，這種模型稱為HE 算法，本文將該算法應用于針對卷曲短纖維隨機場的離散化模型。HE 算法首先將一個空間隨機場H(x，y，z，δ)分解為一個確定域和一個隨機域：

式中：為空間場平均值(它又可以分為α和β角度方向的空間場(見圖3)；δ為初始隨機性；ξ為常變量；λi和φi分別為解析協(xié)方差內核的本征值和本征函數；n為HE 的截斷階。通過求解下面的Fredholm 積分方程[26-30]解出λi和φi的值：

式中：x1(x1,y1,z1)和x2(x2,y2,z2)分別為Ofxfyfzf坐標系內隨機選擇的兩個點；VRVE為RVE 的體積；C(x1,x2)為假設協(xié)方差函數的指數形式，

UC 纖維的離散通過兩個旋轉張量角度α和β通過旋轉變換將場量由主坐標系旋轉至離散UC 坐標系，如圖3 所示。本文采用Galerkin 有限元法[31]對特征函數進行離散。

1.3 短切卷曲纖維強化樹脂基三維隨機排列的概率模型

在卷曲纖維外層一定厚度范圍內存在一層界面相，其屬性是樹脂基體材料，并且界面相承擔了更密集的應力傳遞，它比遠離纖維的基體形變量更大，因此在直徑為dCCF的纖維外部增加一層外徑為di的界面層?？紤]含基體相的獨立短切纖維復合結構UC，如圖4所示。圖4 中，LCCF為CCF的長度。

圖4 包含單根波動碳纖維的單胞結構 Fig. 4 Unit cell with single wavy CCF

CCF 的隨機卷曲度[6]可以定義為

式中：x為微尺度域上的位置向量；LCCF為CCF 的長度。CCF UC 的無限小片段看成是離軸橫斷各向同性層(見圖4)。通過把Chamis 的簡化微觀力學模型[24]應用于短切纖維隨機場，可以得到5 個獨立的彈性屬性E11、E22、G12、G23、V12。

對UC 的E1和υ12使用混合等效應變準則時[30]，可以得到

式中：VCCF(x,δ)為碳纖維的體積；ECCF1為碳纖維長度方向模量；Vm為基體體積；Em為基體彈性模量；Vi(x,δ)為隨機界面的體積；Ei(x,δ)為隨機界面彈性模量；v12(x,δ)為基體泊松比；vCCF12為碳纖維泊松比；νi為隨機界面的泊松比。

界面材料的體積分數可以表示為

式中：di包含了纖維和界面層；dCCF為不包含界面層的單根纖維外經。橫向有效模量可以通過復合材料基礎力學推導而來[32-37]：

同理，隨機剪切模量(G23(x，δ)，G12(x，δ)=G13(x，δ))也可以通過將式(19)中的Em替換為Gm，ECCF2替換為GCCF12得到。隨后的計算可以通過第2 節(jié)所述短切纖維卷曲度坐標變換公式計算。同時，短切碳纖維增強相的長徑比影響可以間接通過局部體積分數和纖維卷曲度來考慮。

本文提出的構建空間隨機CCF 排列模型的方法可以通過兩個取向角來表征(αCCF和βCCF)，分別表示UC 在RVE 中的α和β的離散角。這里將CCF 看成是隨空間變化的隨機場，并且通過HE 方法將其表示為

式中： ＜·＞表示平均值。x軸、y軸、z軸方向的角度均值和相關長度是決定CCF 隨機取向的主要統(tǒng)計參數因素。在實際的聚合物基體中，每根CCF不可能完全對齊，如圖5 所示，圖中〈θCCF〉、θird分別為平均波動和隨機波動。因此，取向的CCF復合材料具有空間變化的各向異性。

圖5 RVE 空間變化的隨機排列角度Fig. 5 Random arrangement angle of RVE space variation

應力張量和應變張量在RVE 坐標系Orxryrzr和UC 坐標系Ouxuyuzu之間來回變換，如式(21)所示：

隨機轉動張量根據式(20)中的角度可以表 示為

式中：Rkl表示隨機旋轉張量變換。

1.4 平均應力，應變及其后處理方法

在三維RVE 結構中，每個離散的UC 可以看成是包含在其中的一個子域，如果要從子域中提取應力和應變，則在對結果進行后處理時不能使用關鍵點自由度(Kdofs)，因為 Kdofs 是針對完整的RVE，而不是子域。應力和應變必須以某種方式平均，涉及到應力和應變在子域上的積分。由于將應力和應變的積分過程作為計算機代碼來實現對于大多數用戶是一項要求很高的工作，下面的簡化可以保持數學上的嚴謹性。高斯定理(散度定理)給出了下列積分公式：

式中：[fx,fy,fz]為三維空間中的連續(xù)可微向量場；?Ω為三維空間中感興趣域的邊界；[nx,ny,nz]為?Ω的單位向外法線；S為?Ω的表面積。特別地，如果在三維空間中引入[fx,fy,fz]=[gxhx,gyhy,gzhz]，這里g和h是兩個單獨的連續(xù)可微向量場，類似于f，高斯定理給出了下列分段積分：

利用這些數學公式，可以導出三維情況下的平均應變和應力，為此使用高斯定理得出以下平均應變的表達式：

1.5 數值模擬算法實現和材料參數

1.5.1 算法實現

本文提出的算法實現過程如圖6 所示。首先對DIGIMAT 腳本進行RVE 參數，材料模型等相關數據輸入，包括HE 算法、RVE 設計尺寸、邊界條件、纖維目標體積分數、纖維半徑、縱橫比、纖維填充數目、纖維最小間距、界面相直徑與纖維直徑之比(di/dCCF)等參數。然后程序根據HE 算法的周期邊界條件對RVE 進行填充，得到相應的RVE FE 結構模型。在MATLAB 軟件腳本對RVE FE 結構模型進行Halpin-Tsai 方程坐標變換，將宏觀RVE 尺度隨機場量離散到細觀的UC 尺度下，并得到以UC 為最小單位的離散化的空間隨機分布模量場，將離散化數據在Abaqus CAE 軟件中基于空間隨機分布模量場進行相應的有限元計算，從而得到相應的數值解。

圖6 HE 算法的流程圖Fig. 6 Flowchart of space direction-sensitive hoteling expansion algorithm

1.5.2 材料屬性和實驗材料制備

表1 為RVE 中3 種材料結構的材料參數和統(tǒng)計參數。需要說明的是，基體材料和界面材料是同種物質，但是界面材料由于緊鄰CCF 相，受到其強化作用的影響，其材料力學行為與基體材料在微力學模型中存在一些區(qū)別，但是這種細微的調整在微力學模型中是必要的，也是可以接受的。雖然CCF 填充樹脂基復合材料復合強度隨著CF 含量的增加遞增，但是短纖維強化與長纖維不同，短纖維增強樹脂基體隨著纖維含量的增加，復合結構脆性斷裂的風險呈指數遞增。因此為兼顧韌性和強度，短纖維增強環(huán)氧樹脂基復合材料的纖維體積分數應該低于3%。本文選取最大體積分數為1.24%。

基體材料考慮塑性變形，采用J2 塑性模型[30]，并且考慮的是指數和線性強化準則，CCF 只考慮彈性效應，相應參數如表1 所示。

表1 RVE 結構各相材料屬性Table 1 Material properties of each phase in RVE structure

RVE 結構微力學實驗模型如圖7 所示。選取平均長度為1 mm 的短切碳纖維作為增強相，均勻分散于環(huán)氧樹脂基體中，并用聚四氟乙烯模具將其定型固化為正方體形狀，如圖7(a)所示。力學實驗選取尺寸為20 mm×20 mm×20 mm 的立方體，填充纖維體積分數分別選取0.2%、0.41%、0.61%、0.82%、1.03%和1.24% 6 組進行實驗。其中，圖7(b)所示為纖維體積分數為1.24%的樣本中短切碳纖維在基體內隨機分布的細觀顯微圖像。然后將固化的立方體試樣再次置于模具中，加入無增強相的環(huán)氧樹脂固化后，成為拉伸力學試樣，如圖7(c)所示。

圖7 復合結構實驗樣本及其細觀顯微圖像Fig. 7 Composite structure test samples and their mesoscopic images

2 實驗結果與討論

2.1 實驗樣本的纖維長度及2 階張量取向角的概率分布測定

為確保數值RVE 模型參數與實驗樣本的一致性，采用樣品基體的化學焚化法提取纖維，通過圖像處理系統(tǒng)對它們進行長度掃描，統(tǒng)計圖7 所示的樣本中短纖維的長度概率分布情況。圖8(a)所示為測試樣本纖維長度的統(tǒng)計數據。圖8 中，f(LCCF)為短纖維長度的概率密度分布函數，而F(LCCF)為纖維長度歸一化概率分布函數，它們之間的關系為f(LCCF)=dF(LCCF)/dx。

測試樣本中，短纖維長度分布為0~4.8 mm，纖維長度分布概率整體呈線性遞增，同時存在一些概率跳躍分布點，這些跳躍點在較短纖維區(qū)間 (< 1 mm)分布較明顯。

同時還確定了纖維取向分布，它需要有關微觀結構幾何形狀的完整的三維信息。如果纖維和基體的密度存在較大不同，例如對于大多數玻璃纖維，碳纖維增強的聚合物，可以通過X 射線計算機斷層掃描以無損方式研究微觀結構。為此，在大量連續(xù)的旋轉角度下，通過X 射線源對樣本進行射線照相。對拍攝的圖像進行計算，可以重建微觀結構；可獲得基于體素的顯微組織三維圖像結果并跟蹤每根纖維，從而提供必要的纖維方向數據，以根據 2 階張量取向法確定纖維取向張量。通過X 射線計算機斷層掃描獲取纖維取向矢量p，并得出樣本中每根纖維的纖維角度。圖8(b)、圖8(c)給出了實驗樣本中根據圖3 所示角度確定的纖維取向分布曲線，該樣本取向纖維以仰角α和投影角β確定。其中投影角α和β的概率分布都是呈高斯正態(tài)分布，β的概率分布統(tǒng)計的周期為π，即[-90°，90°]，而α的概率分布統(tǒng)計的周期為π/2，即[-45°，45°]，從而導致β和α分別呈現單峰和雙峰。

圖8 實驗樣本中隨機取向纖維的長度和方向角的統(tǒng)計分布Fig. 8 Statistical distribution of length and orientation angle of randomly oriented fibers in experimental samples

2.2 RVE 尺寸效應及邊界影響區(qū)長度

2.2.1 RVE 尺寸效應

RVE 通常在下部長度尺內引入，來描述其上部長度尺度的性質，其定義基于代表性這個關鍵詞。如果RVE 能夠再現從有關長度尺度的無限大體積材料中獲得的特性，則認為RVE 具有代表性，這應該作為確定它的尺寸的一個標準。代表性通常是根據感興趣材料上部長尺度的有效性質來定義的。只要材料在其上部長度尺度上是均勻的，就可以在其下部長度尺度上引入RVE 作為材料的適當體積。在多尺度模型中，在上部尺度中的任何有限維在下部尺度上都被認為是無限的。換言之，RVE 的維數必須在下部長度尺度上是有限的，而它在其上部長度尺度上可以被視為數學上的無窮小。在較低長度尺度上的無限維體積總是具有代表性，但不是理想的RVE 候選體。

根據Hill[9]的理論，RVE 尺寸對數值模擬的精度影響很大程度上取決于邊界條件影響的區(qū)域深度。它是指從曲線或曲面到邊界效應減小點的距離，是一個特征度量。雖然它沿RVE 的表面逐點變化，但它是有界的，其上限稱為衰減長度t，t是給定材料的確定值。受邊界條件影響的區(qū)域的體積Vs影響。Vs可以等效地看作是衰減長度厚度的殼，是RVE 的最外層。由于對于給定的材料，t通常是一個固定值，通過衰減長度到RVE 的距離，應力和應變場不再受強加邊界條件的影響，例如，是否通過施加均勻位移或均勻牽引來完成。在將具有衰減長度厚度的殼體從RVE 上剝離后，留下的芯中的應力和應變場應與通過分析無限大體積材料獲得的應力和應變場相同。例如，考慮具有邊長a、衰減長度t和體積為VRVE(RVE 的體積)的立方RVE。隨著RVE 的大小增加，a增加，而t保持不變。當t相對于a變得可以忽略時，兩個體積的比值r趨于消失(見式(31))，此時RVE 預測的有效性無限接近其上部尺度。

目前實驗測定衰減長度采用數字散斑全場應變儀測定試樣在加載前后的局部應變量的不同來計算衰減長度，而數值方法可以直接通過FEA 計算 獲得。

在此模擬中，采用空間方向敏感性Hotelling膨脹算法對CCF 鑲嵌環(huán)氧樹脂的RVE 進行填充。為研究RVE 尺寸對力學精度的影響，考慮的RVE尺 寸 分 別 為 0.25 mm×0.25 mm×0.25 mm 、 0.5 mm×0.5 mm×0.5mm 和1 mm×1 mm×1 mm，并對3 種尺寸的RVE 立方體基于HE 算法進行填充。力學實驗則選擇圖7(b)所示實驗樣本。界面剛度假設為0 Pa，并且在RVE 中實現了彈性剛度隨機離散分布。RVE 結構的網格劃分選取的是8 節(jié)點磚塊單元，單元尺寸均為5 μm×5 μm×5 μm。模型采用的是不考慮孔隙問題的理想化情況。

2.2.2 RVE 邊界影響區(qū)長度

圖9 顯示了3 種尺寸的RVE 結構中，采用周期 邊界條件HE 填充碳纖維的具體過程，體積分數由0.2%遞增至1.24%，圖中每列所對應的體積分數是一致的。其中VRVE=13mm 的最大纖維填充數量為 330根，VRVE=0.53mm的最大纖維填充數量為44根，而VRVE=0.253mm的最大纖維填充數量為8根。

圖9 3 種尺寸RVE 的纖維填充過程(纖維橫縱比：145DCCF=6.9 μm，纖維扭曲度=0.3)Fig. 9 Fiber filling process of CCF for RVE of three sizes(Fiber aspect ratio: 145DCCF=6.9 μm, Torsion of fiber= 0.3)

結構邊界條件的約束方法是通過有限元節(jié)點施加線性約束方程，并將周期邊界條件施加到帶有周期性的RVE 結構中。對所屬面(見圖10)施加約束方程，對于面(也可以是邊或者頂點)的每個節(jié)點上的任意自由度i，邊界約束方程如下：

圖10 RVE 面邊界約束方程Fig. 10 Surface boundary conditions of RVE

段落 true="1">式中：ui為第i個分量的位移，i=1, 2, 3；Lx、Ly、Lz分別為x軸、y軸、z軸方向的初始長度。

同時，對3 種不同尺寸的RVE 在x′軸方向上施加0.03 mm 的位移邊界條件，通過數值計算，獲取了對應的邊界影響區(qū)長度(見圖11)。

由于纖維與基體材料存在非常大的平均應力差，圖11 中把短纖維和界面材料從基體中去除，從而可以準確獲取邊界影響區(qū)長度的細節(jié)信息。從圖11 中可以看到：基體材料在此邊界條件的影響下，應力分布為20~55 MPa 之間；通過對比不同x軸方向深度的切面應力分布圖可知，由于存在邊界影響區(qū)，RVE 結構由外到內存在不同的應力分布梯度；VCCF=1.24%(短切碳纖維體積分數)的RVE結構，應力影響區(qū)長度t取值為15~20 μm 之間，它受到RVE 尺寸的影響并不明顯。在邊界影響區(qū)范圍內，外層平均應力值明顯大于內層。當t值超過20 μm 時，平均應力值不隨x深度變化而改變。

圖11 衰減長度t對RVE 計算精度的影響，VCCF=1.24%Fig. 11 Influence of attenuation lengthton the accuracy of RVE calculation,VCCF=1.24%

為比較隨機纖維體積分數，UC尺寸(見圖3)，以及纖維扭曲度在3種尺寸RVE中對t的影響，對RVE結構模型進行了系統(tǒng)的研究，如圖12所示。首先研究不同VCCF下3種尺寸RVE的衰減長度與實驗測定的值進行對比(見圖11(a))。由于短切纖維的不斷填充，RVE結構內兩相界面的指數式增加以及相界面無序性的增大，使得應力影響區(qū)效應變得明顯，整體的t值隨Vccf呈指數變大，其對t值的影響范圍約為17 μm?？墒峭ㄟ^圖12(a)對比，在其他條件(各相組成、VCCF、扭曲度等)相同前提下，t值與RVE尺寸無關，而實驗測定的t值比數值模擬的結果略大。

此外，填充纖維的彎曲度對t值也存在一定的影響(見圖12(b))，扭曲度的影響與t值基本呈線性關系，但是弱于VCCF的影響，其影響范圍約為4 μm。圖12(c)所示為UC的尺寸對t值的影響。對于UC結構(見圖3)，考慮它是由纖維長度LCCF和基體長度Lm共同組成，并且LCCF為確定的值，通常是通過控制Lm來調節(jié)UC的尺寸，因此這里采用Lm/LCCF的值(UC中基體長度與纖維長度的比值)來確定UC的尺寸。結果顯示，UC中基體相越大，t值越小，其對t值的影響范圍約為1.6 μm。

圖12 RVE 體積，纖維扭曲度和UC 結構與RVE 應力影響區(qū)長度關系Fig. 12 Relationship among RVE volume, fiber tortuosity, UC structure and the length of RVE

對比圖3 可知，對RVE 衰減長度t的影響因素有很多，但是影響程度為VCCF>扭曲度 >Lm/LCCF>VRVE。

2.3 HE 算法的模量均一化響應

由于CCF 在基體中的隨機分散，復合結構中細觀層面的結構模量存在差異性，這主要體現在纖維和基體界面層的模量擾動。過大的擾動也給結構模量均一化過程帶來了相應的誤差。本文研究6 組方向模量，包括3 個方向楊氏模量、3 個剪切模量和3 個泊松比。

圖13 所示為尺寸VRVE=0.5 mm3，VCCF=1.24 以及旋轉角度α=45°、β=90° RVE 中的6 個方向模量和3 個泊松比的計算結果。從圖13 中可以看出：E1′、E2′的均值在270~330 GPa 之間，由于受到旋轉角度的影響，E3′均值模量略小，其范圍在230~300 GPa之間；G1′2′、G2′3′和G3′1′均值在500 GPa 以上；泊松比ν1′2′、ν2′3′和ν3′1′的均值在0.22~0.27 之間，并且基體相分布略大于纖維相，與實際情況相符。

圖13 VCCF=1.24%以及旋轉張量角度α=45°，β=90°時的模量空間隨機分布Fig. 13 Spatial random distribution of modulus withVCCF=1.24% and rotation tensor angleα=45° andβ=90°

RVE結構中方向模量的平均響應主要受到旋轉角度(見圖3)的影響，因此通過固定一個旋轉角度(α或β)來研究方向模量受到另一個角度旋轉的變化，旋轉角度α和β在RVE中的取值范圍均為0°~90°之間。圖14和圖15分別顯示了將旋轉角度α和β固定(45°)，方向模量隨另一個角度的變化情況。結果顯示，復合結構中平均剪切模量變化范圍為550~600 GPa之間，平均楊氏模量變化范圍為280~320 GPa之間，泊松比的平均變化范圍為0.231~0.257。

圖14 隨角β改變的平均方向模量的變化Fig. 14 Average directional modulus varying with rotation angleβ

圖15 隨角α改變的平均方向模量的變化Fig. 15 Average directional modulus varying with rotation angleα

由圖14和圖15可見：G1′2′、G2′3′和E2’受角度β的變化較大，G2′3′、E1′和E2′受角度α的影較大；而泊松比均一化結果是隨機分布的(見圖14(b))，與旋轉角度變化沒有明顯關系；當α=45°時，6個方向模量隨β的0-90°變 化 值 的 大 小 關 系 為G2'3'>G1′2′>G3′1′>E2′>E1′>E3'；β=45°時，6個方向模量隨α的0°~90°變化值的大小關系為G2′3′>E2′>E1′>G1′2′>G3′1′>E3'；受兩個旋轉角度影響最大的方向模量均為G2′3′，它隨兩個角度的變化值分別為50 GPa和72 GPa。表明G2′3′對于旋轉角度更敏感，無論是哪個角度的旋轉，在2′3′平面上，有更多的隨機短纖維排布由緯向轉變?yōu)閺较?，或者由徑向轉變成緯向。

除此之外，為驗證本文中方向模量均一化的準確性，與拉伸實驗樣本(見圖7)測定的結果進行了對比，對比結果如表2 所示。

表2 當α=β=45°時，均一化模量與實驗對比Table 2 Comparison of the homogenization modulus with the experiment when α=β=45°

根據表2 中所示，實驗測定的和有限元模量均一化計算出的結果差異在2~7 GPa 之間，并且實驗測定的數值比有限元計算出的數值偏大，可以根據圖11(a)中的結果來說明，因為實驗測定的應力影響區(qū)厚度較數值分析的大，而應力影響區(qū)存在較大的集中應力，因此會對模量的均一化評估產生偏大的影響。

3 結論

本文提出一種采用Hotelling expansion 模型進行CCF 的三維隨機填充的均一化方法，并對該方法生成的RVE 結構進行了研究。首先通過數字散斑全場應變測量儀測定了CCF 增強環(huán)氧樹脂基復合結構的衰減長度，并利用數值模擬的方法研究衰減長度的影響因素。然后利用所提出的均一化模型研究了短纖維空間分布引起的方向模量各向異性響應行為。根據數值實驗結果，討論了旋轉角度對短纖維的空間隨機分布對復合材料方向模量的影響。得出主要結論如下：

1)提出一種有效預測隨機分布短纖維復合材料方向模量均一化的解決方案。在研究復合材料均一化方向模量時，由于不同方向模量與旋轉角度之間存在的敏感度各異，導致各方向模量均一化結果隨旋轉角度變化存在大小差異。

2)短纖維在空間上的隨機分布狀態(tài)的不同導致了RVE 衰減長度尺寸的變化。受實驗所需成本和時間的限制，僅對RVE 體積VRVE、纖維扭曲度，UC 結構與衰減長度t的測試，并得出其相應的變化情況。數值計算結果對比實驗數據的結果較小主要是由于除去內嵌的纖維，真實的材料中存在氣泡，夾雜物等微小孔洞，而數值分析并沒有將其考慮進去。這一干擾效應的屏蔽及其對有效數值計算精度的修正還有待于進一步的研究。