正交鋪設復合材料層合圓柱殼自由振動分析的辛空間波方法

韓少燕, 李榆銀, 高汝鑫

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

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

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

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

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

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

1 層合圓柱殼的基本方程

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

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

圓柱殼中面沿3個坐標方向的位移分別記為u、v和w,中面繞x和θ軸的旋轉分別記為φx和φθ。復合材料層合圓柱殼由若干單層圓柱殼組成,且認為各層間界面是剛結的。

1.1 幾何方程

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

(1)

(2)

其中,

(3)

(4)

第k層上任意一點處的應變-位移關系可以寫為

(5)

式中:hk+1

1.2 應力-應變關系

根據廣義胡克定律,第k層的應力-應變關系可以寫為

(6)

(7)

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

1.3 運動方程

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

(8)

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

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

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

F=Dε

(9)

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

(10)

其中,

(11)

(12)

復合材料層合圓柱殼的Lagrange密度函數可以表示為

(13)

由Legendre變換可得

(14)

(15)

將式(15)移項可得

(16)

其中,

(17)

(18)

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

(19)

(20)

(21)

則式(16)被修正為

(22)

其中,

(23)

(24)

(25)

此時,將式(22)移項可得

(26)

對式(8)進行移項,并結合式(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)中矩陣各元素的具體表達式見附錄A。

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

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

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

(37)

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

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

(38)

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

(39)

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

(40)

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

(41)

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

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

(42)

其中,

λ=ω2

(43)

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

根據駐波的定義,圓柱殼內產生駐波時,入射波與返回波的相位應該相同,即

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

(44)

式中,N為整數。所以發(fā)生駐波時的波傳播參數μN應滿足

(45)

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

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

ψm,N=φmeimθ(eμNx-e-μNx)

(46)

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

本文方法可進一步擴展到任意邊界條件下層合圓柱殼的強迫振動分析,相關工作正在開展中。

4 數值算例

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

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

4.1 兩端簡支邊界條件

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

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

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

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

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

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

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

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

4.2 其他邊界條件

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

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

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

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

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

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

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

5 結 論

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

附錄A 矩陣H中元素的表達式

由式(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)