基于三階剪切計(jì)算連續(xù)性方法的編織復(fù)合材料板三尺度分析模型

萬傲霜, 王振名,李頂河

中國(guó)民航大學(xué) 航空工程學(xué)院，天津 300300

與傳統(tǒng)的多向?qū)雍蠌?fù)合材料相比，編織復(fù)合材料具有優(yōu)良的面外性能，在航空航天、汽車和船舶等領(lǐng)域的應(yīng)用日益廣泛[1-2]。美國(guó)航空航天局于1988年啟動(dòng)先進(jìn)復(fù)合材料技術(shù)(ACT)計(jì)劃，將編織復(fù)合材料應(yīng)用于飛機(jī)結(jié)構(gòu)(見圖1)[3-4]。編織復(fù)合材料也被大量應(yīng)用于直升機(jī)，包括機(jī)身蒙皮、主旋翼槳葉、尾梁和尾槳等[5]。

編織復(fù)合材料性能不僅取決于組分材料的力學(xué)性能和纖維體積分?jǐn)?shù)，還取決于介觀尺度上編織纖維束的幾何構(gòu)型，需要采用介觀力學(xué)方法精確模擬計(jì)算編織復(fù)合材料的力學(xué)響應(yīng)，這對(duì)計(jì)算能力提出了很高的要求。為了兼顧計(jì)算成本和計(jì)算精度，基于均勻化方法的復(fù)合材料多尺度分析得到了大量研究。Terada和Kikuchi[6]、Matsui等[7]采用相關(guān)變分表述，將均質(zhì)材料的兩尺度建模方法用于具有精細(xì)周期性微觀結(jié)構(gòu)的非均質(zhì)材料。Smit等[8]提出了一種在宏觀/微觀尺度上均考慮材料黏彈性和大變形的均勻化方法。Feyel和Chaboche[9]采用多級(jí)有限元方法建立了新的復(fù)合材料多尺度模型，考慮了纖維和基體力學(xué)性能之間的差異?；跐u近均勻化理論[10-11]，Ghosh等[12]建立了非均質(zhì)材料(包括多孔材料和復(fù)合材料)的多尺度彈塑性分析的有限元模型。Miehe[13]基于精細(xì)尺度的位移場(chǎng)結(jié)果，分析了微觀結(jié)構(gòu)線性位移、周期性波動(dòng)和邊界常應(yīng)力的全局變分問題。Cao[14]計(jì)算得到了復(fù)合材料彈性靜力學(xué)問題解的迭代兩尺度漸近展開式。Zhang等[15]利用多尺度有限元方法，建立了宏觀尺度位移與微觀尺度應(yīng)力應(yīng)變之間的關(guān)系。袁欣和孫慧玉[16]利用均勻化方法計(jì)算分析了三維編織復(fù)合材料在不同溫度下的恒溫松弛模量，并進(jìn)一步建立了由恒溫松弛模量確定變溫松弛模量的一般方法。楊強(qiáng)等[17]基于多尺度漸進(jìn)展開理論，對(duì)復(fù)合材料彈性問題控制方程進(jìn)行尺度分解，建立了復(fù)合材料宏觀和細(xì)觀尺度響應(yīng)之間的關(guān)聯(lián)。胡殿印等[18]采用通用單胞(GMC)方法，建立了平紋編織復(fù)合材料的宏細(xì)觀多尺度模型，通過與傳統(tǒng)有限元單胞模型的宏細(xì)觀力學(xué)響應(yīng)分析計(jì)算結(jié)果對(duì)比，驗(yàn)證了模型的正確性。

圖1 編織復(fù)合材料飛機(jī)結(jié)構(gòu)件[3-4]

現(xiàn)有的均勻化方法研究大多集中于一階均勻化理論，該理論假設(shè)施加于計(jì)算胞元(CUC)或代表性體積元(RVE)的宏觀尺度變形梯度是常值，導(dǎo)致大尺寸CUC或大應(yīng)變梯度情況下的數(shù)值模擬結(jié)果不夠理想。為了捕捉復(fù)合材料場(chǎng)梯度和微觀結(jié)構(gòu)引起的局部效應(yīng)，通過引入適當(dāng)?shù)母唠A校正量，進(jìn)一步發(fā)展了二階均勻化理論，假設(shè)宏觀變形梯度在CUC求解域上線性變化[19-23]。

由于上述兩尺度分析方法難以捕捉編織復(fù)合材料這種具有復(fù)雜內(nèi)部結(jié)構(gòu)的復(fù)合材料在各尺度上的大量振蕩信息，有必要對(duì)編織復(fù)合材料進(jìn)行三尺度分析，包括對(duì)纖維和基體的微觀尺度分析、對(duì)纖維束和基體的介觀尺度分析、對(duì)復(fù)合材料結(jié)構(gòu)的宏觀尺度分析，大量學(xué)者對(duì)此開展了研究。Trucu等[24]、Abdulle和Bai[25]系統(tǒng)研究了采用不同偏微分方程三尺度問題的計(jì)算收斂性。Ma等[26]采用兩尺度漸近展開方法，得到了介觀尺度問題的特征函數(shù)和特征值，并進(jìn)一步推導(dǎo)了均勻化單元函數(shù)，建立了3個(gè)尺度下均勻化系數(shù)和材料系數(shù)之間的關(guān)系?；诰鶆蚧碚摵投喑叨葷u近展開方法，Yang等[27-31]提出了一種研究非均質(zhì)材料在微觀、介觀和宏觀尺度上耦合傳導(dǎo)輻射和熱力學(xué)問題的高階三尺度分析方法。但是，這種方法必須結(jié)合高階單元解和一致解來構(gòu)造三尺度漸近解，限制了該方法的應(yīng)用。

為了克服均勻化理論和廣義連續(xù)理論的局限性，文獻(xiàn)[32-35]提出了計(jì)算連續(xù)性(C2)方法，可以對(duì)有限大小非均質(zhì)體進(jìn)行粗尺度連續(xù)性描述，該方法的粗尺度控制方程在不相交胞元組成的計(jì)算連續(xù)域上表示，確定了胞元位置后，在細(xì)尺度上確定控制方程的弱形式。這種不進(jìn)行尺度分離的多尺度分析方法可以在不假設(shè)細(xì)尺度胞元無限小的條件下提供細(xì)尺度信息，并且不需要高階連續(xù)性、高階邊界條件和引入新的自由度。最近，C2方法被應(yīng)用于等效單層理論框架下的復(fù)合材料層合板多尺度分析[36-37]以及三階剪切變形理論(TSDT)框架下的混凝土結(jié)構(gòu)和復(fù)合材料曲梁多尺度分析[38-40]。

因此考慮到編織復(fù)合材料的宏觀力學(xué)性能取決于細(xì)觀組分材料性能和介觀纖維束編織構(gòu)型，為了在可接受的計(jì)算成本下準(zhǔn)確計(jì)算編織復(fù)合材料的位移場(chǎng)和應(yīng)力場(chǎng)，有必要建立三尺度有限元模型。因此，本文提出一種基于三階剪切計(jì)算連續(xù)性方法的編織復(fù)合材料板三尺度分析模型。采用基于TSDT的八節(jié)點(diǎn)四邊形板單元對(duì)宏觀尺度模型進(jìn)行離散，推導(dǎo)三尺度三階剪切計(jì)算連續(xù)性(TOS-C2)方法及其有限元方程；構(gòu)建編織復(fù)合材料板三尺度TOS-C2方法的計(jì)算框架，利用含夾雜三維立方體的數(shù)值算例，檢驗(yàn)該方法的有效性；利用三尺度TOS-C2方法，模擬計(jì)算平面編織復(fù)合材料板的微觀、介觀、宏觀尺度位移場(chǎng)和應(yīng)力場(chǎng)。

1 理論模型

1.1 三階剪切變形理論

首先三階剪切變形理論的位移假設(shè)為

(1)

其中：

(2)

式中：(x,y,z)為板上某一點(diǎn)的坐標(biāo)；h為板的厚度；u0、v0、w0分別為中心面上一點(diǎn)(x,y)在x、y、z方向的位移；u1、v1分別為(x,y)的剪切轉(zhuǎn)角；?w/?x、?w/?y分別為(x,y)的彎曲轉(zhuǎn)角。根據(jù)位移假設(shè)，應(yīng)變分量可推導(dǎo)為

(3)

其中：

(4)

由式(3)可以看出，橫向剪切應(yīng)變沿板厚呈拋物線變化，并在板的上表面和下表面達(dá)到0，不需要使用剪切校正因子。外載荷的應(yīng)變能和外力功可表示為

(5)

其中:

(6)

式中：C為剛度矩陣；A為面積；w為位移；εx、εy、εxy、εyz、εxz為應(yīng)變分量；σx、σy、σxy、σyz、σxz為應(yīng)力分量；q為面積力。式(5)為推導(dǎo)有限元方程的基礎(chǔ)。

層合板的本構(gòu)關(guān)系為

(7)

式中：Cij為剛度系數(shù)。

(8)

圖2 八節(jié)點(diǎn)四邊形板單元(坐標(biāo)系原點(diǎn)位于板單元中心面上)

根據(jù)式(3)，應(yīng)變-位移的關(guān)系可寫為

(9)

式中：εp和εt分別為面內(nèi)應(yīng)變和彎曲應(yīng)變；Lp和Lt分別為面內(nèi)應(yīng)變和彎曲應(yīng)變的位移關(guān)系矩陣。

任意一點(diǎn)的位移矢量u可表示為

u=Nde

(10)

式中：[N]5×48為形函數(shù)矩陣，下標(biāo)表示該矩陣的維度，具體表達(dá)式見附錄A；節(jié)點(diǎn)位移矢量[de]48×1為

(11)

式中：上標(biāo)(i)為節(jié)點(diǎn)編號(hào)。

那么，面內(nèi)應(yīng)變和彎曲應(yīng)變可以用節(jié)點(diǎn)位移矢量來表示

(12)

式中：Bp和Bt分別面內(nèi)和彎曲為應(yīng)變-位移矩陣，表達(dá)式為

(13)

總應(yīng)變可寫為

(14)

[Bp]3×48、[Bt]2×48、應(yīng)變-位移矩陣[B]5×48的表達(dá)式見附錄A。

1.2 三尺度計(jì)算連續(xù)性方法

基于C2多尺度方法的位移和應(yīng)變分解，三尺度計(jì)算連續(xù)性模型的位移場(chǎng)和應(yīng)變場(chǎng)分別可表示為

(15)

(16)

α,β,…,δ=1,2,3

(17)

1.3 有限元方程

1.3.1 微觀尺度

由單絲纖維和基體組成的微觀尺度問題的弱形式定義為

(18)

將微觀權(quán)函數(shù)umi和試探函數(shù)wmi用形函數(shù)Nmi(ψ)進(jìn)行離散:

(19)

式中：dmi、cmi分別為權(quán)函數(shù)和試探函數(shù)的節(jié)點(diǎn)值。

將式(14)～式(16)、式(19)代入微觀尺度胞元問題的弱形式中，并且考慮到試探函數(shù)cmi的任意性，微觀尺度胞元問題權(quán)函數(shù)在節(jié)點(diǎn)處解的表達(dá)式為

(20)

式(20)可以簡(jiǎn)寫為

(21)

1.3.2 介觀尺度

由纖維束和基體組成的介觀尺度問題的弱形式定義為

(22)

同樣地，將介觀權(quán)函數(shù)ume和試探函數(shù)wme用形函數(shù)Nme(χ)進(jìn)行離散，表達(dá)式為

(23)

式中：dme、cme分別為權(quán)函數(shù)、試探函數(shù)的節(jié)點(diǎn)值。根據(jù)介觀應(yīng)變的泰勒展開式，用權(quán)函數(shù)和試探函數(shù)對(duì)其進(jìn)行表示為

(24)

其中:

α=1,2,3

(25)

將式(16)、式(21)、式(24)代入介觀尺度胞元問題的弱形式中，并且考慮到試探函數(shù)cme的任意性，介觀尺度胞元問題權(quán)函數(shù)在節(jié)點(diǎn)處解的表達(dá)式為

(26)

1.3.3 宏觀尺度

宏觀尺度問題的弱形式可定義為

?wma∈WΩx

(27)

同樣地，將宏觀權(quán)函數(shù)uma和試探函數(shù)wma用形函數(shù)Nma(χ)進(jìn)行離散，表達(dá)式為

(28)

式中：dma、cma分別為權(quán)函數(shù)和試探函數(shù)的節(jié)點(diǎn)值。根據(jù)宏觀應(yīng)變的泰勒展開式，用權(quán)函數(shù)和試探函數(shù)對(duì)其進(jìn)行表示為

(29)

其中

(30)

式(21)、式(26)可以寫為

(31)

(32)

式中：

(33)

將式(31)、式(32)代入宏觀尺度問題的弱形式中，并且考慮到試探函數(shù)cma的任意性，宏觀尺度問題權(quán)函數(shù)在節(jié)點(diǎn)處解的表達(dá)式為

dma=(Kma)-1fma

(34)

其中：

(35)

2 計(jì)算框架

圖3為三尺度三階剪切計(jì)算連續(xù)性(TOS-C2)

圖3 三尺度TOS-C2方法計(jì)算框架

方法的計(jì)算框架，具體流程如下：

1) 建立宏觀問題、介觀胞元和微觀胞元的有限元模型。

4) 由宏觀尺度模型的有限元方程，計(jì)算宏觀位移dma(見式(34))。

7) 后處理得所需尺度模型位移和應(yīng)力分布。

基于Microsoft Visual Studio軟件，通過編寫程序?qū)崿F(xiàn)上述計(jì)算框架。

3 數(shù)值算例

3.1 驗(yàn)證算例

利用含夾雜的三維立方體算例來驗(yàn)證所提出的三尺度TOS-C2方法。分別采用直接數(shù)值模擬(DNS)方法和三尺度TOS-C2方法建立有限元模型(見圖4、圖5)，并進(jìn)行對(duì)比分析。該驗(yàn)證模型中的夾雜是指彈性模量較高的區(qū)域，類似微觀尺度上纖維束胞元模型的的單絲纖維。對(duì)于三尺度TOS-C2模型，微觀尺度胞元如圖4(a)所示；圖4(b)表示8個(gè)微觀尺度胞元組成的介觀尺度胞元模型夾雜域(見圖4(c))中的1個(gè)單元；整個(gè)介觀尺度胞元模型如圖4(d)所示。微觀和介觀尺度胞元模型均采用Hex8單元，其中介觀模型中采用非局部積分。由于1個(gè)宏觀尺度單元由4個(gè)介觀尺度胞元組成，宏觀尺度模型(見圖4(e))共包含256個(gè)介觀尺度胞元，分別在x、y方向有8個(gè)單元。對(duì)于DNS模型，單元尺寸則由微觀夾雜決定，整個(gè)DNS模型(見圖5(a))由夾雜(見圖5(b)、圖5(c))和基體組成。夾雜的材料屬性取彈性模量E=10.0 GPa，泊松比μ=0.3，基體的材料屬性取E=1.0 GPa，μ=0.3。宏觀模型的邊界條件為

(36)

將由DNS模型和三尺度TOS-C2模型計(jì)算得到的應(yīng)力和位移結(jié)果進(jìn)行對(duì)比分析，如圖6～圖8所示。從圖6中可以看出，三尺度TOS-C2模型計(jì)算得到的宏觀位移結(jié)果與DNS模型吻合良好，但宏觀應(yīng)力計(jì)算結(jié)果存在一定偏差。為了將三尺度TOS-C2模型的介觀位移和應(yīng)力結(jié)果與DNS模型進(jìn)行對(duì)比，選取如圖4(e)所示一條線上的A、B、C作為分析對(duì)象。圖7、圖8、表1對(duì)比了介觀尺度胞元模型和DNS模型的位移及應(yīng)力結(jié)果，可以看出2種模型的介觀尺度位移計(jì)算結(jié)果吻合良好，但應(yīng)力計(jì)算結(jié)果存在一定偏差。需要指出的是，因?yàn)榭紤]到需要將DNS模型的計(jì)算量控制在可求解范圍內(nèi)，網(wǎng)格劃分較為稀疏(TOS-C2介觀和微觀胞元模型分別為3×3×3和9×9×9單元數(shù)，而DNS模型中與介觀和微觀尺度相對(duì)應(yīng)均采用3×3×3單元數(shù))，導(dǎo)致應(yīng)力計(jì)算結(jié)果與DNS模型存在一定偏差，考慮到位移計(jì)算結(jié)果吻合良好，認(rèn)為該模型的計(jì)算誤差在可接受范圍內(nèi)。DNS模型和三尺度TOS-C2模型的單元數(shù)、節(jié)點(diǎn)數(shù)、計(jì)算時(shí)間和CPU數(shù)量見表2，可以看出三尺度TOS-C2模型在計(jì)算效率上有顯著優(yōu)勢(shì)。此外，三尺度TOS-C2方法不僅可以有效分析含夾雜三維立方體的應(yīng)力應(yīng)變響應(yīng)，由于該方法采用三階剪切變形理論，不存在剪切閉鎖問題，所以也能有效分析編織復(fù)合材料板的應(yīng)力應(yīng)變響應(yīng)。

圖4 含夾雜三維立方體的三尺度TOS-C2模型

圖5 含夾雜三維立方體的DNS模型

圖6 含夾雜三維立方體位移和應(yīng)力計(jì)算結(jié)果

表1 含夾雜三維立方體A、B、C 3點(diǎn)處DNS模型和TOS-C2介觀胞元模型的位移和應(yīng)力計(jì)算結(jié)果

表2 含夾雜三維立方體DNS模型和三尺度TOS-C2模型的計(jì)算量

3.2 平面編織復(fù)合材料板三尺度TOS-C2模型

圖9～圖11為采用提出的三尺度TOS-C2方法對(duì)平面編織復(fù)合材料板進(jìn)行三尺度分析的過程和結(jié)果。圖9(b)為宏觀尺度模型(見圖9(a))中第25個(gè)單元的第4個(gè)積分點(diǎn)對(duì)應(yīng)的介觀尺度胞元模型，分為2個(gè)部分，即介觀基體(見圖9(c))和編織纖維束(見圖9(d))。圖9(e)為放大后的介觀尺度纖維束，其中1個(gè)單元的求解域包含8個(gè)微觀尺度胞元(見圖9(f)、圖9(h))，圖9(g)、圖9(i)分別為8個(gè)微觀尺度胞元以及微觀尺度胞元中的單絲纖維。該三尺度模型中，介觀尺度和微觀尺度基體的材料屬性都是E=2.1 GPa，μ=0.3,微觀模型中單絲纖維的材料屬性為E=210 GPa，μ=0.3。宏觀尺度模型的邊界條件為四邊固支，并在底部施加均勻分布載荷p=1 Pa。

圖9 平面編織復(fù)合材料板的三尺度TOS-C2模型

圖10 平面編織復(fù)合材料板三尺度TOS-C2模型的位移計(jì)算結(jié)果

圖11 平面編織復(fù)合材料板TOS-C2微觀模型的應(yīng)力計(jì)算結(jié)果

4 結(jié) 論

基于三階剪切變形理論(TSDT)和計(jì)算連續(xù)性(C2)方法，提出了新的編織復(fù)合材料板三尺度分析模型，得到如下主要結(jié)論：

1) 采用基于TSDT的八節(jié)點(diǎn)四邊形板單元對(duì)宏觀尺度模型進(jìn)行離散，可以考慮剪切變形和沿厚度方向非線性剪切應(yīng)變變化引起的轉(zhuǎn)動(dòng)慣量。

2) 利用含夾雜三維立方體的數(shù)值算例，驗(yàn)證了提出的三尺度三階剪切計(jì)算連續(xù)性(TOS-C2)方法的有效性，三尺度TOS-C2模型的位移和應(yīng)力計(jì)算結(jié)果與直接數(shù)值模擬(DNS)模型吻合良好。

3) 提出的三尺度TOS-C2方法可以詳細(xì)描述編織復(fù)合材料板的微觀和介觀尺度結(jié)構(gòu)以及三尺度位移和應(yīng)力分布，為編織復(fù)合材料板航空航天結(jié)構(gòu)設(shè)計(jì)和力學(xué)性能分析提供技術(shù)支持，并且，由于三尺度TOS-C2方法中宏觀和介觀尺度應(yīng)變均采用泰勒級(jí)數(shù)展開，該方法可以適用于應(yīng)力變化劇烈的復(fù)雜單元問題。