翟鋼軍，鮑建宇，馬 哲，孫家文，趙西增

(1. 大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024; 2. 大連理工大學 深海工程研究中心，遼寧 大連 116024; 3. 國家海洋環(huán)境監(jiān)測中心 海域管理技術(shù)重點實驗室，遼寧 大連 116023; 4. 浙江大學 海洋學院，浙江 舟山 316021)

海嘯作為全球最具破壞性的災(zāi)難之一，不僅會造成大量的人員傷亡，還會導致近岸結(jié)構(gòu)物的嚴重破壞，造成財產(chǎn)上的巨大損失。而海堤作為重要的防波堤工程，在近海海域得到了廣泛的應(yīng)用[1]。陳杰等研究了海堤破碎區(qū)至沖瀉區(qū)的水動力特性[2]，此外還關(guān)注了海嘯波作用下的泥沙運動[3]。海嘯波具有巨大的破壞性和強非線性，如何精準地模擬海嘯波與海堤等結(jié)構(gòu)物的相互作用，成為了水動力學領(lǐng)域?qū)W者們的研究熱點。

目前，基于解析解和物理模型試驗的研究方法已相對成熟，如房克照等[4]建立了基于高階Boussinesq水波方程的波浪傳播數(shù)值模型，可以有效處理波浪自由面間斷流動的問題；任冰等[5]開展物理模型試驗，研究了規(guī)則波在島礁地形上的傳播變化；趙西增等[6]開展了潰壩波與海堤相互作用的模型試驗，高速相機能夠較好地捕捉到自由面大變形現(xiàn)象。

在數(shù)值模擬方法的研究方面，郭立棟等[7]基于BEM與VOF波浪耦合模型開展了斜坡堤越浪研究，只用VOF模型很難模擬不規(guī)則波面的非線性流動；王東旭等[8]基于OpenFOAM網(wǎng)格法程序，數(shù)值模擬了孤立波在潛礁上的全傳播過程，但缺少了對液面湍動變形處的細節(jié)模擬?？梢钥闯觯谀M波浪局部液面大變形方面，傳統(tǒng)的網(wǎng)格法(如FDM、VOF法等)仍存在模擬精度不高的問題，液面局部劇烈變形會使網(wǎng)格發(fā)生畸變、扭曲，一定程度上降低該方法的精度[9]，與其他數(shù)值方法耦合或網(wǎng)格的局部加密進一步增加了數(shù)值計算的難度。

近年來，相關(guān)學者在采用潰壩波近似模擬海嘯波的研究中取得了一系列成果。Chanson等[10]通過解析解公式推導和自由波面擬合，驗證了潰壩波與海嘯波涌潮在傳播過程中存在著相似性。此外，由于潰壩波物理模型試驗的可重復性高，技術(shù)條件成熟，隨后越來越多的學者開展了基于潰壩波模擬海嘯波的研究。如Mcgovern等[11]基于SPH法模擬了海嘯波的越浪和爬升過程；Wei等[12]基于SPH法利用潰壩波模擬海嘯波，研究了海嘯波與橋墩的相互作用問題；趙西增等[6]建立潰壩波試驗?zāi)Ｐ湍M海嘯涌浪；張學瑩等[13]應(yīng)用SPH法對二維和三維潰壩波與結(jié)構(gòu)物相互作用過程進行了模擬分析；陳橙等[14]基于潰壩波機制模擬海嘯波沖擊碼頭的模型試驗。SPH法作為新興的無網(wǎng)格粒子法，在提高自由液面大變形計算精度方面取得了長足進展。SPH粒子具有自適應(yīng)性，流體粒子包含速度、壓強、位置等全部信息，無需計算對流項和輸運項，從而消除了自由界面上的數(shù)值耗散，這使得SPH法能夠清晰準確地追蹤自由液面?？偨Y(jié)以上研究成果，下文進一步關(guān)注不同海堤坡角情況下，SPH法邊界條件對數(shù)值求解結(jié)果的影響。

在SPH數(shù)值模型中采用潰壩波數(shù)值模型，基于趙西增等開展的潰壩波物理模型試驗數(shù)據(jù)展開對比驗證。潰壩波的波高與下游水深比值接近于1.64時達到了波浪強非線性條件。通過建立潰壩波數(shù)值模型，對比試驗室中潰壩波沖擊海堤的物理模型試驗結(jié)果，精準地再現(xiàn)了海嘯波的傳播過程，同時再現(xiàn)了堤前液面大變形和堤后射流、崩破波的出現(xiàn)與演變過程。海堤坡角改變時，傳統(tǒng)網(wǎng)格法需進行不同的網(wǎng)格局部加密處理，這里基于無網(wǎng)格法的獨特優(yōu)勢，利用SPH法數(shù)值模擬方法研究了坡角對波浪爬高與衰減以及消波系數(shù)的影響。

1 理論基礎(chǔ)及邊界條件的建立

SPH法是一種無網(wǎng)格粒子法，用一系列任意分布的粒子來離散問題域，通過核函數(shù)積分進行估值近似，一般采用顯式格式對離散的線性方程組進行時間積分[15]，采用標準蛙跳(leapfrog)法，使時間步長Δt滿足CFL條件，保證數(shù)值求解結(jié)果的收斂性。SPH法中構(gòu)造的光滑函數(shù)在支持域上要滿足正則化條件和緊支性條件，同時要保證其對稱性和光滑性。SPH數(shù)值模型基于開源數(shù)值程序SPHERAv.8.0 (RSE SpA)建立[16]，目前常用核函數(shù)有B樣條核函數(shù)、高斯型核函數(shù)、三次樣條核函數(shù)等，下文采用B樣條核函數(shù)進行模擬。

1.1 基本原理

在SPH方法中，函數(shù)f(x)的積分表達式：

(1)

若函數(shù)f(x)在積分域Ω內(nèi)連續(xù)，可以用光滑函數(shù)W(x-x′,h)來代替，此時方程f(x)的積分表達式可寫為：

(2)

式中：f(x)是三維坐標向量x的函數(shù)，Ω為包含x的積分域，即粒子法的整個求解域，δ(x-x′)為狄拉克δ函數(shù)，W(x-x′,h)為核函數(shù)，(x-x′)為粒子間距，h為光滑長度。

1.2 求解N-S方程

拉格朗日形式下流體的連續(xù)性方程：

(3)

簡化后的動量方程：

(4)

(5)

式中：P為粒子壓力，ρ為流體粒子密度，t為時間，g為重力加速度，v為粒子運動速度，Θ為紊動項，v0為水的運動黏性系數(shù)，τ為應(yīng)力張量。

1.3 建立固邊界條件

固邊界的設(shè)置直接決定了邊界附近流場的穩(wěn)定性和計算域的計算精度，一直以來確定合適的邊界條件都是SPH法的研究熱點與難點，主要研究的三種固邊界處理方法如圖1所示。鏡像粒子法最早提出，能夠改善SPH法的邊界缺陷，但每一時間步生成相當數(shù)量的鏡像粒子，極大地增加了計算負擔。邊界排斥力法將固邊界用一系列的粒子來模擬，對靠近固邊界的流體粒子施加作用力，防止流體粒子穿透邊界，但在處理不規(guī)則的邊界條件時邊界力易出現(xiàn)疊加擾動，破壞守恒條件。而海堤坡角改變時，邊界力的形式、大小及方向都會隨之改變，這兩種邊界處理條件都會引起初始擾動，破壞守恒條件。

圖1 三種固邊界處理方法Fig. 1 Three methods of solid boundary treatment

固壁粒子法設(shè)置簡單，適用于復雜幾何固邊界問題，不受作用力大小及方向的限制，文中不同坡角的海堤前坡邊界條件即采用固壁粒子法處理。其中，位于邊界附近的流體粒子i受到邊界粒子k的作用力計算如下，pi是流體粒子的壓力，pk是固邊界粒子的壓力，初始排斥力如式(6)計算。

(6)

基于半解析法邊界條件，對初始作用力求一階偏導，利用粒子的自適應(yīng)性迭代求解各個時刻的邊界排斥力及壓強：

(7)

(8)

(9)

此時，基于半解析固壁粒子邊界條件的SPH連續(xù)性方程推導為式(10)所示，其中Cs為流固耦合項。

(10)

該邊界條件保證了流體粒子不能穿越固邊界，基于Monaco等[17]提出的半解析邊界處理方法，在連續(xù)性方程中引入流固耦合項，設(shè)置為無滑移邊界條件，符合真實流體在水槽中的運動情況。

2 粒子間距敏感度及波面分析

數(shù)值模型水槽尺寸為600 cm×40 cm×60 cm(長×寬×高)，通過擋板將水槽分為長200 cm的蓄水段和長400 cm的試驗段。防波堤模型為斜坡式，下底長45 cm，上底長20 cm，高11 cm，前端距閘門128 cm。上游初始水深20 cm，下游初始水深為與防波堤等高，h=11 cm，試驗布置如圖2所示。

數(shù)值模型參考物理模型試驗布置，物理模型試驗數(shù)據(jù)參考文獻[6]。

圖2 模型試驗布置示意Fig. 2 Layout drawing of model experiment

2.1 粒子間距敏感度分析

本節(jié)設(shè)置三組數(shù)值模擬算例，討論SPH法粒子初始間距對于模擬精度的影響。粒子初始間距分別設(shè)置為0.002 m、0.005 m和0.010 m，對應(yīng)的粒子數(shù)分別為203 015個、32 171個、8 043個。通過分析對比模擬結(jié)果發(fā)現(xiàn)，采用0.010 m的粒子間距時，由于粒子間距過大，初始粒子數(shù)量過少，數(shù)值模擬已無法形成具有清晰輪廓的自由液面圖，故此粒子間距設(shè)置不予采納。

數(shù)據(jù)處理時采用無量綱坐標(x/H,y/H)表示橫、縱坐標，其中H為上游初始水深，取0.20 m。試驗?zāi)M結(jié)果和SPH法數(shù)值模擬波面對比如圖3所示，其中t0時刻代表試驗過程中水頭靠近試驗標記點(x=2.72 m,z=0.11 m)且標記點處液位上升1 cm時刻，經(jīng)試驗記錄結(jié)果確定，時刻t0為0.56 s。圖中圓點代表試驗數(shù)據(jù)，海堤邊界范圍用黑線所畫區(qū)域表示。

在t0+0.56 s時刻，海堤模型上底前端液位開始升高，液體即將進入海堤頂部；在t0+0.85 s時刻，水頭前端運動到海堤頂部右端，此時水頭前端面變形明顯，涌水即將越過海堤右端，隨即越過海堤形成射流，波面演變示意圖如圖3所示。

圖3 試驗與數(shù)值模擬波面結(jié)果對比Fig. 3 Comparison of wavefront change between experiment and numerical simulation

由數(shù)值模擬結(jié)果與試驗數(shù)據(jù)的對比發(fā)現(xiàn)，數(shù)值模型粒子間距設(shè)置為0.002 m時，能更好地擬合液面流態(tài)的劇烈變化；粒子間距設(shè)置為0.005 m時，在潰壩波液面峰值處出現(xiàn)擬合誤差，未能達到較高的模擬精度。由波面結(jié)果對比圖可以看出，粒子間距設(shè)置為0.002 m時，數(shù)值模擬與試驗數(shù)據(jù)的擬合精度最高，數(shù)值模擬7 s海嘯波傳播過程在計算機上需要運行10 h左右，計算成本可接受。在綜合考慮模擬精度和計算成本后，文中數(shù)值模擬采用的粒子間距均為0.002 m。

2.2 海嘯波傳播過程波面分析

通過對比數(shù)值模型與試驗中海嘯波越堤前后的波面形態(tài)，分析海嘯波傳播的波面演變過程，捕捉海嘯波與海堤相互作用時出現(xiàn)的一系列非線性現(xiàn)象。數(shù)值模擬結(jié)果如4左圖所示，高速相機所拍試驗照片如4右圖所示，自上至下分別為t0時刻、t0+0.4 s時刻、t0+0.8 s時刻及t0+1.2 s時刻，試驗及數(shù)值模擬海嘯波翻越海堤過程的波面對比如圖4所示。

圖4 海嘯波翻越海堤過程的數(shù)值模擬及試驗對比圖(試驗照片引自文獻[6])Fig. 4 Comparison of numerical and experimental result during the tsunami wave propagation (from reference [6])

分析試驗數(shù)據(jù)圖像可知，自t0時刻起，潰壩波開始進入高速相機的拍攝范圍內(nèi)，如圖4(a)所示；t0+0.4 s時刻，潰壩涌水峰值開始接近海堤，波峰模擬現(xiàn)象明顯，海堤模型上底前端液位將逐漸升高，涌水將進入海堤頂部，如圖4(b)所示；t0+0.8 s時刻，相對平緩的潰壩水頭前端逐漸升高并變得陡峭，波浪運動至海堤頂部距前端約1/3距離處，如圖4(c)所示；隨后潰壩波沿著海堤頂部繼續(xù)向下游運動，在t0+1.1 s時刻越過海堤頂部下游頂點形成射流，如圖4(d)所示，在形成射流初期，下游近海堤處出現(xiàn)了明顯的液面凹陷，在重力的作用下，凹陷處上部流體隨后閉合，產(chǎn)生了小范圍內(nèi)的崩破波流態(tài)。

通過采用物理模型試驗與數(shù)值模擬相結(jié)合的方法，開展了上述潰壩波模擬海嘯波與海堤相互作用的研究工作。研究結(jié)果表明，SPH法數(shù)值模型不僅可以完整地再現(xiàn)海嘯波翻越海堤的全過程，而且能夠準確地捕獲非線性的自由面大變形現(xiàn)象，驗證了SPH法數(shù)值模型的可靠性和準確性。

3 坡角角度對于海嘯波傳播的影響分析

《海港水文規(guī)范JTS 145-2-2013》中給出了斜坡堤無防浪墻的越浪計算公式：

(11)

式中：Q為單位時間單位堤頂寬度的越浪量，m為海堤前坡坡度的余切值。由式(11)可以看出，在入射波浪有效波高HS、譜峰周期Tp、海堤堤腳處水深d等波浪參數(shù)一致時，海堤前坡度也是越浪量大小的決定因素之一，本節(jié)主要討論海堤坡角對海嘯波消波效果的影響。

3.1 波浪爬高分析

為保證波浪要素輸入條件的一致性，數(shù)值模擬中統(tǒng)一設(shè)置上游潰壩水深為0.2 m，下游水深為0.11 m，在保證海堤堤頂高程相同的前提下，對比建立了7組不同坡角的算例，討論海堤坡角對消波效果的影響?！斗啦ǖ淘O(shè)計與施工規(guī)范JTS 154-2018》中規(guī)定了不同護面型式下最小坡角與最大坡角分別為33.7°及51.3°。本文根據(jù)工程規(guī)范與斜坡堤鋪設(shè)的實際需求，在規(guī)定坡角范圍內(nèi)做加密處理，選取坡角增量為5°，將坡角分別設(shè)置為20°、30°、35°、40°、45°、50°和60°，其中將20°和60°極限坡角作為對比參考值，共計算對比7組不同坡角情況，比較坡角選取范圍的合理性和適用性，模型設(shè)置如圖5所示。

圖5 不同坡角海堤結(jié)構(gòu)模型局部設(shè)置示意Fig. 5 Local layout of submerged breakwater models at different slope angles

圖6 x=3.450 m處波浪爬高變化曲線局部放大圖Fig. 6 Partial enlarged detail of changing curve during the wave run-up at x=3.450 m

在海堤前x=3.450 m處設(shè)置探針，如圖6所示為探針監(jiān)測到的水位變化歷時曲線局部放大圖。從初始時刻開始，液面隨著海嘯波來浪而升高，在t=1.40 s時刻出現(xiàn)最大峰值水位；在t>1.40 s后，由于斜坡面波浪反射的影響，有效越浪量逐漸減少，探針處監(jiān)測的液面高度也將持續(xù)下降。為了分析海堤前波浪爬高現(xiàn)象，獲得探針位置處不同坡角下的液面抬升情況，波浪爬高量通過對比堤前水深h=0.110 m得到，整理坡角角度、波浪爬高峰值和峰值增長率的對應(yīng)關(guān)系如表1所示。

表1 波浪爬高增長率Tab. 1 Growth rate of wave run-up

由監(jiān)測圖像可以看出，坡角對爬高峰值的大小影響顯著，隨著坡角的增大，波浪爬高峰值增大，經(jīng)由海堤斜坡邊界造成的波浪反射不斷增加。由于波浪爬高的大小會對堤后越浪有著顯著的影響，進一步對比了波浪爬高的增長率，統(tǒng)計了平均爬高峰值增長率如表1所示。由表1可知，當坡角增量為5°時，隨著坡角的增大，爬高峰值不斷提升，增長率隨之增大，爬高峰值受坡角變化的影響明顯，坡角是影響波浪爬高和波浪反射大小的主導性因素。

3.2 堤后波浪衰減分析

隨后探討波浪越過海堤后的液面變化情況，這將作為海堤消波性能是否良好的重要指標之一。數(shù)值水槽長6 m，液面監(jiān)測探針2設(shè)置在x=4.22 m處，故不考慮波浪在水槽右側(cè)的反射問題，峰值在t=2.0 s時刻附近出現(xiàn)，通過不同角度海堤的消波系數(shù)對比，分析坡角變化對堤后波浪衰減的影響。

由表2中不同坡角海堤的波浪峰值對比可以看出，坡角對堤后波浪衰減影響顯著，隨著坡角的增大，堤后監(jiān)測點處波浪峰值減小，波浪在越堤后傳播時出現(xiàn)了明顯的衰減現(xiàn)象。在規(guī)定的坡角范圍內(nèi)，當坡角增量為5°時，隨著坡角的增大，峰值差值不斷增大，消波系數(shù)隨之增大。進一步分析圖7中不同角度海堤消波系數(shù)變化對比圖，在20°～40°小坡角范圍內(nèi)，增加坡角對提升消波效果作用明顯，此時坡角是影響消波效果的主導因素之一；在40°～45°范圍內(nèi)，隨著坡角的增加，消波系數(shù)達到0.18左右，此時增長緩慢；在45°～60°大坡角范圍內(nèi)，消波系數(shù)小幅提高，對消波效果提升有限。但考慮到大坡角情況下拋石斜坡堤的施工難度，此時不應(yīng)超出規(guī)范坡角范圍，而要優(yōu)化海堤的堤頂寬度、超高等其他因素，進一步增強消波效果。

圖7 不同坡角海堤的消波系數(shù)對比Fig. 7 Comparison of wave elimination coefficients at different angles

表2 不同坡角海堤的波浪峰值對比
Tab. 2 Comparison of wave-eliminating effects of seawalls with different slope angles

前坡角度/(°)堤前波浪峰值/m堤后波浪峰值/m峰值差值消波系數(shù)200.1700.1510.0190.112300.1760.1500.0260.148350.1780.1480.0300.167400.1800.1480.0320.178450.1810.1480.0330.182500.1820.1470.0350.192600.1840.1450.0380.206

如圖8所示，進一步分析波浪衰減現(xiàn)象逐步加劇的原因?？梢园l(fā)現(xiàn)在波浪越堤傳播到堤后淺水區(qū)時，其波峰變陡，波浪發(fā)生卷曲然后破碎，堤后波浪波陡較大且底坡平緩，此時出現(xiàn)了小范圍內(nèi)的崩破波流態(tài)。淺水沖刷最低位置上方出現(xiàn)波浪回流，回流水體落下后形成裹氣波浪的沖擊波形態(tài)[18]，此時能夠清晰地捕捉到翻轉(zhuǎn)水舌和裹氣現(xiàn)象，崩破波流態(tài)下會產(chǎn)生劇烈能量耗散。隨著波浪爬高的增加，堤前波浪的勢能增大，堤后淺水沖刷最低位置繼續(xù)降低，出現(xiàn)更為顯著的崩破波形態(tài)，造成波浪更加急劇的能量耗散。這也解釋了隨著海堤坡角的遞增，波浪爬升高度增大，消波系數(shù)呈現(xiàn)增大趨勢，堤后波浪衰減現(xiàn)象更加明顯的原因。

圖8 不同時段崩破波演化數(shù)值模擬圖Fig. 8 Numerical simulation of spilling breaker evolution from 1.76 s to 1.84 s

4 結(jié) 語

研究海嘯波與海堤相互作用的水動力學問題，并將數(shù)值模擬結(jié)果與模型試驗數(shù)據(jù)對比，得到結(jié)論如下：

1)根據(jù)數(shù)值模擬精度和計算總時長，粒子間距設(shè)置為0.002 m能夠有效捕捉自由液面，同時引入固壁邊界條件，減少坡角邊界條件改變對計算域精度的影響；

2)準確地模擬了海堤對海嘯波傳播過程的影響，捕捉到射流、崩破波等一系列波浪強非線性現(xiàn)象，對比物理模型試驗中的液位變化圖，驗證了數(shù)值模型的可靠性與準確性；

3)在波浪要素、堤頂水深等要素相同時，海堤坡角對海嘯波的爬高與衰減影響顯著，在20°～60°的海堤坡角范圍內(nèi)，隨著坡度的增加，堤前波浪爬高增大，捕捉到劇烈的崩破波現(xiàn)象，波浪衰減進一步加劇。

4)當坡角較小時，隨著坡角的增大，消波系數(shù)增長明顯，防災(zāi)減災(zāi)效果提升；繼續(xù)增大坡角時，消波系數(shù)增長緩慢，此時應(yīng)考慮優(yōu)化海堤結(jié)構(gòu)，進一步提升消波效果。