復(fù)合材料周期結(jié)構(gòu)數(shù)學(xué)均勻化方法的一種新型單胞邊界條件

朱曉鵬， 陳 磊, 黃 俊

(1.安徽華電工程咨詢(xún)?cè)O(shè)計(jì)有限公司，合肥 230022； 2.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191；3.北京航空航天大學(xué)合肥創(chuàng)新研究院，合肥 230012)

1 引 言

復(fù)合材料具有比強(qiáng)度高和比剛度大等優(yōu)點(diǎn)，廣泛應(yīng)用于航天和航空工業(yè)領(lǐng)域。眾所周知，對(duì)于很多復(fù)合材料的宏觀(guān)解，如低階頻率和模態(tài)，可以使用等應(yīng)變模型或等應(yīng)力模型[1]以及其他均勻化方法[2]，但相對(duì)于宏觀(guān)應(yīng)力分析，細(xì)觀(guān)分析要復(fù)雜很多。為了在計(jì)算精度和效率之間達(dá)到平衡，各種多尺度方法相繼提出，如數(shù)學(xué)均勻化方法MHM[3,4]、廣義有限單元方法GFEM[5,6]、多尺度有限單元方法MsFEM[7,8]、異類(lèi)多尺度方法HMM[9,10]以及多尺度特征單元方法MEM[11,12]，在這些方法中，數(shù)學(xué)均勻化方法是具有代表性的方法，文獻(xiàn)[13-21]都作了詳細(xì)的闡述；文獻(xiàn)[22-28]討論了數(shù)學(xué)均勻化方法的攝動(dòng)階次和使用的單元階次對(duì)計(jì)算精度的影響，在單胞上施加四邊固支邊界條件求解其影響函數(shù)控制方程，結(jié)果顯示對(duì)于一維結(jié)構(gòu)是精確的，但對(duì)于二維結(jié)構(gòu)計(jì)算精度不盡人意。

本文在文獻(xiàn)[22,24]的基礎(chǔ)上詳細(xì)研究了單胞影響函數(shù)控制方程在周期復(fù)合材料結(jié)構(gòu)中的真實(shí)邊界條件，并以此為標(biāo)準(zhǔn)檢驗(yàn)數(shù)值計(jì)算中施加的周期邊界條件計(jì)算精度，提出更為合理的二維結(jié)構(gòu)單胞周期邊界條件。

2 數(shù)學(xué)均勻化方法矩陣列式

基于微觀(guān)周期性和單胞域的一致性假設(shè)，均勻化理論將異質(zhì)邊界值問(wèn)題分解為單胞(微觀(guān))問(wèn)題和結(jié)構(gòu)問(wèn)題(宏觀(guān))，二維周期性復(fù)合材料結(jié)構(gòu)的控制方程為橢圓型偏微分方程：

(1)

(2)

文獻(xiàn)[22,24]強(qiáng)調(diào)了二階攝動(dòng)對(duì)數(shù)學(xué)均勻化方法計(jì)算精度的重要性，所以本文同樣考慮了二階攝動(dòng)項(xiàng)的計(jì)算精度。

表1給出了一階和二階攝動(dòng)位移的解耦形式以及影響函數(shù)的控制方程，其中EH為均勻化彈性模量，控制方程中的算子yj·和yn分別為散度和梯度。影響函數(shù)是單胞尺度(微觀(guān)尺度)y的函數(shù)，而均勻化位移及其各階導(dǎo)數(shù)只依賴(lài)于宏觀(guān)尺度x。從解耦形式和影響函數(shù)控制方程可以得到很多信息，詳見(jiàn)文獻(xiàn)[24]。

表1 攝動(dòng)位移和影響函數(shù)控制方程Tab.1 Perturbed displacements and governing equations of influence functions

數(shù)學(xué)均勻化方法一般是通過(guò)有限元方法得以實(shí)現(xiàn)，所以需要將包含二階攝動(dòng)項(xiàng)的數(shù)學(xué)均勻化方法的數(shù)學(xué)表達(dá)式通過(guò)加權(quán)殘量方法轉(zhuǎn)化為易于編程計(jì)算的矩陣列式，解耦的單胞影響函數(shù)控制方程矩陣列式為

(3)

(4，5)

(6)

宏觀(guān)均勻化位移控制方程為

KHu0=F0

(7)

(8)

式中D?Ω表示周期復(fù)合材料結(jié)構(gòu)特定單胞區(qū)域，g為結(jié)構(gòu)中單胞的數(shù)量。

求解結(jié)構(gòu)宏觀(guān)位移導(dǎo)數(shù)和各階影響函數(shù)之后，宏觀(guān)結(jié)構(gòu)單元節(jié)點(diǎn)的真實(shí)位移計(jì)算公式為

(9)

值得注意的是，對(duì)于周期復(fù)合材料結(jié)構(gòu)，虛擬載荷和虛擬位移是以單胞為基本單位，且具有周期性，所以可以直接通過(guò)單胞結(jié)構(gòu)計(jì)算其虛擬位移，然后復(fù)制到整個(gè)結(jié)構(gòu)上，從而得到宏觀(guān)結(jié)構(gòu)的虛擬載荷和虛擬位移，進(jìn)而計(jì)算真實(shí)位移，此舉主要是減少計(jì)算量以提高計(jì)算效率，這也是數(shù)學(xué)均勻化方法在實(shí)際工程中得以廣泛應(yīng)用的關(guān)鍵優(yōu)勢(shì)所在。另一方面，為了求解方程(3)得到虛擬位移，需要施加周期性邊界條件，所以虛擬位移的計(jì)算精度很大程度上取決于周期邊界條件選擇的合理性，最理想的情況是施加的周期邊界條件能夠反映單胞在結(jié)構(gòu)中的真實(shí)邊界變形情況，即可不因邊界條件的施加而引入誤差。

3 二維周期復(fù)合材料結(jié)構(gòu)

為了更加詳細(xì)地研究周期復(fù)合材料單胞邊界條件問(wèn)題，考慮一個(gè)四邊固支的矩陣平面周期復(fù)合材料結(jié)構(gòu)，其結(jié)構(gòu)尺寸為45 mm×45 mm，如圖1所示，每個(gè)單胞尺寸為9 mm×9 mm，內(nèi)部包含一個(gè)矩陣夾雜，基體和夾雜的彈性模量分別為E1=2.97 GPa和E2=90.585 GPa，泊松比均為0.33。

圖1 二維周期復(fù)合材料及其單胞結(jié)構(gòu)

3.1 單胞真實(shí)位移邊界條件

為了施加更為合理的單胞周期邊界條件，首先必須弄清求解影響函數(shù)控制方程(3)時(shí)，單胞在宏觀(guān)結(jié)構(gòu)中真實(shí)的邊界情況，有效的方法是將整個(gè)宏觀(guān)結(jié)構(gòu)作為一個(gè)大單胞來(lái)處理，即可避免額外施加單胞邊界條件，直接使用結(jié)構(gòu)邊界條件求解虛擬位移(方法1)，這樣得到的虛擬位移在結(jié)構(gòu)內(nèi)部每個(gè)單胞邊界上的情況即為單胞真實(shí)位移邊界條件。

一階虛擬載荷F1由材料彈性參數(shù)矩陣決定，與單胞邊界條件無(wú)關(guān)，如方程(6)所示，F(xiàn)1矩陣的 3列 作為3個(gè)相互獨(dú)立的虛擬載荷作用在結(jié)構(gòu)上，其矢量如圖2所示，分別對(duì)應(yīng)F1矩陣的第1至第3列。可以看出，與F1相關(guān)的 3個(gè) 獨(dú)立虛擬載荷只有在基體和夾雜的交界線(xiàn)上具有非0的分段線(xiàn)性載荷。

圖2 宏觀(guān)結(jié)構(gòu)F1矩陣矢量

(1) 宏觀(guān)結(jié)構(gòu)內(nèi)部單胞的真實(shí)邊界條件具有周期性。

(2) 宏觀(guān)結(jié)構(gòu)內(nèi)部單胞邊界的虛擬位移是非0的，圖3(c)表現(xiàn)尤為明顯。

圖3 宏觀(guān)結(jié)構(gòu)真實(shí)邊界條件一階虛擬位1變形圖

(3) 宏觀(guān)結(jié)構(gòu)單胞內(nèi)部虛擬位移和邊界虛擬位移相差較大。

圖4 宏觀(guān)結(jié)構(gòu)F2矩陣矢量

(1) 宏觀(guān)結(jié)構(gòu)內(nèi)部單胞的真實(shí)邊界條件仍然具有周期性。

(2) 宏觀(guān)結(jié)構(gòu)內(nèi)部單胞邊界二階虛擬位移相對(duì)一階虛擬位移變形更為明顯。

3.2 單胞影響函數(shù)控制方程邊界條件

3.2.1 四邊固支邊界條件

盡管四邊固支邊界條件廣泛應(yīng)用于數(shù)學(xué)均勻化方法處理單胞邊界問(wèn)題，但其對(duì)各階影響函數(shù)求解精度的影響尚未得到充分的研究。為了弄清該邊界條件的合理性，從結(jié)構(gòu)中取出一個(gè)單胞施加四邊固支周期邊界條件分別計(jì)算一階和二階影響函數(shù)，然后利用單胞結(jié)構(gòu)的周期性復(fù)制到整個(gè)結(jié)構(gòu)中去，得到結(jié)構(gòu)的虛擬變形圖(方法2)，與3.1節(jié)(方法1)得到的虛擬位移結(jié)構(gòu)變形圖進(jìn)行對(duì)比。

結(jié)構(gòu)的一階和二階虛擬位移結(jié)構(gòu)變形圖分別如圖6和圖7所示，分別與3.1節(jié)的圖3和圖5進(jìn)行對(duì)比可以得到，

圖6 使用單胞四邊固支周期邊界條件的一階虛擬位移1結(jié)構(gòu)變形圖

圖7 使用單胞四邊固支周期邊界條件的二階虛擬位移2結(jié)構(gòu)變形圖

(1) 對(duì)于一階虛擬位移，方法1和方法2的虛擬位移結(jié)構(gòu)變形圖的差別主要集中在單胞邊界上，但并不明顯。四邊固支邊界條件決定了單胞邊界上節(jié)點(diǎn)虛擬位移為0，而真實(shí)單胞邊界條件中邊界虛擬位移較小，所以使用四邊固支邊界條件求解一階虛擬位移是合理的。

(2) 對(duì)于二階虛擬位移，方法1和方法2得到的虛擬位移結(jié)構(gòu)變形圖的差別十分明顯，同時(shí)表現(xiàn)在單胞邊界上的節(jié)點(diǎn)和單胞內(nèi)部的節(jié)點(diǎn)，所以使用四邊固支邊界條件求解二階虛擬位移是不合理的，會(huì)造成很大的計(jì)算誤差。

3.2.2 超單胞周期邊界條件

從3.2.1節(jié)可以看出，盡管使用四邊固支周期邊界條件求解單胞虛擬位移較為方便，但對(duì)于二階虛擬位移的求解精度卻難以接受。實(shí)際上，單胞作為結(jié)構(gòu)基本周期單元包含了所有的微觀(guān)信息，而且單胞的虛擬位移在邊界上取決于其在結(jié)構(gòu)內(nèi)部周?chē)鷨伟南嗷プ饔?，由此提出超單胞周期邊界條件(方法3)，即在9個(gè)單胞(3×3)組成的一個(gè)超單胞上施加四邊固支邊界條件，如圖8(a)所示，具體實(shí)現(xiàn)過(guò)程如下。

(1) 使用公式(3)計(jì)算超單胞的一階和二階虛擬位移。

(2) 將超單胞中9個(gè)單胞根據(jù)其在超單胞中的位置分為9種單胞，分別為A，B，C，D，E，F(xiàn)，G，H和I，如圖8(b)所示，每個(gè)單胞均有獨(dú)立的一階和二階虛擬位移。

(3) 將單胞E的一階和二階虛擬位移復(fù)制到結(jié)構(gòu)內(nèi)部的每一個(gè)單胞上，A，B，C，D，F(xiàn)，G，H和I共8個(gè)單胞的一階和二階虛擬位移根據(jù)其在超單胞中的位置分別對(duì)應(yīng)復(fù)制到結(jié)構(gòu)邊界上的單胞，如圖8(c)所示。

圖8 超單胞邊界條件實(shí)現(xiàn)過(guò)程

由超單胞周期邊界條件計(jì)算得到的一階和二階虛擬位移結(jié)構(gòu)變形分別如圖9和圖10所示，將其分別與圖3和圖5對(duì)比，可以得到，

圖9 使用超單胞周期邊界條件得到的一階虛擬位移1結(jié)構(gòu)變形

圖10 使用超單胞周期邊界條件得到的二階虛擬位移2結(jié)構(gòu)變形

(1) 使用超單胞周期邊界條件得到的虛擬位移結(jié)構(gòu)變形圖和真實(shí)邊界條件得到的虛擬位移結(jié)構(gòu)變形圖差異較小。

(2) 超單胞周期邊界條件與真實(shí)單胞邊界條件相似，可以有效處理數(shù)學(xué)均勻化方法單胞問(wèn)題。

3.3 虛擬勢(shì)能泛函

為了從數(shù)值上直觀(guān)反映超單胞周期邊界條件的計(jì)算精度，提出了與虛擬位移或虛擬載荷相對(duì)應(yīng)的虛擬勢(shì)能泛函Π的概念，其計(jì)算公式如下，

Πk l p=Uk l p+Wk l p

(10)

表2和表3分別給出一階和二階虛擬勢(shì)能泛函的數(shù)值，其中ΠFEM表示使用真實(shí)單胞邊界條件(方法1)計(jì)算得到的虛擬勢(shì)能泛函，在此作為標(biāo)準(zhǔn)檢驗(yàn)四邊固支周期邊界條件和超單胞周期邊界條件的計(jì)算精度，Π1為使用四邊固支周期邊界條件(方法2)計(jì)算得到的虛擬勢(shì)能泛函，Π2為使用超單胞周期邊界條件(方法3)計(jì)算得到的虛擬勢(shì)能泛函，從表2和表3的數(shù)值結(jié)果可以得到，

表2 對(duì)應(yīng)一階虛擬位移的虛擬勢(shì)能泛函(107J)

表3 對(duì)應(yīng)二階虛擬位移/載荷的虛擬勢(shì)能泛函(107J)

(1) 二階虛擬位移計(jì)算誤差要遠(yuǎn)大于一階虛擬位移，即隨著虛擬位移階次的提高，無(wú)論施加何種單胞周期邊界條件，計(jì)算誤差都隨之增大。

(2) 使用四邊固支邊界條件計(jì)算一階虛擬位移的誤差較小，而二階虛擬位移計(jì)算誤差過(guò)大。

(3) 使用超單胞周期邊界計(jì)算虛擬位移的誤差大幅降低，說(shuō)明超單胞周期邊界條件可以有效提高數(shù)學(xué)均勻化方法高階虛擬位移的計(jì)算精度。

3.4 數(shù)值算例分析

為了驗(yàn)證上文分析結(jié)果的正確性，使用雙線(xiàn)性平面單元，在圖1結(jié)構(gòu)中沿坐標(biāo)x1方向施加面載荷f=0.1 MPa，每個(gè)單胞分為9×9個(gè)單元。為了方便研究單胞邊界條件對(duì)數(shù)學(xué)均勻化方法計(jì)算精度的影響，均勻化宏觀(guān)結(jié)構(gòu)單元尺寸與單胞結(jié)構(gòu)單元尺寸相同，均為1 mm×1 mm。

對(duì)四邊固支的二維周期復(fù)合材料結(jié)構(gòu)分別使用上文中的三種邊界條件進(jìn)行計(jì)算。

(1) 將整個(gè)結(jié)構(gòu)作為一個(gè)大單胞進(jìn)行處理，此時(shí)的單胞邊界條件與結(jié)構(gòu)邊界條件相同，小參數(shù)ε=1，沒(méi)有引入周期邊界條件所帶來(lái)的誤差。

(2) 對(duì)單胞施加四邊固支周期邊界條件求解虛擬位移，值得注意的是，此時(shí)單胞單元的尺寸擴(kuò)大5倍為5 mm×5 mm，小參數(shù)ε=1/5。

(3) 使用超單胞周期邊界條件求解虛擬位移，此時(shí)超單胞內(nèi)的單元尺寸仍然擴(kuò)大5倍為5 mm×5 mm，小參數(shù)ε=1/5。

使用這三種方法計(jì)算分別得到圖1中4條縱線(xiàn)A(穿過(guò)夾雜內(nèi)部)、B(沿著單胞邊界)、C(沿著基體和夾雜交界線(xiàn))和D(穿過(guò)基體)上的節(jié)點(diǎn)真實(shí)位移曲線(xiàn)，結(jié)果列入表4，其中FEM為有限元細(xì)網(wǎng)格計(jì)算結(jié)果，作為標(biāo)準(zhǔn)檢驗(yàn)三種單胞邊界條件的計(jì)算精度，MsAEM1表示包含一階攝動(dòng)項(xiàng)的節(jié)點(diǎn)真實(shí)位移，MsAEM2表示包含二階攝動(dòng)項(xiàng)的節(jié)點(diǎn)真實(shí)位移，表5給出了勢(shì)能泛函的計(jì)算結(jié)果，計(jì)算公式為

(11)

由表4和表5的數(shù)值計(jì)算結(jié)果對(duì)比可以得到,

表4 沿著縱線(xiàn)A，B，C和D節(jié)點(diǎn)的真實(shí)位移曲線(xiàn)Tab.4 Actual displacements comparisons of line A,B,C and D

(1) 使用方法2，即四邊固支邊界條件得到的節(jié)點(diǎn)宏觀(guān)位移精度較差，根本原因在于四邊固支導(dǎo)致單胞邊界上的節(jié)點(diǎn)虛擬位移為0，即攝動(dòng)位移為0，此時(shí)的節(jié)點(diǎn)宏觀(guān)位移就是均勻化位移；MsAEM2的精度在單胞內(nèi)部較低，主要是由于二階虛擬位移計(jì)算精度較差。

(續(xù)表)

(2) 使用方法1和方法3得到的包含二階攝動(dòng)位移的節(jié)點(diǎn)宏觀(guān)位移計(jì)算精度較高，也驗(yàn)證了文獻(xiàn)[22]二階攝動(dòng)項(xiàng)必要性的結(jié)論。

(3) 方法1和方法3的計(jì)算結(jié)果十分相近，說(shuō)明超單胞周期邊界條件很好地反映了單胞在結(jié)構(gòu)中的真實(shí)邊界條件，可以有效地處理數(shù)學(xué)均勻化方法中的單胞邊界問(wèn)題。

4 結(jié)論

詳細(xì)研究了數(shù)學(xué)均勻化方法單胞邊界條件計(jì)算精度的問(wèn)題，將影響函數(shù)作為虛擬位移得到了周期結(jié)構(gòu)內(nèi)單胞的真實(shí)邊界條件，并提出了與之相應(yīng)的虛擬載荷和虛擬勢(shì)能泛函的概念，得到結(jié)論如下。

(1) 固支周期邊界條件不適合作為二維結(jié)構(gòu)單胞邊界條件。

(2) 提出超單胞周期邊界條件，大幅提高了二維周期復(fù)合材料結(jié)構(gòu)數(shù)學(xué)均勻化方法的計(jì)算精度。

(3) 使用虛擬勢(shì)能泛函驗(yàn)證了不同邊界條件和攝動(dòng)階次對(duì)結(jié)構(gòu)宏觀(guān)位移計(jì)算精度的影響，進(jìn)一步確定了二階攝動(dòng)項(xiàng)對(duì)數(shù)學(xué)均勻化方法計(jì)算精度的必要性。