有限元／間接邊界元法求解浸水板振動特性

劉 城，洪 明，2，劉曉冰

（1.大連理工大學船舶工程學院，遼寧大連116024；2.大連理工大學工業裝備結構分析國家重點實驗室，遼寧大連116024）

有限元／間接邊界元法求解浸水板振動特性

劉 城1，洪 明1，2，劉曉冰1

（1.大連理工大學船舶工程學院，遼寧大連116024；2.大連理工大學工業裝備結構分析國家重點實驗室，遼寧大連116024）

針對浸水板結構振動，考慮流體可壓縮特性，利用間接邊界元法計算附加質量矩陣，包含計算中奇異邊界單元的處理。結合結構有限元法求解無限域浸水板的振動特性，通過編制的相應計算機程序，對浸水懸臂矩形平板進行了流體間接邊界元與結構有限元耦合數值模擬；并對水下自由板進行了模態識別實驗，實驗結果驗證了數值計算結果的可靠性。典型算例表明了流體介質改變了結構的動力特性，在低頻段流體的可壓縮性對結構的動力特性影響甚微；推導的流固耦合數值方法對研究浸水結構動力和聲輻射研究有很好的參考價值。

有限元法；間接邊界元法；浸水板；附加質量矩陣；矩形平板；固有頻率；流體可壓縮性；模態識別

當振動結構周圍流體對結構振動的影響不可忽略時，稱為強流固耦合問題，諸如關系到艦船結構的減振降噪性能。流固耦合工程上廣泛遇到，但由于問題的復雜性，比如輕流體與重流體、可壓與不可壓流體以及系統振動頻率的不同等，都會對耦合的物理分析模型有很大的差異。計算方法上也是包含解析、數值和實驗，本文僅針對低頻段強耦合流固耦合系統振動特性的邊界元／有限元數值分析研究。

在低頻階段，流體介質對結構的作用可通過在結構運動方程中加入附加質量項來實現流體固體耦合［1-9］。附加質量的數值計算方法有很多種，如文獻［5］利用有限元法對結構的附連水質量進行了研究，利用有限元法計算無界輻射聲場時，必須人為的截斷流體域，且計算量大；文獻［6］和文獻［7］利用直接邊界元法（direct boundary element method，DBEM）求解浸水結構的附連水質量，但此方法只適用于水下封閉結構，對水下開口結構無法進行計算；對于無限大剛性障板中的板或者加筋板結構，可以通過Rayleigh積分來計算流體對結構的影響，但實際工程中水下結構多處于無限域或者具有自由液面的半無限域水域中，而非Rayleigh積分要求的具有剛性壁面的半無限域，所以該方法具有很大局限性；文獻［8］通過經驗公式并求解水下懸臂板的附加質量，從而求解結構的自振頻率和振型，該方法只適合水下懸臂板結構，且對于不同的長寬比及厚寬比，懸臂板的附加質量計算公式還有所不同，因此該方法也具有很大的局限性。

間接邊界元法（IBEM）是從直接邊界元法（DBEM）推導而來，它以結構表面的壓力差（雙層勢）和壓力梯度差（單層勢）為未知變量，能同時計算內場問題和外場問題，同時適用于無限域中的封閉和開口結構。間接邊界元法的另一個優點是它是通過變分原理推導出系統方程，得到的系統矩陣是對稱的，這種對稱性使其能夠更有效的與結構有限元法相結合。本文根據文獻［9］利用FEM／DBEM計算水下封閉結構的固有頻率的原理上，運用間接邊界元法（IBEM）［10-11］推導在低頻段任意浸水結構的附加質量陣的數值算法，并考慮流體的可壓縮性對水下結構動力特性的影響［12］，通過編制 FORTRAN計算機程序，得到結構的附加質量陣，然后與基于有限元結構質量矩陣相疊加，實現了浸水結構系統的振動特性數值模擬，并進行了相關模態識別實驗，數值計算結果與實驗結果的對比驗證了本文算法的準確性。

1 理論基礎

1.1 間接邊界元方程推導

設如圖1所示結構處無限流體域中，S為結構浸水表面，S兩側均有流體域。若結構浸水表面存在振動時，結構的振動會引起周圍流體介質的擾動，從而在流體域中將產生輻射壓力場，而流體的擾動也會反過來影響結構的振動。對于無粘且可壓縮流體，在線性小擾動的情況下，流體域中各點的壓力p滿足Helmholtz方程：

式中：p為計算場點聲壓，k＝ω／c定義為波數，其中ω為流體介質運動圓頻率，c為聲波在流體介質中的傳播速度。在浸水振動結構表面，流體壓力法向梯度與法向振動速度滿足：

式中：n為結構浸水表面法向（如圖1所示），vn為結構表面法向振動速度，ρ為流體介質密度。

圖1 IBEM聲場域及邊界示意圖Fig.1 Representation of acoustic domain and boundary of IBEM

利用格林公式，壓力輻射域中Helmholtz微分方程（1）可轉為振動結構邊界上Helmholtz積分方程：

式中：S、Ω1和Ω2分別表示結構表面、外域和內域，Y為振動結構表面（源點），X為流體域中計算點（場點），對于無限域，G X，Y( )＝e-ikr／4πr為Y點處的基本解，r為源點到場點距離，C X()為影響系數，它與Y及結構表面光滑度相關。

方程（3）為直接邊界元對應的外場問題邊界積分方程，內場問題邊界積分方程也有類似形式［9］。間接邊界元是從直接邊界元方程推導而來，它以結構表面的聲壓差（雙層勢）和聲壓梯度差（單層勢）為未知變量，在邊界表面兩側分別對內場問題和外場問題應用上述直接邊界元內外場的Helmholtz邊界積分方程，然后將兩方程相加，即可得聲場域內任意點的聲壓：

式中：σ為結構表面的法向壓力梯度差（單層勢），μ為結構表面的聲壓差（雙層勢），其表達式為

假設結構表面滿足Neumann邊界條件，由上述關系式，可得邊界條件與未知量的關系為

定義如下泛函［10］：

式中：泛函F（μ）具有超強奇異性，可通過下式改善其奇異性：

在邊界上對方程（10）進行邊界元離散，利用變分原理，即可得到系統方程，求解此系統方程可求出結構表面的聲壓差μ，然后利用式（4）即可求出流場中任意點的壓力。

1.2 間接邊界元方程的離散

本文以三結點三角形線性邊界單元在流體固體耦合邊界上對方程（9）進行數值離散，最終可得到如下形式的表達式：

式中：Q（ω）為間接邊界元的對稱影響矩陣；μ為結構表面結點上的未知變量壓力差列向量；A為結構表面流體單元面積矩陣。矩陣Q和矩陣A的單元塊表達式如下：

利用變分原理，根據方程（15）可以得到間接邊界元的表達式：

由此可以看出，影響矩陣Q（ω）為振動頻率ω的函數，根據式（13）求出結構表面聲壓差后，任意點的壓力由式（4）即可求得。

而對無粘、不可壓流體，流體中聲度c→∞，則k＝0，聲場中任意點的壓力p滿足Laplace方程［12］：

與可壓縮流體中的推導類似，流體的控制方程最終可寫成

從上述推導中可以看出，可壓縮流體與不可壓縮流體中，最終得到的流體方程均可表達成式（16）的形式。需要注意的是，可壓縮流體中的影響矩陣Q（ω）與頻率有關，而不可壓縮流體中的影響矩陣Q′與頻率無關。

1.3 結構有限元與流體間接邊界元耦合描述

考慮流體與結構振動的耦合影響時，結構振動方程的有限元形式為

式中：M為結構的質量矩陣；Cd為結構的阻尼矩陣；K為結構的剛度矩陣；u為結構有限元結點的位移向量，Fs為結構上作用外激勵向量，Fa為結構與流體耦合面上流體作用在結構上流體動壓力向量。

結構表面兩側均有流體介質時（如浸水非封閉結構），流體動壓力可表示為

式中：T為方向余弦轉換陣，A為結構表面流體單元面積矩陣，μ為結構表面結點壓力差向量。

假設結構為簡諧振動，則滿足關系a＝iωv，v和a為結構振動的速度和加速度。流體變量和結構變量之間的幾何關系如下：

方程（18）中的流體動壓力項可以轉化為

低頻時忽略流體阻尼效應，流體動壓力最終轉化為結構的附加質量效應項：

根據前面的推導可知，可壓縮流體中附加質量陣Ma與頻率有關，而不可壓縮流體中附加質量陣Ma與頻率無關，在計算過程中需注意兩者的區別。

計算出流體的附加質量陣后，水下結構流固耦合作用下的動力學方程的有限元形式變為

對應的廣義特征值問題為

即

式中：Mc＝M+Ma為耦合的質量矩陣。

求解上述廣義特征問題，即可求得浸水結構的固有頻率和振型。

基于上述有限元和間接邊界元理論，本文基于FORTRAN語言編寫了相關的計算程序，其具體的程序流程如圖2和圖3所示。

圖2 FEM／IBEM流固耦合動力分析主程序流程圖Fig.2 Flow chart of the fluid-structure interaction dynamic analysis program

圖3 IBEM計算附加質量陣子程序流程圖Fig.3 Flow chart of the added mass matrix program

從流程圖中可以看出，首先利用有限元知識構造結構的質量陣和剛度陣；然后利用間接邊界元理論，根據式（21）得到結構的附加質量矩陣；最后將結構的附加質量陣和結構本身的質量陣相疊加，得到耦合的質量陣，然后求解廣義特征值問題，計算出浸水結構的模態。

由間接邊界元推導出的附加質量矩陣如式（21）所示，從式（21）的求解形式可以看出，附加質量矩陣是一個實數對稱矩陣，且結構本身的質量陣也是一個實數對稱矩陣，故耦合的質量矩陣為實數對稱矩陣。為了與間接邊界元耦合方便，結構有限元也是采用三角形三結點殼單元，在求得結構的剛度陣和耦合的質量陣后，即可求得結構的自振頻率和振型，此類問題都可歸結為求解形如式（24）所示的高階廣義特征值問題，本文應用子空間迭代法來求解此類廣義特征值問題。

2 數值算例分析

對如圖4所示的浸沒于無限水域中的懸臂矩形彈性平板結構，平板長為0.406 4 m，寬為0.203 2 m，厚為2.67 mm，材料泊松比為0.3，楊氏模量為1.95× 1011N／m2，密度為7 700 kg／m3。分別將流體視為可壓縮流體和不可壓縮流體，比較流體的可壓縮性對結構振動特性的影響。

圖4 矩形懸臂板單元剖分示意圖Fig.4 Representation of mesh discretization for the cantilever plate

2.1 結構有限元與流體間接邊界元耦合描述

假設流體不可壓縮，流體的密度為1 000 kg／m3，浸水懸臂板前四階固有頻率如表1所示。

表1 空氣與不可壓流體中懸臂板固有頻率計算結果Table 1 Natural frequency of the cantilever plate in theair and incompressible fluid Hz

從表1的計算結果可以看出，本文的計算結果與文獻［13］的計算結果是相吻合的。

2.2 可壓縮流體中懸臂板振動計算

從前面的理論推導中可知，在可壓流體中結構的附加質量矩陣是振動頻率的函數，在不可壓流體中與計算頻率無關。假設流體的密度1 000 kg／m3，聲速1500 m／s，分別計算振動頻率為10、100、500和1 000 Hz時結構的前四階固有頻率，計算結果如表2所示。

表2 可壓縮流體中附加質量計算頻率對固有頻率的影響Table 2 Impact of the added mass matrix frequency on natural frequency of the plate in the incompressible fluid Hz

從計算結果可以看出，附加質量計算振動頻率在低頻范圍內變化對結構的固有頻率影響非常小，僅當增大到高頻范圍時才對結構的固有頻率有較明顯的影響。因此計算可壓縮流體中板的固有頻率時，低頻段附加質量計算頻率的大小對結果影響不大。

2.3 流體可壓縮性對結構動力特性的影響

將流體分別視作可壓縮與不可壓縮2種模型下計算出的板的各階固有頻率進行對比，考慮流體的可壓縮性對結構動力特性的影響，表3分別給出了結構在空氣中、不可壓流體和可壓流體中的前四階固有頻率對比，其中可壓流體中附加質量陣計算頻率為10 Hz，計算結果如表3所示。

表3 流體可壓縮性對固有頻率的影響Table 3 Impact of fluid compressibility on the structure natural frequency Hz

從表3結果可以看出，結構在不可壓縮與可壓縮流體中的計算結果相差非常小，流體的可壓縮性對結構的固有頻率影響極其微小。因此可以說，在分析結構流固耦合自由振動頻率時，可以忽略流體的可壓縮性對結構固有頻率的影響。

3 實驗研究

為驗證本文所推導的流固耦合數值方法與所編寫程序的正確性，本文對水下全自由矩形鋼板進行了模態識別實驗，實驗證明了數值計算結果的可靠性。

3.1 實驗模型與儀器簡介

本實驗采用的是一塊長為0.6 m、寬為0.25 m、厚為2.0 mm的矩形鋼板。為模擬出無限大水域的自由場環境，該實驗在消聲水池中進行，消聲水池采用長470 mm、底徑直徑50 mm的橡膠尖劈密排，可用空間為7.0 m×4.5 m×4.0 m。實驗測量的主要儀器有東華DH5922動態信號測試分析儀、YE1311掃頻信號發生器、JZK系列激振器和ICP加速度傳感器。在實驗過程中，用4根彈簧吊起矩形鋼板放入消聲水池中，這樣可以模擬鋼板的全自由邊界條件，實驗示意圖如圖5所示。

圖5 水下自由板模態識別實驗示意圖Fig.5 Representation of modal identification experiment for the free submerged plate

3.2 實驗模型的數值分析

運用本文編寫的流固耦合動力分析程序來計算矩形鋼板的自由模態，計算中取材料泊松比為0.3，楊氏模量為2.1×1011N／m2，密度為7 850 kg／m3；消聲水池中水介質視為可壓縮流體，密度為1 000kg／m3，水中聲音傳播速度為1 500 m／s；附加質量陣計算頻率取為10 Hz，鋼板的單元剖分示意圖如圖6所示。

圖6 矩形平板單元剖分示意圖Fig.6 Representation of mesh discretization for the rectangular plate

3.3 模型的實驗測量

本實驗在消聲水池中進行。實驗測量中，用4根彈簧吊起鋼板放入消聲水池中，采用激振器作為激勵源，鋼板上布置了15個加速度傳感器，傳感器的布置示意圖如圖7所示。

圖7 水下自由板模態識別實驗模型實物圖Fig.7 Model of modal identification experiment

3.4 結果對比

對實驗模型進行了耦合模態數值分析和模態識別實驗后，兩者得到的結果如表4所示。

表4 空氣中與水中自由板前三階固有頻率實驗值與數值解對比Table 4 Natural frequency of the free plate in the air and in the water

從表4結果可以看出，運用本文所推導的方法計算得到的結構振動模態結果與實驗結果是相吻合的。說明本文推導的計算方法與所編寫程序的正確性。圖8表示的是空氣中和水中全自由鋼板的前3階振型，從圖中可以看出，由于流體附加質量的影響，減小了結構的固有頻率，但對結構的振型影響很小。

圖8 空氣中和水中全自由矩形鋼板的前三階振型Fig.8 The first three mode shapes in the air and water

4 結論

本文在間接邊界元理論下，推導出低頻段水下結構的附加質量矩陣，編寫相應的FEM／IBEM流固耦合動力分析程序，考慮流體的可壓縮性對結構動力特性的影響，分析了水下懸臂矩形平板的動力特性；并針對水下全自由矩形鋼板進行了模態識別實驗，得到如下結論：

1）本文的間接邊界元求解附加質量的方法簡單有效且具有可靠的精度，并為求解水下結構物的振動模態提供了一定的依據。本文方法的最大優點是能同時計算內場和外場問題，適合無限域中任意形狀的復雜結構；

2）在低頻段，流體對結構的影響可以忽略流體阻尼的影響，僅考慮流體的附加質量效應，及將流體作用簡化為結構的附加質量項；而隨著頻率的增大，流體的附加質量效應減小，流體的阻尼作用增大，此時不能忽略流體阻尼作用的影響；

3）在低頻段，流體的可壓縮性對結構振動模態的影響甚微，在計算中可以忽略流體可壓縮性對結構振動模態的影響。

［1］錢勤，黃玉盈，劉忠族.求附連水質量的一種直接方法［J］.力學與實踐，1996，18（5）：19-21.

QIAN Qin，HUANG Yuyin，LIU Zhongzu.A direct solution method for the add water matrix［J］.Mechanics and Practice，1996，18（5）：19-21.

［2］王基盛，楊慶山.流體環境中結構附加質量的計算［J］.北京交通大學學報，2003，27（1）：40-43.

WANG Jisheng，YANG Qingshan.Calculation of add mass matrix of submerged structure［J］.Journal of Beijin Jiaotong University，2003，27（1）：40-43.

［3］王基盛，楊慶山，武岳，等.薄膜結構附加質量的計算方法［J］.空間結構，2003，9（3）.38-41.

WANG Jisheng，YANG Qingshan，WU Yue，et al.Calculation of add mass matrix of thin plate structure［J］.Space Structure，2003，9（3）：38-41.

［4］蘇海東，黃玉盈.求半無限域流場中附連水質量的一種簡便解法［J］.華中科技大學學報：城市科學版，2003，20（4）：14-16.

SU Haidong HUANG Yuying.Novel numerical method for solving added mass of structure embedded in semi-infinite liquid field［J］.Journal of Huazhong University of Science and Technology：Urban Science Edition，2003，20（4）：14-16.

［5］金占禮，王宗利，李紅云，等.結構在無限流體域中振動時附連水質量的數值計算方法［J］.上海交通大學學報，2000，34（8）：1079-1082.

JIN Zhanli，WANG Zongli，LI Hongyun，et al.Numerical calculation method of the add water mass of vibratory structure in infinite fluid domain［J］.Journal of Shanghai Jiao Tong University，2000，34（8）：1079-1082.

［6］安小同，洪明，李艮田.基于FEM／BEM計算浸水結構振動特性［J］.船海工程，2011，40（6）：55-58.

AN Xiaotong，HONG Ming，LI Yintian.Calculation of vibration characteristics of submerged plates based on FEM／BEM［J］.Naval Architecture and Ocean Engineering，2011，40（6）：55-58.

［7］李華東，朱錫，羅忠，等.附連水質量的邊界元法求解［J］.海軍工程大學學報，2009，21（2）：45-49.

LI Huadong，ZHU Xi，LUO Zhong，et al.Calculation of the added water mass based on boundary element method［J］.Journal of Naval University of engineering，2009，21（2）：45-49.

［8］LIANG Chochung，LIAO Chingchao，TAI Yuhshiou，et al.The free vibration analysis of submerged cantilever plates［J］.Ocean Engineering，2001，28（9）：1225-1245.

［9］EVERSTINE G C.Prediction of low frequency vibrational frequencies of submerged structures［J］.Journal of Vibration and Acoustics，1991，113（2）：187-191.

［10］ALIA A，SOULI M，ERCHIQUI F.Variational boundary element acoustic modelling over mixed quadrilateral-triangular element meshes［J］.Communications in Numerical Methods in Engineerging，2006，22（7）：767-780.

［11］WANG Weiping，ATALLA N.A numerical algorithm for double surface integrals over quadrilaterals with a 1／R singularity［J］.Communications in Numerical Methods in Engineerging，1997，13（11）：885-890.

［12］張升明.流體的可壓縮性對彈性結構振動的影響［J］.水動力學研究與進展，1994，9（4）：429-436.

ZHANG Shengming.The impacts on vibration characteristics of submerged structures of fluid compressibility［J］.Journal of Hydrodynamics，1994，9（4）：429-436.

［13］JEANS R A，MATHEWS I C.Solution of fluid-structure interaction problems using a coupled finite and variational boundary element technique［J］.J Acoust Soc Am，1990，88（5）：2459-2466.

The solution for vibration characteristics of submerged plates by applying FEM／IBEM

LIU Cheng1，HONG Ming1，2，LIU Xiaobing1

（1.School of Naval Architecture Engineering，Dalian University of Technology，Dalian 116024，China；2.State Key Laboratory of Structure Analysis for Industrial Equipment，Dalian University of Technology，Dalian 116024，China）

With a focus on the structural vibration of a submerged plate and consideration for the compressibility of fluid，the added mass matrix was calculated by utilizing the indirect boundary element method（IBEM）in this paper，including the treatment for the singularity boundary element in the calculation.In combination with the finite element method（FEM）for structure，the vibration characteristics of a submerged plate in the infinite fluid field were calculated.Through the use of corresponding programming，a numerical analysis of the coupling between the indirect boundary element of the fluid and the finite element of the structure was conducted for the submerged cantilever rectangular plate；in addition，a mode identification experiment was conducted on the underwater free plate.The experiment results verify the reliability of the numerical calculation results.The typical numerical example shows that the fluid medium changed the dynamic characteristics of the structure.At a low frequency band，the compressibility of the fluid had little impact on the dynamic characteristics of the fluid.The deduced numerical method for the fluid-structure coupling has excellent reference value for the dynamic research on a submerged structure and the research on acoustic radiation.

finite element method；indirect boundary element method；submerged plate；added mass matrix；rectangular plate；natural frequency；fluid compressibility；modal identification

10.3969／j.issn.1006-7043.201303052

O326

A

1006-7043（2014）04-0395-06

http：／／www.cnki.net／kcms／doi／10.3969／j.issn.1006-7043.201303052.html

2013-03-21. 網絡出版時間：2014-03-15 20：44：29.

國家自然科學基金資助項目（51079027）.

劉城（1987-），男，碩士研究生；洪明（1959-），男，教授，博士生導師.

洪明，E-mail：mhong＠dlut.edu.cn