基于Bezier函數方法的球殼振動特性分析

池旭帆，束永平

(東華大學 機械工程學院，上海 201620)

彈性開口球殼結構具有良好的物理特性，因此被廣泛應用于航天返回艙、潛水器、管道封頭、高壓容器等裝置。在實際工作環境中，由彈性球殼構成的機械結構需要面臨復雜多變的外界環境對結構的振動沖擊和負載，當激勵頻率接近球殼結構固有頻率時，可能會導致結構發生共振甚至被破壞，因此研究彈性開口球殼結構的振動特性對保障相關結構的安全具有重要意義。

針對板殼類結構的振動特性，國內外學者分別采用Ritz法[1-3]、有限元法、區域能量分解法[4-6]、動力剛度法[7]等方法相繼開展研究。Gautham等[8]采用半解析殼單元對層合球殼結構的軸對稱和非軸對稱自由振動特性進行了研究。Jin等[9]針對圓柱、圓錐、球殼類結構建立了基于一階剪切變形理論的改進傅里葉解，采用Rayleigh-Ritz方法對經典約束和彈性約束下結構的頻率和模態振型進行求解。Ye等[10-11]采用切比雪夫多項式作為位移容許函數，對任意邊界條件下層合開口圓柱殼、錐殼和球殼結構的自由振動特性進行了研究。Qu等[12]采用區域能量分解法對回轉錐殼、柱殼和球殼結構在任意邊界條件下的振動特性進行了研究。高晟耀等[13]利用改進傅里葉級數表示位移函數以消除邊界的不連續性，并利用彈簧參數法來模擬一般邊界條件，采用Ritz法求解分析一般邊界條件下中等厚度功能梯度球環結構的自由振動特性，探討了相關參數的影響。李善傾等[14]采用準Green函數方法，將夾支任意形狀底扁球殼自由振動問題化為第二類Fredholm積分方程，通過離散化方程對扁球殼的固有頻率進行了求解。

現有文獻針對殼體結構振動特性展開研究時，多采用引入輔助函數的方法來避免殼體端點位置的不連續問題，同時布置邊界彈簧，通過對不同彈簧的適當賦值來實現一般邊界條件的模擬[15]，這種方式雖然保證了位移函數能夠滿足連續性要求，但是計算量較大，且在計算前必須先開展彈簧剛度值與一般邊界條件的對應關系研究。為避免端點處的不連續問題，同時能夠有效模擬一般邊界條件，本文采用分段Bezier函數構建的位移函數進行開口球殼結構受迫振動特性分析，旨在為一般邊界條件下開口球殼結構的振動問題提供一種新型有效的思路，進一步豐富球殼受迫振動的半解析求解方法。

1 理論方法

1.1 開口球殼結構幾何模型

開口球殼的幾何模型及坐標系如圖1所示。以球殼中面為參考面建立坐標系，其中φ、θ、z分別表示母線方向、周向和法向坐標。設球殼中面曲率半徑為R，開口角度為φ0，在z軸方向的厚度為h。

圖1 開口球殼結構幾何模型Fig.1 Geometric model of open spherical shell structure

1.2 運動學關系與應力合力

根據一階剪切變形理論，開口球殼上任意一點的位移分量可表示為

(1)

式中：u、v和w分別為結構中面上任意點在φ、θ和z方向的位移；Ψφ和Ψθ分別為中面法線繞φ和θ方向的轉角；t為時間。基于球殼中面應變-位移關系得到如下關系式：

(2)

式中：εφ、εθ為中面應變；γφθ、γφz、γθz為中面切應變；χφ、χθ、和χφθ為中面曲率變化；R1=Rsinφ。對于各向同性材料，球殼應變能UV可表示為

(3)

式中：E為材料彈性模量，μ為泊松比。球殼的動能TV可表示為

(4)

式中：ρ為材料密度。當球殼結構受簡諧力作用時，假設外部載荷全部作用在球殼中面位置處，沿φ、θ、z方向的簡諧力分別為Fφ、Fθ和Fz，此時外力做功We可表示為

(5)

1.3 位移函數

在利用能量變分法對結構的振動特性進行求解過程中，結構位移函數的構造對計算結果的精度具有重要意義。為便于適應不同邊界條件的情況，選用Bezier曲線來構造位移函數。Bezier基函數的數學表達式如下：

(6)

式中：M為Bezier函數截斷階數，m=0, 1, …,M；x∈[0, 1]。Bezier曲線在端點處的n階導數僅與相鄰的n+1個控制點有關, 滿足：

Bi, M(0)=Bi, M(1)=0，i=1, 2, …,M-1

B0, M(0)=BM, M(1)=1

B0, M(1)=BM, M(0)=0

B-1, M-1(x)=BM, M-1(x)=0

B′m, M(x)=M[Bm-1, M-1(x)-Bm, M-1(x)]

利用Bezier函數這種局部性質，在改變邊界條件時，只需對個別控制點進行相應的條件約束就可以適應不同的情況，因此邊界條件的設置較為方便，并且能有效避免球殼頂點處的不連續問題。相較于引入輔助函數來避免頂點不連續問題的方法而言，選用Bezier函數大大減小了計算量。同時，考慮高階項需要耗費更多的計算時間，為了在控制點數量不變的情況下，降低函數的最高階數，從而減少計算量。采用三段Bezier函數表示球殼位移分量沿母線方向的變化，用傅里葉級數表示球殼位移分量沿周向方向的變化。不同方向的位移可表示為

(7)

式中：k為球殼沿φ方向的分解子殼段；n為周向波數，當n=0時，可以得到軸對稱模態；akmn、bkmn、ckmn、dkmn、ekmn為5項位移函數的展開系數；j為虛數單位。

如使三段Bezier函數滿足應變位移關系式，位移曲線需滿足一階導數連續, 即令

(8)

式中：ι=a,b,c,d,e；k=0, 1。

1.4 振動特性方程的求解

球殼結構的整體能量泛函可表示為

(9)

將式(3)～(5)代入式(9)進行變分運算，并沿一個周期積分，可得到球殼的振動微分方程。當球殼自由振動時，We=0，球殼結構的自由振動微分方程可表示為

(K-ω2M)G=0

(10)

式中：K、M和G分別為剛度矩陣、質量矩陣及未知系數矩陣，G=[akmn，bkmn，ckmn，dkmn，ekmn]T；ω為圓頻率。通過求解式(10)即可求得球殼結構的固有頻率。

當球殼受簡諧激勵產生受迫振動時，強迫振動微分方程可表示為

(11)

式中：F為外界簡諧激勵經變分后得到的載荷向量；K*為阻尼力做功時的剛度矩陣，在考慮復阻尼影響的情況下，K*=(1+jη)K；ωc為激勵圓頻率。在激勵頻率給定時，可求得位移展開項的系數向量為

(12)

將位移展開項的系數向量代回各位移函數，即得到在給定激勵下球殼上任一點的位移響應幅值。

1.5 邊界條件

利用Bezier曲線的局部性，要控制曲線在端點處滿足邊界條件，只要對少數幾個系數進行適當賦值即可。根據Bezier曲線特性可知：

球殼底端固支時，在φ=φ0處，u=v=w=Ψφ=Ψθ=0，可得：

a2Mn=b2Mn=c2Mn=d2Mn=e2Mn=0

球殼底端簡支時，在φ=φ0處，u=v=w=Ψθ=Mφ=0，Mφ為φ方向的力矩，可得：

2 自由振動特性分析

2.1 收斂性和有效性驗證

2.1.1 收斂性驗證

圖2 開口球殼軸對稱模態的收斂性Fig.2 Convergence of axisymmetric modes of open spherical shell

2.1.2 有效性驗證

開口球殼在不同邊界條件下的軸對稱無量綱頻率的計算結果與已有計算結果的對比如表1所示，球殼幾何尺寸同上文算例。由表1可知，前5階軸對稱無量綱頻率計算結果與文獻[16]吻合良好，最大相對誤差為7.09%。由此說明本文所用方法適用性強，在不同邊界條件下均具有良好的精度。

表1 不同邊界條件下球殼軸對稱無量綱頻率對比

開口球殼除了軸對稱模態以外還存在非軸對稱模態，在承受非軸對稱載荷情況下，非軸對稱模態也可能被激發出來。底端固支條件下球殼的非軸對稱(n≠0)第一階頻率參數對比如表2所示，其中球殼開口角度為φ0=15°，λ為球型扁殼的無量綱幾何參數，λ4=12(1-μ2)R2φ04/h2。由表2可知，隨著周向波數n的增加，無量綱頻率增大，本文計算結果與文獻[16]結果的最大相對誤差為13.95%，與ABAQUS計算結果的最大相對誤差為15.62%。隨著球殼幾何參數λ與周向波數n的增大，無量綱頻率相對誤差增大，這可能與周向波數增大后模態振型的復雜變化以及模型簡化方式的不同有關。由表1和2可知，本文方法得到的開口球殼軸對稱和非軸對稱無量綱頻率計算結果與參考數據均具有較高一致性，這說明基于Bezier曲線計算開口球殼結構的振動特性是有效的。

表2 底端固支條件下球殼非軸對稱無量綱頻率對比

2.2 結構參數對頻率影響

固支邊界條件下半球殼(φ0=90°)的前4階軸對稱模態對應的無量綱頻率與厚度的關系如圖3所示，其中，實線為本文方法計算數據，虛線為有限元計算得到的參考數值。由圖3可知，固支邊界條件下半球殼固有頻率隨厚度增加而增大，且球殼的高階固有頻率受厚度影響較大。

圖3 半球殼前4階頻率隨厚度變化曲線圖 Fig.3 First four frequencies of hemispheric at different thickness

固支邊界條件下開口球殼無量綱頻率與開口角度的關系如圖4所示，其中，實線為本文方法計算數據，虛線為有限元計算得到的參考數值。由圖4可知，開口球殼固有頻率隨開口角度增大而減小，且固有頻率在開口角度較小時變化較大。

圖4 開口球殼前4階頻率隨角度變化曲線Fig.4 First four frequencies of open spherical shells at different angles

3 受迫振動響應

球殼的幾何參數和材料性質同第2節算例，邊界條件為周邊固支，研究受簡諧力作用下的球殼受迫振動。為研究結構內摩擦對于響應損耗的影響，本文基于復阻尼理論，采用復數彈性模量來表示阻尼力的大小以及阻尼力與應變之間的相位差，E*=(1+jη)E，其中η表示復阻尼，本小節取η=0.001，掃頻范圍f=10～2 000 Hz，步長Δf=10 Hz，響應點為球殼頂點。選取兩種常見簡諧激勵形式，分別為作用在頂點處的法向點載荷以及作用在中面上的法向面載荷。

3.1 模型驗證

在兩種載荷作用下，本文方法和有限元方法的球殼結構頂點位置的法向位移響應曲線對比如圖5所示。在考慮復阻尼影響時，球殼位移響應曲線在固有頻率值附近出現多個大小不一的波峰，這是由復阻尼的特性所導致的，在研究球殼共振特性時，主要考慮峰值較大位置處的共振情況。由圖5可以看出，當開口球殼受徑向載荷時，兩種方法計算得到的振動位移響應值峰值產生一定偏移，在面載荷作用下兩者最大偏移相對誤差小于5%，在點載荷作用下兩者最大偏移相對誤差小于2%，因此本文所用方法計算得到的振動位移響應值與有限元方法計算得到的振動位移響應值擬合良好。同時，從圖5還可以看出，在低頻段位移響應波峰較大，隨著激勵頻率的增大，位移響應峰值逐漸減小，這也是在實際工程中總是關注結構基頻的原因。

(a) 面載荷作用下

3.2 結構參數對受迫振動影響

在兩種載荷作用下不同開口角度球殼結構頂點的法向位移響應曲線如圖6所示。從圖6可以看出：在結構參數不變的情況下，球殼受點載荷作用時頂點位置處產生的法向位移響應大于受面載荷作用時；隨球殼開口角度的增大，球殼位移響應曲線向著低頻方向發生移動；在面載荷作用下，隨著開口角度的增大，最大共振峰值增大；而在點載荷作用下，開口角度小的球殼最大共振峰值較大。

(a) 面載荷作用下

在兩種載荷作用下不同厚度半球殼結構頂點的法向位移響應曲線如圖7所示。從圖7可以看出：隨著球殼結構厚度的減小，位移響應曲線向著低頻方向發生移動，同時響應峰值增大；在低頻段(10～1 000 Hz)響應波峰較分散，而在高頻段(1 000～2 000 Hz)波峰較為密集。在兩種載荷作用下，不同阻尼系數的半球殼結構頂點的法向位移響應曲線如圖8所示。由圖8可知，當結構阻尼增大時，其共振峰位置幾乎不發生改變，而響應峰值明顯降低。這是由于隨著阻尼力的增大，結構內摩擦所損耗的能量更多，從而造成響應值的下降。

(a) 面載荷作用下

(a) 面載荷作用下

4 結 語

將Bezier函數作為開口球殼結構的位移試解函數建立了開口球殼振動分析模型，采用Ritz法研究了邊界條件、結構參數的變化對開口球殼結構振動特性的影響，并通過與已有計算結果進行對比，驗證了本文方法的正確性。研究主要得到的結論如下：

(1) 本文方法具有較好的收斂性和較高的計算精度。在收斂性方面，Bezier函數階數M≥10 時，計算結果收斂，收斂速度快，且具有良好的穩定性。在計算精度方面，本文結果與已有公開發表的文獻結果一致性較好。

(2) 球殼結構固有頻率不僅與邊界條件有關，而且與其自身結構屬性有關。開口球殼不同階次的頻率參數均隨著厚度的增加而增加，其中高階頻率受到厚度的顯著影響。開口角度對球殼振動響應有顯著影響，無量綱頻率隨著開口角度的增大而減小。

(3) 球殼結構受迫振動響應峰值位置受到結構參數和激勵形式的影響，其中低頻段的響應位置對球殼開口角度、厚度的變化更為敏感。