韓少燕, 李榆銀, 高汝鑫

(1. 西安交通大學(xué)城市學(xué)院 機(jī)械工程系,西安 710018;2. 株洲中車時(shí)代電氣股份有限公司,湖南 株洲 412001;3. 北京理工大學(xué) 先進(jìn)結(jié)構(gòu)技術(shù)研究院,北京 100081;4. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室,大連 116024)

復(fù)合材料層合圓柱殼因其加工簡單、比強(qiáng)度和比剛度高、力學(xué)性能優(yōu)異等特點(diǎn),被廣泛應(yīng)用于航空航天、艦艇、油氣輸送等各個領(lǐng)域。在工作環(huán)境中,圓柱殼結(jié)構(gòu)會受到復(fù)雜的動力載荷作用,從而產(chǎn)生有害振動甚至共振行為,對結(jié)構(gòu)的安全以及內(nèi)部設(shè)備的工作狀態(tài)產(chǎn)生不利的影響。因此,準(zhǔn)確快速地分析復(fù)合材料層合圓柱殼的動力學(xué)特性是十分必要的。

在過去的幾十年,學(xué)者們提出了不同的殼體理論以及大量的計(jì)算方法,這些理論可以在Leissa等[1-3]的書中找到詳細(xì)的介紹。殼體的理論主要分為3種:經(jīng)典殼體理論(classical shell theory, CST)、一階剪切變形理論(first-order shear deformation theory, FSDT)和高階剪切變形理論(high-order shear deformation theory, HSDT)。CST理論是針對薄壁均質(zhì)殼結(jié)構(gòu)提出的,忽略了橫向剪切變形,其要求殼體的厚度遠(yuǎn)小于其他兩個維度的尺寸。CST理論根據(jù)滿足Kirchhoff假設(shè)的薄板擴(kuò)展而來,根據(jù)不同的假設(shè)和簡化,發(fā)展出了各種薄殼理論,例如Love殼理論、Flügge殼理論、Sanders殼理論和Donnell殼理論等。隨著復(fù)合材料的發(fā)展,CST理論也被擴(kuò)展到薄壁層合殼的分析中。Lam等[4]通過引入?yún)?shù),將上述4種薄殼理論納入到統(tǒng)一的框架下,通過設(shè)定參數(shù)的值,可退化到不同的薄殼理論,用于旋轉(zhuǎn)薄壁層合殼的動力學(xué)分析中。過去的幾十年,基于薄壁層合殼體理論的數(shù)值方法也得到了充分的發(fā)展[5-6]。

一般情況下,相比于各向同性殼,層合殼的厚度更大,剪切剛度更小,考慮橫向剪切的影響是必要的。FSDT考慮了橫向剪切變形的影響,其假定橫向剪切應(yīng)變沿厚度方向?yàn)槌?shù),學(xué)者們基于FSDT發(fā)展了很多適用于中厚殼的動力學(xué)分析方法[7-12];HSDT放松了厚度方向上橫向剪切應(yīng)變?yōu)槌?shù)的假定,進(jìn)一步發(fā)展了殼體理論;此外,為考慮橫向正應(yīng)力的影響,基于HSDT的三維(3D)殼體理論[13-14]也被提出和發(fā)展。需要指出的是,雖然學(xué)者們基于FSDT和HSDT發(fā)展了眾多殼體動力分析方法,但是這些方法的推導(dǎo)和使用往往比較復(fù)雜,且一般情況下對于薄壁層合殼的分析,CST理論的精度完全滿足工程需求,故對于基于CST的薄壁層合圓柱殼的動力分析方法的研究仍然受到學(xué)者們的廣泛關(guān)注。

Yao等[15]將分析力學(xué)的Hamilton理論體系引入到彈性力學(xué)問題的求解中,其將彈性力學(xué)問題從Lagrange體系轉(zhuǎn)化到辛空間下的Hamilton體系,打破了傳統(tǒng)求解方法的局限,得到了適用于多種彈性力學(xué)問題的理性解析解。許多學(xué)者在該理論體系下進(jìn)行了深入研究。Li等[16]提出了矩形薄板自由振動的辛疊加方法,給出了有/沒有彈性基礎(chǔ)矩形薄板自由振動的解析解。Zhou等[17]利用Hamilton體系研究了圓板和環(huán)板的自由振動問題。Tong等[18]基于Donnell和Reissner殼體理論利用辛方法給出了任意邊界條件下圓柱殼自由振動的精確解。Gao等[19]針對正交各向異性圓柱殼的自由和強(qiáng)迫振動問題,提出了辛空間波傳播方法。Pan等[20-21]針對靜水壓下加環(huán)肋和雙層加筋加環(huán)肋圓柱殼的強(qiáng)迫振動問題,發(fā)展了辛空間波傳播方法。

本文針對正交鋪設(shè)復(fù)合材料層合圓柱殼的自由振動分析,提出了辛空間波方法。首先通過選取合適的狀態(tài)向量,將物理空間中的彈性力學(xué)三大方程轉(zhuǎn)換成辛空間中的振動控制方程,實(shí)現(xiàn)了問題的降維和變量的分離;其次利用分離變量法將層合圓柱殼的自由振動問題轉(zhuǎn)化為Hamilton體系下的辛本征值問題,根據(jù)辛空間波形與圓柱殼模態(tài)的對應(yīng)關(guān)系,得到不同邊界條件下的自由振動波數(shù);最后求解圓柱殼自由振動問題的代數(shù)方程組,即可得到正交鋪設(shè)復(fù)合材料層合圓柱殼的固有頻率和模態(tài)。數(shù)值算例對比了本文方法和其他方法的計(jì)算得到的圓柱殼固有頻率,驗(yàn)證了本文方法的有效性,指出了本文方法在處理非兩端簡支邊界條件時(shí)的誤差來源。

1 層合圓柱殼的基本方程

考慮如圖1所示的復(fù)合材料層合圓柱殼,圓柱殼的長度、中面半徑和厚度分別以L、R和h表示。使用固定在圓柱殼上的柱坐標(biāo)系(x、θ和r)描述圓柱殼的運(yùn)動,見圖1。

圖1 正交鋪設(shè)復(fù)合材料層合圓柱殼示意圖Fig.1 Schematic diagram of a thin composite laminated circular cylindrical shell

圓柱殼中面沿3個坐標(biāo)方向的位移分別記為u、v和w,中面繞x和θ軸的旋轉(zhuǎn)分別記為φx和φθ。復(fù)合材料層合圓柱殼由若干單層圓柱殼組成,且認(rèn)為各層間界面是剛結(jié)的。

1.1 幾何方程

根據(jù)Kirchhoff-Love薄殼理論,圓柱殼中面上的膜應(yīng)變分量和曲率應(yīng)變分量可以分別寫為

(1)

(2)

其中,

(3)

(4)

第k層上任意一點(diǎn)處的應(yīng)變-位移關(guān)系可以寫為

(5)

式中:hk+1

1.2 應(yīng)力-應(yīng)變關(guān)系

根據(jù)廣義胡克定律,第k層的應(yīng)力-應(yīng)變關(guān)系可以寫為

(6)

(7)

式中,Aij、Bij和Dij,i,j=1,2,6分別為第k層的膜剛度、膜彎耦合剛度和彎曲剛度。對于正交鋪設(shè)復(fù)合材料層合殼體,有A16=A26=D16=D26=B16=B26=0。

1.3 運(yùn)動方程

根據(jù)Kirchhoff-Love薄殼理論,層合圓柱殼自由振動的運(yùn)動方程可以寫為

(8)

式中:ω為圓頻率;ρ為圓柱殼的質(zhì)量密度。

2 辛體系下復(fù)合材料層合圓柱殼的控制方程

將式(7)寫為矩陣的形式有

F=Dε

(9)

取位移基本變量為q={uvwφx}T,可將ε重新寫為以下形式

(10)

其中,

(11)

(12)

復(fù)合材料層合圓柱殼的Lagrange密度函數(shù)可以表示為

(13)

由Legendre變換可得

(14)

(15)

將式(15)移項(xiàng)可得

(16)

其中,

(17)

(18)

將式(17)和式(18)展開可寫為

(19)

(20)

(21)

則式(16)被修正為

(22)

其中,

(23)

(24)

(25)

此時(shí),將式(22)移項(xiàng)可得

(26)

對式(8)進(jìn)行移項(xiàng),并結(jié)合式(3)和式(4),可得到

(27)

其中,

(28)

(29)

(30)

將式(26)和(27)寫成統(tǒng)一的形式,有

(31)

其中,

(32)

(33)

(34)

(35)

HC=Cγ

(36)

從式(31)～式(36)容易看出,矩陣HB和HC是對稱矩陣,所以矩陣H為Hamilton算子矩陣。式(33)中矩陣各元素的具體表達(dá)式見附錄A。

3 復(fù)合材料層合圓柱殼的自由振動分析

式(31)為齊次微分方程,其可利用分離變量法來求解,狀態(tài)向量z可以表示為以下形式

z(x,θ)=η(θ)eμx

(37)

式中,μ為軸向的波傳播參數(shù)。式(37)實(shí)現(xiàn)了變量x與θ的分離。圓柱殼在周向具有周期性,所以η(θ)可以寫為

η(θ)=φmeimθ,m=0,±1,±2,…

(38)

將式(37)～式(38)代入式(31)得到

(39)

將式(39)移項(xiàng),可以得到

(40)

式(40)有非平凡解的條件為

(41)

式(41)可以進(jìn)一步簡化為

λ3+c1(μ)λ2+c2(μ)λ+c3(μ)=0

(42)

其中,

λ=ω2

(43)

式(42)是一個一元三次方程,其可以被解析地求解,其中c1(μ)、c2(μ)和c3(μ)是由軸向波數(shù)μ決定的系數(shù)。式(42)存在3個根,分別對應(yīng)圓柱殼在3個方向上的自振頻率,其中值最小的根對應(yīng)于圓柱殼的彎曲振動。

根據(jù)駐波的定義,圓柱殼內(nèi)產(chǎn)生駐波時(shí),入射波與返回波的相位應(yīng)該相同,即

2μL+i(φr+φl)=2iNπ

(44)

式中,N為整數(shù)。所以發(fā)生駐波時(shí)的波傳播參數(shù)μN(yùn)應(yīng)滿足

(45)

式(45)即為自然模態(tài)對應(yīng)的波傳播參數(shù)的值,給定一個整數(shù)N,即可得到一個相應(yīng)的波傳播參數(shù)μN(yùn),然后將該波傳播參數(shù)μN(yùn)代入到式(42)并求解,即可得到圓柱殼的自然頻率。

將入射波與返回波的波幅疊加后即可得到自然模態(tài),以兩端簡支的圓柱殼為例,自然模態(tài)可以寫為

ψm,N=φmeimθ(eμN(yùn)x-e-μN(yùn)x)

(46)

由式(31)和式(32)可以看出,本文方法可以解析地得到圓柱殼的自然頻率,具有較高的精度和效率;另外,利用式(38)、式(39)和式(46)可以快速地得到圓柱殼的自然模態(tài)。

本文方法可進(jìn)一步擴(kuò)展到任意邊界條件下層合圓柱殼的強(qiáng)迫振動分析,相關(guān)工作正在開展中。

4 數(shù)值算例

根據(jù)Reddy的研究,Lévy方法可以解析地求解兩端簡支的層合圓柱殼的振動問題,然而Lévy方法無法求解其他(非兩端簡支)邊界條件下的層合圓柱殼的振動問題。FEM(finite element method)可以容易地處理各種復(fù)雜的邊界條件,然而FEM的計(jì)算成本較高,尤其對于求解高階模態(tài),FEM需要劃分精細(xì)的網(wǎng)格才能準(zhǔn)確地描述高階模態(tài),進(jìn)一步加劇了計(jì)算負(fù)擔(dān)。

本章首先考慮兩個兩端簡支邊界條件下的層合圓柱殼,對比本文方法、Lévy解、FEM計(jì)算得到的結(jié)果以及已有文獻(xiàn)中的結(jié)果,證明本文方法的正確性。其次依次考慮兩端固支-固支(CC)、固支-簡支(CS)、固支-自由(CF)和自由-自由(FF)邊界條件下的正交鋪設(shè)復(fù)合材料層合圓柱殼,對比利用本文方法的和FEM計(jì)算得到的自振頻率,展示本文方法處理非兩端簡支邊界條件的能力,解釋本文方法結(jié)果的誤差來源。本章使用FEM計(jì)算圓柱殼的自振頻率時(shí)采用Lanczos算法作為矩陣特征值分解方法。

4.1 兩端簡支邊界條件

考慮一個正交鋪設(shè)的復(fù)合材料層合圓柱殼,圓柱殼兩端簡支。其幾何尺寸如下:R=1.0 m、L/R=5、h/R=0.1,兩層,鋪設(shè)方案為[0°/90°]。每層的材料屬性為:質(zhì)量密度ρ=1 532 kg/m3、兩個主方向上的彈性模量E1=124.5 GPa和E2=10.2 GPa、剪切模量G12=6.3 GPa、泊松比ν12=0.34。

根據(jù)Qatu和Reddy的研究,兩端簡支邊界條件下存在Lévy解。分別利用本文方法、Lévy解和FEM求解復(fù)合材料層合圓柱殼的固有頻率,結(jié)果如表1所示,其中m和n分別表示圓柱殼在θ方向上的波數(shù)和x方向上的半波數(shù)。

表1 兩端簡支邊界條件下不用方法計(jì)算的自振頻率

由表1可知,本文方法計(jì)算的結(jié)果與Lévy解一致,與FEM計(jì)算結(jié)果的誤差也非常小,具有較好的一致性。驗(yàn)證了本文方法對兩端簡支復(fù)合材料層合圓柱殼自由振動分析的正確性。

以m=3,n=2為例,利用本文方法和FEM分別計(jì)算圓柱殼的模態(tài),其中本文方法采用和FEM相同的網(wǎng)格進(jìn)行渲染,如圖2所示?？梢钥吹?本文方法和FEM得到的圓柱殼模態(tài)形狀吻合很好,進(jìn)一步驗(yàn)證了本文方法對兩端簡支復(fù)合材料層合圓柱殼模態(tài)計(jì)算的正確性。

圖2 兩種方法計(jì)算得到的模態(tài)對比Fig.2 Comparison of the natural frequency calculated using FEM and the present method

接下來,考慮文獻(xiàn)[4]中給出的例子,兩端簡支的正交鋪設(shè)復(fù)合材料層合圓柱殼,幾何參數(shù)為:R=1.0 m、L/R=5或10、h/R=0.002,三層,鋪設(shè)方案為[0°/90°/0°],材料屬性為:質(zhì)量密度ρ=1 643 kg/m3、兩個主方向上的彈性模量E2=7.6 GPa和E1=2.5E2、剪切模量G12=4.1 GPa、泊松比ν12=0.26。

表2 兩端簡支正交鋪設(shè)復(fù)合材料層合圓柱殼的無量綱頻率參數(shù)

4.2 其他邊界條件

為了進(jìn)一步驗(yàn)證本文方法的正確性及對不同邊界條件的適用性,本節(jié)考慮不同邊界條件下的正交鋪設(shè)復(fù)合材料層合圓柱殼。圓柱殼尺寸為R=1.0 m、L/R=5、h/R=0.05,由三層單殼正交鋪設(shè)而成,鋪設(shè)方案為[0°/90°/0°]。每層的材料屬性為:質(zhì)量密度ρ=1 700 kg/m3、兩個主方向上的彈性模量E1=25.0 GPa和E2=1.0 GPa、剪切模量G12=0.5 GPa、泊松比ν12=0.25。圓柱殼的兩端邊界條件分別考慮為:CC、CS、CF和FF。不失一般性,圓柱殼的周向波數(shù)取m=1～10,軸向半波數(shù)取n=1～2。此處不對比剛體模態(tài),故FF邊界條件下,圓柱殼的周向波數(shù)取m=2～10,軸向半波數(shù)取n=0～1。

由于在非兩端簡支的情況下,不存在理論的Lévy解,所以此處將FEM計(jì)算得到的結(jié)果作為參考解。分別利用本文方法和FEM計(jì)算復(fù)合材料層合圓柱殼的自振頻率。4種不同邊界條件下兩種方法計(jì)算得到的自振頻率的對比,如表3～表6所示。其中m和n分別為圓柱殼在θ方向上的波數(shù)和x方向上的半波數(shù),表中給出的誤差為本文方法與FEM結(jié)果的相對誤差百分比。

表3 CC邊界條件下不同方法計(jì)算的自振頻率

表4 CS邊界條件下不同方法計(jì)算的自振頻率

表5 CF邊界條件下不同方法計(jì)算的自振頻率

表6 FF邊界條件下不同方法計(jì)算的自振頻率

由表3～表6可知,當(dāng)周向波數(shù)較小時(shí),本文方法計(jì)算的自振頻率結(jié)果誤差較大,這是由非簡支邊界條件下近場波的影響導(dǎo)致的。近場波是衰減波,正如表中數(shù)據(jù)所展示的,隨著周向波數(shù)的增加,近場波急劇衰減,從而使得本文方法計(jì)算的結(jié)果與FEM得到的結(jié)果吻合的非常好。

5 結(jié) 論

本文針對正交鋪設(shè)復(fù)合材料層合圓柱殼的自由振動分析提出了辛空間波方法。該方法在辛空間中建立了復(fù)合材料層合圓柱殼的振動控制方程,實(shí)現(xiàn)了問題的降維和變量的分離,進(jìn)一步利用辛空間波形與模態(tài)的對應(yīng)關(guān)系求解自由振動問題的多項(xiàng)式方程,得到復(fù)合材料層合圓柱殼的固有頻率和模態(tài)。本文方法可以處理不同的邊界條件,數(shù)值算例表明,對于兩端簡支邊界條件,本文方法總是能給出高精度的結(jié)果;對于其他邊界條件,近場波的影響導(dǎo)致本文方法存在一定的誤差,然而隨著軸向波數(shù)的增加,近場波的影響減小,此時(shí)本文方法同樣可以給出精度較高的結(jié)果。相比于Lévy解,本文方法可以處理非兩端簡支的邊界條件;相比于有限單元法,本文方法計(jì)算成本較小,尤其對于高階模態(tài)的求解,本文方法具有更顯著的效率優(yōu)勢。另外,本文方法可擴(kuò)展到任意邊界條件下層合圓柱殼的強(qiáng)迫振動分析,相關(guān)工作正在開展中。

附錄A 矩陣H中元素的表達(dá)式

由式(33)知矩陣H由3個子矩陣HA、HB和HC組成,下面分別給出上述3個矩陣的非零元素,由于HB和HC均為對稱矩陣,此處只給出兩者的下三角部分。

首先給出矩陣HA

(A.1)

(A.2)

(A.3)

(A.4)

HA,34=-1

(A.5)

(A.6)

(A.7)

其次給出矩陣HB

(A.8)

(A.9)

(A.10)

(A.11)

最后給出矩陣HC

(A.12)

(A.13)

(A.14)

(A.15)

(A.16)

(A.17)

其中,

(A.18)

α2=A12D11-B11B12

(A.19)

α3=B11D12-B12D11

(A.20)

α4=A11B12-A12B11

(A.21)

α5=A11D12-B11B12

(A.22)

α6=A66R2+2B66R+D66

(A.23)

α7=A66R+B66

(A.24)

α8=B66R+D66

(A.25)

(A.26)

(A.27)

(A.28)

(A.29)

(A.30)

α14=A12B12-A22B11

(A.31)

α15=B12D12-B22D11

(A.32)

(A.33)

(A.34)

α18=B12D12-B11D22

(A.35)