平紋機織碳纖維復合材料的多尺度隨機力學性能預測研究1)

許 燦 朱 平,2) 劉 釗 陶 威

?(上海交通大學機械系統與振動國家重點實驗室,上海 200240)

?(上海市復雜薄板結構數字化制造重點實驗室,上海 200240)

??(上海交通大學設計學院,上海 200240)

引言

碳纖維增強聚合物基 (carbon fiber reinforced polymer,CFRP) 復合材料,簡稱碳纖維復合材料,具有密度低、比剛度高、比強度高、可設計性強等優點,能夠很好地實現結構的輕量化和高性能設計[1].作為一種新型輕質復合材料,碳纖維復合材料因其優異的力學性能在航空航天、汽車、船舶等行業得到了廣泛的應用[2-3].平紋機織CFRP 的制備過程包括編織、鋪層、浸潤和固化等,得到的材料在結構上具有多尺度特征.在CFRP 的制備過程中,由于存儲條件和組成相成分、批次的不同導致組成相材料性能有所差異,同時由于織物和預浸織布的生產、存儲和搬運等過程導致紡織結構的變動.這種材料和結構的不確定性反映到CFRP 性質上表現為各尺度力學性能均有一定的隨機性,并最終影響材料的宏觀性能[4-7].因此,有必要發展可靠的CFRP 隨機力學性能預測方法,充分考慮各尺度材料和結構參數的不確定性,研究其對材料宏觀力學性能預測結果的影響作用.

為了實現對考慮隨機性的復合材料力學性能的準確預測,國內外學者開展了廣泛的研究工作.與有限元技術相結合的計算細觀力學方法可以在各個尺度上建立考慮隨機性的精細仿真模型來進行力學性能預測[8-11].在纖維束尺度上,目前研究主要是通過纖維角度和位移等參數的擾動來生成滿足給定參數組合和約束條件下不規則排列結構.Zhu 等[12]提出了基于序列隨機擾動算法生成纖維隨機分布結構,并基于有限元模型研究了纖維隨機分布對纖維束力學性能的影響.類似的研究還包括Yang 等[13]提出的隨機序列展開方法,Zhang 等[14]提出的結合隨機分布和完全彈性碰撞的生成算法.在介觀尺度上,主要是基于圖像分析與幾何參數統計學描述方法相結合的方法生成具有結構隨機性的代表性體積單元(representative volume element,RVE) 進行仿真分析.Bale等[15]基于微計算機斷層掃描技術得到三維圖像并對纖維束位置和形狀進行了統計學描述.Blacklock等[16]則在此基礎上重構了具有隨機性的RVE 模型.Zhu 等[17]則進一步考慮了參數之間的相關性,提出了基于相關性混合高斯隨機序列的平紋機織CFRP的RVE 重構算法.

基于計算細觀力學的方法本質上是建立多個能夠滿足統計學表征的有限元模型并統計仿真分析的結果,因此計算代價過高.而解析細觀力學方法由于具有較高的計算效率,在進行考慮不確定性的預測研究時具有很大的優勢.Chakraborty 等[18]采用多項式相關函數展開實現了復合材料層合板的隨機自由振動分析.考慮到微觀結構的不確定性,隨機有限元[19]被采用進行有效的預測.Cui 等[20]采用Copula 函數表征復合材料中的纖維鋪設角度和材料參數之間的相關性,并基于攝動法進行不確定性分析.Mehrez等[21]提出了基于混沌多項式展開(polynomial chaos expansion,PCE) 的層級式傳遞方法得到宏觀材料力學性能的不確定性.Xu 等[22]則在層級傳遞的過程中進一步考慮了分布的真實形式和分布之間的高維相關性,并應用于三維正交機織復合材料的力學性能預測中.

現有的面向平紋機織CFRP 的研究主要存在以下幾個問題:一是所考慮的隨機參數不夠全面,大多研究僅考慮部分尺度的隨機參數.且采用的不確定性傳遞方法大多具有局限性:隨機有限元法,作為一種嵌入式的不確定性算法,對問題不具備很好的適應性;攝動法則對概率分布形式具有要求;二是在力學性能預測中對隨機參數之間的相關性考慮不夠充分.針對以上這些問題,本文提出一種基于PCE 和Vine Copula 的平紋機織CFRP 隨機力學性能預測方法,該方法考慮了微觀和介觀尺度中不確定性材料、結構參數,基于自下而上層級傳遞的策略逐層研究平紋機織CFRP 的各尺度力學性能參數的不確定性.所采用的非嵌入式PCE 方法能夠高效準確地實現不確定性地傳遞.此外,Vine Copula 函數能準確構造相關隨機參數之間的高維聯合概率分布,從而有效掌握相關性對性能響應的影響,提高平紋機織CFRP 力學性能參數預測的精度和可靠性.

1 確定性模型構建

采用解析細觀力學方法實現平紋機織CFRP 的多尺度確定性隨機彈性力學性能預測,該方法能充分考慮影響平紋機織CFRP 力學性能的因素,包括組分材料的力學性能、各尺度的幾何參數、體積分數-等[23].尺度包括纖維絲尺度(微觀尺度)、纖維束尺度(介觀尺度)和單胞尺度(宏觀尺度),如圖1 所示.該方法基于連續介質力學及均勻化理論,采用自下而上多尺度式串行的策略逐級預測力學性能.

圖1 平紋機織CFRP 多尺度特征Fig.1 Multiscale characteristics of plain woven CFRP

1.1 纖維束模型

平紋機織CFRP 中,每一根浸潤過的纖維束都包含了嵌入在基體中的許多單向纖維,浸潤纖維束被認為是單向復合材料.若纖維束的纖維絲體積分數為Vf,則力學性能可由下述Chamis 方程得到[24]

式中,下標f 和m 分別表示纖維絲和基體,V表示體積分數,E表示彈性模量,G表示剪切模量,u表示泊松比.

纖維絲的體積分數可以通過下式計算

其中,Sf表示纖維束中纖維絲的總截面面積,d表示纖維絲直徑,Nfiber表示纖維束中纖維絲數量,a和b分別表示纖維束的長軸和短軸長度,如圖2 所示.

圖2 單胞截面Fig.2 Cross section of unit cell

由于纖維束通常被假設為橫向各向同性的,故其剛度矩陣可以定義為

其中,剛度矩陣各元素可通過下式獲得

由于平紋機織CFRP 的纖維束存在卷曲變形,對于卷曲的纖維束,如圖3 所示.

圖3 纖維束的正弦函數表達Fig.3 Sinusoidal expression for yarn

本文采用正弦函數來描述其構型

式中,T表示周期,A表示幅度.綜合圖2 和圖3 可以得到T=2(a+l),A=(h?b)/2.將卷曲纖維束劃分為無數的微段,各微段沿纖維主方向均可視作單向纖維復合材料,如圖4 所示,xL方向即為主方向.

圖4 纖維束材料主方向Fig.4 Material principal orientation in yarn

定義第N微段在其局部坐標系下的剛度矩陣為

式中,上標L表示微段纖維局部坐標系(xL,yL),上標F表示纖維束坐標系(xF,yF),上標T 表示轉置,H為轉換矩陣,其求解公式為

式中

式中

式中θ1代表微段纖維局部坐標系與纖維束坐標系的夾角,如圖4 所示.在纖維束坐標系下,纖維束的等效剛度矩陣可通過對各微段沿x方向進行積分來獲得

1.2 單胞模型

對于CFRP,在宏觀尺度上選取單胞來表征其宏觀性能.由于平面機織CFRP 是由經紗和緯紗機織而成,單胞整體的剛度矩陣可根據各纖維束的平均性能計算得到,由于經紗和緯紗的方向不同,需通過坐標轉換,將纖維束剛度矩陣轉化為單胞坐標系下的纖維束剛度矩陣

轉換矩陣Hunit如下

式中

式中,θ2代表纖維束坐標系(xF,yF) 與單胞坐標系(xunit,yunit) 的夾角,如圖2 所示.從而可得平紋機織CFRP 的宏觀剛度陣

式中,Vunit表示單胞的體積,下標warp 表示經紗,下標weft 表示緯紗,表示基體占單胞的體積分數.由于紗線是卷曲的,需要對其沿x方向的長度進行積分獲得伸展的纖維束長度,如圖3 所示.經緯紗和基體的體積分數可通過式(15)來獲得

通過對宏觀剛度矩陣進行求解,可獲得單胞的彈性模量、剪切模量和泊松比.

1.3 模型驗證

本文所制備的平紋機織CFRP 碳纖維絲和基體材料力學性能參數如表1 所示,碳纖維絲選用臺麗碳纖維TC33,碳纖維以3K 的形式集合成碳纖維束.基體選擇Huntsman 公司的環氧樹脂,牌號為LY1564 SP,對應的固化劑牌號為Aradur 3486,兩者的體積比為5:2.采用真空導入工藝進行力學性能試驗所用的平紋機織CFRP 的成型.通過上述工藝制備得到的CFRP 的密度為1.47 g/cm3.

表1 碳纖維和基體基本力學性能參數Table 1 Basic mechanical properties of carbon fiber and matrix

為了從宏觀上探究CFRP 的力學性能,分別根據試驗標準ASTM D638[25]和ASTM D7078[26]進行平紋機織CFRP 的軸向拉伸和面內剪切試驗,試驗樣件如圖5 所示,尺寸單位為mm.準靜態軸向拉伸力學性能試驗使用島津5 t EHF-UM 電液伺服疲勞試驗機進行,準靜態面內剪切試驗采用SUNS 電子萬能試驗機.所有試驗均采用位移控制方式,并根據樣件的尺寸設定加載速率,保證對應的工程應變率在0.001 s?1范圍.每種試驗均保證測量得到5 個有效試驗數據并取平均值.最終試驗得到的宏觀軸向拉伸彈性模量為60.89 GPa,面內剪切模量為3.63 GPa.

纖維絲和纖維束的幾何參數通過對纖維束界面顯微圖像和Micro-CT 圖像(如圖6 和圖7 所示) 進行數據統計可獲得,纖維絲的直徑為d=6.23μm,纖維束的幾何參數如表2 所示.表中平均值為統計得到的均值,上下限值為統計得到的±3σ 置信區域數值.纖維絲占纖維束體積分數的均值為65.32%,纖維束占單胞體積分數的均值為67.90%.

圖5 準靜態力學性能試驗樣件Fig.5 Specimens for quasi-static mechanical tests

圖6 纖維束截面顯微圖像Fig.6 Microscopic image of cross section for yarns

圖7 Micro-CT 圖像Fig.7 Micro-CT images

表2 幾何參數Table 2 Geometric parameters

以各變量的均值為輸入參數,經過前述的解析模型預測得到的單向纖維和單胞的確定性力學性能如表3 所示,通過與有限元仿真結果的對比可看出解析法的相對誤差較小.特別地,解析法得到的結果與試驗結果的相對誤差分別為0.82%和6.33%,預測結果與試驗結果有較好的一致性.

表3 力學性能預測結果Table 3 Prediction results of mechanical properties

2 纖維束隨機力學性能預測

本文主要考慮組分材料力學性能和幾何參數不確定性對單向纖維力學性能預測結果的影響,組分材料的力學性能參數和幾何參數的概率密度函數(probability density function,PDF)均設定為具有1%變異系數的高斯分布.

圖8 平紋機織CFRP 層級結構Fig.8 Hierarchical structure of plain woven CFRP

平紋機織CFRP 的層級結構如圖8 所示,由于纖維束長軸和短軸既參與體積分數Vf的計算,又參與了體積分數Vwarp的計算,屬于跨尺度的共享變量.單向纖維力學性能輸入包括7 個組分材料參數和幾何參數a,b,各參數之間相互獨立,故直接采用混沌多項式方法來進行預測.考慮到單向纖維具有5 個獨立的材料參數,分別選擇C11yarn,C22yarn,C12yarn,C23yarn,C66yarn作為纖維束尺度模型力學性能響應輸出.

2.1 混沌多項式展開

PCE 實質上是對隨機變量構造一個具有隨機性的代理模型[27-29].為了描述的方便,設定X=(E11f,E22f,G11f,G23f,v12f,Em,vm,a,b),其任一元素設定為X.采用PCE,隨機響應Q的截斷的PCE 可以簡潔地表達如下

式中,δαβ是克羅內克符號,E(·) 表示期望算子.式(16)所對應的截斷集定義如下

利用PCE 進行不確定性傳遞時,最關鍵的步驟是求解式(16)中待定的多項式系數qα,本文采用回歸法求解待定系數,其回歸表達式為

求解上式得到多項式中的待定系數.基于正交性質,前四階矩,包括均值μ、標準差σ、偏度系數Cs、峰度系數Ck可以通過下式結合多項式系數獲得.

基于前四階并采用概率分布擬合方法可以獲得響應的概率分布.常用的概率分布方法包括最大熵原理,鞍點逼近等[30-32].若偏度系數約為0,峰度系數約為3,則可直接假設該分布服從正態分布從而得到概率分布.

本文采用留一交叉驗證方法來驗證所構造的PCE 的預測精度.假設原始數據有N個樣本,那么每個樣本單獨作為驗證集,其余的N?1 個樣本作為驗證集.留一交叉驗證能充分利用所有樣本,并得到N個模型,用這N個模型驗證集的相對誤差平均數作為驗證誤差以評估模型的精度.

2.2 預測結果分析

將PCE 方法應用到平紋機織CFRP 單向纖維的力學性能預測中,設定多項式階數為3,基于拉丁超立方采樣方法采取原始數據集,包含50 個樣本.預測得到的結果及模型驗證誤差如表4 所示.

表4 單向纖維力學性能隨機預測結果Table 4 Stochastic prediction results of mechanical properties for unidirectional fiber

所構造的PCE 模型采用留一交叉驗證方法來驗證模型的精度,由表4 可知,力學性能參數預測模型的交叉驗證誤差最高為0.470,模型精度很高.偏度系數均在0 左右,峰度系數均在3 左右,說明單向纖維力學參數基本服從正態分布.

3 單胞隨機力學性能預測

單胞尺度下,纖維束預測得到的力學性能響應屬于同一模型的不同響應,故不可避免地存在變量相關性.此外,幾何參數a,b由于既參與纖維絲占纖維束體積分數Vf的計算,又參與纖維束占單胞體積分數Vwarp和Vweft的計算;而Vf又與C11yarn等剛度矩陣元素存在因果關系,故幾何參數a,b與C11yarn等剛度矩陣元素之間也存在相關性.隨機響應預測過程中若不考慮變量之間的相關性,會使得預測結果出現較大的偏離[33].因此,單胞尺度的力學性能預測需要充分考慮隨機參數之間的相關性,本文基于Vine Copula 理論首先構造相關隨機變量的聯合概率分布,然后基于Rosenblatt 轉換實現相關樣本至獨立樣本的轉換,最后基于相互獨立的樣本并采用PCE 來實現單胞尺度的拉伸彈性模量Ex=Ey,拉伸彈性模量Ez和面內剪切模量Gxy,面外剪切模量Gxz=Gyz,泊松比vxy,及vxz=vyz的隨機預測,如圖8 所示.

3.1 基于Vine Copula 的相關性模型構建

Copula 函數能夠建立一元邊緣累積概率分布與多元聯合分布之間的關系[34-35],為了描述的方便,設定Y=(C11yarn,C22yarn,C12yarn,C23yarn,C66yarn,a,b,h,l),其聯合累積概率分布(cumulative distribution function,CDF)表達式為

式中,NU為變量個數,Pr(·) 為累積概率算子.基于Sklar’s 理論,聯合CDF 可以重構成一系列邊緣CDF的函數.

式中,C(·) 表示Copula 函數,當所有的邊緣CDF 都是連續函數時,Copula 函數是唯一的; ? 表示Copula參數向量.式(22)可進一步寫成

式中,u為累積概率值向量.Copula 函數有其對應的Copula 密度函數

基于式(24)可得聯合PDF

對于相關性測度,本文采用一種非線性測度Kendall 系數來進行描述.相關模型的構建是基于累積概率樣本集,因此原始樣本集需要轉換到對應空間來進行相關模型的構造,即Y→U.Copula 理論為多元相關建模提供了一種合理的方法,但是目前主要針對二元Copula 函數.隨著維度的增加,構造一個合適的多元Copula 會變得越來越困難.為了克服上述局限性,采用Vine Copula 來將多元Copula 函數分解為一系列二元Copula 函數的組合來構造多元相關模型.分解方法如下

式中,fj|1,2,···,j?1,j=2,3,···NU是邊緣PDF,其可以通過Copula 密度函數的偏導求得

觀察式(26)不難發現通過改變分解變量的順序,多元Copula 函數的分解方式多樣.本文采用基于圖形模型的R-vine 來進行結構構造[36].該方法通過分析變量之間的相關系數并基于最大生成樹方法來確定R-Vine 的結構[37].二元Copula 函數存在多種形式,如Gumbel Copula,Gaussian Copula,Frank Copula等,且不同函數存在對應的Copula 參數.分別引入最大似然估計方法和赤池信息準則(Akaike information criterion,AIC)[38]來確定最優Copula 參數和函數.

式中,k?是需要確定的Copula 參數數目.通過上述方法的組合可逐步構造二元Copula 函數直到獲得所需要的多元Copula 函數.

3.2 獨立性轉換

基于Vine Copula 理論構造相關性模型后,可通過該模型生成相關累積概率樣本集u,該樣本集無法直接用于后續的不確定性分析,故須先通過獨立性轉換方法得到相互獨立的樣本集.由于最大生成樹方法已經確定了Vine Copula 的結構,即確定了隨機變量的轉換次序,故此時采用Rosenblatt 轉換會得到唯一一組獨立樣本集.轉換方法如下

Rosenblatt 轉換實現了樣本從相關到獨立的轉換,即進一步地采用可得到用于不確定性傳遞的樣本集,其轉換表達式為

3.3 預測結果分析

將Vine Copula 理論應用到平紋機織CFRP 單胞尺度力學性能預測的相關模型建模中.具有相關性的參數為(C11yarn,C22yarn,C12yarn,C23yarn,C66yarn,a,b),幾何參數a和b的PDF 已設定為高斯分布,力學性能參數(C11yarn,C22yarn,C12yarn,C23yarn,C66yarn)的PDFs 則已通過纖維束尺度的隨機預測獲得.得到的R-Vine樹結構1 如圖9 所示.

圖9 R-Vine 樹結構1Fig.9 R-Vine tree structure 1

為了驗證所構造的相關性模型的準確性,基于該模型產生隨機樣本并求解得到相關性系數,該相關性系數與原樣本集相關性系數進行對比,對比結果如圖10 中相關性矩陣所示.相關性矩陣上對角矩陣的相關性系數為原樣本集計算得到的結果,下三角矩陣的相關系數為Vine Copula 模型計算得到的結果.對比矩陣里面的元素可發現相關性系數幾乎一致,絕對誤差最大僅為0.03,說明構造的相關性模型精度很高.

圖10 相關性矩陣對比Fig.10 Comparison of correlation matrix

基于獨立性轉換方法將原樣本集Y轉換為用于不確定性傳遞的樣本集,圖11 給出了(C11yarn,C22yarn)在空間Y→U→→轉換過程.可以看出,最終獲得了相互獨立的樣本集.基于該樣本集并采用前述的PCE 方法來實現平紋機織CFRP 單胞的力學性能預測.設定多項式階數為5,預測得到的結果如表5 所示.

圖11 樣本集轉換過程Fig.11 Transformation of samples

表5 單胞力學性能隨機預測結果Table 5 Stochastic prediction results of mechanical properties of unit cell

所構造的PCE 模型的驗證誤差最大為0.080 6,模型精度很高.由于偏度系數值在0 左右,峰度系數值在3 左右,說明單胞彈性力學性能參數基本服從正態分布.將單胞力學性能隨機預測結果與表1 中的確定性結果進行對比可以發現,隨機預測結果的均值與確定性結果基本一致,這說明了平紋機織CFRP 軸向拉伸彈性模量和面內剪切模量的多尺度模型非線性程度較低,使得隨機預測結果的均值相對于確定性結果波動較小; 也從側面證明了傳統確定性方法實際上就是基于隨機變量的均值來進行分析.此外,通過計算得到隨機預測結果的變異系數分別為1.64%,2.02%,2.32%,2.63%,1.38%,1.13%相比較于輸入參數1%的變異系數來說均有所提高,表明不確定性相互耦合且傳遞之后會對隨機預測結果的方差產生放大的影響.

4 結論

本文針對現有平紋機織CFRP 力學性能預測研究中的考慮隨機參數不夠全面,不確定性傳遞方法大多具有局限性以及對隨機參數之間的相關性考慮不夠充分等問題,提出一種基于PCE 和Vine Copula的平紋機織CFRP 隨機力學性能多尺度預測方法.具體工作與結論如下:

(1)該方法綜合考慮了平紋機織CFRP 微觀及介觀尺度的材料、結構隨機參數,基于層級傳遞的策略逐尺度地研究了力學性能參數的不確定性.本方法從根源上量化了不確定性,并逐尺度地揭示了參數的不確定性對各尺度力學性能的影響作用.

(2)隨機力學性能預測結果表明,所構造的各尺度PCE 模型的交叉驗證誤差均很低,說明非嵌入式PCE 方法能夠準確地實現各尺度力學性能的隨機預測.同時,由于PCE 模型能夠解析地給出預測結果的均值、標準差、偏度系數和峰度系數,能夠進一步得到參數的PDF.

(3)該方法在預測過程中充分考慮了隨機變量之間的高維相關性,基于Vine Copula 方法和Rosenblatt轉換實現了對相關性的建模與轉換.相關性矩陣中系數之間的絕對誤差最大為0.03,表明所構造的相關性模型能夠很好地反映原數據之間的相關性;而Rosenblatt 轉換與所構造的相關性模型相結合則能夠很好地實現樣本集的獨立性轉換.