楔形體勻速和自由落體水彈性砰擊的半解析法求解

馮 松， 高 健， 陳毓珍， 孫 哲， 張桂勇,4

(1.大連理工大學 船舶工程學院 工業裝備結構分析國家重點實驗室,大連 116024；2.陸軍裝備部駐大連地區軍代室,大連 116001；3.中國船舶集團有限公司 上海船舶研究設計院，上海 201203；4.高新船舶與深海開發裝備協同創新中心,上海 200240)

入水現象廣泛存在于船舶與海洋工程(如惡劣海況中高速船船艏在與波浪之間循環砰擊)以及航天工程(如水上飛機在海平面上的降落)。入水過程中濕表面迅速增大，附連水質量急速增加，從而在固液交界面形成強烈的非線性沖擊載荷，容易造成局部結構屈服失效或者疲勞損傷，因此有必要在設計階段根據特定物理模型預報結構的響應從而評價結構抗砰擊強度。

彈性體入水過程中流體壓力作用于結構上，改變了固液交界面幾何形狀和速度；反之交界面形狀和速度變化會影響流場的邊界條件，進而改變流體的壓力。因此彈性體入水是典型的流固耦合問題，稱之為水彈性砰擊。由于楔形體在高速船中廣泛應用，因此經常作為分析對象。為準確模擬彈性楔形體和流體之間相互作用，以往學者采用了半解析法[1-2]、數值耦合[3-8]和試驗[9-12]3種方法來研究水彈性砰擊。半解析解中流場壓力的理論基礎大部分源于Wagner[13]提出的相當平板擬合理論。相比于數值方法，半解析法重點關注于造成結構響應的主要壓力成分，忽略了流體的壓縮性和黏性、重力和非線性自由液面等次要因素和射流等物理細節。因此半解析法具有計算效率高和求解穩定的優勢，并能提供滿足工程應用的結構響應值，同時也為驗證數值方法提供了重要參考數據，所以半解析法得到了諸多學者的關注并取得長足的發展。

為考慮彈性效應對速度勢分布的影響，不同學者采用了不同的處理方法。一種方法是在相當平板理論基礎上對砰擊速度直接修正，如平均濕表面彈性速度修正[14-16]和直接彈性速度修正[17]；另一種是基于陣型函數展開[18]，將濕表面速度勢表達為陣型函數關于正則坐標速度的疊加形式，這類方法需要求解一個邊值問題，對于簡單的邊界條件該邊值問題可以由解析方法獲取。借助上述基礎，流體的壓力與結構的彈性響應得以相關聯，從而為后續耦合方程的建立打下基礎。Korobkin等耦合Wagner理論和基于梁理論的有限元法(finite element method，FEM)，求解了恒速入水的二維剖面在砰擊階段(結構完全浸濕之前)的水彈性問題。Khabakhpasheva等耦合Wagner理論和模態疊加法，并將求解時間拓展至分離階段(結構全部浸濕之后)，其中分離階段的附加質量為分離瞬間附加質量，等效載荷設定為0。于鵬垚等[19]通過引入虛擬物面的方法，模擬了分離階段水彈性響應，壓力求解基于MLM(modified logvinovich model)從而考慮二階非線性壓力的影響。

目前半解析法求解水彈性砰擊主要局限于恒速入水，自由落體入水研究較少，而真實的物理問題中結構整體加速度、彈性變形和流場壓力之間相互影響。為考慮結構在自由入水狀態下的水彈性問題，Panciroli等[20]提出一種半解析顯示求解格式；但這種方法需要試驗測量加速度和流體壓力數據的導入，不具備理論獨立性。Shams等[21]考慮了結構的自由入水砰擊，結構固有陣型通過施密特正交化獲取，耦合方程通過隱式的Newmark-β法求解；但Shams等忽略了彈性加速度對整體加速度的影響并使用線性化壓力結果，在與Panciroli等的試驗(底升角30°)對比中發現計算應變偏大，Shams等將其原因歸結為非線性壓力的忽略和三維效應的影響。

本文根據平均濕表面彈性速度修正模型，建立了水彈性砰擊過程中剖面加速度、結構彈性響應和流體壓力之間耦合關系，基于MATLAB語言開發出一種適用于求解楔形體垂直恒速和自由入水的流固耦合問題的算法。為了提高大底升角楔形體的預報精度，本文考慮了伯努利方程中非線性項。另外，在求解浸濕長度過程中，本文基于Wagner條件考慮了剖面彈性變形的影響，并采用顯式求解，避免了傳統隱式求解浸濕長度的復雜性。

1 理論基礎

如圖1所示，二維對稱結構剖面垂直向下運動，剖面底部外板簡化為龍骨和舷端之間的兩個單跨梁。外板與龍骨相連，但在舷端可能為自由狀態，以后面具體算例定義為準。梁的厚度為h，初始底升角為β，龍骨與舷端距離為L。以t表示時間變量，在t=0時刻龍骨與靜止狀態的自由液面接觸，入水過程開始。流體假設為無黏、不可壓縮且運動無旋，即滿足勢流理論。由于一般工程問題中入水過程持續時間較短，因此重力對流體和結構的影響均忽略不計。定義整體坐標系xy用來描述流體運動，其中x軸指向右側為正，位于初始靜止液面，y軸向上為正，位于結構對稱線。定義局部坐標系ηζ用來描述結構的彈性運動，其中η軸與初始時刻無變形的結構重合并指向舷端，ζ軸垂直η軸向上并通過龍骨點，ηζ坐標系隨剖面垂直向下運動并保持原點與龍骨的重合。在t時刻，自由液面與結構相交于點Pwet。Pwet在xy坐標系橫坐標c(t)定義為橫向浸濕長度，在ηζ坐標系橫坐標s(t)定義為斜向浸濕長度。此時結構整體的運動速度為v(t)，向下為正；加速度為a(t)，向上為正；彈性變形為w(η,t)，向剖面內凹為正。

圖1 二維彈性剖面入水模型定義

1.1 流體壓力計算

基于勢流理論，引入速度勢φ，流體速度由?φ表示，其中?為梯度算子?？紤]到結構的彈性變形，結構濕表面的速度勢分布采用平均濕表面彈性速度修正模型，即

(1)

(2)

式中，ρf為流體密度，本文均取為1 000 kg/m3。將式(1)代入式(2)有

(3)

1.2 整體加速度計算

假設結構整體質量為M，結構下落過程中受到流體壓力，在y方向由牛頓定律有

(4)

(5)

通過式(5)整體加速度與彈性變形的關系得以建立。

1.3 水彈性方程的建立

本文中固體求解采用歐拉伯努利梁理論，作用于梁上外力包括水壓力以及由于整體結構非慣性運動產生的慣性力。梁的運動控制方程如下

(6)

式中，Es為結構的彈性模量。梁上端面的應變可由下式求解

(7)

引入模態疊加法，將結構變形表達為有限個陣型函數的線性疊加形式

(8)

將式(8)代入式(6)，同時將式(8)代入式(3)并將結果代入式(6)，兩端乘以Wi并在[0,L]區間積分，并利用形函數的正交性和質量歸一性，最終得到正則坐標控制方程如下

(9)

式中,ωi為i階固有頻率。針對i=1～Nmod重寫式(9)并按照矩陣形式給出正則坐標向量的控制方程如下

(10)

(11)

(12)

(13)

式(10)表明在入水過程中，流體對結構一部分作用力等效于附連水慣性力和阻尼力。由于浸濕長度的變化，因此附連水質量和阻尼矩陣均為時間的函數。f的右端第二項涉及結構彈性速度平方項，若放置于左端則出現非線性項，給求解帶來困難。本文將其處理為一外力項，由上一時間步結構響應獲取。在某一時間步內，Mf，Cf和f認為不變，式(10)變為一線性有阻尼強迫振動方程，可由隱式Newmark-β算法求解。

1.4 濕表面長度求解

流體壓力、剖面加速度和結構響應通過式(3)、式(5)和式(10)得以建立耦合關系，但在對3個方程求解之前，濕表面長度和其關于時間的導數仍然未知。首先根據圖 1中的幾何關系，結構節點在xy坐標系的位置可表達為ηζ坐標系中的參數函數，即

x(η)=ηcosβ-w(η,t)sinβ

(14)

(15)

橫向浸濕長度c(t)和斜向浸濕長度s(t)關系類似可以寫作

c=scosβ-w(s,t)sinβ

(16)

(17)

式中,y=fb(x,t)為考慮整體位移和彈性變形的剖面形狀方程。傳統對浸濕半寬的求解大多基于式(17)采用隱式求根方法，需迭代多次才能得到滿足精度要求的結果。本文將式(17)轉化為s的初值問題，避免了求根過程。對式(17)兩端關于時間求導有

(18)

(19)

(20)

式中，wη=?w/?η為剖面在局部坐標系斜率。將式(19)和式(20)代入式(18)有

(21)

(22)

1.5 耦合計算過程

根據式(3)、式(5)和式(10)以及浸濕長度的計算方法，可以建立起耦合計算的遞進過程，計算步驟簡略示意如圖 2，其詳細計算過程如下。

圖2 耦合計算示意圖

步驟1準備階段

步驟1.2求解結構固有陣型，固有頻率。

步驟2迭代計算階段

步驟2.1根據s(tn)按照式(16)計算c(tn)。

步驟2.2按照s(tn)將濕表面均等劃分為Nmesh個網格，并計算各個節點對應的η、x、y、θ=arcsin[x/c(tn)]、Wi和fp。

步驟2.5基于式(8)計算結構變形w及其導數。

步驟2.7按照式(11)～式(13)計算附連水質量矩陣Mf、阻尼矩陣Cf和外力向量f。

步驟2.9基于式(5)計算結構整體加速度a(tn+1)，對于恒速入水模型，強制另a(tn+1)=0。

步驟2.10根據加速度結果積分預估計算下一時刻速度：v(tn+1)=v(tn)-0.5[a(tn)+a(tn+1)]dt(負號是由于速度與加速度正方向規定不同)。

步驟2.12若浸濕長度超過L，退出循環，否則進入下一步，n=n+1。

上述計算過程中，式(11)～式(13)和式(21)涉及的積分求解，可利用步驟2.2離散化的結果由數值積分完成。

2 算例驗證

為驗證第1章算法并討論整體加速度、流體壓力和彈性響應耦合，本章給出具體的幾個算例，其模型參數匯總如表1所示。

表1 計算模型參數

表1中:BC一欄為邊界條件,第一個字母為龍骨的約束條件,第二個字母為舷端約束條件,字母C為固支,S為簡支;F為自由;W1和W2以恒速運動狀態入水；W3和W5自由落體入水；W4以恒速或者自由落體入水。材料Al的彈性模量為Es=67.5 GPa，密度為2 700 kg/m3；材料Fe的彈性模量為Es=200 GPa，密度為7 800 kg/m3。

2.1 收斂性分析

首先基于模型W1給出算法的收斂性分析，并為后續計算提供網格密度、模態數和時間步長的選擇依據。每組收斂性分析均選取結構中點變形的時程曲線作為對比目標。

圖 3(a)給出的是網格密度收斂性分析結果，其中模態數取為Nmod=10，時間步長取為dt=1×10-5s，網格數目Nmesh依次取為50，100和200。從圖 3(a)中可以看出，不同網格數目對應的計算結果非常接近，網格數100和200對應的最大變形相差為0.02%，可以認為網格數200達到收斂要求。

圖 3(b)給出的是模態收斂性分析結果，其中網格數目取為Nmesh=200，時間步長取為dt=1×10-5s，模態數Nmod依次取為1，2，5和10。從圖 3(b)可以看出，不同模態數對應的計算結果非常接近，前5階與前10階計算最大變形相差0.02%，可以認為前10階模態滿足收斂要求。

圖 3(c)給出的是時間步長收斂性分析結果，其中網格數目取為Nmesh=200，模態數取為Nmod=10，時間步長dt依次取為1×10-4s，5×10-5s和1×10-5s(對應的網格和時間步長無量綱比依次為L/(Nmeshvdt)=1.02，2.03和10.15)。從圖 3(c)可以看出，不同時間步長對應的計算結果非常接近，時間步長5×10-5s與1×10-5s前10階計算的最大變形相差為1.5%，可以認為dt=1×10-5s滿足收斂要求。

(a) 網格收斂性

基于上述收斂性結果，在后續流固耦合分析中，將選取以下計算參數：Nmesh=200，Nmod=10，dt= 0.1×L/(Nmeshv)。

2.2 恒速入水

圖4和圖5分別給出W1模型0.25L監測點處壓力和總體砰擊力時程曲線，其中壓力考慮了非線性項。計算結果與基于移動粒子半隱式法(moving particle semi-implicit method，MPS)的耦合求解器MPS-MPS[26]和基于光滑粒子流體動力學(smooth particle hydrodynamics，SPH)與FEM的分區耦合求解器SPH-FEM[27]進行了對比。為了比較剛體與彈性體的水壓力不同，剛體的計算結果也參與對照(剛體的解可由本文拓展得出，即不考慮彈性變形對壓力和浸濕長度變化的影響)。本文模擬的0.25L監測點砰擊壓力和砰擊力在整體變化趨勢與文獻吻合較好，壓力峰值計算結果為25 MPa，與MPS-MPS結果相對誤差為4%，最大砰擊力與SPH-FEM預估的誤差為3%。從圖 4和圖 5可以看出，剛體和彈性體監測點壓力在接觸流體一瞬間達到峰值而后迅速減小，之后剛體與彈性體的壓力演化有著明顯不同。由于剛體結構剖面形狀不隨時間改變，因此壓力最終收斂，整體砰擊力也關于時間呈線性變化。彈性體由于剖面的變形，砰擊力不再隨時間呈線性變化。按照砰擊力增長速率變化，彈性體入水可以分為3個階段：第一階段為0～0.5 ms，此過程結構彈性變形很小，彈性體與剛體砰擊壓力值相差較小，兩者砰擊力幾乎吻合；第二階段為0.5～1.2 ms，此過程中結構內凹變形，龍骨處局部底升角變大，舷端底升角變小，由于自由液面未到達舷端附近，固液交點所處位置局部底升角大于初始底升角且逐漸變大，因此該過程砰擊力相對于剛體變弱；第三階段為1.2 ms舷端浸濕，此過程中固液交界所處位置局部底升角逐漸減小至初始10°并繼續減小，導致砰擊力繼續增加，最后的砰擊載荷甚至超過了剛體。

圖4 W1模型0.25L監測點壓力時程曲線

圖6給出了W1模型中點彈性變形非線性耦合和線性耦合計算結果，并與MPS-MPS和SPH-SPH[28]的結果進行了對比。從圖6可以看出，非線性解與文獻結果預估的最大變形相對誤差小于4%。相比之下，由于線性解忽略速度平方項，導致壓力被高估，線性解得到的結構變形略有提高。相比于非線性解，線性解的中點最大變形值提高了13%。

由于線性解的理論基礎為相當平板擬合，因此在小底升角楔形體中線性解與非線性解相差較小。但在較大底升角情況下，非線性項影響變大，若忽略非線性項勢必會過高估算結構響應。為了進一步驗證本文算法對大底升角楔形體模擬，我們采用線性解和非線性解對模型W2(底升角為30°)的楔形體進行了計算，并給出中點變形時程曲線，見圖7?；谶吔缭?boundary element method，BEM)耦合FEM的數值方法BEM-FEM和基于開源代碼的計算流體力學(computational fluid dynamics，CFD)方法[29]的結果也一并繪出作為對比。從圖7可以看出，非線性解相比于文獻結果仍具有較高吻合度，而線性解則過高估算了結構的響應，兩者計算的中點最大變形相差達到了50%。

為比較線性解和非線性解對不同底升角楔形體應變預報的差別，表 2匯總了模型W2在兩種解計算下得到最大應變。其中最大應變為整個梁在入水過程所有時間中出現的最大應變的絕對值，εL,max為線性解的最大應變，εN,max為非線性解的最大應變。從表2可以看出，在不同底升角下線性解得到的最大應變均高于非線性解。10°底升角下兩個解差別最小，最大應變相差為12%。隨著底升角的提高，線性解與非線性解之間差別增大。當底升角為30°時，兩個解得到的最大應變相差達到48%。以上分析表明在大底升角情況下非線性項的影響不可忽略。

表2 線性解與非線性解計算W2模型最大應變對比

通過恒速入水算例可以看出，結構的彈性變形和流體互相影響，因此彈性楔形體入水需要采用耦合的方法來研究，特別是結構變形較大的情況(如較大的砰擊速度，較小的底升角和較弱的結構形式)。另外，線性解在大底升角情況下會嚴重高估結構響應，因此考慮非線性壓力對大底升角結構輕量化設計具有重要意義。

2.3 自由入水

相比于恒速入水問題，自由入水問題中整體垂向加速度對結構變形和流體壓力均有影響，反之流體壓力和結構彈性加速度也影響整體垂向加速度，因此彈性體自由入水中需要同時考慮整體加速度、流體壓力和彈性響應的變化。

為驗證本文算法對自由下落彈性楔形體響應預報的可靠性，本節利用非線性耦合給出模型W3在監測點η=0.1L的應變預報，并與Panciroli等的試驗結果和Shams等的線性半解析結果進行對比，如圖8所示。Panciroli等的試驗中給出監測點3次測量結果，本文中取3次測量平均值為基準值，每次測量結果的不同以誤差棒形式體現。從圖8中可以看出，在0～20 ms時間內，本文與Shams等的研究均能很好地與試驗結果相吻合，在20～45 ms時間內，本文與Shams等的研究均高估了該測點應變。試驗中入水初期浸濕長度較小，模型寬度(為0.25 m)遠大于浸濕長度，三維效應可以忽略。在入水后期，隨著浸濕長度提高，三維效應影響逐步提高。模型完全浸濕時濕表面長度與模型寬度比例為6/5。由于三維效應的存在，流體從模型前后方向逃逸，導致實際模型中流體壓力降低。因此，基于2D理論的本文和Shams等的研究在后半段時間內均高估了結構應變。相比于Shams等的研究，由于本文考慮了非線性影響，因此得到的結果更為接近試驗值，應變的預測得到明顯改善。

(a) v(0)=3 m/s

為揭示加速度的引入對耦合系統的影響，首先比較模型W4在不同剖面質量下入水速度和中點應變，見圖9。其中恒速入水作為非耦合解也一并給出，即不考慮加速度對流體壓力和結構變形的影響。圖9(a)比較了各個工況的入水速度。從圖9中可以看出，當剖面質量為20 kg，整體加速度較大，速度變化較大，最終入水速度僅為1 m/s。隨著剖面質量的增大，剖面整體加速度變小，且入水整體時間較短，因此入水速度趨近于初始入水速度3 m/s。圖9(b)給出各工況下中點應變。當剖面質量為20 kg，由于砰擊速度的急速下降，砰擊載荷降低，因此結構的應變相對于恒速入水明顯降低。隨著剖面質量增大，自由落體砰擊載荷逐漸趨向于恒速入水，因此中點應變趨近于恒速入水結果。當剖面質量為640 kg，自由落體中點應變幾乎與恒速入水結果相同。

(a) 入水速度

下面，為揭示結構彈性對整體加速度的影響，我們比較了模型W5在不同板厚下舷端變形和加速度，見圖10?；趧傮w模型的非耦合解也在圖10中給出，即不考慮彈性變形對流體壓力和加速度的影響。圖10(a)給出舷端的變形時程曲線。從圖10(a)可以看出，當板厚為14 mm，結構變形較小，中點最大變形僅為0.4 mm，結構可近似為剛體。當板厚為3 mm時，結構剛度較低，結構發生較為明顯的變形，端部最大變形達到了30 mm。圖10(b)給出整體加速度時程曲線。結構的剛度降低對整體加速度影響體現在兩個方面:一方面，由于結構的內凹變形，結構砰擊力相對于剛體降低，因此彈性體加速度小于剛體加速度；另一方面，整體砰擊力的一部分由彈性變形產生的慣性力吸收。且由于結構彈性變形，砰擊力和彈性慣性力出現周期性變化，因此加速度出現波動。當板厚為3 mm，結構剛度較低，這種現象最為明顯。隨著結構剛性的提高，結構彈性變形對整體加速度的影響變弱。當板厚達到14 mm時，彈性體下落加速度與剛體幾乎相同。

(a) 舷端變形

通過彈性體自由入水的分析可以看出，整體加速度對結構的應變發展和流體壓力演化具有重要作用，而結構的彈性對加速度也有很大的影響。在整體結構質量較輕、入水過程中砰擊速度產生明顯變化的情況下，若想精確捕捉整體結構非慣性運動下結構響應，必須考慮整體加速度在整個耦合體系中的作用。

3 結 論

本文提出基于耦合Wagner理論和模態疊加法的半解析方法，并計算楔形體以垂直恒速和自由落體兩種狀態入水，考慮了整體加速度、結構彈性響應和流體壓力三者之間耦合關系，并與文獻進行了對比，得到如下結論：

(1) 本文所提出的半解析法能夠精確預報恒速入水過程中結構響應和流場壓力。通過對不同底升角楔形體采用線性解和非線性解計算結果對比可以看出，線性解高估結構響應，兩組解的差異隨底升角增加而變大。對于本文所采用的算例，在30°底升角下線性解高估結構應變達到48%。

(2) 針對自由落體楔形體入水問題，本文提出的方法能夠總體上預報結構應變響應變化趨勢，且能提供準確的最大響應結果。在入水后期，由于未考慮到三維效應影響，本文計算結果略高于試驗結果。彈性體自由入水過程中彈性響應和整體加速度通過流體壓力相互影響。通過對彈性體自由落體和恒速入水對比可以看出，自由落體由于砰擊速度降低，結構應變幅值減??；通過對剛體和彈性體自由入水對比可以看出，彈性體由于結構剖面的內凹變形以及彈性慣性力的影響，其整體加速度小于剛體，且出現周期性波動。因此在自由落體入水過程中需要同時考慮整體加速度、流體壓力和結構彈性三者的影響關系并按照耦合的方法求解。