楊立洲，卿光輝

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

在工程結(jié)構(gòu)有限元仿真領(lǐng)域，眾多學(xué)者致力于探索高效簡潔且精確可靠的數(shù)值分析模型，例如，F(xiàn)RP(Fiber-reinforced Polymers)層合板的降維處理模型[1]、分層殼單元高效的非線性分析模型[2]和用于變剛度層合板的靜力學(xué)與動力學(xué)分析的先進(jìn)模型等。

20 世紀(jì)90 年代初，HYER 等[3]首先提出變剛度層合板(Variable Stiffness Laminates)的概念。變剛度層合板的纖維鋪放角度沿某一方向不斷變化，這極大增加了變剛度層合板的可設(shè)計性。但是，其應(yīng)力-應(yīng)變關(guān)系復(fù)雜，理論分析與數(shù)值仿真變得困難[4]。

部分學(xué)者采用解析法[5-9]對變剛度層合板結(jié)構(gòu)的力學(xué)問題求解，雖然計算效率高，但需要對模型有較多的簡化。從彈性力學(xué)的基本方程出發(fā)，附加若干假設(shè)(如KIRCHHOFF-LOVE 假設(shè)、KIRCHHOFF 假設(shè)等)的板殼理論對變剛度層合板的分析較為多見。如GUPTA 等[10]基于改進(jìn)的三階剪切變形理論，建立了用于變剛度層合板殼線性和非線性靜力分析的有限元公式，研究了曲率比、層數(shù)和邊界條件等因素對薄殼結(jié)構(gòu)非線性靜力性能的影響。譚萍等[11]基于板殼理論，采用加權(quán)殘值法對纖維角度沿徑向變化的變剛度圓環(huán)板的撓度變化進(jìn)行相關(guān)研究。陳曉東等[12]基于經(jīng)典層合理論，研究了不同邊界條件下纖維角度變化時變剛度層合板的撓度及平面內(nèi)節(jié)點位移受到的影響。板殼理論中的經(jīng)典層合理論和高階剪切變形理論，都假定板的撓度與z坐標(biāo)無關(guān)，忽略了擠壓變形的影響。由于中厚板剪切剛度小，存在較大的剪切變形，會導(dǎo)致經(jīng)典層合理論的預(yù)測結(jié)果精度較低。當(dāng)然逐層理論可以應(yīng)用于中厚板分析，但基于該理論的模型變量個數(shù)依賴層合板的層數(shù)[13]。

YAN 等[14]基于改進(jìn)的分層Legendre 展開式和微分求積有限元法，采用Carrera Unified Formulation模型，提出了一種新的統(tǒng)一的準(zhǔn)三維解，對變剛度層合板進(jìn)行振動分析。HASIM 等[15]基于精化鋸齒理論(Refined Zig-zag Theory, RZT)建立了一種新的等幾何公式，并應(yīng)用于變剛度層合板和夾芯板的靜力分析。該方法將RZT 與非均勻有理B 樣條公式融合，避免了單元模型出現(xiàn)剪切自鎖等問題。

廣義混合元法[16-17]同時引入了應(yīng)力邊界條件，使得應(yīng)力精度明顯提高。避免了經(jīng)典混合元法的整體剛度矩陣包含零子矩陣的問題，提高了計算過程的穩(wěn)定性。QING 等[18-19]建立了非協(xié)調(diào)廣義混合元，劉艷紅等[20]將該方法擴(kuò)展至熱彈性復(fù)合材料，建立了一種含參數(shù)的非協(xié)調(diào)辛元。王燮等[21]建立了層合板的非協(xié)調(diào)廣義部分混合元模型，分析了層合板的自由邊界效應(yīng)問題。王聿航等[22]擴(kuò)展了非協(xié)調(diào)廣義部分混合元，用于壓電材料靜力學(xué)分析，使用較少的單元數(shù)便可獲得了和精確解一致性很高的結(jié)果。

考慮到解析法、板殼理論以及位移有限元法的局限性和變剛度層合板的特點，本文為變剛度層合板的靜力學(xué)分析提出一種修正的非協(xié)調(diào)廣義部分混合元模型。

1 非協(xié)調(diào)廣義部分混合元列式

1.1 常規(guī)的非協(xié)調(diào)廣義部分混合元列式

含參數(shù)的廣義混合變分原理：

式中：V、S分別為層合板模型的體積和表面積；σ =[σxzσyzσzσxσyσxy]T為應(yīng)力向量；u=[ux uy uz]T為位移向量；C為材料的剛度矩陣；?為微分操作算子矩陣；T=[TxTyTz]T為給定的應(yīng)力邊界條件。式(1)中參數(shù)α 的取值范圍為0～1，在后文的算例分析中α 取0.75[23]。

將位移和應(yīng)力分別表示為：

式中：Nd和Ns分別為位移和應(yīng)力的形函數(shù)矩陣；Nr為非協(xié)調(diào)項的形函數(shù)矩陣；qe和pe分別為單元節(jié)點的位移向量和應(yīng)力向量；re為單元內(nèi)部節(jié)點的非協(xié)調(diào)位移向量。

將式(2)和式(3)代入式(1)后進(jìn)行變分處理有：

式中：

其中：

采用文獻(xiàn)[21]中方法消去式(4)的第一個方程中的面內(nèi)應(yīng)力向量 [σxσyσxy]T后，可得到只包含面外應(yīng)力和位移的部分混合列式：

式中的相關(guān)符號表達(dá)式見參考文獻(xiàn)[21]。

1.2 修正的非協(xié)調(diào)廣義部分混合元列式

對傳統(tǒng)定向纖維層合板的離散處理，因為其單層板內(nèi)纖維角度不變，可以直接將單層板的材料性質(zhì)賦予每個單元，即每個單元具有統(tǒng)一的材料性質(zhì)。因此對于傳統(tǒng)的層合板結(jié)構(gòu)，單元角度變化的偏軸應(yīng)力-應(yīng)變關(guān)系為[24]：

式中： σ′為單層材料主方向坐標(biāo)系下的應(yīng)力分量；T為坐標(biāo)轉(zhuǎn)換矩陣。

而變剛度層合板與傳統(tǒng)定向纖維層合板不同，如圖1(a)中的模型以及纖維路徑函數(shù)所示，變剛度層合板的鋪層為連續(xù)曲線纖維，即纖維角度θ 連續(xù)變化，且θ 是關(guān)于x的函數(shù)。對變剛度層合板作如圖1(b)所示離散處理后，變剛度層合板的偏軸應(yīng)力-應(yīng)變關(guān)系式如下：

圖1 變剛度層合板及有限元模型Fig.1 Variable stiffness laminate and its finite element model

其中：

x為幾何模型離散后第n列單元中間處的橫坐標(biāo)，坐標(biāo)原點為板的中心，x軸方向為長邊a的方向，L為總的單元列數(shù)。

考慮變剛度層合板中單元材料參數(shù)剛度的變化情況，在推導(dǎo)列式(4)或式(5)時，按文獻(xiàn)[15，25]中通用且可靠的方法，取單元中間位置的纖維角度作為該單元的纖維角度，如圖1(b)所示。將式(7)中的 (T(θ(x)))-1C((T(θ(x)))-1)T作為單元的材料剛度參數(shù)代入列式(4)或式(5)中。因此，在本文用于數(shù)值分析的Wolfram 語言代碼編寫中，對常規(guī)的非協(xié)調(diào)廣義混合元列式(4)和非協(xié)調(diào)部分廣義混合元列式代碼進(jìn)行了修改，即根據(jù)單元的幾何位置修改單元的材料剛度參數(shù)。

為了區(qū)分常規(guī)的非協(xié)調(diào)廣義部分混合元，本文稱通過以上處理后的有限元列式為修正的非協(xié)調(diào)廣義部分混合元列式。

2 算例分析

例1：如圖2 所示，考慮由兩層等厚石墨-環(huán)氧樹脂材料制成的正方形變剛度層合板[26]，尺寸為a=b=1.0 m，板厚h=0.1 m，鋪層方式為[＜45°/90°＞/＜45°/0°＞]，四邊固支，上表面施加法向載荷Pz=-10 kN/m2。材料屬性如下：E1=137.9 GPa，E2=8.96 GPa，E3=8.96 GPa；G12=7.1 GPa，G13=7.1 GPa，G23=6.21 GPa；μ12=0.3，μ13=0.3，μ23=0.49。

圖2 變剛度層合板模型及坐標(biāo)系 /mFig.2 Variable stiffness laminate and coordinate system

以變剛度層合板的中心為原點，建立如圖2所示坐標(biāo)系。纖維路徑函數(shù)為：

式中：β1和β0分別為層合板邊界處(x=-0.5a)和中心處的曲線切向與x軸的夾角。單層板的路徑曲線標(biāo)記為＜β1/β0＞。

為了評估本文非協(xié)調(diào)廣義部分混合元法的適用性以及準(zhǔn)確性，對比ABAQUS 的參考解，在ABAQUS 中采用二十節(jié)點六面體單元網(wǎng)格計算分析，網(wǎng)格劃分為80×80×20，共128 000 個單元。

本文的方法采用八節(jié)點六面體單元，網(wǎng)格劃分為32×32×20，共20 480 個單元。分析沿板厚方向點Q(-0.25a, -0.25b)(參見圖3(a))處的位移和應(yīng)力。從圖4 可以看出，本文方法得到的沿板厚方向的位移uy和uz結(jié)果與ABAQUS 結(jié)果曲線吻合較好，盡管uz結(jié)果存在差異，但是最大誤差在1%以內(nèi)。由圖5 可以看出，盡管本文方法的低階單元數(shù)量小于ABAQUS 中20 節(jié)點位移單元數(shù)量。但是，本文方法得到的應(yīng)力曲線與ABAQUS得到的應(yīng)力曲線基本吻合。例如，在界面處的σxy結(jié)果與ABAQUS 的相對誤差僅為1.9700%，其它 應(yīng) 力 σy、 σz和 σxz包 括 本 文 未 給 出 應(yīng) 力 σx和σyz與ABAQUS 的相對誤差均小于1.9700%。另一方面，由于本文考慮了應(yīng)力邊界條件，因此在層合板上下表面處的 σxz結(jié)果為0(如圖5(d)所示)，符合真實情況。

圖3 變剛度層合板的鋪層方式Fig.3 Lamination method of the variable stiffness laminate

圖4 四邊固支變剛度層合板位移沿厚度方向的分布Fig.4 Distribution of displacement along the thickness of variable stiffness plate with sides clamped

圖5 四邊固支變剛度層合板應(yīng)力沿厚度方向的分布Fig.5 Distribution of stresses along the thickness of variable stiffness with sides clamped

例2：考慮TORNABENE 和BACCIOCCHI[25]提出的三層變剛度夾芯板的彎曲問題，尺寸為a=b=1.5 m、h=0.1 m，四邊固支，上表面處(z=h/2)受均布載荷qz=-10 kPa 作用。三層夾芯板由下到上每層厚度依次為0.02 m、0.06 m、0.02 m。夾芯板第一、第三層由曲線玻璃纖維環(huán)氧樹脂制成。第二層的芯為泡沫橡膠，視為各向同性材料。相關(guān)材料屬性見表1。以夾芯板的中心為原點，板長方向為x軸方向建立坐標(biāo)系。第一、第三層變剛度板的纖維路徑函數(shù)分別為βlayer1、βlayer2，表達(dá)式如下：

表1 玻璃纖維環(huán)氧樹脂和泡沫橡膠的材料屬性Table 1 Material parameters of glass fiber reinforced epoxy resin and foam rubber

式中：β 為x=-0.5a處纖維切向與x軸的夾角，稱為初始角度，由表2 給定； θ′的表達(dá)式為：

表2 鋪層方式Table 2 Laying methods

考慮如表2 所示的四種不同的纖維初始角度(β1)的鋪層方式(a、b、c、d)。其中β1、β0分別為x=-0.5a、x=0 處的纖維角度。四個案例的網(wǎng)格劃分方式均為28×28×30。選擇P(-0.25a,-0.25b)點沿厚度方向節(jié)點作為研究對象，本文方法的位移結(jié)果與文獻(xiàn)[15]的位移結(jié)果曲線列于圖6，與文獻(xiàn)[25]的應(yīng)力結(jié)果曲線列于圖7。

圖6 變剛度夾芯板沿厚度方向的ux 和uy 的分布Fig.6 ux and uy along the thickness of the variable stiffness sandwich plate

圖7 變剛度夾芯板沿厚度方向的應(yīng)力分布Fig.7 Stresses distribution along the thickness of variable stiffness sandwich plate

從圖6 可以清楚地看出，本文方法的位移曲線與文獻(xiàn)[15]位移曲線變化趨勢吻合較好。同時可以看出，隨著纖維鋪層角度不斷增加，即由案例a～案例d，ux曲線逐漸左移，uy曲線逐漸右移?；诎咐齞 的纖維鋪設(shè)可解釋位移曲線如此變化的原因：案例d 的纖維鋪設(shè)導(dǎo)致夾芯板第一層在x方向和第三層在y方向的剛度顯著降低。故導(dǎo)致此時的夾芯板的ux在第一層最小，uy第三層最大。

圖7 為本文方法沿夾芯板厚度方向的應(yīng)力解與文獻(xiàn)[25]結(jié)果的對比。由圖7(a)和圖7(b)可以清楚地看出，本文方法得到的面內(nèi)應(yīng)力結(jié)果與文獻(xiàn)[25]的結(jié)果吻合較好。從圖7(c)和圖7(d)可以推斷，本文方法與文獻(xiàn)[25]的橫向剪切應(yīng)力在四種鋪層方式下都有較好吻合。然而在案例c 和案例d 的鋪層方式下，文獻(xiàn)[25]的 σxz(圖7(d))應(yīng)力曲線在界面處出現(xiàn)突變，這是不真實的情況，而本文方法在界面處的應(yīng)力曲線無異常突變。

文獻(xiàn)[25]以板殼理論為基礎(chǔ)，對該變剛度層合板進(jìn)行分析。鑒于分析對象屬于厚板類結(jié)構(gòu)，因此，本文從最基礎(chǔ)的線彈性理論出發(fā)，采用非協(xié)調(diào)廣義部分混合三維實體元仿真該變剛度層合板的靜力學(xué)行為，可以避免板殼理論中的諸多假設(shè)。

3 結(jié)論

考慮到解析法、板殼理論以及位移有限元法的局限性和變剛度層合板的特點，本文在常規(guī)的非協(xié)調(diào)廣義部分混合元理論的基礎(chǔ)上，為變剛度層合板的靜力學(xué)分析提出一種修正的非協(xié)調(diào)廣義部分混合元模型。主要結(jié)論如下：

(1) 在四邊固支和均布載荷作用下，對兩層的變剛度層合板進(jìn)行分析，數(shù)值結(jié)果表明，相對于ABAQUS 高階單元、高網(wǎng)格密度的模型，本文的修正的非協(xié)調(diào)廣義部分混合元模型在網(wǎng)格比較稀疏的情況下表現(xiàn)出了良好的適用性和數(shù)值結(jié)果精度高等特點。

(2) 有關(guān)四邊固支的三層變剛度夾芯板的分析表明，平面內(nèi)的位移受纖維角度的變化影響明顯。主要原因是單層板內(nèi)纖維角度變化可顯著影響該層的剛度大小。例如，單層板內(nèi)纖維角度越接近90 度，相應(yīng)的其沿x方向的剛度越小，則x方向的位移ux越大。