劉瑞杰，王雪仁,2，賈地

(1.中國人民解放軍92578 部隊，北京 100161；2.哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱，150001)

梁結(jié)構(gòu)廣泛應(yīng)用于船舶與海洋工程、建筑工程和航空航天等領(lǐng)域中，對于此類結(jié)構(gòu)的振動特性及其參數(shù)影響規(guī)律研究受到科研工作者的廣泛關(guān)注[1-11]。隨著材料技術(shù)和加工工藝的進步，復(fù)合材料由于自身比剛度大、可設(shè)計性強、機械性能突出等優(yōu)點，成為了工程應(yīng)用中的首選材料。復(fù)合材料往往呈現(xiàn)出較為明顯的正交各向異性，對各向異性梁結(jié)構(gòu)動力學(xué)特性的研究具有十分重要的意義。經(jīng)典的梁理論模型主要包含Euler-Bernoulli和Timoshenko 2種理論。然而，歐拉-伯努利梁理論忽略了結(jié)構(gòu)橫向拉伸和剪切變形以及轉(zhuǎn)動慣量的影響，其計算結(jié)果與實際情況相比存在較大誤差。鐵木申科梁理論雖然考慮了橫向剪切變形和轉(zhuǎn)動慣量的影響，但其引入的剪切修正因子對梁結(jié)構(gòu)的幾何截面形狀和材料參數(shù)均具有較高的依賴性，故而不適用于復(fù)合材料梁結(jié)構(gòu)的研究。另外，經(jīng)典梁理論更多地關(guān)注于梁結(jié)構(gòu)沿著特定方向上的形變或振動問題。當(dāng)梁結(jié)構(gòu)的長徑比較小或截面為不規(guī)則形狀時，梁結(jié)構(gòu)在外界載荷作用下不僅會出現(xiàn)彎曲變形，亦會產(chǎn)生扭轉(zhuǎn)、拉伸或彎扭組合變形。此時需尋求具有較高計算精度和適用性的三維理論。Carrera等[12]基于卡諾統(tǒng)一定理提出了一種改進的高階梁理論模型。該理論將梁結(jié)構(gòu)的位移容許函數(shù)分解為截面和軸向2部分，利用高階泰勒公式擬合梁結(jié)構(gòu)的截面位移，從而以泰勒展開的階次控制理論的計算精度。該理論模型也被成功地應(yīng)用于求解一些復(fù)雜問題，如薄壁結(jié)構(gòu)、復(fù)合材料問題、旋轉(zhuǎn)葉片等[13-16]，并展現(xiàn)出較高的計算精度。

本文結(jié)合卡雷拉理論(Carrera unified formulation,CUF)和有限元理論建立正交各向異性梁結(jié)構(gòu)的三維振動分析模型。基于梁結(jié)構(gòu)的應(yīng)變能和動能，由最小勢能原理推導(dǎo)梁單元的質(zhì)量矩陣和剛度矩陣。分別從質(zhì)量和剛度矩陣中提取出3×3階的質(zhì)量、剛度核心矩陣，利用循環(huán)得到整個梁的質(zhì)量矩陣和剛度矩陣。隨后對多種邊界條件下各向異性梁結(jié)構(gòu)進行了自由振動分析，并將計算結(jié)果與有限元分析軟件結(jié)果對比分析。

1 正交各向異性梁三維理論模型

1.1 梁結(jié)構(gòu)位移場離散

如圖1所示，考慮一矩形正交各向異性梁結(jié)構(gòu)，建立直角坐標(biāo)系(x，y和z)。正切軸x與生長環(huán)相切，縱向軸y與纖維方向平行,徑向軸z與生長環(huán)垂直。梁形結(jié)構(gòu)上任意一點沿坐標(biāo)x，y和z3個方向上的位移變形分別用Ux，Uy和Uz表示。梁結(jié)構(gòu)的長度為L，寬度和高度分別為b和h。

圖1 正交各向異性梁結(jié)構(gòu)示意Fig.1 The dimensions of an orthotropic beam

為了簡化三維實體模型，本文通過分離變量將梁形結(jié)構(gòu)的位移函數(shù)分為截面位移函數(shù)和軸向位移函數(shù)2個部分。其中，利用卡諾高階理論對梁結(jié)構(gòu)截面面內(nèi)位移進行擬合，卡諾高階理論中用到的二維泰勒展開分量可以簡寫為統(tǒng)一的形式。泰勒展開的階次可直接控制該方法的計算精度。梁結(jié)構(gòu)的位移函數(shù)為：

(1)

式中：N為二維泰勒展開的階次；Γ為二維泰勒展開的總項數(shù)；Fτ表示泰勒展開中的第τ項；變量uxτ、uyτ和uzτ為梁形結(jié)構(gòu)的軸向位移變量，是關(guān)于軸向y的函數(shù)。

軸向位移變量uxτ、uyτ和uzτ利用有限元方法進行軸向解域離散，可表示為：

(2)

以2節(jié)點線性梁單元(M=2)為例，函數(shù)為：

(3)

將式(2)代入式(1)中，梁結(jié)構(gòu)的位移函數(shù)為：

(4)

通過式(4)可以看出，對于任意一梁單元，每一節(jié)點有3×Γ個未知變量，用向量表示為：

q=(A11,B11,C11,…,AΓ1,BΓ1,CΓ1)T

(5)

值得注意的是，通過增加分析中所使用的線性單元的數(shù)目，或使用高階的插值函數(shù)，都能夠提高有限元結(jié)果的精度。

1.2 梁結(jié)構(gòu)自由振動控制方程

考慮結(jié)構(gòu)發(fā)生微小變形，根據(jù)應(yīng)變位移關(guān)系可得：

(6)

根據(jù)廣義胡克定律，正交各向異性結(jié)構(gòu)的應(yīng)力應(yīng)變關(guān)系可表示為：

(7)

正交各向異性材料的彈性系數(shù)可表示為[17]：

(8)

(9)

(10)

另外，υij和υji存在以下關(guān)系：

υijEj=υjiEi

(11)

如此，需要由9個獨立參數(shù)來確定正交各向異性材料的物性參數(shù)：E1、E2、E3、G12、G13、G23、υ12、υ13、υ23。

基于能量法推導(dǎo)梁單元的剛度矩陣梁單元的應(yīng)變能：

(12)

式中：ε為應(yīng)變向量；σ為應(yīng)力向量；K即為梁單元剛度矩陣。仔細觀察方陣K的內(nèi)部構(gòu)造可以發(fā)現(xiàn)，該方陣是由一個3×3階方陣Kijτs堆積而成，其中的9個組成部分如下：

(13)

(14)

(15)

式中：〈FτFs〉Ω表示截面內(nèi)積分；〈NiNj〉l表示軸向單元內(nèi)積分。當(dāng)梁結(jié)構(gòu)的厚度較小時，此處的有限元方法會產(chǎn)生剪切自鎖現(xiàn)象，本文用低階的高斯-勒讓德數(shù)值積分消除此處的剪切自鎖現(xiàn)象[18]。

梁單元的動能由各連續(xù)質(zhì)點的動能累加而成：

(16)

式中：ρ為梁的密度；u′為梁上任意一點的速度，是關(guān)于時間t的函數(shù)。將位移函數(shù)代入上式中，則：

(17)

與剛度矩陣相似，Mijτs為質(zhì)量矩陣的核心，質(zhì)量矩陣可由3×3階方陣Mijτs堆積而成，其中的9個組成部分如下：

(18)

采用劃行化列、賦大數(shù)等有限元邊界施加方法處理正交各向異性梁結(jié)構(gòu)的邊界條件，正交各向異性梁自由振動控制方程可表示為：

(Kz-ω2Mz)q=0

(19)

式中：Kz為梁的總剛度矩陣；Mz為總質(zhì)量矩陣。

2 正交各向異性梁自由振動分析

本文通過對一正交各向異性梁結(jié)構(gòu)的自由振動分析來驗證本文模型的合理性。考慮梁結(jié)構(gòu)的材料為正交各向異性材料，其材料屬性可表示為：E1=145 GPa，E2=9.6 GPa，G12=G13=4.1 GPa，G23=3.4 GPa，υ12=0.3，密度ρ=1 389 kg/m3。為了方便與文獻[19]對比，纖維方向角取為90°，梁結(jié)構(gòu)的長高比取為10，考慮兩端固支邊界條件，并引入無量綱頻率參數(shù)：

(20)

表1給出了正交各向異性梁結(jié)構(gòu)的前6階彎曲模態(tài)的無量綱頻率參數(shù)對比結(jié)果。本文計算結(jié)果與文獻[19]結(jié)果具有較好的一致性。值得注意的是，本文所建立的梁結(jié)構(gòu)理論模型為三維模型，不僅可以預(yù)測梁結(jié)構(gòu)的橫向彎曲振動，也可以獲得梁結(jié)構(gòu)的扭轉(zhuǎn)模態(tài)和縱振模態(tài)。

表1 正交各向異性梁前6階彎曲模態(tài)的無量綱頻率參數(shù)Table 1 First six non-dimensional frequencies of bending modes of an orthotropic beam

為了進一步驗證本文三維振動分析結(jié)果的收斂性和正確性，本文將與有限元商業(yè)軟件ANSYS的三維仿真結(jié)果進行對比。選取梁結(jié)構(gòu)材料(graphite fabric-carbon matrix layers，GFCML)，其材料屬性可表示為：E1=173.1 GPa，E2=33.1 GPa，E3=5.17 GPa，G12=9.38 GPa，G13=8.27 GPa，G23=3.24 GPa，υ12=0.036，υ13=0.25，υ23=0.171，密度ρ=1 650 kg/m3。令梁結(jié)構(gòu)的寬度和高度相等b=h=0.2 m，長度L為2 m。為了較為準(zhǔn)確地預(yù)測梁結(jié)構(gòu)的扭轉(zhuǎn)模態(tài)，本文選取截面擬合多項式的階次N為4。

表2為正交各向異性梁結(jié)構(gòu)在簡支邊界條件下，前10階模態(tài)的收斂情況，每一行代表結(jié)構(gòu)不同的模態(tài)階次，每一列代表著不同的單元類型和單元數(shù)。B2表示2節(jié)點線性單元，B3表示3節(jié)點拋物線單元。從表中可以較為清楚地看到隨著單元數(shù)量的增加，本文方法的計算結(jié)果快速收斂至穩(wěn)定值。選用較高的單元類型會進一步加快計算方法的收斂速度。同時，表2也列出了利用ANSYS軟件計算出的數(shù)據(jù)，對于該仿真結(jié)果，本文選用了SOLID186單元類型，劃分640個單元，3 665個節(jié)點。其有限元網(wǎng)格劃分結(jié)果如圖2所示。另外，2種計算方法所需的自由度數(shù)(DOFs)均在表2中給出。通過對比可以發(fā)現(xiàn)，本文方法所需較少的自由度數(shù)而獲得與ANSYS相近的計算結(jié)果。

表2 簡支邊界下正交各向異梁結(jié)構(gòu)前10階模態(tài)收斂性分析Table 2 Convergence of the first ten modes of an orthotropic beam with simple boundary condition Hz

圖2 有限元網(wǎng)格劃分結(jié)果Fig.2 Mesh dividing results of the finite element method

表3為3種經(jīng)典邊界條件(簡支、懸臂和固支邊界)下正交各向異性梁結(jié)構(gòu)的前10階固有頻率。通過與ANSYS仿真結(jié)果的對比，可以發(fā)現(xiàn)，對于不同的約束條件，本文方法均具有較高的計算精度。另外圖3給出了簡支邊界條件下正交各向異性梁結(jié)構(gòu)的振型圖，圖3(a)～(c)為ANSYS仿真結(jié)果，圖3(d)～(f)為本文方法繪制的結(jié)果。(a)和(d)為梁結(jié)構(gòu)第1階彎曲振型，(b)和(e)為結(jié)構(gòu)1階扭轉(zhuǎn)振型，(c)和(f)為結(jié)構(gòu)第1階縱振振型。可以看出本文方法不僅可以用于預(yù)測正交各向異性梁結(jié)構(gòu)的彎曲模態(tài)，也可以較為準(zhǔn)確地預(yù)測梁結(jié)構(gòu)的扭轉(zhuǎn)模態(tài)和縱振模態(tài)。

表3 正交各向異性梁結(jié)構(gòu)不同邊界條件下的前10階固有頻率分析結(jié)果Table 3 Frequency of first ten modes of an orthotropic beam with different boundary conditions Hz

圖3 簡支邊界條件下正交各向異性結(jié)構(gòu)的三維振型圖Fig.3 Three-dimensional mode shapes of an orthotropic beam with a simple boundary condition

3 結(jié)論

1)本文方法不僅可以提供高精度的解，還可以較為全面地展現(xiàn)其振動特性(彎曲模態(tài)、扭轉(zhuǎn)模態(tài)、縱振模態(tài))。

2)本文方法所需較少的自由度數(shù)而獲得與商業(yè)有限元相同的計算精度。

3)對于正交各向異性梁結(jié)構(gòu)的幾何參數(shù)設(shè)計及動力學(xué)特性優(yōu)化具有指導(dǎo)意義。