熱沖擊下的復(fù)合材料殼剛?cè)狁詈蟿?dòng)力學(xué)研究

潘科琪，劉錦陽(yáng)

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240)

在熱沖擊載荷作用下，復(fù)合材料的各向異性性質(zhì)將引起大范圍運(yùn)動(dòng)和拉伸變形、彎曲變形、扭轉(zhuǎn)變形、剪切變形和翹曲變形的相互耦合，呈現(xiàn)高度非線(xiàn)性特征。在軌航天器的太陽(yáng)翼展開(kāi)，或者航天器變軌的過(guò)程中，受到太陽(yáng)輻射的影響，其熱流條件的突變，其溫度場(chǎng)會(huì)發(fā)生變化，從而容易誘發(fā)天線(xiàn)的熱振動(dòng)，影響航天器的姿態(tài)運(yùn)動(dòng)。復(fù)合材料在航天器中的廣泛應(yīng)用，彈性變形對(duì)航天器姿態(tài)運(yùn)動(dòng)的影響愈加明顯，因此對(duì)復(fù)合材料薄壁柔性附件精確的變形位移描述是提高復(fù)雜結(jié)構(gòu)航天器仿真精度的前提，也是航天器安全性的保障。為了保證剛－柔耦合動(dòng)力學(xué)仿真精度和動(dòng)應(yīng)力的計(jì)算精度，建立復(fù)合材料層合殼多體系統(tǒng)的動(dòng)力學(xué)模型是很有必要的。

前人工作表明，復(fù)合材料層合殼的橫向剪切變形影響比各向同性材料顯著。對(duì)于高模量復(fù)合材料(例如碳硼纖維增強(qiáng)機(jī)體復(fù)合材料)，在厚跨比h/L小于1/30時(shí)，橫向剪切變形的影響不可忽略。基于Mindlin假設(shè)，假定垂直于中面的法線(xiàn)在變形后仍為直線(xiàn)，以3個(gè)中面變形位移u，v，w和2個(gè)法線(xiàn)轉(zhuǎn)角φx，φy為變形位移坐標(biāo)，非中面上任意一點(diǎn)的變形位移的表達(dá)式中包含u，v，w，φx，φy與z的一次項(xiàng)，建立復(fù)合材料層合板的動(dòng)力學(xué)模型。基于復(fù)合材料板理論，Ding等［1］研究了閉合復(fù)合材料殼的熱彈性耦合特性。Krejaa等［2］在應(yīng)變－位移關(guān)系中考慮高次項(xiàng)，研究了大變形復(fù)合材料層合殼的動(dòng)力學(xué)特性，Abea等［3］研究了一邊固支的復(fù)合材料層合殼的共振特性，Zhou等［4］基于幾何非線(xiàn)性理論對(duì)復(fù)合材料層合殼的運(yùn)動(dòng)穩(wěn)定性進(jìn)行了分析。

隨著復(fù)合材料在工程中的大量應(yīng)用，近幾年來(lái)，國(guó)內(nèi)外學(xué)者對(duì)復(fù)合材料多體系統(tǒng)的動(dòng)力學(xué)開(kāi)展了一系列工作。Yoo等［5］等考慮幾何非線(xiàn)性和橫向剪切變形，研究了作旋轉(zhuǎn)運(yùn)動(dòng)復(fù)合材料梁的頻率特性。劉錦陽(yáng)等［6］建立了復(fù)合材料梁系統(tǒng)的動(dòng)力學(xué)模型，研究了溫度變形對(duì)復(fù)合材料梁動(dòng)力學(xué)特性的影響。Neto等［7］基于線(xiàn)彈性理論，建立了復(fù)合材料層合板多體系統(tǒng)的動(dòng)力學(xué)模型，在建模過(guò)程中考慮了橫向剪切變形，但是沒(méi)有考慮幾何非線(xiàn)性。Yoo等［8］研究了作旋轉(zhuǎn)運(yùn)動(dòng)復(fù)合材料層合板的頻率特性，考慮幾何非線(xiàn)性的影響。同樣，在熱載荷下的復(fù)合材料曲梁也是力學(xué)界研究的重點(diǎn)，早在1992年，Huang等［9］研究了受瞬時(shí)熱沖擊的殼的動(dòng)力學(xué)特性，Khdeir等［10］使用模態(tài)方法研究復(fù)合材料殼的瞬時(shí)受熱問(wèn)題。Oguamanam等［11］關(guān)于復(fù)合材料圓柱殼在溫度場(chǎng)的大范圍運(yùn)動(dòng)的動(dòng)力學(xué)性能的分析采用16節(jié)點(diǎn)的等參拉格朗日差值的單元，研究對(duì)象為常曲率的殼結(jié)構(gòu)，使其在工程中的適應(yīng)性受到很大的限制。

基于Mindlin假設(shè)，考慮了幾何非線(xiàn)性和材料的各向異性，建立了任意形狀復(fù)合材料殼的剛?cè)狁詈蟿?dòng)力學(xué)模型，研究了在熱沖擊的條件下，中心剛體－復(fù)合材料殼動(dòng)力學(xué)響應(yīng)。為了提高計(jì)算的精度，采用了16節(jié)點(diǎn)的等參單元進(jìn)行差值，由于節(jié)點(diǎn)數(shù)目較多，廣義質(zhì)量陣、力陣較為復(fù)雜，采用高斯積分法對(duì)常值陣進(jìn)行積分。殼的動(dòng)力學(xué)方程為高度非線(xiàn)性的微分代數(shù)方程，為提高計(jì)算效率，將廣義－α法結(jié)合Newton-Raphson迭代法對(duì)剛－柔耦合的動(dòng)力學(xué)方程進(jìn)行積分，求出了動(dòng)力學(xué)響應(yīng)的時(shí)間歷程。最后給出了仿真算例，對(duì)復(fù)合材料殼的動(dòng)力學(xué)特性進(jìn)行了研究。

1 復(fù)合材料殼的數(shù)學(xué)模型

圖1表示一個(gè)中心剛體和一個(gè)任意形狀柱狀殼結(jié)構(gòu)，殼的一端固支與中心剛體。為了方便描述，引入了慣性坐標(biāo)系，連體基坐標(biāo)系0(bbb)以及局部坐標(biāo)系c。圖2表示殼的中線(xiàn)面在連體基坐標(biāo)系bb面的形狀，z=f(x)為形函數(shù)，對(duì)于只有一個(gè)方向曲率的柱狀殼，柱狀殼在連體基上的方程為以及連體基內(nèi)任意點(diǎn)的曲率半徑么分別為:

其中:H表示殼寬度方向的長(zhǎng)度。

設(shè)熱沖擊ΔT沿殼的厚度方向均勻分布，具體表示為:

圖1 復(fù)合材料層合殼Fig.1 Laminated composite shell

其中:ΔTss為溫度載荷穩(wěn)態(tài)時(shí)的幅值，t為仿真時(shí)間，Tf與熱沖擊時(shí)間相關(guān)的溫度常數(shù)。

圖2 任意形狀殼在連體基面的投影函數(shù)bbFig.2 The functions of shell center plane in body coordinate with arbitrary shape

根據(jù)Mindlin理論，殼上非中面上任意點(diǎn)的變形可以表示為:

2 復(fù)合材料殼的運(yùn)動(dòng)學(xué)描述

式(6)的變分表達(dá)式以及相對(duì)于時(shí)間的二階導(dǎo)數(shù)分別為:

其中:

3 有限元離散

等參單元離散，上任意點(diǎn)變形和轉(zhuǎn)角可以分別表示為:

其中:－1≤ξ，η≤1，ξ，η均為無(wú)量綱變量，pe為單元坐標(biāo)陣為 Ni，i=1，2，3，4，5 為型函數(shù)，具體表達(dá)式請(qǐng)參見(jiàn)附件。

單元上任意點(diǎn)的在連體基上的坐標(biāo)可以差值為:

根據(jù)式(11)，殼體微元表示為:

其中:h為殼的厚度，ls為曲線(xiàn)的弧長(zhǎng)。Jxe=

根據(jù)式(3)～式(5)以及式(10)，任意點(diǎn)在e→c上的變形為:

4 考慮熱效應(yīng)的應(yīng)力應(yīng)變關(guān)系

結(jié)合式(9)和(10)，考慮幾何非線(xiàn)性的殼結(jié)構(gòu)的應(yīng)變－位移關(guān)系［11］可以表示為:

考慮熱載荷及橫向剪切變形的復(fù)合材料板第k層應(yīng)力應(yīng)變的本構(gòu)關(guān)系為:

其中:

5 動(dòng)力學(xué)方程

采用虛功原理來(lái)建立系統(tǒng)的動(dòng)力學(xué)方程，對(duì)于任意單個(gè)殼體，不計(jì)體力以及其它外載荷情況，虛功原理的形式如下:其中:δWf為慣性力的虛功，δWin為彈性力的虛功:

設(shè)中心剛體－復(fù)合材料殼系統(tǒng)的廣義坐標(biāo)q=［θpT］T，根據(jù)式(10)和式(11)，并結(jié)合有限元離散化方法，慣性力的虛功表示為:

其中:

各項(xiàng)具體表示為:

其中，常值陣為:

將式(25)和式(26)代入上式，得到復(fù)合材料殼彈性力的虛功為:

其中:N為復(fù)合材料的層數(shù)。

需要特殊說(shuō)明的是因?yàn)棣腃Θ=CδΘ，所以非線(xiàn)性應(yīng)變的變分可以表示為:δεN=CGδpe，據(jù)式(33)，復(fù)合材料殼的彈性力的虛功可以具體表示為:

其中:fb，fs，fT具體表達(dá)式為:

其中:矩陣 A，B，D，Aα，Bα的表達(dá)式請(qǐng)參見(jiàn)文獻(xiàn)［12］。

幾何非線(xiàn)性變形與熱載荷耦合引起的彈性力表示為:

若鉸約束方程為Φ(q，t)=0，系統(tǒng)封閉的第一類(lèi)拉格朗日動(dòng)力學(xué)方程為:

其中:λ為拉格朗日乘子坐標(biāo)陣。

6 數(shù)值計(jì)算方法

本文將廣義－α［13］法和Newton-Raphson迭代法結(jié)合求解動(dòng)力方程(24)，為了方便表述，將動(dòng)力學(xué)方程表示為:

主要步驟為:

(1)基于前一時(shí)刻t的動(dòng)力學(xué)響應(yīng)，根據(jù)廣義－α法的預(yù)估關(guān)系式估計(jì)下一時(shí)刻t+Δt的位移、速度、加速度第一次(n=1)的迭代值為:

(3)根據(jù)式(27)，第 n+1次的迭代值可以表示為:

(4)第n+1次的速度、加速度以及拉格朗日乘子的迭代值可以表示為:

其中:α'=(1 －αm)/［Δt2α(1 －αf)］，β'=β/(Δtα)

(5)將位移、速度、加速度代入式(25)中的Φ—1，判斷是否滿(mǎn)足允許誤差:

若在某次的迭代值滿(mǎn)足上式關(guān)系，則把其作為t+Δt時(shí)刻動(dòng)力學(xué)響應(yīng)精確值，返回第一步進(jìn)行下一時(shí)刻的動(dòng)力學(xué)響應(yīng)的求解。

7 仿真算例

其中，本算例中 H=0.1 m，L=0.8 m，d分別取為0.05 m，0.10 m，0.15 m 以及 d=0。d=0 時(shí)殼變?yōu)槠矫姘濉Ｓ蓤D3可知，若L保持不變，d越大曲率越大，本算例中H和L都保持不變。

圖3 仿真算例復(fù)合材料殼在連體基坐標(biāo)bb面的函數(shù)Fig.3 The projection functions in simulating example

復(fù)合材料的材料參數(shù)為:彈性模量E1=181 GPa，E2=103 GPa，G12=G13=7.7 GPa，G23=2.87 GPa，復(fù)合材料的鋪層為對(duì)稱(chēng)形式［0°/90°/0°/90°/0°/0°/90°/0°/90°/0°］，泊松比 v12=0.28，主方向的熱膨脹系數(shù)分別為 α11=2.08 × 10－8，α22=22.5 × 10－6，密度 ρ=1 560 kg/m3。本算例中熱沖擊時(shí)間即式(3)中的Tf=0.2 s，熱沖擊在穩(wěn)態(tài)值即式(3)中的ΔTss=2 K。中心剛體的轉(zhuǎn)動(dòng)慣量J0=0.01 kg·m2。初始時(shí)刻中心剛體－殼系統(tǒng)的角速度均為零。熱沖擊作用下引起殼結(jié)構(gòu)變形，并且由于剛－柔耦合的作用，殼結(jié)構(gòu)變形的同時(shí)也引起了中心剛體的轉(zhuǎn)動(dòng)。

7.1 本文模型正確性和收斂性的驗(yàn)證

為了驗(yàn)證本文模型的正確性，在中心剛體被約束不動(dòng)的情況下，將本文提出的等參單元計(jì)算的懸臂殼的頻率與ANSYS軟件用殼單元離散得到的頻率進(jìn)行對(duì)比，對(duì)比結(jié)果如表1所示。從表中可以看出，無(wú)論曲率較大d=0.15m還是較小d=0.05m時(shí)，通過(guò)本文模型計(jì)算得到的結(jié)果與ANSYS軟件計(jì)算得到的結(jié)果誤差都在5%，說(shuō)明現(xiàn)有模型的正確性。

為了驗(yàn)證本文單元的收斂性，在中心剛體不受約束、初速度為零、有熱沖擊的情況下，對(duì)曲率較大，d=0.15 m的情況下進(jìn)行了動(dòng)力學(xué)仿真，如圖4所示為非線(xiàn)性模型計(jì)算得到的中心剛體的轉(zhuǎn)角，可以看到最多15個(gè)單元就可以收斂證明了本文單元的收斂性。

表1 現(xiàn)有方法與ANSYS軟件的頻率比較Tab.1 Comparison of present frequency results with those obtained by ANSYS Software

7.2 線(xiàn)性和非線(xiàn)性模型動(dòng)力學(xué)響應(yīng)比較

如圖5～7中所示為 d=0.15 m，d=0.05 m，d=0 m的工況下，線(xiàn)性模型和非線(xiàn)性模型在熱沖擊作用下的縱向變形對(duì)比圖。圖5所示為d=0.15 m時(shí)，線(xiàn)性模型和非線(xiàn)性模型縱向變形的對(duì)比。由圖5可見(jiàn)線(xiàn)性模型在熱沖擊載荷作用下，縱向變形明顯小于非線(xiàn)性模型，熱沖擊穩(wěn)態(tài)時(shí)縱向變形有高頻振蕩，而非線(xiàn)性模型的時(shí)間歷程曲線(xiàn)較為平滑。線(xiàn)性模型與非線(xiàn)性模型的仿真結(jié)果差異顯著，主要是由于線(xiàn)性模型沒(méi)有考慮幾何非線(xiàn)應(yīng)變項(xiàng)εN，熱沖擊作用下的仿真結(jié)果誤差較大。如圖6和圖7所示分別為d=0.05 m和d=0 m時(shí)縱向變形的對(duì)比，可以看到當(dāng)d=0.05 m時(shí)，線(xiàn)性和非線(xiàn)性模型的結(jié)果差異比d=0.15 m時(shí)二者的差異有顯著地減少。當(dāng)d=0 m時(shí)，即沒(méi)有曲率時(shí)，線(xiàn)性和非線(xiàn)性模型的結(jié)果幾乎沒(méi)有差別，說(shuō)明殼結(jié)構(gòu)比平面板結(jié)構(gòu)在受到熱沖擊時(shí)更容易產(chǎn)生非線(xiàn)性效應(yīng)，并且曲率越大非線(xiàn)性效應(yīng)越明顯。

圖4 非線(xiàn)性模型仿真結(jié)果Fig.4 Results of nonlinear model

圖5 縱向變形Fig.5 Longitudinal deformation

圖6 縱向變形Fig.6 Longitudinal deformation

7.3 曲率、材料屬性對(duì)動(dòng)力學(xué)響應(yīng)的影響

本小節(jié)的結(jié)果均為非線(xiàn)性模型計(jì)算得到。圖8和圖9為不同曲率的復(fù)合材料殼在熱沖擊下的動(dòng)力學(xué)響應(yīng)。從圖8可以看出，曲率越大，也是d越大，縱向變形越明顯。同樣，從圖9可以看，曲率越大，中心的剛體的轉(zhuǎn)角就越大。另外，從圖8還可以看出，d=0和d=0.05時(shí)的縱向變形差異較小，熱載荷穩(wěn)態(tài)時(shí)，兩者幅值的差異值不到d=0和d=0.1時(shí)的差異值的1/3倍，但是如圖9所示，d=0和d=0.05時(shí)的中心剛體的轉(zhuǎn)動(dòng)角度差異值已經(jīng)接近d=0和d=0.1時(shí)轉(zhuǎn)角差異值的1/2，說(shuō)明由于剛?cè)狁詈系淖饔茫⑿〉目v向變形就能引起較大的中心剛體的轉(zhuǎn)動(dòng)。

圖8 縱向變形Fig.8 Longitudinal deformation

圖9 中心剛體的轉(zhuǎn)角Fig.9 Rotation angle of the center hub

圖10 中心剛體的轉(zhuǎn)角Fig.10 Rotation angle of the center hub

為了研究熱沖擊工況下異性材料和同性材料對(duì)系統(tǒng)的動(dòng)力學(xué)響應(yīng)的影響，圖10對(duì)比了復(fù)合材料和同性材料殼的仿真結(jié)果。本文定義的同性材料是將算例中的復(fù)合材料參數(shù)設(shè)置為E1=E2=181 GPa，G12=G13=G23=7.7 GPa，并將鋪層設(shè)為對(duì)稱(chēng)形式［0°/0°/0°/0°/0°/0°/0°/0°/0°/0°］，熱膨脹系數(shù) α11= α22=22.5 ×10－6，其它參數(shù)都變，此時(shí)可以視為同性材料，以便減少其它因素的引入對(duì)計(jì)算結(jié)果的影響。同性材料和復(fù)合材料殼的曲率參數(shù)d=0.05 m。從圖10中可以看出，復(fù)合材料殼的中心剛體的轉(zhuǎn)角明顯大于同性材料的，這是因?yàn)閺?fù)合材料的各向異性，使得熱膨脹引起的熱應(yīng)力分布不均，熱應(yīng)力梯度引起了中心剛體的轉(zhuǎn)動(dòng)，而同性材料殼，即使有曲率的存在，由于各處熱膨脹系數(shù)相同，各處的熱應(yīng)變相同，不會(huì)對(duì)中心剛體產(chǎn)生力偶矩。

8 結(jié)論

本文建立了考慮熱載沖擊的中心剛體－變曲率復(fù)合材料殼的動(dòng)力學(xué)模型，并與ANSYS軟件計(jì)算得到的頻率對(duì)比說(shuō)明本文模型的正確性。通過(guò)將廣義－α法與Newton－Raphson迭代法結(jié)合的方式對(duì)動(dòng)力學(xué)方程進(jìn)行積分。仿真結(jié)果表明，在熱沖擊載荷下，復(fù)合材料殼線(xiàn)性模型仿真結(jié)果誤差較大，并且隨著曲率的增加，線(xiàn)性和非線(xiàn)性模型的結(jié)果差異越明顯，需要考慮幾何非線(xiàn)性項(xiàng)才能滿(mǎn)足精度要求。復(fù)合材料殼的曲率越大，在熱沖擊作用下的縱向變形越大，變形和大范圍運(yùn)動(dòng)的耦合作用就越明顯。同性材料的殼在熱沖擊作用下不會(huì)引起中心剛體的大范圍運(yùn)動(dòng)。

［1］Ding K W.The thermoelastic dynamic response of thick closed laminated shell［J］.Shock and Vibration，2005，12(4):283－291.

［2］Krejaa I，Schmidt R.Large rotations in first-order shear deformation FE analysis of laminated shells［J］.International Journal of Non-linear Mechanics，2006，41(1):101 －123.

［3］Abea A，Kobayashib Y，Yamada G.Nonlinear dynamic behaviors of clamped laminated shallow shells with one-to-one internal resonance［J］.Journal of Sound and Vibration，2007，304(3－5):957－968.

［4］Zhou C T，Wang L D.Nonlinear theory of dynamic stability for laminated composite cylindrical shells［J］.Applied Mathematics and Mechanics，2001，22(1):53 －62.

［5］Yoo H H，Lee S H，Shin S H.Flapwise bending vibration analysis of rotating multi-layered composite beams［J］.Journal of Sound and Vibration，2005，286(4－5):745－761.

［6］劉錦陽(yáng)，潘科琪.考慮熱效應(yīng)的復(fù)合材料多體系統(tǒng)動(dòng)力學(xué)研究［J］.動(dòng)力學(xué)與控制學(xué)報(bào)，2009，7(1):9－13.LIU Jin-yang， PAN Ke-qi. Dynamic invetigation on composite flexible multi－body system considering thermal effect［J］.Journal of Dynamics and Control，2009，7(1):9－13.

［7］Neto M A，Ambrósio J A C，Leal R P.Composite materials in flexible multibody systems［J］.Computer Methods in Applied Mechanics and Engineering，2006，195(50－51):6860－6873.

［8］Yoo H H，Kim S K，Inman D J.Modal analysis of rotating composite cantilever plates［J］.Journal of Sound and Vibration，2002，258(2):233－246.

［9］Huang N N，Tauchert T R.Thermally induced vibration of doubly curved cross-ply laminated panels［J］.Journal of Sound Vibration，1992，154(3):485－494.

［10］Khdeir A A. Thermallyinduced vibration ofcross-ply laminated shallow shells［J］.Acta Mechanica，2001，151(3－4):135－147.

［11］Oguamanam D C D，Hansen J S，Heppler G.R.Nonlinear transient response of thermally loaded laminated panels［J］.Journal of Applied Mechanics，2004，71:49 －56.

［12］Pan K Q，Liu J Y.Dynamic investigation on composite flexible multi-body system considering thermal effect［J］.Journal of Shanghai Jiaotong University，2010，15(4):414 －422.

［13］Arnold M，Brüls O.Convergence of the generalized － α scheme for constrained mechanical systems［J］.Multibody System Dynamics.2007，18:185 －202.

［14］王勖成.有限單元法［M］.北京:清華大學(xué)出版社，2006.

附 錄

16 節(jié)點(diǎn)的拉格朗日型單元的型函數(shù) Si，i=1，2，…，16，可參見(jiàn)文獻(xiàn)［14］，對(duì)角陣 Si=diag(Si)16×16，S=［S1S2S3… S16］，式(9)和式(10)中的差值函數(shù)為:Ni=S(i，∶)(i=1，2，3，4，5)，S(i，∶)表示矩陣 S 中第 i行的所有列，(ξ，η)=(S1S2… S16)。