基于能量變分原理的梁-圓柱殼耦合系統振動特性分析

聶 睿， 李天勻， 朱 翔， 陳 旭

(1.華中科技大學 船舶與海洋工程學院，武漢 430074； 2.高新船舶與深海開發裝備協同創新中心，上海 200240；3.船舶與海洋水動力湖北省重點實驗室，武漢 430074)

以圓柱殼為主體的結構廣泛存在于船舶與海洋工程、航空航天工程等領域中。其振動特性被國內外學者進行了廣泛研究。受螺旋槳脈動力激勵而產生的推進軸系-艇體整體結構振動是近年來水下航行器振動特性研究重點。推進軸系和艇體會相互傳遞振動，存在著明顯耦合效應，建立兩者耦合模型進行振動特性計算十分必要。但是在以往工作中，水下航行器推進軸系與殼體往往不考慮彼此的耦合作用，僅僅只探討了隔離后軸系或殼體的振動問題[1-5]。如欲對軸-艇耦合系統的振動特性進行更全面分析，對組合結構振動這類問題開展更深入的研究勢在必行。

對于組合結構，通常先將整個結構劃分成多個子結構，再利用結構之間的連續協調條件體現結構間耦合作用。周海軍等[6]利用改進傅里葉級數方法和波傳播法分別求取了梁和圓柱殼的導納，再利用阻抗綜合法，建立了耦合計算模型。Qu等[7]采用廣義變分法模擬梁、殼等結構單元，再基于Hamilton原理，計算了梁-球殼-柱殼組合結構振動響應。郭亮等[8]采用解析法對錐柱組合殼結構及多跨梁結構分別建模，并在軸承處建立橫向連續條件，形成了有效的系統橫向振動解析求解模型。李攀碩[9]分別利用解析法和有限元法計算子結構導納，再通過頻響函數綜合法建立了軸與加筋殼體的耦合系統模型。Merz等[10]采用有限元和邊界元方法建立水下加筋殼體模型，并通過簡化的彈簧阻尼系統與用桿結構模擬的推進軸系耦合，計算了系統振動與聲特性。

總的來說，已有文獻中軸-艇耦合系統振動分析方法主要有數值法和解析法兩種。有限元等數值方法可以用于分析各類耦合系統結構振動問題，但它們難以從機理上對系統振動行為提供合理解釋。而解析法能給出較精準的解析結果，易于從機理上解釋某些研究現象，但其往往又存在求解范圍有限的缺點。而實際應用中存在大量復雜的邊界約束，解析法針對不同邊界需要配置不同位移函數和邊值條件，難以建立統一的動力學分析模型。

Rayleigh-Ritz法作為一種基于能量變分原理的(近似)解析法，由于其假定振型函數滿足結構位移邊界條件和振動方程，由此得到的函數組合往往接近真實振型，因此有很高精度；同時，由于假定振型函數之間具有正交性，實際運算成本極低。此外，Rayleigh-Ritz法還十分容易處理附加質量或附加剛度。只需添加附加動能項或附加應變能項，即可在不改變結構原振型函數的前提下，求取新系統振動特性。但傳統Rayleigh-Ritz 法的限制是：對復雜邊界條件，難以選取恰當的假設振型函數。Li[11]提出了一種改進傅里葉級數，該級數通過在傳統傅里葉級數增加輔助多項式函數，保證了梁結構在邊界上的連續性，結合兩類彈簧模擬一般邊界的思想，成功建立了梁的自由振動計算模型。隨后，國內外學者借鑒該研究思路，結合Rayleigh-Ritz法，研究了梁、板、殼等各類結構的振動特性[12-15]。這些研究證明，引入改進傅里葉級數的Rayleigh-Ritz法，在面對復雜邊界條件時，可構建準確統一的分析模型。

當對艇體中段進行截斷處理后，軸-艇耦合結構可簡化為梁-圓柱殼耦合模型[16]。在對軸-艇系統進行合理簡化后，本文將以梁-圓柱殼結構為研究對象，基于能量變分原理，采用改進傅里葉級數統一描述梁結構和圓柱殼結構的振動位移；通過人工約束彈簧組模擬各類邊界條件，耦合彈簧組來滿足梁殼之間的連續條件，建立了梁-圓柱殼耦合系統振動計算模型。如前文所述，本文方法是對Rayleigh-Ritz 法的改進，故所得計算模型不僅可用于圓柱殼取各類經典邊界條件下的梁-圓柱殼耦合結構，同時由于級數良好的收斂性亦只需極低運算成本即獲得準確結果。通過與有限元軟件的計算結果對比，驗證了本文計算模型正確性。最后，分析了系統在部分參數改變時，振動響應性能的變化規律。

1 梁-圓柱殼耦合系統振動計算模型

1.1 模型描述

當進行截斷處理后，軸-艇耦合系統可合理簡化為梁-圓柱殼結構。即推進軸系可簡化為帶集中質量點的彈性梁，艇體可簡化為兩端受彈簧約束的圓柱殼，中間軸承等結構可簡化為連接梁結構與殼結構的耦合彈簧組。

梁-圓柱殼耦合模型如圖1所示，全局坐標系為直角坐標系O-xyz，坐標原點O位于梁結構左端點處；圓柱殼的長度和截面半徑分別為L和R，殼體中面軸向位移、周向位移、徑向位移及繞y軸轉角分別為ucy、vcy、cy、βcy(βcy=?cy/?x)；梁通過兩處節點與圓柱殼連接，每段長度分別為l1和l2，梁的自由端有一附加集中質量點，其質量為m；kL和kR分別為左連接點與右連接點的空間點對點耦合彈簧組(由三向線位移彈簧kx、ky、kz和兩向轉角彈簧Ky、Kz構成)；圓柱殼兩端引入兩組連續線分布的約束彈簧用于模擬各類邊界(由三種位移彈簧ku、kv、k和一種轉角彈簧Kβ構成)；fx、fy、fz為作用在質量點處的簡諧激勵力。

圖1 梁-圓柱殼耦合模型

1.2 理論推導

如前文所述，與直接求解結構控制微分方程的方法不同，本文將基于能量變分原理建立梁-圓柱殼耦合模型。耦合結構能量泛函Π的形式如下

Π=Tb-Ub+b+Tcy-Ucy-Uc

(1)

式中:Tb為梁結構的總動能;Ub為梁結構的總勢能;b為梁結構上外部激勵力的外力功；Tcy為圓柱殼結構的總動能;Ucy為圓柱殼結構的總勢能；Uc為結構耦合彈簧儲存的彈性勢能。

本文計算模型中梁結構為含集中質量點的彈性歐拉梁，其自身局部坐標系為xbybzb，坐標原點Ob位于梁左端，如圖2所示。

圖2 含集中質量點的梁結構

梁結構的總動能Tb由梁結構自身動能Tb,s和集中質量點的動能Tb,m構成

Tb=Tb，s+Tb，m

(2)

式中:ub，vb和b分別為梁結構在縱向(x)，橫向(y)和垂向(z)的位移;ρb和Ab分別為梁的密度和截面面積;m為集中質量點的質量。集中質量點附加在梁結構左端面上，即梁結構坐標系中原點處。

由于本文計算模型中梁結構無外部約束，其勢能Ub僅由自身應變勢能構成，形式如下

(3)

式中:Eb為梁材料的楊氏模量;Ib,y和Ib,z是梁的截面慣性矩。

集中質量點處作用有縱向、橫向及垂向三個方向的簡諧點激勵力，由此產生的外力功b如下

b=fxub(xb,t)|xb=0+fyvb(xb,t)|xb=0+

fzb(xb,t)|xb=0

(4)

圓柱殼結構采用柱坐標系(R,θ,xcy)，xcy和θ分別代表軸向和周向坐標,坐標原點Ocy位于圓柱殼左端面圓心處。圓柱殼結構如圖3所示。

圖3 圓柱殼結構示意圖

采用Ressiner薄殼理論計算圓柱殼的動能與應變能。動能Tcy具體形式如下

(5)

式中:ρcy為殼體密度;hcy為殼體厚度。

在圓柱殼兩端引入兩組沿端面圓周方向均勻分布的線性彈簧組用于模擬各類邊界條件。因此，圓柱殼的總勢能Ucy由自身應變能Ucy,p和存儲在邊界彈簧的彈性勢能Ucy,e組成。具體計算形式如下

Ucy=Ucy,p+Ucy,e

(6)

式中:Ecy和μcy分別為圓柱殼的楊氏模量和泊松比;ku0，kv0，k0和Kβ0為圓柱殼左端約束彈簧;ku1，kv1，k1和Kβ1為圓柱殼右端約束彈簧。連續線分布的位移彈簧和轉角彈簧的單位分別為N/m2和N/(rad·m)。

梁結構和圓柱殼結構將通過分布彈簧系統實現耦合。每對耦合節點之間的耦合彈簧組均包含縱(x)向、橫(y)向和垂(z)向三種線位移彈簧與y向及z向兩種轉角彈簧。圓柱殼的振動位移分量將從局部柱坐標轉換到全局笛卡爾坐標下，以便構造耦合彈簧組的能量泛函。梁結構和圓柱殼結構的振動位移如圖4所示。

圖4 梁和圓柱殼的振動位移

vcy=vcycosθ-cysinθ,cy=cycosθ+vcysinθ

(7)

式中，vcy和cy為圓柱殼結構轉化到全局坐標系下的振動線位移分量。

通常，薄殼理論在解決殼體問題時，殼體模型僅考慮三向線位移和繞y軸旋轉位移。但本文工作中為保證梁、殼結構的耦合振動分量數目相當，需要再加入繞z軸旋轉位移αcy。引入額外旋轉位移并進行轉化后的形式如下

Rotcy,y=βcycosθ-αcysinθ

Rotcy,z=αcycosθ+βcysinθ

(8)

式中，Rotcy,y和Rotcy,z為圓柱殼結構轉化到全局坐標系下的旋轉位移分量。

如前文所述，本文計算模型含kL、kR兩組梁殼耦合彈簧組。耦合彈簧存儲的彈性勢能Uc具體形式如下：

Uc=Uc,L+Uc,R

(9)

式中:Uc,L為左耦合彈簧組儲存的彈性勢能;kL,x，kL,y和kL,z分別代表x向，y向和z向的線位移耦合彈簧;KL,y和KL,z分別代表y向和z向的旋轉位移耦合彈簧;(xL，θL)為左耦合彈簧組上圓柱殼耦合節點在柱坐標系上的位置坐標。下標R代表右耦合彈簧組對應的變量，故不贅述。

本文將引入改進傅里葉級數作為梁結構振動位移函數及圓柱殼結構軸向振動位移函數。同時，由于完整的圓柱殼擁有軸對稱外形，其周向振動位移函數可由傳統的兩組單傅里葉級數描述。具體形式如下。

梁結構：

(10)

式中:Mb為位移函數的截斷項數;Am,u和Bl,u是梁軸向位移容許函數中的未知傅里葉系數;eiωt被作為時間簡諧因子引入來描述梁在不同時刻的位移;Am,v、Bl,v、Am,和Bl,是梁橫向及垂向位移容許函數中的未知傅里葉系數。式(10)中的位移容許函數即為一種改進傅里葉級數，其中，λmb=mπx/(l1+l2)，λlb=lπx/(l1+l2)。它在傳統的單余弦傅里葉級數基礎上引入輔助正弦函數，使得容許函數滿足軸向振動位移一階導數或彎曲振動位移三階導數在梁上的任意一點存在且連續，保證滿足位移邊界條件，同時加速了級數自身的收斂速度。

殼結構：

(11)

式中:Ncy、Mcy為位移函數截斷項數;Amn等為未知傅里葉系數，λmcy=mπx/L，λlcy=lπx/L。

將式(1)～(11)聯立，并根據最小勢能原理，對能量泛函求取極值

(12)

式中，q為未知系數向量，q=[Am,u，Bl,u，…，Amn，al，…，Cmn，cl]。

最終得到的梁-圓柱殼耦合系統振動方程組矩陣化形式如下

([K]-ω2[M]){q}={F}

(13)

式中:[K]和[M]分別為耦合系統的剛度矩陣和質量矩陣;ω為圓頻率;{F}為集中質量點處激勵力的廣義力向量。廣義力向量取0時，亦可由式(13)求取系統的固有振動特性。

2 計算模型驗證

為驗證本文方法收斂性和計算模型準確性，下面將從傅里葉級數的截斷項數出發，分析方法收斂性。同時，將本文計算模型計算結果與有限元軟件計算結果進行比對，驗證本文計算模型正確性。有限元模型共有31 566個單元。其中，殼單元31 400個，梁單元165個，質量點單元1個。

本文計算模型如前文所述，具體參數如表1所示。

表1 耦合模型參數

2.1 傅里葉級數截斷項數收斂性

改進傅里葉級數的截斷項數對計算結果精度有直接影響，其展開的數目越多，計算結果越趨于真實值，但過大的截斷項數會降低求解效率。一般來說，截斷項數達到一定值后，計算結果即可穩定(收斂)。本文以圓柱殼為自由邊界條件(殼體兩端約束彈簧組取值為0)的梁-圓柱殼模型為例，探討計算結果收斂性隨Mb、Mcy的變化。而圓柱殼周向截斷數Ncy并不影響計算結果的收斂，僅僅會影響固定頻率范圍內結果是否會缺少某些模態。

驗證模型參數見表1，Ncy=6，圓柱殼上的連接點為(xL，θL)=(0.5，π), (xR，θR)=(1.5，π)，梁的連接點為xb=l1和xb=l1+l2處。前8階固有頻率隨截斷項數Mb、Mcy的變化結果如表2所示。

表2 固有頻率與截斷項數關系

從表2可以看出，隨著截斷項數的增加，本文方法的計算結果很快就趨于收斂了，而與有限元法(finite element method,FEM)軟件計算結果的對比證明了本文方法有很高精度。綜合計算成本和精度，本文其余計算工作均取截斷項數為10。

2.2 彈簧模擬經典邊界條件

如前所述，本文采用兩類彈簧模擬圓柱殼邊界。彈簧剛度增大，儲存的彈性勢能也越大，圓柱殼邊界處約束越強。但由剛度增大而產生的邊界約束效應必定會趨于完全約束(位移或轉角)。因此，取值達到一定大小的彈簧組可以用于模擬各類經典邊界，例如固支邊界即應為所有彈簧組取足夠大值。下面將驗證該假設。

仍取前文計算模型，但假定ku0、kv0、k0、Kβ0值取0，ku1、kv1、k1、Kβ1值取Q。顯然，此時耦合模型圓柱殼左端邊界條件仍為自由。圖5將顯示耦合系統計算模型固有頻率隨彈簧剛度取值Q的變化關系。

圖5 邊界彈簧組剛度取值Q對左端自由系統固有頻率影響

從圖5可知，耦合計算模型的前10階固有頻率隨彈簧剛度增加都有著明顯的“緩慢增加—快速增加—穩定不變”趨勢，這證明了前文彈簧剛度收斂性假設正確可靠。此外，雖然不同階固有頻率的曲線變化規律略有不同，但在剛度值取1×1012左右均可以達到穩定。因此，本文其余工作均認為彈簧組剛度值取1×1012時，其對應方向上的殼體位移或轉動為完全約束。由于本文計算模型在y(橫)和z(垂)向是基本對稱的，因此前10階固有頻率為兩兩一組，成對出現。

下面假定ku0、kv0、k0取值為Q，Kβ0為0，ku1、kv1、k1、Kβ1取值為1×1012。此時，圓柱殼右端所有位移均被完全約束，即對應固支條件。圖6為此條件下耦合系統計算模型固有頻率隨彈簧剛度取值Q的變化關系。

圖6 邊界彈簧組剛度取值Q對右端固支系統固有頻率影響

從圖6可以看出，在右端為固支條件時，增加圓柱殼左端線彈簧組的剛度，耦合系統固有頻率也有類似的規律，彈簧剛度取值達到1×1012左右也趨于收斂。對應地，ku0、kv0、k0取1×1012且Kβ0取0，這一條件下，圓柱殼左端應視為簡支邊界條件。為進一步證明該假設的正確性，下面將這兩種邊界組合下的耦合系統固有頻率計算值與有限元計算結果進行對比。結果如表3所示。

表3 本文方法計算結果與有限元結果對比

從表3的結果可以看出，本文方法和FEM計算結果十分接近，使用彈簧模擬邊界的計算方法是正確可靠的。

2.3 計算模型正確性

在截斷項數收斂性和彈簧剛度收斂性的探討中，本文計算模型正確性已經從固有頻率的比對中得到了初步驗證。下面將進一步從對應模態振型及受迫振動響應的對比角度出發，驗證本文計算模型正確性。

這里的驗證模型取圓柱殼邊界條件為簡支-簡支組合，其余參數保持不變。首先，對比梁結構z(垂)向振動的前5階模態振型，結果如表4所示。

從表4可以看出，本文計算模型的梁結構z向振動前2階模態，圓柱殼基本可視為剛體結構，不參與系統的振動。進一步觀察得知，這2階模態中，耦合點處于梁結構振型的不動點處，因此梁結構和圓柱殼結構難以通過耦合點進行振動傳遞。這也可以說明圖6中，為何圓柱殼約束增強，但耦合系統前4階(包括2階垂向和2階橫向)振動的固有頻率幾乎沒有變化。因為此時梁結構與圓柱殼結構可以視為解耦狀態，圓柱殼可等效為剛性基礎，其約束的增強并不能對系統產生明顯影響。而后面的模態中，圓柱殼結構參與了整個系統的振動，因此其約束的增強會使得整個系統勢能增加，固有頻率也會隨之增大。

然后，對比本文方法和FEM的受迫振動響應。在集中質量點處加z向單位幅值激勵力，掃頻范圍為1～200 Hz，提取梁結構在質量點、左耦合點和右耦合點處的z向位移響應及殼結構在左耦合點和右耦合點處的z向位移響應。取參考位移A0=1×10-12m，A為提取點的z向位移響應幅值，則振動位移級響應為P=20lg(A/A0)。

從圖7可以看出，本文計算模型與FEM軟件計算結果十分吻合，這進一步驗證本文方法的正確性。在190 Hz附近兩種方法得到的位移級響應峰值位置存在微小偏差，表明兩者計算得到的系統固有頻率存在差別，這也與表3中內容相符。

(a) 梁結構

3 參數對耦合系統振動特性影響

3.1 耦合彈簧組剛度對系統振動特性影響

由前文理論推導部分可知，耦合彈簧組儲存的彈性勢能是構成整個系統的重要部分。顯然，耦合彈簧剛度值將直接影響其儲存的彈性勢能大小，進一步影響耦合系統的振動特性。

在前文的計算中，耦合彈簧組選取了較大的耦合剛度值1×108。參照邊界彈簧收斂性討論部分的內容，這種情況下梁結構和圓柱殼結構耦合點之間應該為一種“強耦合”的狀態，即梁節點與圓柱殼節點的對應位移嚴格相等。當耦合彈簧的剛度減弱時，耦合點之間可能會轉化為一種“弱耦合”的狀態。當這種狀態發生的時候，耦合系統振動特性應該會產生明顯的變化。下面將保持系統其余參數不變且圓柱殼取簡支-簡支邊界條件，通過改變耦合彈簧組的線彈簧剛度取值，驗證這一猜想。令kL,x=0，其余各耦合線彈簧剛度采用相同取值。

從圖8中可以看出，與邊界彈簧隨剛度值增加的變化規律類似，耦合彈簧剛度值增大到一定量之后，系統固有振動特性將會趨于穩定。此時，耦合剛度的改變對系統振動性能影響可以忽略不計，具體表現為剛度取值為1×108和1×1010的響應曲線有著高度重合性。而前文中關于“強弱”耦合狀態轉變時耦合系統將產生振動特性變化的猜想也在圖中得到了明顯體現。系統振動特性的變化可以分為兩個方面。首先，從圖8和圖9(a)、9(b)可以看出，在“弱耦合”階段，梁結構與殼結構在對應耦合點的位移級響應曲線雖然變化趨勢基本吻合但存在明顯“落差”，說明由梁結構質量點處激勵力產生的振動在耦合點處產生了明顯衰減；隨著耦合剛度增加，達到“強耦合”階段后，對應耦合點處兩者位移級響應曲線基本重合了，表明此狀態下兩者振動傳遞幾乎沒有損失。其次，在耦合剛度較小的“弱耦合”階段，隨著耦合彈簧剛度增加，系統的勢能增大，固有頻率整體有增大的趨勢，在圖8中表現為受迫振動位移響應共振峰整體右移，直至達到“強耦合”階段后，共振峰不再產生明顯移動。

圖8 耦合彈簧剛度取值變化時梁及殼結構左耦合點處

圖9 不同耦合彈簧剛度值下梁及殼結構耦合點處位移級響應對比

3.2 旋轉耦合彈簧對系統振動特性影響

以上工作僅考慮梁-殼結構之間的線位移耦合，但實際上兩者之間還存在耦合彎矩作用，即還應考慮本文耦合彈簧組中旋轉彈簧的影響。仍然取圓柱殼為簡支-簡支條件下的計算模型，而各耦合彈簧如存在則剛度值取1×108。表5為各耦合彈簧對系統固有振動特性影響。

如表5所示，如考慮梁-殼結構之間的耦合彎矩作用，即添加旋轉耦合彈簧后，由于耦合彈性勢能的增加，系統的固有頻率將會增大。此外，由于本文模型在垂(z)向及橫(y)向上無耦合，因此耦合旋轉彈簧也僅會對對應方向上的振動產生影響。

表5 不同耦合彈簧組合下系統前10階固有頻率

3.3 集中質量點質量對系統振動特性影響

由本文計算模型推導部分可知，附加集中質量的大小將直接影響梁結構的附加動能項，進而引起系統振動性能變化。下面分別取耦合彈簧剛度值為1×104和1×108且圓柱殼為簡支-簡支條件下的計算模型，研究僅改變附加質量大小時耦合系統振動性能變化。

從圖10(a)及10(b)中可以看出，在本文探討的變化范圍內，集中質量點質量的增加對梁結構的振動位移響應有一定的抑制作用，響應曲線有下移的趨勢；此外，當集中點質量增大時，附加動能項能量增加，系統固有頻率減小，在圖中反映為共振峰的左移。實際上，本文質量點為螺旋槳的簡化，其質量大小的改變可視為由螺旋槳與水耦合作用產生的附連水質量效應，而本文計算結果也與附連水導致含槳系統固有頻率降低的實際情況相符。

圖10 不同耦合剛度值下集中質量點質量變化對結構左耦合點處位移響應影響

4 結 論

本文基于Rayleigh-Ritz法，引入改進傅里葉級數作為結構位移函數，結合邊界彈簧組和耦合彈簧組，建立了梁-圓柱殼耦合動力學計算模型。通過改變邊界彈簧組剛度值即可模擬各種經典圓柱殼邊界條件，構建了適用范圍廣的動力學分析模型。與有限元法(FEM)軟件計算結果的對比，驗證了本文方法在計算固有特性及受迫振動響應時的正確性。隨后，進一步分析了耦合剛度值和集中點質量等參數對梁-圓柱殼耦合系統振動特性的影響，得出了以下結論：

(1) 針對本文研究的梁-圓柱殼耦合系統，耦合系統的前4階模態可視為支撐于剛體上的梁結構彎曲振動，梁結構與圓柱殼結構可視為解耦狀態。此時，殼體邊界約束條件變化對系統固有頻率幾乎沒有影響。

(2) 耦合彈簧剛度值處于較小的范圍時，由于對耦合節點處對應方向的連續性約束較弱，梁殼耦合結構將在耦合節點處出現位移不連續的現象。此時，兩者的振動位移響應曲線雖然趨勢相同，但響應幅值出現落差，表明振動在結構中的傳遞產生了衰減。

(3) 此外，由于集中質量點附加質量對結構振動動能項的影響，梁端點處的集中質量點質量對耦合系統振動特性有一定的影響，系統會出現一定的固有頻率改變(響應曲線共振峰偏移)，對應測點處的振動位移響應結果也會產生一定的變化。