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

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

0 引 言

砰擊現(xiàn)象廣泛存在于在洶涌的海面上航行的船舶與海洋結(jié)構(gòu)物中，例如小型的高速艇，大型的商船和豪華郵輪等。當(dāng)砰擊發(fā)生時(shí)，結(jié)構(gòu)物在短時(shí)間內(nèi)受到較大的沖擊力，這種沖擊力有可能會(huì)破壞結(jié)構(gòu)甚至造成人員的傷亡。因此，在船舶設(shè)計(jì)中，準(zhǔn)確預(yù)報(bào)結(jié)構(gòu)物砰擊載荷對(duì)結(jié)構(gòu)設(shè)計(jì)至關(guān)重要。

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

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

1 數(shù)值模型

1.1 流體控制方程

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

質(zhì)量守恒

動(dòng)量守恒

標(biāo)量形式的一般輸運(yùn)方程

空間守恒定律

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

1.2 結(jié)構(gòu)響應(yīng)方程

在Abaqus 中計(jì)算水動(dòng)力載荷下結(jié)構(gòu)變形時(shí)，需求解瞬態(tài)動(dòng)力學(xué)方程[11]：

式中：M 為質(zhì)量矩陣；Mh為附加質(zhì)量矩陣；K 為剛度矩陣；D 為阻尼矩陣，由M 和K 決定；Dh為附加阻尼矩陣；Fce為旋轉(zhuǎn)所產(chǎn)生的離心力；Fco為慣性力，當(dāng)小變形時(shí)可忽略不計(jì)；Fh為水動(dòng)力載荷；u 為位移。

1.3 流固耦合

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

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

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

在STAR 中設(shè)置運(yùn)動(dòng)規(guī)范為“變形”，在Abaqus

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

圖2 流體流動(dòng)與流動(dòng)誘導(dǎo)的結(jié)構(gòu)響應(yīng)求解流程圖Fig. 2 Flow chart of a coupled simulation of fluid flow and flow induced motion in a floating body

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

2 計(jì)算模型及模型設(shè)置

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

2.1 計(jì)算模型

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

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

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

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

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

2.2 測(cè)點(diǎn)分布

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

2.3 模型驗(yàn)證

圖5 計(jì)算區(qū)域網(wǎng)格劃分Fig. 5 Computational domain mesh

圖6 楔形體底部測(cè)點(diǎn)布置圖Fig. 6 Point layout at the bottom of the wedge

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

3 計(jì)算結(jié)果分析

板厚的變化會(huì)影響結(jié)構(gòu)彈性，為了研究厚度參數(shù)對(duì)彈性結(jié)構(gòu)響應(yīng)的影響，本文針對(duì)不同板厚的楔形體進(jìn)行了計(jì)算，通過對(duì)比結(jié)構(gòu)的砰擊壓力、變形位移和應(yīng)力分布等參量來分析結(jié)構(gòu)彈性變化對(duì)結(jié)構(gòu)動(dòng)態(tài)響應(yīng)的影響。

圖7 測(cè)點(diǎn)實(shí)驗(yàn)值和數(shù)值仿真結(jié)果對(duì)比Fig. 7 Comparison of experimental and numerical simulation results

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

圖9 給出了彈性模量（E=2.1E11 Pa）和入水速度（v=10 m/s）保持不變時(shí)，板厚依次增加，楔形體的最大壓力峰值對(duì)比圖。由圖可知，隨著板厚的增加，楔形體的最大壓力極值是逐漸增大的。這是因?yàn)榘搴裨酱笫沟媒Y(jié)構(gòu)剛性越大，質(zhì)量越大，砰擊壓力極值越高。這也從側(cè)面表明了在結(jié)構(gòu)入水砰擊響應(yīng)計(jì)算中，為準(zhǔn)確預(yù)報(bào)砰擊載荷，需要采用流固耦合的數(shù)值模型并計(jì)及結(jié)構(gòu)變形效應(yīng)的影響。

圖8 不同板厚各測(cè)點(diǎn)砰擊壓力曲線對(duì)比Fig. 8 Comparison of slamming pressure of each point with different thickness of plate

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

圖10 不同板厚時(shí)各測(cè)點(diǎn)變形位移Fig. 10 Deformation displacement of each measuring point at different plate thickness

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

圖11 不同板厚P1，P4 測(cè)點(diǎn)變形位移對(duì)比Fig. 11 Comparison of deformation displacement of P1 points and P4 points with different thickness of plate

圖11 為楔形體在不同板厚時(shí)P1，P4測(cè)點(diǎn)的變形位移對(duì)比。由圖可知，不同板厚時(shí)各測(cè)點(diǎn)的位移曲線均為隨時(shí)間的增大而增大；板厚12 mm 時(shí)，位移近似一條直線，而板厚2 mm 時(shí)，位移近似一條曲線，顯示出較為復(fù)雜的變形位移。這是因?yàn)榘搴裨胶?，結(jié)構(gòu)的剛度越大，結(jié)構(gòu)越不易發(fā)生變形。

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

4 結(jié) 語

本文通過基于VOF 和FEM 方法的流固耦合數(shù)值模型，對(duì)三維彈性楔形體進(jìn)行入水仿真，對(duì)比分析了結(jié)構(gòu)彈性變形對(duì)砰擊壓力、結(jié)構(gòu)動(dòng)態(tài)響應(yīng)和結(jié)構(gòu)應(yīng)力分布的影響。

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

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

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

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