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

汪志強， 李學斌， 黃利華

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

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

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

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

1 理論分析

1.1 頻率求解

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

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

設圓柱殼的材料主軸和坐標系主軸重合，則滿足如下關系式

vxEθ=vθEx

(1)

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

(2)

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

(3)

根據波傳播方法，可假設位移函數形式為

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為軸向波數；n為環向波數；t為時間；ω為圓頻率(1/s)；i2=-1。將這個位移函數代入式(2)，可以得到如下矩陣形式的方程

(5)

矩陣A的元素如下

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

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

(6)

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

表1 梁彎曲振動特征函數和波數解

1.2 多變量分析

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

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

(7)

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

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

(8)

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

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

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

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

步驟1 初始化 權向量和學習率初始化，建立初始優勝鄰域；

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

步驟3 合作 競爭優勝神經元確定獲勝區域，確定鄰域節點；

步驟4 調整 對優勝神經元獲勝區域內的神經元進行節點權值調整；

步驟5 檢查判斷 檢查學習率是否達到預設要求。否則轉到步驟2。

SOM技術通過這兩層網絡形式能夠把多維數據(輸入層)轉換到二維平面圖(輸出層)，它既是一種數據分析工具，也是一種可視化技術[14]。通過形成的二維圖形，能夠獲得更多數據之間的內在關系。本文擬采用這個技術分析結構幾何參數以及材料參數對于頻率的影響，以及輸入參數之間的關聯關系。

2 數值算例和討論

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

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

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

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

2.2 多元統計和數據分析

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

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

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

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

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

表4 最低頻率Ω1和參數的相關性指標以及p值

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

表5 最低頻率多元方差分析結果

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

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

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

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

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

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

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

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

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

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

3 結 論

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

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

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

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

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

(5) 對應各向異性材料，對于動態特性設計，Ex和Eθ的配合需要協調，可提高制造經濟性。

[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] 李學斌. 正交各向異性圓柱殼的穩態動力響應分析[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] 馬成良, 張海軍, 李素平. 現代試驗設計優化方法及應用[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] 高慧璇. 應用多元統計分析[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] 李學斌. 正交各向異性圓柱殼靜動態特性分析及比較研究[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.