盧 翔，屈鶴鵬，劉彥波

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

隨著復(fù)合材料在航空航天、汽車工業(yè)、交通運(yùn)輸、化工紡織等領(lǐng)域得到越來越多的應(yīng)用，復(fù)合材料層合板的熱彈性問題逐漸成為材料界和力學(xué)界的研究重點(diǎn)。研究復(fù)合材料層合板的理論主要有以下幾種：經(jīng)典層合板理論、一階剪切變形理論、高階剪切變形理論、高階層合板理論、分層理論以及三維彈性理論[1-3]。上述層合板理論沒有考慮復(fù)合材料承載能力遠(yuǎn)遠(yuǎn)大于普通材料的特點(diǎn)，在研究復(fù)合材料層合板問題時均存在不足，不能兼顧精確性和方便性。楊正林等[4]以Hamilton 正則方程理論為基礎(chǔ)，給出了層合板熱彈性問題體的半解析法；鄒貴平[5]在一階層合板理論的基礎(chǔ)上，采用共軛辛正交的方法給出了層合板熱應(yīng)力分析的精確解法；隨后，相關(guān)學(xué)者通過增維積分法建立了熱彈性體的齊次狀態(tài)方程，并給出了相應(yīng)的精確解法[6-9]。傳統(tǒng)的半解析法在求解層合板位移和應(yīng)力時有一定的震蕩性，且隨著網(wǎng)格密度增大，高維指數(shù)運(yùn)算需要耗費(fèi)大量機(jī)時，不適應(yīng)于大規(guī)模有限元問題。

以彈性力學(xué)最小勢能原理和修正的H-R 變分原理為基礎(chǔ)，將熱彈性體的位移單元[10]與Hamilton 半解析混合單元[11]相結(jié)合，構(gòu)建了新的半解析混合單元，提出了新的半解析法，結(jié)合位移法和Hamilton 正則方程半解析法的優(yōu)點(diǎn)，大大簡化了問題的計算。通過實(shí)例計算，研究了梯度溫度載荷作用下層合板的熱彈性問題，通過與不同方法的計算結(jié)果[12]相對比，驗(yàn)證了新方法的可靠性。

1 基本理論推導(dǎo)

1.1 基本理論

彈性體在變溫情況下的應(yīng)變分量包括自由熱膨脹引起的應(yīng)變分量和熱膨脹時彈性體內(nèi)各部分之間相互約束引起的應(yīng)變分量，其與熱應(yīng)力之間同樣服從胡克定律。因此，參照彈性力學(xué)的基本方程，可得各向正交異性復(fù)合材料熱彈性問題的本構(gòu)方程為

其中：σ = [σxxσyyσzzσyzσxzσxy]T為應(yīng)力變量向量；T 為溫度載荷向量；ε=[εxxεyyεzzεyzεxzεxy]T為應(yīng)變向量；S = C-1為熱彈性體的柔度系數(shù)矩陣；α =[α1α2α30 0 0]T為線膨脹系數(shù)矩陣；

為剛度系數(shù)矩陣。

1.2 Hamilton 正則方程和最小勢能原理

1.2.1 Hamilton 正則方程

根據(jù)Hellinger-Reissner 廣義變分原理，層合板熱彈性問題的泛函表達(dá)式為

其中：Su為應(yīng)力邊界；Sσ為位移邊界；Q=[u v w]T為彈性體的位移向量為邊界面上已知的位移向量；p=[pxxpyypzz]T為彈性體受到的表面力；為邊界面上已知的應(yīng)力向量。

根據(jù)彈性力學(xué)的Hamilton 正則方程理論，可得熱彈性體修正后的H-R 變分原理為

其中：p = [ pxzpyzpzz]T為平面外應(yīng)力向量；H 為Hamilton 函數(shù)，其表達(dá)式為

其中：G1和G2為偏微分算子矩陣；Φij為本構(gòu)關(guān)系式（1）經(jīng)轉(zhuǎn)置后應(yīng)力向量與應(yīng)變向量之間的關(guān)系矩陣[6]。

按照文獻(xiàn)[4]中的方法，分別以P 和Q 為變量對Hamilton 函數(shù)進(jìn)行變分，可得熱彈性體的Hamilton 正則方程為

其中：A，B，C 和D 分別為Hamilton 函數(shù)經(jīng)變分后P和Q 的系數(shù)矩陣，即

1.2.2 最小勢能原理

彈性力學(xué)中的傳統(tǒng)有限元位移法是以最小勢能原理[13]為基礎(chǔ)的，不考慮體積力，則熱彈性體最小勢能原理的泛函可表示為

熱彈性體的幾何方程可表示為

其中，X 為熱彈性體應(yīng)變向量ε 與位移量Q 之間的關(guān)系矩陣，即

將式（8）代入式（7），可得用位移變量表示的最小勢能原理為

2 半解析混合單元模型

將文獻(xiàn)[11]中的Hamilton 等參元法運(yùn)用到求解層合板的熱彈性問題中。對于任意形狀的四邊形單元，以平面四節(jié)點(diǎn)Hamilton 等參元為例，將插值后的變量[11]代入式（9）即可得到經(jīng)Hamilton 等參元離散后任意單元的最小勢能原理的泛函表達(dá)式，即

由于δΠP=0，因此以Qe為變量，對式（10）進(jìn)行變分，可以得到

其中

按照同樣方法對熱彈性體的Hamilton 正則方程式（6）進(jìn)行Hamilton 等參元離散，即可得到離散后的Hamilton 正則方程為

其中：N 為插值函數(shù)矩陣；J 為雅各比矩陣式[11]，選取式（12）中的第1 個方程

將其與式（11）合并，可建立熱彈性體經(jīng)Hamilton等參元離散后的半解析混合單元模型（不計體積力），即

3 層合板熱彈性問題混合方程

層合板的任一子層均可用上述方法進(jìn)行離散并建立半解析混合單元。經(jīng)過單元組裝，可以得到任一子層的混合狀態(tài)方程。具體對式（15）求和后可得

其中：m 為層合板第m 子層；Am、Bm和Km分別為經(jīng)單元組裝后的系數(shù)矩陣；f1m、Γ1m和Γ2m分別為f1e、Γ1e和Γ2e進(jìn)行單元求和后的結(jié)果，即

在具體求解時，首先對式（16）第2 個方程進(jìn)行位移求解，然后將所得位移結(jié)果一一代入第1 個方程進(jìn)行平面內(nèi)應(yīng)力的求解。因此，式（16）中第1 個方程就變成更為簡潔的一階微分方程，即

其中：Φm（z）=BmQ（z）m+f（z）1m

式（17）的解為

其中

對于n 層的層合板，根據(jù)層間的連續(xù)關(guān)系，則有

其中：m = 2，3，…，n；hm為層合板第m 子層的厚度；Pm（0）為第m 層下表面的平面外應(yīng)力向量；Pm（hm）為第m 層上表面的平面外應(yīng)力向量。

引入傳遞矩陣并結(jié)合式（18）～式（19），可以得到層合板上下表面應(yīng)力之間的關(guān)系表達(dá)式。然后根據(jù)層合板上下表面的應(yīng)力位移初始值即可求得層合板任意位置處的應(yīng)力和位移。

以上求解方法實(shí)際上就是先求位移，然后求平面外應(yīng)力的微分方程。很明顯，對比式（6）和式（17）的系數(shù)矩陣可以發(fā)現(xiàn)，后者的指數(shù)維數(shù)是前者的1/4。由于傳統(tǒng)位移法在求解位移時是準(zhǔn)確的，因此從理論上講，新半解析法在保證計算精度的同時，可以有效提高運(yùn)算速度，減少運(yùn)算時間。

4 實(shí)例分析

選取1 塊等厚度的矩形薄板，如圖1 所示（a，b 遠(yuǎn)遠(yuǎn)大于厚度h），溫度載荷T 為y 的函數(shù)。無量綱化后的幾何性質(zhì)：a=100，b=20；材料特性：彈性模量E=210，泊松比v=0.3；應(yīng)力溫度系數(shù)：α11=0.000 1，T0=100；板的上下表面和側(cè)面均為應(yīng)力自由表面。分別考慮兩種溫度載荷條件下矩形薄板的熱彈性問題：①T（y）=T0（1-y2/b2）；②T（y）=T0（1-y3/b3）。

圖1 自由邊界條件下的矩形薄板Fig.1 Rectangular plate with free edges

由文獻(xiàn)[13]可知矩形薄板內(nèi)熱應(yīng)力的精確解為

離散過程中，由于網(wǎng)格的疏密程度對計算結(jié)果影響較大，因此只選用32 × 12 網(wǎng)格密度下新半解析解的計算結(jié)果進(jìn)行誤差分析，即

A（0，0）和B（0，-20）兩點(diǎn)在不同溫度載荷條件下的應(yīng)力計算結(jié)果如表1～表4 所示。

表1 和表2 分別給出了施加溫度載荷1 時不同網(wǎng)格密度條件下A（0，0）點(diǎn)和B（0，-20）點(diǎn)熱應(yīng)力σxx的計算結(jié)果。圖2 和圖3 分別是溫度載荷條件1 作用下A（0，0）點(diǎn)和B（0，-20）點(diǎn)的熱應(yīng)力σxx隨網(wǎng)格密度的變化。由圖2 和圖3 可以直觀地看出，相對于ABAQUS解，新半解析解的計算結(jié)果受網(wǎng)格疏密程度的影響更小。通過橫向?qū)Ρ瓤梢园l(fā)現(xiàn)，當(dāng)網(wǎng)格疏密程度一樣時，新半解析解比ABAQUS 解更接近精確解，具有更高的計算精度。即使在網(wǎng)格密度較小的情況下，新半解析解的應(yīng)力計算結(jié)果也和精確解十分接近。

表1 溫度載荷條件1 作用下點(diǎn)A（0，0）處熱應(yīng)力σxx的3 種解Tab.1 σxxat Point A（0，0）under Thermal Case One

表2 溫度載荷條件1 作用下點(diǎn)B（0，-20）處熱應(yīng)力σxx的3 種解Tab.2 σxxat Point B（0，-20）under Thermal Case One

表3 溫度載荷條件2 作用下點(diǎn)A（0，0）處熱應(yīng)力σxx的3 種解Tab.3 σxxat Point A（0，0）under Thermal Case Two

表4 溫度載荷條件2 作用下點(diǎn)B（0，-20）處熱應(yīng)力σxx的3 種解Tab.4 σxxat Point B（0，-20）under Thermal Case Two

表3 和表4 分別給出了施加溫度載荷2 時不同網(wǎng)格密度條件下A（0，0）點(diǎn)和B（0，-20）點(diǎn)熱應(yīng)力σxx的計算結(jié)果。由表3 可以看出，在溫度載荷條件2 作用下，矩形薄板中心處的熱應(yīng)力σxx=0，無論是新半解析解還是ABAQUS 有限元解，都很好地符合了這一結(jié)果。由表4 可以看出，在溫度載荷條件2 作用下，隨著網(wǎng)格密度增大，兩者的計算結(jié)果均越來越接近精確解。不同的是，與ABAQUS 有限元解相比，新半解析解的計算結(jié)果誤差更小、精確度更高，并在較小的網(wǎng)格密度下可保證計算的精確性。網(wǎng)格密度決定了計算機(jī)的工作量，因此新半解析解在保證較高計算精度的同時，可有效縮短計算時間。同時對計算機(jī)的運(yùn)行內(nèi)存和運(yùn)行速度要求不高，適用于大規(guī)模有限元計算。

圖2 A（0，0）點(diǎn)σxx隨網(wǎng)格密度的變化Fig.2 σxxvariation at Point A along with mesh density change

圖3 B（0，-20）點(diǎn)σxx隨網(wǎng)格密度的變化Fig.3 σxxvariation at Point B along with mesh density change

在合適的網(wǎng)格密度下，ABAQUS 解也可得到較為精確的結(jié)果。由于研究溫度梯度對層合板層間熱應(yīng)力影響的文獻(xiàn)較少，因此在沒有精確解做參考的情況下，也可以將新半解析解和ABAQUS 解進(jìn)行對比分析。

5 結(jié)語

新半解析法承襲了傳統(tǒng)位移法和Hamilton 半解析法的優(yōu)點(diǎn)，在保證計算精度的同時，可以大大降低計算量，有效減少計算時間，相對于原有的半解析混合法和傳統(tǒng)有限元法，更適應(yīng)于大規(guī)模有限元問題的分析。分別采用新半解析法和ABAQUS 有限元仿真計算不同溫度載荷條件下矩形薄板的熱彈性參數(shù)，并與精確解進(jìn)行對比，驗(yàn)證了該方法的優(yōu)越性。