于鵬垚， 趙 勇， 王天霖， 甄春博， 蘇紹娟

(大連海事大學 船舶與海洋工程學院， 遼寧 大連 116026)

船舶砰擊載荷與結構響應在船舶與海洋工程領域具有重要的研究意義.高頻的砰擊載荷可能會引起結構的塑性變形和疲勞損傷，而且伴隨著輕質復合材料的應用，沖擊過程中的結構變形與流場間的耦合作用也愈發(fā)突出.V形剖面是高速船舶的典型剖面形式，因此，在砰擊問題研究中，楔形體入水沖擊模型被廣泛地采用.

Faltinsen[1]結合正交板理論和Wagner解析砰擊模型[2]分析了雙體船濕甲板的砰擊響應，其中，沖擊速度考慮了濕甲板彈性變形的影響.Khabakhpasheva等[3]基于Wagner理論和梁理論求解二維楔形體的入水沖擊響應，并給出3種近似求解方法，討論了不同耦合程度下各方法的求解精度.采用不同的結構響應求解方法，Shams等[4]和Datta等[5]也利用解析Wagner砰擊理論分析了楔形體的入水沖擊響應.陳震等[6]利用Dytran軟件對二維楔形體的入水砰擊問題進行仿真分析.Stenius等[7]利用LS-DYNA軟件分析了二維彈性楔形體入水沖擊的流固耦合機理，曹正林等[8]也采用類似方法進行三體船砰擊載荷的分析.Piro等[9]采用CFD軟件Openfoam求解流場，采用模態(tài)疊加法求解結構響應，進而分析了二維彈性楔形體入水和出水過程的振動響應.莫立新等[10]開展了不同剛度楔形體的入水砰擊試驗.王文華等[11]采用一種新型的CFD方法分析了各狀態(tài)參數對彈性楔形體入水性能的影響.

按照流場的求解方法不同，彈性楔形體入水沖擊的流固耦合分析方法可分為基于解析砰擊模型的分析方法和基于計算流體力學的分析方法，與基于計算流體力學的分析方法相比，基于解析砰擊模型的分析方法計算效率更高[12].現有的基于解析砰擊模型的分析方法普遍是基于Wagner模型進行流場求解，Wagner模型僅考慮了砰擊壓力的一階線性項，而且未考慮流動分離后的水動力.本文基于半解析砰擊模型[13-15](引入虛擬物面的 Modified Logvinovich Model[12]，MLM)建立彈性楔形體入水沖擊響應的流固耦合分析方法，其中，流場的求解是在MLM的基礎上實現的，結構響應采用模態(tài)疊加法來計算，并計入結構彈性振動對沖擊速度的影響來實現耦合作用的分析.與基于Wagner模型的水彈性分析方法相比，該方法不僅計及砰擊壓力的二階項，而且能夠實現流動分離后楔形體動力響應的分析.

1 半解析砰擊模型

如圖1所示，初始時刻t=0，流體靜止，OX軸與初始靜水面重合，楔形體為對稱剖面，其與靜水面的交點位于XOZ平面的原點O處.O′ξ軸沿著楔形體的斜邊方向，初始時刻，O′與原點O重合.n為楔形體斜邊的法向方向，遠離流體的方向為正.c(t)為水平方向的浸濕半寬，s(t)為沿著斜邊的浸濕長度.假定流體為不可壓縮的理想流體且忽略重力的影響，利用伯努利方程并考慮速度勢空間導數對應的二階項，可得流場內任意一點的壓力

(1)

式中：ρ為流體的密度；φ(x,y,t)為流場內的速度勢；φt為速度勢的時間導數.

圖1 入水砰擊的坐標系與參數Fig.1 Coordinates and parameters of water entry

在砰擊剖面為鈍形剖面的前提下，Korobkin[12]采用漸近分析方法詳細地推導了MLM砰擊模型，得到勻速入水時，楔形剖面表面壓力沿OX軸水平方向的分布為

(2)

式中：V為剖面勻速入水的速度；β為楔形剖面的斜升角.

c(t)可由下式的Wagner條件進行求解:

(3)

對于斜升角為β的楔形剖面，則有

(4)

(5)

當β已知時，τ即為一個定值.此時利用τ將 -a(t)

(6)

可見，當x=c(t)時，P(x,t)=0，避免了不合理的負壓.

當楔形剖面為有限寬度時，沖擊過程中水流會沿著楔形體的表面發(fā)生分離.對于流動分離后的砰擊載荷，可以采用引入虛擬物面的方法進行分析[13-15].如圖2所示，在流動分離邊界處，引入斜升角為α的虛擬物面，假設在流動分離發(fā)生后，液面沿虛擬物面上升，對應斜邊長度為L的楔形剖面，其剖面形狀表達式為

f(x)=

(7)

圖2 流動分離后的參數Fig.2 Parameters after flow seperation

(8)

當發(fā)生流動分離后，引入θs，假定當θ=θs時，有c(t)sinθs=Lcosβ，那么

(9)

(10)

為便于計算廣義力，利用c(t)=s(t)cosβ，x=ξcosβ，可得砰擊壓力沿O′ξ軸的分布

(11)

相應的楔形體斜邊浸濕速度為

(12)

2 彈性體周圍的壓力場與廣義力

在彈性結構的入水沖擊過程中，結構的自身振動會與流場間產生耦合作用，文獻[1,7]中的處理方法，是利用楔形體浸濕范圍內的平均振動速度來修正剛性楔形體法向n方向的沖擊速度，進而得到計及結構彈性效應的壓力分布

(13)

同樣，考慮彈性振動對沖擊速度的影響，可得楔形體斜邊浸濕速度

(14)

為便于后續(xù)推導，將彈性楔形體周圍的壓力場分解為剛體沖擊對應的壓力項Pr(ξ,t)和彈性效應對應的壓力項Pe(ξ,t)，即P(ξ,t)=Pr(ξ,t)+Pe(ξ,t)，則有

(15)

(16)

本文采用模態(tài)疊加法計算楔形體入水沖擊過程中的結構響應，需給出各階模態(tài)對應的廣義力.若第i階振動模態(tài)下，沿著斜邊法向方向的位移振型為ψi(ξ)，則在斜邊浸濕長度為s(t)時，對應的廣義砰擊力為

(17)

令fi(s(t))=fexc,i(s(t))+fela,i(s(t))，則有

(18)

(19)

式中：fexc,i(s(t))為楔形體剖面的廣義激勵力，其僅與剖面的振型和剛性剖面所受的砰擊壓力有關，而與剖面的彈性振動無關；fela,i(s(t))為彈性效應引起的廣義水動力.

若楔形體j階振動模態(tài)對應的主坐標為pj(t)，可得濕表面上的瞬時平均振動位移、瞬時平均振動速度和瞬時平均振動加速度的表達式

(20)

(21)

將式(16)和(20)代入式(19)，彈性效應引起的廣義力可以進一步表達為

(22)

式中：Aij(s(t))為楔形體的廣義流體附加質量；Bij(s(t))為振動速度線性項對應的廣義流體附加阻尼；Dij(s(t))為振動速度平方項對應的廣義流體附加阻尼.其表達式應為

ψi(ξ)dξ

(23)

(24)

(25)

那么，在第i階振動模態(tài)下，楔形體剖面對應的廣義砰擊力為

(26)

3 流固耦合方程的建立與求解

在不計及結構阻尼時，基于模態(tài)疊加法建立的楔形體動力學方程為

(27)

式中：a為結構廣義質量矩陣；c為結構廣義剛度矩陣；F為廣義力矩陣.

將式(26)代入式(27)，可得彈性楔形體入水沖擊的流固耦合方程

(28)

式中:

m為計入的彈性模態(tài)數目.

(29)

在具體求解中，本文采用適用范圍更廣的有限元方法來進行結構模型的模態(tài)分析，得到楔形體的離散振型ψi(ξ).對于式 (23)～(25)中的積分值，可借助有限元模型的網格節(jié)點通過梯形積分來計算.如圖3所示，O′ξ線上的短線代表有限元模型中離散節(jié)點的位置，由于楔形剖面上的網格是預先劃分，當計算中s(t)從零開始逐漸增大，會存在s(t)的端點恰好位于兩節(jié)點之間.若計算該時刻廣義力，需要重新插值得到s(t)端點位置對應的振型，并求解式(23)～(26)中的積分.而在每個時間步上都進行振型的插值和求積分需要消耗大量時間，也將導致計算更加復雜.

圖3 網格與浸濕長度的關系Fig.3 Relationship between mesh and wetted length

為解決上述問題，以第i階彈性模態(tài)為例，預先計算一系列s值下的下列各式的函數值

(30)

(31)

(32)

(33)

式中：sk取值恰好位于有限元模型中的節(jié)點處.當t時刻s(t)位于線元的中間時，采用線性插值的方法確定瞬時時刻的上述積分T1,i(s(t))、T2,i(s(t))和T3,i(s(t))，則t時刻Aij(s(t))、Bij(s(t))和fexc,i(s(t))可以分別表達如下

可以看出，利用上述方法更新水彈性方程式(28)中的相應系數，避免了每個時間步上振型的插值和積分的求解，加快了計算速度.

4 算法分析

4.1 計算模型

以文獻 [9] 中的彈性楔形體為例進行分析，其中，楔形體底升角為10°，斜邊長度L=0.5 m，板條梁的厚度h=18 mm，結構材料為鋼材，密度為 7 850 kg/m3，彈性模量E=210 GPa，流體的密度為 1 000 kg/m3，板條梁兩端取為簡支邊界.楔形體的沖擊速度為4 m/s.本文計算中，楔形體的模態(tài)分析利用商業(yè)有限元軟件Nastran來實現.采用殼單元進行板條梁的建模，建模中板條梁的寬度取值為 0.04 m，并且沿著寬度方向僅劃分一個網格.通過模態(tài)分析，得到各階模態(tài)的振型、固有頻率、廣義質量和廣義剛度信息，如圖4和表1所示.本文僅列出前5階振型.

圖4 模態(tài)振型Fig.4 Mode shapes

表1 模態(tài)信息Tab.1 Modal information

4.2 收斂性分析

對本文方法中涉及到的楔形體網格離散數目和計算模態(tài)數目進行收斂性分析.將楔形體的斜邊分別離散為100、200和400等份進行計算，用于計算的模態(tài)數目為前5階模態(tài)，時間步長為 4.0 μs，得到楔形體斜邊中心位置處沿著斜邊法線方向的振動位移如圖5所示.圖中：m為楔形體的網格離散數目；w為振動位移.可以看出，不同網格密度下的振動響應十分接近，200與400等份之間的最大位移響應偏差為 0.077%.將楔形體的斜邊離散為200份，計算的模態(tài)數目為前5階模態(tài)，時間步長Δt分別取值為 4.0、2.0 和 1.0 μs，得到楔形體斜邊中心位置處沿著斜邊法線方向的振動位移如圖6所示.可以看出，不同時間步長下的振動響應十分接近，4.0 μs與 1.0 μs對應的最大位移響應偏差小于 0.001%.將楔形體的斜邊離散為200份，時間步長取值為 4.0 μs，用于計算的模態(tài)數目分別為前1、3和5階模態(tài)，得到楔形體斜邊中心位置處沿著斜邊法線方向的振動位移如圖7所示.可以看出，不同計算模態(tài)下的振動響應十分接近，前3階與前5階計算模態(tài)的最大位移響應偏差為 0.083%.

圖5 網格密度的收斂性Fig.5 Convergence of mesh densities

圖6 時間步長的收斂性Fig.6 Convergence of time steps

圖7 模態(tài)數目的收斂性Fig.7 Convergence of mode numbers

總體看來，本文算法具有較好的網格收斂性和模態(tài)收斂性.后續(xù)將斜邊劃分為200等份，并取前5階振動模態(tài)用于計算.

4.3 算法驗證

圖8 解耦結果的比較Fig.8 Comparison of decoupling results

圖9所示為本文方法與文獻 [9] 的計算結果的對比.可以看出，本文方法不僅可以求解流動分離前的結構變形，而且也可獲得流動分離后的振動響應.

與基于Wanger砰擊模型的計算結果相比，本文方法與利用CFD方法求解流場的計算結果更為接近.總體看來，本文方法可以應用于彈性楔形體入水沖擊響應的分析.

圖9 本文方法與文獻結果的對比Fig.9 Comparison of this method with the literature results

5 流固耦合的影響分析

為分析流固耦合作用的影響，進行不同沖擊速度下彈性楔形體入水沖擊響應的解耦方法和耦合方法的對比計算，結果如圖10所示.圖中，楔形體參數為4.1中的計算模型，比較對象為楔形體中間位置沿著法向方向的位移.可以看出，在低速狀態(tài)下，解耦方法和耦合方法在流動分離前的結構響應比較接近(見圖10(a))，而流動分離后兩者存在明顯不同，主要由于流動分離后結構的自身振動響應其主要作用，解耦方法未考慮流體的附加質量，以干模態(tài)下的固有頻率振動，而耦合方法考慮了附加質量的影響，將以濕模態(tài)固有頻率進行振動.隨著沖擊速度的增加，在流動分離前的階段，2種方法計算結果也將出現明顯差別.定義wDC和wC分別為解耦和耦合方法的位移響應峰值，如圖10(d)和表2，在某些沖擊速度范圍內，解耦方法的最大位移響應小于耦合方法；而在某些沖擊速度范圍內，解耦方法的最大位移響應大于耦合方法，表明十分有必要采用合理的考慮彈性效應的耦合方法來分析結構的入水沖擊響應.

圖10 不同沖擊速度下解耦與耦合方法的比較Fig.10 Comparison between decoupling and coupling methods under different impact velocities

表2解耦方法和耦合方法位移響應峰值的比較

Tab.2Comparisonofdisplacementresponsepeakvaluebetweendecouplingandcouplingmethods

V/(m·s-1)wDC/mmwC/mm|wC-wDC|wC%20.363470.361590.5241.55631.51782.5452.41112.31754.0463.53943.41483.6574.98885.21694.3786.67097.495411.00

6 結論

本文基于對半解析砰擊模型、彈性結構的廣義力和耦合動力方程的研究，系統(tǒng)地建立了彈性楔形體入水沖擊的流固耦合分析方法.研究得到結論如下：

(1) 通過對比不同網格密度、時間步長和模態(tài)數目下的計算結果，驗證了該方法的收斂性.

(2) 通過與商業(yè)軟件結果、文獻結果的對比，表明該方法可用于沖擊過程中解耦和耦合響應的預報.而且與基于Wagner模型的水彈性分析方法相比，該方法不僅能夠合理地預報流動分離前的結構響應，而且能夠反映出流動分離后的振蕩響應.

(3) 采用耦合方法和解耦方法對比分析不同沖擊速度下楔形體的結構響應，表明在彈性結構入水沖擊分析時應合理地考慮流固耦合的影響，而本文方法可為該問題的分析提供參考.