劉艷紅，張 瑩

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

傳統(tǒng)混合有限元法（簡稱傳統(tǒng)混合法）[1-5]的主要優(yōu)點是：一方面，可同時求得應(yīng)力和位移結(jié)果，而且應(yīng)力結(jié)果的精度較高，一些算例表明其應(yīng)力結(jié)果可接近精確解[1-2]；另一方面，不同單元同一結(jié)點的應(yīng)力結(jié)果是連續(xù)的。但一般情況下，傳統(tǒng)混合元系數(shù)矩陣的主對角線上有0元素。因此，不經(jīng)特殊處理的傳統(tǒng)混合元的數(shù)值結(jié)果不穩(wěn)定。傳統(tǒng)混合法中各種非常規(guī)的穩(wěn)定元素技術(shù)方案所涉及的數(shù)學(xué)理論較復(fù)雜，不利于工程師掌握和應(yīng)用[4]。另外，基于傳統(tǒng)混合法的模型主要涉及具體的二維板殼類問題。有關(guān)三維工程問題的文獻[5]表明，相應(yīng)的穩(wěn)定元素技術(shù)的理論過程更為復(fù)雜，因而很少有針對具體三維工程問題進行數(shù)值分析的文獻。

無需任何穩(wěn)定混合元素技術(shù)，基于兩個變分原理（即最小勢能原理和H-R變分原理）或廣義混合變分原理可得到本質(zhì)上穩(wěn)定的廣義混合元[6-7]，即系數(shù)矩陣主對角線上沒有0元素。更進一步，文獻[8]根據(jù)修正的H-R變分原理直接建立了廣義的非協(xié)調(diào)辛元。

由于所采用的變分原理不同，廣義混合元[6-8]的形式當(dāng)然也不同。多個實例的數(shù)值分析表明：廣義混合法[6-8]和被廣泛應(yīng)用的位移法或雜交應(yīng)力法一樣，用于分析線彈性靜力學(xué)問題沒有原則上的困難。不必諱言，廣義混合法和傳統(tǒng)混合法一樣，在有限元控制方程中增加了應(yīng)力變量。但正因如此，不需要從單元級依據(jù)材料的本構(gòu)關(guān)系求解應(yīng)力結(jié)果，可避免應(yīng)力數(shù)值結(jié)果不準(zhǔn)確的問題。非協(xié)調(diào)廣義混合元的位移結(jié)果和應(yīng)力結(jié)果收斂速度快。位移和應(yīng)力變量的數(shù)值結(jié)果精確度高是文獻[8]中非協(xié)調(diào)辛元最突出的特點。

超級單元的形式多種多樣，構(gòu)造方法種類繁多，如模態(tài)綜合子結(jié)構(gòu)法就是一種很常見的超級單元方法。在桁架或鋼架結(jié)構(gòu)分析中，人們常用“組合單元”這個名稱來表述超級單元。一般情況下，超級單元法或組合單元法可統(tǒng)稱為子結(jié)構(gòu)法。鐘萬勰[9]提出的區(qū)段方法也可理解為子結(jié)構(gòu)法。

研究以非協(xié)調(diào)辛元為基礎(chǔ)，針對復(fù)合材料層合板的結(jié)構(gòu)特點，給出了構(gòu)建復(fù)合材料層合板超級混合單元的基本過程和步驟，并用實例證明其可靠性。

1 理論公式

1.1 修正的H-R變分原理與廣義部分混合元列式

為了清楚地表示復(fù)合材料層合板超級部分混合元的基本理論，首先簡要地給出非協(xié)調(diào)廣義部分混合元的基本步驟。

式（2）中：σo=[σxzσyzσzz]T為平面外應(yīng)力；Φ11、Φ21和Φ22是與材料參數(shù)相關(guān)的子矩陣；Δ

1、Δ

2和Δ

3是微分算子矩陣，即

平面內(nèi)應(yīng)力可由位移變量和平面外應(yīng)力變量表示為

根據(jù)非協(xié)調(diào)位移元[13-14]理論，可將位移插值表達式分為協(xié)調(diào)部分和非協(xié)調(diào)部分，即

其中：diag[N]=[NeNeNe]，Ne為對應(yīng)8結(jié)點六面體元等參形函數(shù)；位移參數(shù)向量qe=[uevewe]T；diag[Nr]=[NreNreNre]為對應(yīng)單元內(nèi)非協(xié)調(diào)位移參數(shù)向量re的形函數(shù)。

平面外應(yīng)力σo變量可與位移u變量一樣在相同的空間域內(nèi)被離散，設(shè)其表達式為

平面外應(yīng)力的結(jié)點參數(shù)向量pe=[σexzσeyzσezz]T。

根據(jù)文獻[8]修正的H-R變分原理（1）離散泛函的變分結(jié)果，可得如下非協(xié)調(diào)廣義部分混合元

1.2 廣義部分混合元的辛對稱性質(zhì)

事實上，式（6）可改寫為如下形式

為了方便，不妨設(shè)

根據(jù)辛代數(shù)的理論，設(shè)J為單位辛陣，即

式中：I24×24為一般意義下的單位陣，其下標(biāo)表示維數(shù)。

不難證明以下關(guān)系

從式（7）可看出，其系數(shù)矩陣本身是不對稱的，而在辛代數(shù)中，其系數(shù)矩陣辛對稱。

1.3 超級混合元

假設(shè)層合板結(jié)構(gòu)如圖1（a）所示，其第i層中結(jié)點的排列如圖1（b）所示，則可方便地分離式（7）中上下表面的結(jié)點，因而得到如下表達式

圖1 層合板結(jié)構(gòu)和第i層的結(jié)點排列Fig.1 Laminated plate structure and joint array on Level i

對式（11）交換第2行（列）和第3行（列）有

上式可簡寫為

式中

對式（13）求和，可得第i層的控制方程為

假設(shè)考慮如圖1（a）所示的4層結(jié)構(gòu)，根據(jù)層合板中各層間上下表面的平面外應(yīng)力和位移的關(guān)系（即層間上下表面的平面外應(yīng)力和位移相等的關(guān)系），對式（14）求和，可得總體結(jié)構(gòu)的控制方程為

根據(jù)文獻[6-8]中的求解方法，自然也可對控制方程（15）同時引入平面外應(yīng)力和3個位移變量的邊界條件。求得位移和平面外應(yīng)力結(jié)果后，平面內(nèi)應(yīng)力的求解方法與文獻[6-8]相同。

2 數(shù)值算例

考慮如圖 2（a）所示的 2 層復(fù)合材料板[0°/90°]，已知板的長和寬a=b=1.0，厚度h=2h1=0.1。材料參數(shù)[15]為：E11=10E22=10E33，G12=G13=0.6E33，G23=0.5E33，μ12=μ13=μ23=0.25。板上表面（z=0.1處）受垂直向下的均布載荷q=0.1作用。四邊簡支（SSSS）：在x=0和x=a處，σx=w=v=0；在y=0和y=b處，σy=w=u=0。

圖2 復(fù)合材料板、坐標(biāo)系與層數(shù)劃分Fig.2 Composite laminates,coordinate system and layer division

厚度方向的層數(shù)劃分如圖2（b）所示。上下2層中各4薄層，每層厚度s1=0.01125h；接近中性面有上下2個薄層，每層厚度s2=0.005h。所以，總體上板被劃分為10層?？紤]到結(jié)構(gòu)的對稱性，有限元的數(shù)值分析采用1/4模型。網(wǎng)格劃分為12×12×10，分別采用了廣義部分混合元及超級混合元方法計算應(yīng)力及位移，兩種算法所消耗的時間如表1所示。

表1 兩種方法對比分析Tab.1 Comparative analysis between two methods

因為每層4個薄層的單元及材料相同，另一薄層的單元及材料也相同，因此，上下2層的4個薄層分別只需分析1層，并考慮接近中性面的2個薄層，所以超級混合單元法僅對576個單元進行數(shù)值積分。就此實例而言，從表1中可看出，超級混合單元法大約可節(jié)省40%的時間?？梢姴捎贸壔旌显治鰧雍辖Y(jié)構(gòu)靜力學(xué)問題的成效明顯。

圖3～圖8為SME（代表超級混合元）中位移變量和部分應(yīng)力變量與精確解（EXACT）[15]的對比結(jié)果，其中最大誤差發(fā)生于σxz（a/2，b/8，z）沿厚度的分布這一變量上，如圖5所示，不超過0.8%。由于材料不同，圖7表明平面內(nèi)應(yīng)力σxx在2層材料的界面之間不連續(xù)，由此可見所研究的復(fù)合材料層合板的超級混合元同樣能反映這一客觀特性。圖3～圖8中，x3指z方向，u1指位移 u，u3指位移 w，數(shù)字 1、2、3 分別指 x、y、z方向。

圖3 位移u（a/8，b/2，z）沿厚度方向的分布Fig.3 Distribution of displacement u（a/8，b/2，z）along thickness direction

圖4 位移w（a/2，b/2，z）沿厚度方向的分布Fig.4 Distribution of displacement w（a/2，b/2，z）along thickness direction

圖5 應(yīng)力σxz（a/2，b/8，z）沿厚度方向的分布Fig.5Distribution of stress σxz（a/2，b/8，z）along thickness direction

3 結(jié)語

圖6 應(yīng)力σzz（a/2，b/8，z）沿厚度方向的分布Fig.6Distribution of stress σzz（a/2，b/8，z）along thickness direction

圖7 應(yīng)力σxx（a/2，b/2，z）沿厚度方向的分布Fig.7Distribution of stress σxx（a/2，b/2，z）along thickness direction

圖8 應(yīng)力σxy（a/8，b/8，z）沿厚度方向的分布Fig.8Distribution of stress σxy（a/8，b/8，z）along thickness direction

考慮到復(fù)合材料層合板的特點，提出了基于非協(xié)調(diào)辛元的超級混合單元。超級混合單元的組裝僅涉及矩陣的加法，因而運算快捷。算例表明，用超級混合元分析層合結(jié)構(gòu)靜力學(xué)問題的效果明顯。但必需說明：對于不同的問題，節(jié)省時間不僅與結(jié)構(gòu)本身的形式相關(guān)，還與所劃分的相同層數(shù)多少有關(guān)。對于復(fù)雜結(jié)構(gòu)，可能沒有簡單結(jié)構(gòu)的效果明顯。但對于大型的復(fù)雜結(jié)構(gòu)而言，計算模型所用時間較長，只要網(wǎng)格劃分得當(dāng)，速度提高的效果也可能非常有意義。

如果要獲得壓電材料類復(fù)雜層合結(jié)構(gòu)靜力問題的高精度結(jié)果，那么廣義部分混合法是可選的方法之一，而所提出的超級混合單元法同樣可用于這類問題的數(shù)值分析。