基于CIP方法對波浪沖擊過程的數值模擬

孫宏月，孫昭晨，梁書秀，任喜峰

(大連理工大學 海岸與近海工程國家重點實驗室，大連 116024)

波浪對結構物的沖擊是造成海洋結構受損的重要原因之一，因此波浪沖擊荷載的預測是結構設計的重要參考指標。影響波浪沖擊荷載的主要因素有[1-2]：入射波的規模及入射波與結構物之間的角度；入射波浪的破碎與結構物的相對位置關系；波浪沖擊時卷入空氣的方式及卷入的空氣量。以往對沖擊問題的研究主要集中前兩點，相應的研究重點在于結構的整體受力分析，而通常采用的預測方法為半經驗公式法或物理模型實驗法，如李玉成等[3-4]分別研究了不規則遠破波和近破波對直墻的波浪力；周益人等[5]研究了透空式水平板受波浪上托力的沖擊；以及駱俊彬等[6]對基于不規則波對直立式建筑物沖擊的物理模型實驗結果的概率統計分析，擬合出相應的經驗公式。

近年來，隨著計算機運算速度的發展，數值模擬方法正成為計算波浪沖擊荷載的有效工具。尤其為了對波浪沖擊中的細節問題有更深入的認識，開展了大量的基于數值波浪水槽對波浪沖擊問題的研究工作。如VOF方法：1999年，王永學等[7]建立了基于VOF方法的數值波浪水槽，并研究了直墻上波浪荷載的分布；任冰等[8]將此數值波浪水槽應用于非線性波對水平板的沖擊作用研究；2004，任冰等[9]驗證了VOF方法在計算不規則波沖擊問題中的可行性及準確性；隨后，丁兆強[10]將VOF方法應用于計算波浪對透空式結構的沖擊。無網格光滑粒子法(SPH)：Sun等[11]在2010年及Gao等[12]在2012年分別將SPH方法應用于計算波浪對水平板的沖擊問題；2014年，Didier等[13]應用SPH方法研究了波浪對直墻的沖擊；2018年，任喜峰等[14]基于改進的CDP-SPH方法模擬了規則波對開孔沉箱的沖擊?；陂_源CFD類庫OpenFOAM開發的數值模型：2012年，Jacobsen等[15]開發了基于函數法造波及消波的數值波浪水槽；2014年，Chen等[16]驗證了OpenFOAM類庫在研究波浪與結構物相互作用問題中的準確性；Hu等[17]在2017年應用基于OpenFOAM的數值波浪水槽研究了破碎波對直墻結構的沖擊作用。CIP方法：Hu等[18]2008年應用CIP方法計算了波浪與結構物的相互作用問題，討論了CIP方法在求解非線性問題中的優越性；2012年Zhao等[19]對比實驗結果，驗證了CIP方法在求解畸形波與浮體相互作用問題中的準確性；王佳東[20]基于CIP方法建立的數值水槽模型充分討論了孤立波與水平板及垂向板相互作用問題。

以上所述數值波浪水槽模型都基于不可壓縮假設，忽略了氣體可壓縮性的影響。單就波浪的一次沖擊來看，當沖擊過程中有大量氣體卷入時，卷入的空氣會降低波浪沖擊荷載的極值[21]。由此看來，以往忽略空氣相的影響對沖擊荷載的預測可能偏安全，但Takahashi等[22]以及Peregrine等[23]提出空氣泡在降低沖擊極值的同時也增加了沖擊荷載的持續時間。Lugni等[24]通過實驗指出，當波浪對直墻沖擊過程中卷入氣體時，由于氣體的可壓縮性，造成沖擊壓強出現劇烈的震蕩，這與不可壓縮流沖擊的結果有很大不同，有可能造成結構更大的損傷。因此，在研究波浪沖擊問題時，有必要考慮空氣相的存在對結構的影響。

Hu等[25]最先提出了CIP方法可應用于求解可壓縮兩相流沖擊問題；Yang等[26]將其用于楔形體入水的砰擊模擬；Sun等[27]則充分討論了基于CIP方法的可壓縮兩相流模型在計算兩相流沖擊問題的準確性，指出該方法可以較好的模擬出受氣泡影響產生的壓強震蕩等現象。因此，本文旨在基于前人工作，將CIP方法用于構建可壓縮兩相流數值波浪水槽模型，并將本模型應用于計算波浪對直墻及水平板結構的沖擊作用，分析空氣層對沖擊荷載的影響。

1 數值方法

1.1 控制方程

假定計算區域為包含水和空氣的可壓縮粘性兩相流，為簡化計算，忽略溫度的影響，僅考慮連續性方程和Navier-Stokes方程作為控制方程組

(1)

(2)

式中：ρ為流體的密度；u為流體的速度矢量；p為流體壓強；μ為流體的動力黏性系數；g為流體所受重力。

(3)

1.2 基于CIP的分步計算法

對于由方程(1)，(2)和(3)組成的控住方程組應用分步法進行求解計算

(1)對流項求解

(4)

(5)

(6)

對流項的求解采用CIP方法，具體步驟詳見Yabe等[28]。由此可得中間值u*，ρ*，p*。

(2)擴散項求解

(7)

式(7)的求解采用中心差分法，得到中間速度u**。

(3)壓力場求解

(8)

(9)

(10)

對式(9)求散度，與式(10)相結合，得到如下壓力泊松方程

(11)

對方程(11)采用超松弛迭代法求解，得到pn+1，進一步通過方程(8)與方程(9)求得速度un+1和密度ρn+1。

1.3 體積函數

選取體積函數φm作為各相流體在計算網格內所占的體積分數，其中對于水的體積函數φ1求解如下輸運方程

(12)

求解方法采用類VOF方法的高精度插值曲線THINC方法[18]。本文中固體作為剛體處理，固體相φ3的求解通過固體邊界在計算網格所占比例自動求解。氣體相的計算如下

φ2=1-φ1-φ3

(13)

各物理量λ可以由各相的體積分數來確定：λ=φ1λ1+φ2λ2+φ3λ3，式中λ可以是流體密度ρ，動力黏性系數μ和流體中的聲速Cs。

2 兩相流數值水槽的建立與驗證

圖1 數值水槽造波及消波邊界示意圖Fig.1 Sketch of the relaxation and absorbing zone in the numerical wave tank

2.1 造波及消波方法

為了使數值水槽的造波邊界具有消除二次反射的功能，本文選取了松弛區(Relaxation Zone)造波法。在松弛區內，通過波浪的解析解對數值解進行修正以得到目標波要素。如圖 1所示，左側松弛區內生成波浪并吸收遇到結構物后反射回的波浪。松弛區內的修正采用如下松弛函數

(14)

式中：χR為松弛因子。在松弛區內

φ=αRφC+(1-αR)φT

(15)

式中：φ為速度u和體積分數φm。其中下標C和T分別為計算值和解析解的目標值。

同理為了使右側邊界具有消波功能，選取吸收區(Absorbing Zone)進行消波，具體方法與松弛區造波法類似，只是最終的目標值為靜止水面。造波及消波方法的具體步驟詳見Jacobsen等[15]。

2-a 波面抬高的y向網格收斂性驗證2-b 波面抬高的x向網格收斂性驗證圖2 波面抬高的無量綱時間過程線Fig.2 Dimensionless time history of free surface elevation

2.2 數值水槽的驗證

為驗證數值水槽的計算性能，本文首先對規則波的網格收斂性進行了充分的驗證。受篇幅所限，本文僅選取其中一組結果呈現。水槽結構如圖1所示，選取的數值試驗參數分別為水深D=0.5 m，入射波波高H=0.05 m，周期T=1.2 s。計算網格選用x向均勻，y向不均勻網格(其中加密處為自由液面附近)。本文分別對x向及y向網格的收斂性進行了驗證，圖2給出了距離造波區4 m處的無因次波面抬升過程與二階Stokes波的對比。分別從圖2-a與2-b中可以看出，三種尺寸的x向與y向網格結果均與理論值吻合，表明了該模型的具有較好的網格收斂性。

3 規則波沖擊過程的數值模擬

3.1 波浪對直墻的沖擊

圖3 數值波浪水槽及直墻結構示意圖Fig.3 Schematic of the numerical wave tank and the vertical wall

直墻式建筑是海洋工程中常見的結構形式，Didier等[13]曾在LNEC波浪水槽中進行了多組直墻式防波堤物理模型實驗。本文選取了與Didier等人的物理實驗相同的模型尺寸及工況進行了數值模擬。模型的具體結構形式及參數如圖3所示，實驗中造波方法為推板式造波，斜坡前沿距離造波板3.62 m。其中，試驗水深D=0.3 m，規則波周期T=1.3 s，波高H=0.1 m。如圖3所示，在直墻上布置一排點壓力計，點壓力計a高出上平臺0.055 m，從下至上，每個點壓力計(b，c，d)間距均為0.055 m。數值水槽的造波方法選用上節所訴的松弛區造波法，其余參數與實驗保持一致。其中x向網格在直墻結構附近加密，最小網格△x=0.2 m，y方向選用均勻網格，△y=0.4 m。

圖4給出了本文基于CIP方法的數值結果與實驗數據及SPH模擬結果[13]關于沖擊壓強歷時曲線的對比。波浪經過前兩次較小的沖擊后，開始出現明顯的波浪沖擊過程。波浪的沖擊造成壓強在較短時間內達到壓強極值，壓強達到極值隨后快速下降。但從圖中可以看出，沖擊后的荷載依然持續約有0.5 s，最后回落，這一段的壓力主要與水體所受慣性力有關。受慣性力影響，水體沿墻體繼續爬坡，然后回落，接著下一次的沖擊開始。由圖4可以看出，對于a測點，基于CIP方法的數值結果優于SPH的模擬結果，尤其對于沖擊壓強極值及后續沖擊荷載的持續過程的模擬都更接近實驗數據。且由于a測點始終位于水下，不同于其它測點會經歷離水后的0值過程，其壓強歷時始終在變化且有負壓段出現。對于其他測點(b，c，d點)，本文對于沖擊荷載過程的模擬都略大于SPH的結果，與實驗值接近，整體再現了實驗結果中壓強的歷時過程。對于圖中a，b測點在3.5 s前，數值結果略大于實驗值，以及7 s后的沖擊壓強極值小于實驗結果，這主要與本文選取的造波方法不同于實驗有關。

整體來看，本文數值結果與物理模型實驗結果的吻合較好，表明本數值波浪水槽模型可以應用于進行波浪沖擊問題的計算與模擬。

4-a a點4-b b點4-c c點4-d d點圖4 波浪沖擊壓強的CIP數值模擬結果與實驗及SPH方法[13]的對比Fig. 4 Comparison of CIP-based simulated pressure results with experiments and SPH method[13]

3.2 波浪對水平板的沖擊

為研究空氣層在波浪沖擊荷載中的影響，本文模擬了一組規則波對水平板的沖擊物理模型實驗[29]。模型實驗在環境水槽中進行，水槽全長22 m，寬0.45 m，實驗水深為0.35 m。實驗采用的水平板為厚0.6 cm，板寬0.44 m，長1.2 m的有機玻璃板。在水平板的中線上從前向后布置P1～P24共計24個點壓力計，如圖 5所示。其中P1點距離平板前沿2.5 cm，后續點間距均為5 cm。選取的規則波的波長L=2.4 m，波高H=8 cm。數值波浪水槽及結構布置如圖6所示，由于松弛區造波法可消除二次反射，為減少計算區域長度，水平板放置在造波區3～4倍波長處，計算中其余參數均與實驗相同。為提高計算精度，x向及y向網格在自由水面及水平板附近進行加密，其中最小網格分別為△x=0.5 cm，△y=0.25 cm。

圖5 平板底部壓力傳感器布置圖(cm)Fig.5 Sketch of pressure gauge underneath the flat plate圖6 數值波浪水槽及平板結構示意圖Fig.6 Schematic of the numerical wave tank and the flat plate

任喜峰等[29]的物理模型實驗中，對同一組工況共進行了6次重復試驗。由文獻中給出了其中三次沖擊壓強的對比可知，波浪對水平板沖擊壓強的歷時過程具有周期性，但同時也具有一定隨機性。本文選取文獻中提到的第三次(Trial 3)試驗結果與數值結果進行對比分析。圖7～圖9分別為實驗中水平板前(P1點)、中(P13點)、后(P24點)各點沖擊壓強的歷時過程與基于CIP方法的可壓縮及不可壓縮數值結果的對比。由圖可以看出，數值結果很好的模擬了波浪沖擊荷載的發展過程。對于水平板前端(P1點)：波浪的沖擊造成壓強迅速達到較大極值，經歷一個較短的震蕩期(這一過程主要受空氣墊層收縮與擴張的影響)后，進入約為0.3倍波周期的荷載階段；之后，隨著水體的回落，壓強開始下降，直至出現一定的負值；當水體完全脫離平板后，壓力計讀數歸零，約0.5倍波周期。對于平板中段(P13點)：沖擊荷載經歷兩種交替出現的歷時過程，第一種與P1點類似，但負值的歷時過程更長，且受封閉空氣層的影響，壓力計讀數不會有歸零階段；第二種歷時過程較第一種復雜，沖擊壓強震蕩更為劇烈，且出現兩次峰值。對于平板后部(P24點)：沖擊荷載較小且呈現規律的升高降低過程，但受水、氣強烈摻混的影響，會有較大的沖擊壓強極值出現。此處數值模型對于較大沖擊壓強的模擬小于實驗值，這可能與實驗中的湍流現象及水、氣摻混中產生的非常小的氣泡對壓強結果的影響有關。

從圖中對比可以看出，可壓縮模型對于壓強極值的模擬結果更接近實驗值，尤其對于受空氣墊層影響較大的P1點及P13點。由此看來，不可壓縮數值波浪水槽對于卷氣沖擊的模擬可能會高估沖擊壓強的極值。

7-a P1點實驗結果 7-b P1點可壓縮數值模型結果7-c P1點不可壓縮數值模型結果圖7 P1點壓力的CIP數值模擬與實驗結果的對比Fig.7 Comparison of CIP-based simulated pressure results with experimental data at P1

8-a P13點實驗結果 8-b P13點可壓縮數值模型結果8-c P13點不可壓縮數值模型結果圖8 P13點壓力的CIP數值模擬與實驗結果的對比Fig.8 Comparison of CIP-based simulated pressure results with experimental data at P13

9-a P24點實驗結果 9-b P24點可壓縮數值模型結果9-c P24點不可壓縮數值模型結果圖9 P24點壓力的CIP數值模擬與實驗結果的對比Fig.9 Comparison of CIP-based simulated pressure results with experimental data at P24

為了更好的說明空氣墊對沖擊壓強的影響，圖10給出了前、中、后各點沖擊壓強歷時過程的細節圖(基于可壓縮數值波浪水槽模型)，與之對應圖11和圖12給出了造成P13點兩種不同沖擊荷載過程對應的氣墊層的形態發展過程：圖11中氣墊層相對較寬，呈現扁平狀；圖12中氣墊層較窄，呈現渾圓狀。在t=17.04 s時，P1點及P13點同時達到最大沖擊壓強值，此時自由液面分布如圖11所示，波浪剛剛沖擊到水平板前端，此時平板中段產生的沖擊壓強主要受空氣墊層的影響；隨后在平板底部形成一個完整的封閉的扁平氣墊層，如t=17.14 s時刻所示。t=17.4 s時，隨波浪的向后推進，P24點出現一次小幅壓強升高；在t=17.64 s時，水體的回落造成P1點與P13點幾乎同時達到負壓強極值；而在t=17.7 s 時，P13點與P24點出現第二次壓強峰值，這與水體在向后運動中再一次對平板中后部產生沖擊有關。圖12為另一種氣墊層的形態發展過程。在t=18.5 s時，P1和P13點同時達到最大沖擊壓強峰值；t=18.71 s時，P1、P13點的出現第二次的沖擊荷載，此時平板底部形成一個相對渾圓的氣墊層；t=19.13 s時，P1點受水體脫離平板的影響出現負壓強值；隨著水體的運動，約在t=19.42 s后，封閉的氣墊才從平板后部溢出。

10-a P1點沖擊壓強 10-b P13點沖擊壓強10-c P24點沖擊壓強圖10 各點沖擊壓強歷時曲線細節圖Fig.10 Detailed time history of impact pressure

數值結果中模擬出了交替出現的壓強歷時過程及兩種氣墊層的形態發展，這一現象與任喜峰等[29]的實驗結果一致。由此看出氣墊層對沖擊壓強除了有緩沖作用的影響外，也會使平板中后部受到多次沖擊，以致有二次壓強極值的出現。

圖11 水平板底部形成扁平氣墊層時的自由面變化過程 Fig.11 Evolution of the free surface when a flat air cushion formed under the horizontal plate圖12 水平板底部形成渾圓氣墊層時的自由面變化過程Fig.12 Evolution of the free surface when a concentrated air pocket formed under the horizontal plate

4 結論

本文建立了基于CIP方法的可壓縮兩相流數值波浪水槽模型，為消除造波邊界的二次反射，采用松弛區造波法。首先，對數值波浪水槽模型的網格收斂性進行了充分的驗證；然后，在基礎上對波浪與直墻的相互作用進行了模擬。波浪對直墻沖擊荷載的計算結果與實驗值吻合較好；與SPH方法的比較可以看出，本模型在沖擊壓強極值及砰擊荷載持續過程的模擬有一定優勢。

應用基于CIP方法的可壓縮兩相流數值水槽模型模擬了規則波對水平板的沖擊，數值結果準確的再現了波浪沖擊荷載的壓強歷時過程。同時將本模型結果與實驗及不可壓縮數值結果進行了比較：可壓縮模型對波浪沖擊極值的模擬結果比不可壓縮模型更小，與實驗值更接近，尤其在受空氣墊層影響較大的水平板的中前段。這一差異說明在波浪沖擊問題中，應用不可壓縮模型可能會高估有氣墊層沖擊的壓強極值；同時，可壓縮模型對于沖擊中出現的負壓強值的模擬更接近實驗的結果。因此，有必要在波浪沖擊問題的模擬中考慮空氣的可壓縮性。

最后，對水平板底部自由水面發展過程與沖擊壓強的歷時曲線進行了對比，分析認為沖擊荷載不單受到波浪形態的影響，也受到平板底部空氣墊層發展的影響。數值結果很好的模擬了實驗中兩種形態的空氣墊層交替出現的現象。

綜上，基于CIP方法的可壓縮數值波浪水槽為波浪沖擊過程提供了新的、有效的模擬工具。