基于波傳播方法和多元分析的正交各向異性圓柱殼振動特性研究

汪志強， 李學斌， 黃利華

(武漢第二船舶設(shè)計研究所，武漢 430205)

圓柱殼結(jié)構(gòu)在工業(yè)上得到了廣泛的應(yīng)用。隨著復(fù)合材料的使用，很多圓柱殼結(jié)構(gòu)形式都能夠采用正交各向異性圓柱殼的模型進行討論[1-2]。這種圓柱殼在軸向和周向具有不同的楊氏彈性模量和泊松比，并滿足一定的約束。很多學者對這種圓柱殼結(jié)構(gòu)形式的動態(tài)特性進行了研究，如自由振動[3]，沖擊響應(yīng)研究[4]等。

在對圓柱殼的研究中，為了探索物理參數(shù)以及幾何參數(shù)對于動態(tài)特性(如頻率)，通常的做法是，先設(shè)定物理參數(shù)以及幾何參數(shù)的范圍，確定需要研究的因素(如厚度和半徑比值)，通過變化這些因素來求解它對頻率的影響。這樣往往會固定某些參數(shù)以突出該因素的作用進行分析。這樣的方法效率比較低，難以顧及其他參數(shù)的同時變化，因此對于振動特性的全局性分析存在一定的局限性。

本文基于波傳播方法[5]討論了正交各向異性圓柱殼在多種邊界下的自由振動特性。給定一定的物理和幾何參數(shù)之后，通過試驗設(shè)計[6]的方法獲得很多頻率解，這樣就得到了比較大的設(shè)計空間。為了討論該設(shè)計空間的屬性，本文引入了多元變量分析的概念[7]，討論了這些幾何參數(shù)、材料參數(shù)與頻率的關(guān)聯(lián)性。通過方差分析(Analysis of Variance, ANOVA)判斷出這些參數(shù)對頻率影響的重要程度，可以為圓柱殼的動態(tài)特性設(shè)計提供參考。自組織映射[8]是神經(jīng)網(wǎng)絡(luò)技術(shù)的一種，它在數(shù)據(jù)分析上具有重要的地位。本文還采用該技術(shù)對這個空間進行了數(shù)據(jù)分析以及可視化研究。

1 理論分析

1.1 頻率求解

正交各向異性圓柱殼和坐標系示意圖，如圖1所示。圓柱殼半徑R，長度L，厚度h。假定材料是正交各向異性，剪切彈性模量G。圓柱殼中面點在坐標軸方向的位移分別為u，v和w(向內(nèi)為正)。

圖1 正交各向異性圓柱殼及振型示意圖

設(shè)圓柱殼的材料主軸和坐標系主軸重合，則滿足如下關(guān)系式

vxEθ=vθEx

(1)

式中：Ex，vx和Eθ，vθ分別是軸向和切向的楊氏彈性模量與泊松比。應(yīng)用Flügge經(jīng)典殼體理論[9]，軸向壓力P作用下正交各向異性圓柱殼的振動方程為

(2)

Dx,Dθ,Dxθ為拉伸剛度；Kx,Kθ,Kxθ為彎曲剛度

(3)

根據(jù)波傳播方法，可假設(shè)位移函數(shù)形式為

u(x,θ,t)=U0e-iknsxcos(nθ)eiωt，

v(x,θ,t)=V0e-iknsxsin(nθ)eiωt，

w(x,θ,t)=W0e-iknsxcos(nθ)eiωt

(4)

式中：U0、V0和W0為振型幅值；kns為軸向波數(shù)；n為環(huán)向波數(shù)；t為時間；ω為圓頻率(1/s)；i2=-1。將這個位移函數(shù)代入式(2)，可以得到如下矩陣形式的方程

(5)

矩陣A的元素如下

圓柱殼的振動頻率就是這個齊次方程組的非零解，即要求系數(shù)行列式等于零。將行列式展開，得到如下關(guān)于圓頻率ω2的3次方程

g6(ω2)3+g4(ω2)2+g2ω2+g0=0

(6)

式(6)可通過解析公式準確求解。通過式(6)求解圓柱殼的頻率，還需要考慮圓柱殼兩端的邊界條件。本文采用波傳播方法，即利用對應(yīng)梁彎曲振動的邊界條件[10]近似為圓柱殼的邊界條件。這里考慮了3種邊界情況：兩端薄膜簡支，兩端固定以及簡支-固定邊界。這些邊界條件對應(yīng)的梁特征函數(shù)以及解的形式在表1中給出。表中m是軸向形成的振型波數(shù)。

表1 梁彎曲振動特征函數(shù)和波數(shù)解

1.2 多變量分析

在圓柱殼振動的研究和設(shè)計過程中，研究者往往都會通過變換不同的參數(shù)，獲得它們對于頻率的影響，并研究各參數(shù)的重要性?；蛘哒f，頻率對哪個參數(shù)的變化更加敏感一些。通常還可以設(shè)計成圖譜等形式。對于正交各向異性圓柱殼而言，結(jié)構(gòu)參數(shù)比較多，特性變化也比較復(fù)雜。本文擬采用試驗設(shè)計的方法，研究圓柱殼的材料特性以及幾何參數(shù)對于自振頻率的影響。本文主要關(guān)注兩個事項，材料特性以及幾何參數(shù)中，參數(shù)和頻率的相關(guān)性如何，以及參數(shù)改變對頻率變化的影響分析。

從前面的分析可知，可認為圓柱殼頻率ω實際上是材料特性、幾何參數(shù)、邊界條件和外力的函數(shù)，可以表示為

(7)

多元統(tǒng)計分析就是以p個變量，N次觀察數(shù)據(jù)所形成的數(shù)據(jù)矩陣X=[xij]，i=1,2,…,N；j=1,2,…,p為依據(jù)。矩陣X的第i行表示第i個樣品的觀察值，它是一個p維的向量。矩陣X的第j列表示對第j個變量的n次觀測值，它是一個N維變量。

變量之間的相互依賴關(guān)系多元分析的一個重要內(nèi)容，它研究1個或幾個變量的變化是否依賴一些變量的變化。對于本文，就是研究頻率變量和幾何參數(shù)以及材料特性參數(shù)的變化依賴性。X矩陣中，兩個變量(N維)a和b之間的相關(guān)性可以由如下公式獲得[11]

(8)

式中，Cov(a,b)是向量a和b的協(xié)方差，Var(a)和Var(b)分別表示a和b的方差。ρab也稱為Pearson相關(guān)系數(shù)。

了解輸入?yún)?shù)和頻率的相關(guān)性能之后，作為圓柱殼動態(tài)特性的設(shè)計者，還希望知道設(shè)計參數(shù)對于頻率的影響程度，以更好幫助設(shè)計。ANOVA是數(shù)據(jù)分析的一種方法，它解決的基本問題是，通過數(shù)據(jù)分析弄清楚和研究對象有關(guān)的各因素以及各因素之間相互作用對于該對象的影響。ANOVA的研究可以放在一般線性模型的框架內(nèi)，它給出的假設(shè)是：觀測值獨立并且成正態(tài)分布，所有這些觀測值都有同樣的方差，觀測值的均值可以表示成為某些參數(shù)的線性組合[12]。多元方差分析是ANOVA的拓展，它的主要步驟包括：問題確定，選擇變量，假設(shè)，模型估算，結(jié)果分析和驗證等環(huán)節(jié)。這種技術(shù)已經(jīng)是多元統(tǒng)計分析的重要部分，有很多統(tǒng)計類軟件提供這個功能[13]。本文擬采用該技術(shù)研究幾何參數(shù)以及材料參數(shù)各因素對頻率的影響程度。

從上述過程可以看到，經(jīng)過試驗設(shè)計得到的解空間和性能空間是多維的。當我們考慮的系統(tǒng)超過3維之后，其可視化就會遇到很多困難。對于這個多維數(shù)據(jù)構(gòu)成的設(shè)計空間，本文引入自組織神經(jīng)網(wǎng)絡(luò)(Self-Organizing Facture Map, SOM)技術(shù)進行了研究。SOM是一種無導師的神經(jīng)網(wǎng)絡(luò)學習方法。它最重要的特點是通過自動尋找樣本中的內(nèi)在規(guī)律和本質(zhì)屬性，自組織、自適應(yīng)地改變網(wǎng)絡(luò)參數(shù)與結(jié)構(gòu)，廣泛用于數(shù)據(jù)挖掘以及可視化研究。

SOM網(wǎng)絡(luò)共有兩層，即輸入層和輸出層。輸入層的神經(jīng)元和樣本的維數(shù)相等。輸出層也就是競爭層，通常的形式是二維平面陣。輸出神經(jīng)元常常采用矩形或者六角形側(cè)向連接。在網(wǎng)絡(luò)訓練過程中，競爭獲勝神經(jīng)元和其周圍一定半徑內(nèi)的神經(jīng)元可以依據(jù)特定的形式調(diào)整權(quán)向量。SOM的算法過程和簡要描述如下：

步驟1 初始化 權(quán)向量和學習率初始化，建立初始優(yōu)勝鄰域；

步驟2 競爭 樣本輸入以及計算獲勝節(jié)點；

步驟3 合作 競爭優(yōu)勝神經(jīng)元確定獲勝區(qū)域，確定鄰域節(jié)點；

步驟4 調(diào)整 對優(yōu)勝神經(jīng)元獲勝區(qū)域內(nèi)的神經(jīng)元進行節(jié)點權(quán)值調(diào)整；

步驟5 檢查判斷 檢查學習率是否達到預(yù)設(shè)要求。否則轉(zhuǎn)到步驟2。

SOM技術(shù)通過這兩層網(wǎng)絡(luò)形式能夠把多維數(shù)據(jù)(輸入層)轉(zhuǎn)換到二維平面圖(輸出層)，它既是一種數(shù)據(jù)分析工具，也是一種可視化技術(shù)[14]。通過形成的二維圖形，能夠獲得更多數(shù)據(jù)之間的內(nèi)在關(guān)系。本文擬采用這個技術(shù)分析結(jié)構(gòu)幾何參數(shù)以及材料參數(shù)對于頻率的影響，以及輸入?yún)?shù)之間的關(guān)聯(lián)關(guān)系。

2 數(shù)值算例和討論

2.1 不同邊界下圓柱殼的自振頻率

這里首先使用波傳播方法計算了不同邊界和材料特性情況下，受不同軸向力作用時圓柱殼的頻率，見表2。 對于正交各向異性圓柱殼，引入頻率因子Ω2=ρR2ω2[15]。

表2 不同邊界條件下正交各向異性圓柱殼的頻率Ω，Ω×10-5

對于簡支邊界條件，表2還給出使用經(jīng)典的解析解方法的結(jié)果[16]。對于這種邊界，波傳播方法和解析解的結(jié)果非常接近。對于這兩個方法的比較，Gan等[17]針對環(huán)肋圓柱殼結(jié)構(gòu)還給出理論分析，證明了兩者之間的等價性。

2.2 多元統(tǒng)計和數(shù)據(jù)分析

這里討論一個數(shù)值算例，計算圓柱殼在一定的材料特性和幾何參數(shù)范圍的頻率特性。引入試驗設(shè)計技術(shù)研究這些參數(shù)對于自振頻率的影響。參數(shù)范圍以及試驗設(shè)計的水平列在表3。 假設(shè)不考慮外力，這樣共有7個參數(shù)，每個參數(shù)選取3個水平，采用全因子設(shè)計方式。

表3 材料特性、幾何參數(shù)范圍和水平取值

通過這些參數(shù)的組合，并計及式(1)以及vθ<0.2的約束，可以得到1 782個數(shù)據(jù)。這樣X矩陣的大小就是1 782×10。如果僅考慮最低頻率，則X矩陣的大小就是1 782×8。

根據(jù)式(8)可以得到最低頻率和這些參數(shù)的相關(guān)系數(shù)ρ以及p值見表4。計算中采用顯著性水平α=0.05，p值是否定假設(shè)檢驗中原假設(shè)是否適當?shù)膮?shù)。

從表4可知,L/R和最低頻率的相關(guān)系數(shù)為負值，這說明殼體越長，則其最低頻率越低。殼體越厚(h/R增加)以及軸向模態(tài)越高，頻率也越大。兩個坐標軸方向的楊氏彈性模量和頻率的相關(guān)性相差不大。注意到表中vx和n對應(yīng)的p值并不為零，這說明它們在置信度為0.95時并不是顯著相關(guān)的。

表4 最低頻率Ω1和參數(shù)的相關(guān)性指標以及p值

使用一般線性模型[18]，對于最低頻率進行多元方差分析的結(jié)果見表5。表5中自由度等于變量的水平-1；Seq-SS(連續(xù)平方和)為組間平方和(因子)以及組內(nèi)平方和(誤差)，MS-平方和是連續(xù)平和和除以自由度得到的均方。f是將因子MS除以誤差得到。p值是用以確定某個因子是否顯著，通常和α=0.05進行比較。如果p<0.05，這表示該因子是顯著的。

表5 最低頻率多元方差分析結(jié)果

從表5可知，可以得到7個輸入?yún)?shù)對于最低頻率的重要程度。L/R的重要性達到5 055.44/17 916.62×100%=28.22%. 這說明圓柱殼的長度/半徑比值對頻率的影響最大。其次是軸向模態(tài)m以及厚度因素，分別達到了16.34%和7.05%。其他參數(shù)對于頻率的影響重要性均比較小。

針對表3中的數(shù)據(jù)，因子范圍保持不變，加大因子的水平到5，采用全因子設(shè)計共得到59 375個試驗樣本。以圓柱殼兩端簡支，軸向壓力P=0為例，將這個59 375×10維的數(shù)據(jù)(含中間頻率Ω2和最高頻率Ω3)通過SOM技術(shù)映射到二維空間，可以得到圖2。本文利用SOM Toolbox軟件[19]進行分析求解。該映射圖能夠?qū)⒄桓飨虍愋詧A柱殼的7個輸入?yún)?shù)和3個頻率(即共計10維的高維空間數(shù)據(jù))用二維的形式展示出來，可以方便研究者從整體上了解振動特性。和傳統(tǒng)用曲線表示兩兩變量間的變化方法相比，這里能更加直觀給出高維空間里變量之間的關(guān)系。

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

在圖2中，神經(jīng)元以六角形側(cè)向緊密連接，映射圖中網(wǎng)格大小為47(縱向)×26(橫向)，即59 375個數(shù)據(jù)共投影到這1 222個神經(jīng)元內(nèi)。當考察某屬性映射圖里(比如L/R圖)單個神經(jīng)元特性的時候，需要同時計及該神經(jīng)元在其它映射圖中(對應(yīng)坐標網(wǎng)格)的屬性(比如h/R等)。

在SOM映射過程中的權(quán)向量是多維的，有多種表示權(quán)向量差別的度量方法，主要有基于距離的統(tǒng)一距離矩陣(U-matrix)和基于密度估計的密度度量矩陣(P-matrix)方式。圖2中給出U-matrix結(jié)果以及各個變量的變化情況。U-matrix矩陣中存儲的是對應(yīng)位置節(jié)點和相鄰節(jié)點權(quán)向量之間的距離。從這個圖形能夠依據(jù)距離對原始空間的數(shù)據(jù)點進行聚類分析。本文形成SOM輸入數(shù)據(jù)的是通過全因子分析得到的，每個參數(shù)給出的水平值相同，在U-matrix圖像上顯示的聚類信息并不顯著。

觀察7個輸入變量以及3個頻率參數(shù)的圖像，可以分析如下：

1)最低頻率Ω1(圖2(i))在整個范圍內(nèi)變化均比較小，Ω1遠小于1。在右下角有比較集中的大值，這對應(yīng)于軸向波數(shù)m較大、h/R較大以及殼體較短(L/R較小)的情況。此區(qū)域的Ex值要遠大于Eθ。結(jié)合表5的數(shù)據(jù)可知，Ex的重要程度也大于Eθ。比較L/R和Ω1的分布圖，殼體越短，在其他參數(shù)不變的情況下，頻率就越高。分析這兩個圖從左上角變化到右下角的趨勢可以看出，這和表4中L/R與Ω1的相關(guān)系數(shù)-0.48(負值)是對應(yīng)的(大小變化順序相反)。軸向模態(tài)m對最低頻率Ω1的影響，與L/R的影響類似，不過是正向的，即m越大，則Ω1也越大。這與m和Ω1的相關(guān)系數(shù)0.40(正值)也是符合的。

2)周向波數(shù)n對應(yīng)Ω1的影響并不顯著，除了從相關(guān)系數(shù)可以看出外，對比圖2(i)可知，n值的取值范圍內(nèi)可以得到比較小最低頻率。在最左下角處，即使是n值很大，結(jié)合h/R，L/R值，頻率也比較低。

3)對于各向異性材料，為了獲得較高的基礎(chǔ)頻率，在幾何參數(shù)確定之后，比較圖2(i)的左半幅和圖2(f)與圖2(g)，Ex和Eθ的增加導致頻率增加恰恰是比較小的。Ex和Eθ的配合需要協(xié)調(diào)，并不是越高越好。這對于材料設(shè)計和制造經(jīng)濟性而言是有啟示意義的。

4)圓柱殼的最高頻率Ω3和最低頻率Ω1有相似的變化趨勢。最大值出現(xiàn)在右下角，殼體越短，m越大則越高。不過，在殼體比較薄(h/R小)以及低周波數(shù)n時，也會出現(xiàn)比較大的最高Ω3頻率。對應(yīng)于軸向位移為主的中間頻率，其變化規(guī)律和Ω3相近。

3 結(jié) 論

本文根據(jù)Flügge經(jīng)典殼體理論以及波傳播方法討論了圓柱殼在不同邊界下的自振頻率特性。引入試驗設(shè)計的思想，結(jié)合多元統(tǒng)計分析方法以及自組織映射技術(shù)，深入分析了幾何參數(shù)以及物理特性參數(shù)對于受到軸向壓力作用正交各向異性圓柱殼的頻率影響。本文給出的方法能夠有效探知設(shè)計參數(shù)對于頻率的影響，可以為圓柱殼的動力設(shè)計提供有效手段。從上述過程可以得到如下結(jié)論：

(1) 殼體的長度/半徑比值和最低頻率的相關(guān)系數(shù)為負值。殼體越長則其最低頻率越低。殼體越厚(h/R增加)以及軸向模態(tài)越高，頻率也越大。

(2) 長度/半徑比值對頻率的影響最大。其次是軸向模態(tài)m以及厚度因素，其他參數(shù)對于頻率的影響重要性均比較小。

(3) 通過全因子分析得到的數(shù)據(jù)，在U-matrix圖像上顯示的聚類信息并不顯著。

(4) 最低頻率Ω1在整個范圍內(nèi)變化均比較小且遠小于1。比較集中的大值對應(yīng)于m較大、h/R較大以及殼體較短(L/R較小)的情況。Ex的重要程度也大于Eθ。周向波數(shù)n對應(yīng)Ω1的影響并不顯著，即使是n值很大，結(jié)合h/R，L/R值，頻率也可能比較低。

(5) 對應(yīng)各向異性材料，對于動態(tài)特性設(shè)計，Ex和Eθ的配合需要協(xié)調(diào)，可提高制造經(jīng)濟性。

[1] LEISSA A W. Vibration of shells (NASA SP 288), NASA SP 288[R]. Washington DC, US: Government Printing Office, 1993.

[2] 趙鑫, 張博, 李躍明. 不同邊界條件下截頂正交各向異性圓錐殼的振動和聲輻射研究[J]. 振動與沖擊, 2016,35(3):99-106.

ZHAO Xin,ZHANG Bo,LI Yueming.Influences of boundary conditions on the vibration and sound radiation of a truncated orthotropic conical shell[J].Journal of Vibration and Shock,2016,35(3):99-106.

[3] LIU B, XING Y F, QATU M S, et al. Exact characteristic equations for free vibrations of thin orthotropic circular cylindrical shells[J]. Composite Structures, 2012,94:484-493.

[4] 李學斌. 正交各向異性圓柱殼的穩(wěn)態(tài)動力響應(yīng)分析[J]. 船舶力學, 2007,11(1):79-87.

LI Xuebin. Hairmonic response analysis of orthotropic cylindrical shells[J]. Journal of Ship Mechanics, 2007,11(1):79-87.

[5] LI Xuebin. Study on free vibration analysis of circular cylindrical shells using wave propagation[J]. Journal of Sound and Vibration, 2008,311:667-682.

[6] 馬成良, 張海軍, 李素平. 現(xiàn)代試驗設(shè)計優(yōu)化方法及應(yīng)用[M]. 鄭州: 鄭州大學出版社, 2007.

[7] JOHNSON R A, WICHERN D W. Applied multivariate statistical analysis[M]. 6th. New Jersey: Prentice Hall, 2007.

[8] KOHONEN T. Self-organizing maps[M]. 3rd. Berlin: Springer-Verlag, 2001.

[9] FLüGGE. Stresses in shells[M]. Berlin: Springer-Verlag, 1973.

[10] 季文美, 方同, 陳松淇. 機械振動[M]. 北京: 科學出版社, 1985.

[12] DEGROOT M H. Probability and Statistics[M]. Wokingham, England: Addison-Wesley Publishing Company, 1989.

[13] 高慧璇. 應(yīng)用多元統(tǒng)計分析[M]. 北京: 北京大學出版社, 2005.

[14] MWASIAGI J I. Self organizing maps-applications and novel algorithm design[M]. Rijeka, Croatia: InTech, 2011.

[15] WARBURTON G B, SONI S R. Resonant response of orthotropic cylindrical shells[J]. Journal of Sound and Vibration, 1977,53(1):1-23.

[16] 李學斌. 正交各向異性圓柱殼靜動態(tài)特性分析及比較研究[D]. 武漢: 華中科技大學, 2004.

[17] GAN L, LI X, ZHANG Z. Free vibration analysis of ring-stiffened cylindrical shells using wave propagation approach[J]. Journal of Sound and Vibration, 2009,326:633-646.

[18] TIMM N H. Applied multivariate analysis[M]. New York: Springer-Verlag, 2001.

[19] VESANTO J, HIMBERG J, ALHONIEMI E, et al. Self-organizing map in MATLAB:the som toolbox:[C]//Proceedings of the Matlab DSP Conference. Espoo, Finland, 1999.