王加夏，陳 月，彭丹丹，劉 昆

(江蘇科技大學 船舶與海洋工程學院，江蘇 鎮(zhèn)江 212003)

0 引 言

砰擊現象廣泛存在于在洶涌的海面上航行的船舶與海洋結構物中，例如小型的高速艇，大型的商船和豪華郵輪等。當砰擊發(fā)生時，結構物在短時間內受到較大的沖擊力，這種沖擊力有可能會破壞結構甚至造成人員的傷亡。因此，在船舶設計中，準確預報結構物砰擊載荷對結構設計至關重要。

砰擊問題的研究最早是由Von Karman[1]開展的，他將水上飛機的降落過程簡化為二維楔形體砰擊水面的理論模型。隨后，Wagner[2]在Von Karman 的基礎上，考慮了結構沖擊時的液面抬升現象，并且引入了濕表面的概念。Dobrovol'skaya[3]將二維剛性楔形體等速入水的勢流問題轉化為復平面內的自相似流問題，求解了整個流場的流體速度和壓力。Zhao 等[4]采用完全非線性邊界元法研究斜升角為4°～81°二維結構的入水問題，采用非線性物面條件、自由液面條件進行計算。Muzafrija 等[5]使用有限體積法（FVM）和VOF（Volume of Fluid）類型的自由表面建立模型，這種方法通過精細的網格能很好地捕捉射流的發(fā)展，結合自由表面捕獲模型，預測了結構入水的粘性自由表面流。在試驗研究方面，Stavovy[6]通過剛性平底結構及斜升角為1°～15°的剛性楔形體模型的入水沖擊試驗，得到了最大砰擊壓力與斜升角之間的關系。Ochi[7]總結上百個船模在波浪中的砰擊試驗數據，提出了船舶剖面底部的砰擊壓力公式。近年來，有限單元法（FEM）廣泛應用于流固耦合問題中。Peseux[8]應用有限元法和實驗方法研究了三維剛體和彈性結構的砰擊問題，結果表明剛體試驗的規(guī)律性較好，而彈性體結構的結果有待進一步解釋。趙九龍等[9]基于三維邊界元方法和Wagner 自由液面抬升理論，建立“等效自由液面”的新數值計算方法。

當不考慮結構變形時，流場演化特性以及結構砰擊載荷的相互作用已經通過數值模型進行了廣泛研究。然而在實際應用中，流場砰擊載荷和結構砰擊響應的過程是相互耦合、相互影響的，是由流場演化和結構響應特征共同決定的，因此需要建立計及結構變形效應的精確砰擊數值模型。本文利用STAR-CCM+（STAR）和Abaqus 軟件進行外部耦合，建立了考慮結構變形效應的砰擊數值模型，并研究三維楔形體結構入水砰擊時彈性效應對砰擊壓力和結構動態(tài)響應的影響，觀察不同入水階段，結構變形模式的變化，并考慮板厚變化對楔形體結構響應的影響。

1 數值模型

1.1 流體控制方程

STAR 軟件基于FVM 法將連續(xù)控制方程轉化為可以求解的代數方程組，使用VOF 多相流模型追蹤自由液面的形態(tài)變化。對于粘性的三維流動，假定流動由RANS 方程控制，其中湍流效應包括渦度粘性模型。在這種情況下，求解連續(xù)性方程、動量方程和2 個湍流特性方程。首先空間域被細分為有限數量的鄰接控制量（CVs），并且當CVs 隨著結構運動而移動時，必須注意空間守恒定律。控制方程如下[10]：

質量守恒

動量守恒

標量形式的一般輸運方程

空間守恒定律

其中：ρ 為流體速度；v 為流體速度矢量；vb為CV 表面的速度；n 為表面為S 體積為V 的CV 表面的單位向量法線；T 為指應力張量（用速度梯度和渦粘性表示）；p 為壓力；I 為單位張量；φ 指代標量（k，ε 或者ω）；Γ 為擴散系數；b 為單位質量的體積矢量；bφ為φ 的源和匯方向。

1.2 結構響應方程

在Abaqus 中計算水動力載荷下結構變形時，需求解瞬態(tài)動力學方程[11]：

式中：M 為質量矩陣；Mh為附加質量矩陣；K 為剛度矩陣；D 為阻尼矩陣，由M 和K 決定；Dh為附加阻尼矩陣；Fce為旋轉所產生的離心力；Fco為慣性力，當小變形時可忽略不計；Fh為水動力載荷；u 為位移。

1.3 流固耦合

結構入水砰擊問題屬于強非線性的流固耦合問題，本文通過協同耦合功能自動在流固交界面處進行數據的交互傳遞實現雙向的流固耦合計算。其中STAR 求解流體方程，Abaqus 求解結構方程，屬于分區(qū)式流固耦合FSI（Fluid-structure interaction）。本文主要采用隱式耦合算法，并在每個外迭代計算流場載荷及更新結構物的受力、位移及變形。

根據流固耦合所遵循的守恒原則，在流固耦合交界面上，流體域和固體域之間相互傳遞的應力σ 與位移d 等變量應當是守恒的，即

式中：下標f 表示描述流體域的變量；下標s 表示描述固體域的變量。流固耦合中流體求解和固體求解的方程及數據交換如圖1 所示。

在STAR 中設置運動規(guī)范為“變形”，在Abaqus

圖1 流固耦合方程及數據交換示意圖Fig. 1 Sketch map of fluid-structure coupling equation and data switching

圖2 流體流動與流動誘導的結構響應求解流程圖Fig. 2 Flow chart of a coupled simulation of fluid flow and flow induced motion in a floating body

設置模型的彈性模量、板厚、初始速度、重力、時間步長和耦合面等參量。具體計算流程如下：首先STAR 開始流場載荷計算，然后將壓力和剪切力傳遞給FEM 求解器，之后Abaqus 根據耦合界面的載荷進行結構響應分析并將得到位移和變形量傳遞給STAR，STAR 再根據傳遞的位移量通過“變形”功能實現網格的移動，這樣完成一次協同耦合進行計算更新，重復直至達到最大物理時間。

2 計算模型及模型設置

本文以30°傾角的楔形體為例，建立變形楔形體入水砰擊的數值模型，并和模型試驗相對比，驗證了數值模型的有效性。同時增加板厚使得結構偏于剛性，研究楔形體入水砰擊時變形效應的影響規(guī)律，研究砰擊載荷時空分布、結構響應規(guī)律等特性，分析得到在砰擊現象這一強流固耦合問題中不同板厚時結構變形的影響規(guī)律。

2.1 計算模型

由于結構在入水砰擊過程中，涉及到結構-液體-氣體等三相物質，同時還需進一步追蹤結構和自由液面大位移的運動，結構的變形使得每一時間步需要對網格進行重構，增加了計算的復雜度。為解決此問題，本文引入三維重疊網格技術（overset mesh），該技術可考慮結構大變形，適合求解該類復雜的流固耦合問題且無需網格重構。另外，固液交界面會出現流場射流現象，為清晰捕捉射流的演化形態(tài)，本文在交界面處采用網格細化技術得到精細網格。

整體計算背景區(qū)域為0.8 m×0.8 m×0.8 m，重疊網格區(qū)域為0.4 m×0.7 m×0.35 m。坐標系原點在楔形體底部中點處，z 軸豎直向上為正，x 軸水平向右為正，右手螺旋法則，其中水面在z=-0.01 m 處。計算域在Y=0 處剖面如圖3 所示，楔形體尺寸如圖4 所示。

圖3 Y=0 處計算域剖面Fig. 3 Computational domain profile in Y=0

圖4 楔形體模型尺寸Fig. 4 Size of the wedge model

計算區(qū)域內部采用表面重構、自動表面修復、切割體網格單元、棱柱層網格的方法進行網格的生成，為了在最大程度上消除在2 個網格間插入變量時產生的錯誤，對背景網格和重疊網格的重疊區(qū)域使用相同的網格密度數量級。背景區(qū)域除頂部采用壓力出口條件，其余都采用壁面邊界條件，重疊區(qū)域中楔形體采用壁面邊界條件，其余采用重疊網格邊界條件。計算區(qū)域的網格劃分具體如圖5 所示。

2.2 測點分布

通過在底板不同高度位置處布置測點對結構入水時壓力時間變化進行監(jiān)測，考慮到楔形體的對稱性，只在一側取點。測點布置如圖6 所示。底板中心為O 點，在底板一側的橫向中心線距離O 點15 mm 處布置測點P1，之后每隔30 mm 布置1 個測點，即P2，P3，P4，P5。

2.3 模型驗證

圖5 計算區(qū)域網格劃分Fig. 5 Computational domain mesh

圖6 楔形體底部測點布置圖Fig. 6 Point layout at the bottom of the wedge

為了驗證本文所建立的考慮結構變形效應的流固耦合數值模型的有效性，通過STAR 與Abaqus 進行耦合對彈性楔形體入水砰擊過程進行仿真模擬，并和模型試驗[12]進行對比，試驗中模型的板厚為5 mm，從0.4 m 的高度自由下落，模型尺寸見圖4。圖7 給出了P3，P4測點壓力曲線的實驗值和數值仿真結果對比，楔形體在接觸水面的瞬間，砰擊壓力迅速達到峰值，隨著時間的推移，各測點出現峰值后迅速減小至趨于穩(wěn)定，測點P3的峰值要大于測點P4。由圖可知，兩者曲線吻合良好，彈性仿真的最大峰值相比試驗測量值稍小，可以認為彈性體入水耦合仿真的結果是可靠的。

3 計算結果分析

板厚的變化會影響結構彈性，為了研究厚度參數對彈性結構響應的影響，本文針對不同板厚的楔形體進行了計算，通過對比結構的砰擊壓力、變形位移和應力分布等參量來分析結構彈性變化對結構動態(tài)響應的影響。

圖7 測點實驗值和數值仿真結果對比Fig. 7 Comparison of experimental and numerical simulation results

為分析板厚對楔形體入水動態(tài)響應的影響，改變斜升角30°楔形體的板厚，定義板厚分別為2 mm，5 mm，8 mm，12 mm，設定楔形體以10 m/s 的速度沖擊水面。不同板厚時P1～P4測點砰擊壓力曲線如圖8所示。由圖可知，楔形體在接觸水面的瞬間，砰擊壓力迅速達到峰值，隨著時間的推移，各測點相繼出現明顯的峰值后迅速減小，板厚越大砰擊壓力越大。另外，從圖中圓圈處可以清晰地看出，在0.008 s 之后，當板厚減小時，砰擊壓力曲線出現了明顯的波動，板厚為12 mm 時曲線較為平穩(wěn)，板厚5 mm和8 mm 時曲線有微小波動，板厚2 mm 時曲線波動較大甚至出現負壓，這是由于板厚的變化造成結構彈性特征改變引起的，板厚越小，結構彈性越大，流固耦合時彈性效應的影響越大，因此結構的變形越大，同時變形引起的壓力波動越大。

圖9 給出了彈性模量（E=2.1E11 Pa）和入水速度（v=10 m/s）保持不變時，板厚依次增加，楔形體的最大壓力峰值對比圖。由圖可知，隨著板厚的增加，楔形體的最大壓力極值是逐漸增大的。這是因為板厚越大使得結構剛性越大，質量越大，砰擊壓力極值越高。這也從側面表明了在結構入水砰擊響應計算中，為準確預報砰擊載荷，需要采用流固耦合的數值模型并計及結構變形效應的影響。

圖8 不同板厚各測點砰擊壓力曲線對比Fig. 8 Comparison of slamming pressure of each point with different thickness of plate

圖9 不同板厚時楔形體最大砰擊壓力峰值對比Fig. 9 Comparison of maximum slamming pressure peak values of wedges under different thickness of plates

圖10 不同板厚時各測點變形位移Fig. 10 Deformation displacement of each measuring point at different plate thickness

圖10 為不同板厚的楔形體測點P1-P5的變形位移曲線?？梢钥闯?，不同板厚的各測點位移曲線總體趨勢均為隨時間的增大而增大，但當板厚較小時（圖10（a）），各測點的變形幅度和變形模式有所差別，測點P1與其他測點的位移響應具有反相位特性；而在板厚較大時（圖10（b）），各測點的位移值非常相近，位移曲線幾乎重合。這是由于板厚改變了結構的剛度，板厚越厚，結構剛性越大，結構各位置變形幅度越一致。另外，在0.008 s 之前，即楔形體完全入水前，位移曲線斜率要大于完全入水后，這是由于完全入水前結構砰擊壓力很大，結構變形速度相對較大，而完全入水之后結構受力迅速降低，因此變形位移增加趨勢減緩。

圖11 不同板厚P1，P4 測點變形位移對比Fig. 11 Comparison of deformation displacement of P1 points and P4 points with different thickness of plate

圖11 為楔形體在不同板厚時P1，P4測點的變形位移對比。由圖可知，不同板厚時各測點的位移曲線均為隨時間的增大而增大；板厚12 mm 時，位移近似一條直線，而板厚2 mm 時，位移近似一條曲線，顯示出較為復雜的變形位移。這是因為板厚越厚，結構的剛度越大，結構越不易發(fā)生變形。

圖12 給出了楔形體板厚5 mm 時不同時刻的各方向應力云圖，從圖12（a）中可以看出，當楔形體剛接觸水面時，最大應力集中在楔形體底部尖端處；從圖12（b）、圖12（c）可以看出，在入水過程中，應力集中開始由底部尖端處向上轉移，最終出現在楔形體頂板中心；從圖12（d）可以看出，在楔形體完全入水后，楔形體的整體應力迅速降低，這與圖8 中0.008 s之后砰擊壓力迅速降低相符；從圖12（e）中可以看出，當楔形體完全入水一段時間后壓力穩(wěn)定時最大應力出現在楔形體底部尖端處，且頂板中心由于受擠壓也存在大的應力集中現象。

4 結 語

本文通過基于VOF 和FEM 方法的流固耦合數值模型，對三維彈性楔形體進行入水仿真，對比分析了結構彈性變形對砰擊壓力、結構動態(tài)響應和結構應力分布的影響。

1）通過STAR 與Abaqus 耦合對彈性楔形體的入水模擬得到的砰擊壓力值與實驗值吻合較好，驗證了本文建立的數值分析模型的有效性。

2）研究了結構在不同板厚情況下入水砰擊時結構變形效應的影響規(guī)律。結果表明：結構板厚越小，結構彈性效應對流固耦合的影響越大，結構所受砰擊壓力越小且壓力波動越大，同時變形位移增大且變形模式較為復雜。表明了結構設計過程中，計及結構變形效應的必要性。

3）通過分析結構響應觀察到結構在不同入水階段的應力變化規(guī)律為：楔形體完全入水前最大應力由底部尖端處向頂板中心處轉移，完全入水后結構整體應力迅速下降，最大應力出現在楔形體底部尖端處，頂板中心由于受擠壓也存在大的應力集中現象。

圖12 板厚5 mm 的楔形體不同時刻的應力云圖Fig. 12 Stress contours for different time of wedge in thickness 5 mm