呂超凡，趙西增,2，殷銘簡

(1. 浙江大學 海洋學院，浙江 舟山 316021; 2. 浙江海洋大學 海洋工程裝備學院，浙江 舟山 316022)

重力式防波堤雖然能較好地起到擋浪作用，但它在一定程度上阻止港池內(nèi)外水體之間的自然循環(huán)和交換[1-2]，使得港內(nèi)水質(zhì)惡化，產(chǎn)生一系列的環(huán)境問題[3]。為在擋浪的同時，加強港池內(nèi)外水體交換，可在防波堤的水下部分開挖涵洞，形成涵洞式防波堤。利用波能主要集中在水面附近的特點，防波堤的上部結構可以反射波浪，水下涵洞部分可進行水體交換[4]。該措施既可滿足港池的泊穩(wěn)條件，又可提高港池的水質(zhì)；相較于其他常見的水體交換方式，具有結構強度高、經(jīng)濟性好、施工容易等優(yōu)點[2]。

總體而言，針對涵洞式防波堤的消浪及透波特性研究己較為詳盡。但在水體交換更新方面，雖然已有基于淺水方程的大尺度數(shù)值模擬，但其針對的是潮流，波浪作用下的水體交換機理及其規(guī)律尚缺乏相關研究。為此，采用計算流體力學(CFD)數(shù)值模擬[11-13]結合物理模型試驗[14-15]的方式精細化研究涵洞式直立堤在波浪作用下涵洞內(nèi)部水體的交換特性。數(shù)值模擬基于高精度的VPM (volume-average/point-value method)-THINC (tangent of hyperbola for interface capturing)/QQ (quadratic surface representation and Gaussian quadrature)模型[16-17]，該模型可以精確捕捉自由面，重現(xiàn)波浪與涵洞式直立堤的相互作用，為分析其運動機理提供依據(jù)[18]。此外，在其基礎上提出一種能夠標記涵洞內(nèi)外水體變化的雙液相流體體積(VOF)方法[11]。該方法主要是把原有的計算水體分成多個子水體并分別標記，而后計算其各自輸運方程，進而獲取不同水體的運動特征，以此來反映在水動力的作用下，水體內(nèi)部的流動規(guī)律。該方法可以同時獲取涵洞內(nèi)部的水體流動以及整個水體的流動，相對簡單，計算成本低[19]。由此，文中利用上述方法定性以及定量描述涵洞內(nèi)外水體的交換特性，而后結合局部流場揭示水體交換機理，并分析不同涵洞參數(shù)以及波浪參數(shù)對水體交換的影響。

1 數(shù)值方法

1.1 基本方程

CFD數(shù)值模擬采用不可壓縮黏性流模型，其控制方程如下：

(1)

(2)

(3)

式中：U(x,z)為流體質(zhì)點速度，ρ為流體密度，p為相對動壓力，m為動力黏性系數(shù)，U(x,z)為加速度，U(x,z) 為坐標點，F(xiàn)σ(x,z)為表面張力。Fσ表達式如下：

(4)

式中：σ為張力系數(shù)，κ為界面平均曲率，α為流體體積分數(shù)。當網(wǎng)格單元為空氣時，α=0；為水時，α=1；水汽摻混時，0<α<1。另外每個網(wǎng)格的流體特性定義如下：

λ=λ1α+λ2(1-α)

(5)

式中：λ可為網(wǎng)格內(nèi)流體的密度ρ或黏性系數(shù)μ。

其中，Navier-Stokes方程采用高階有限體積法VPM求解，自由面采用高精度的THINC/QQ算法[17]重構。建立基于OpenFOAM底層函數(shù)庫的VPM-THINC/QQ模型，具體可參考文獻[18]。

為實現(xiàn)數(shù)值波浪水槽的造波與消波，采用速度邊界法造波[18]，以及 Jacobsen等[20]提出的松弛法來消除透射波和二次反射波。

1.2 雙液相算法

為直觀反映涵洞式直立堤的水體交換特性，在VPM-THINC/QQ模型的基礎上提出一種能夠模擬港池內(nèi)外水體交換的雙液相VOF方法[19]。該方法假設兩種性質(zhì)相同的液體分別為β1和β2，兩種液體均滿足VOF輸運方程：

(6)

(7)

其中，

β1+β2=α

(8)

(9)

(10)

由此，港池外部以及涵洞內(nèi)部水體標記為β1，港池內(nèi)部水體標記為β2，開展基于雙液相VOF方法的涵洞式直立堤水體交換特性研究。

2 試驗以及CFD數(shù)值模型布置

為定性以及定量分析涵洞式直立堤的水體交換特性，物理模型試驗以及數(shù)值模擬布置如下。

2.1 試驗設備及模型布置

物理模型試驗在波浪水槽內(nèi)進行，水槽尺寸為25 m×0.7 m×0.7 m(長×寬×高)。水槽首端裝有主動吸收式推板造波機，可產(chǎn)生單向二維正弦規(guī)則波，目標波浪生成質(zhì)量不受造波時間影響，水槽尾端配有消波設施來消除波浪反射。模型試驗布置見圖1，試驗水槽中設置6根浪高儀，分別為EG1～EG6，所在位置依次為5 m、10 m、11 m、13 m、16 m、17 m，防波堤模型置于14 m處。模型采用亞克力板制作，防波堤長度為B=0.5 m，涵洞高度及深度可調(diào)節(jié)，試驗水深h=0.4 m。具體試驗工況如表1所示，其中S為涵洞高度，dS為涵洞中軸線位置，Ti為試驗周期，Hi為試驗波高。

圖1 物理模型試驗布置Fig. 1 Sketch of an experimental setup

表1 模型試驗以及CFD模擬工況設置Tab. 1 Conditions of model experiments and CFD simulations

2.2 CFD數(shù)值模擬布置

CFD數(shù)值波浪水槽布置見圖2，水槽布置與模型試驗布置基本一致。但數(shù)值水槽另外在14.25 m處設置流速測點FV1，用來監(jiān)測整個涵洞斷面的水平流速u(x,z)，進而計算出整個斷面的振蕩流量Qn；同時，在0～4 m處設置波浪造波松弛區(qū)，用來吸收從結構物傳來的反射波浪，18～26 m設置為消波區(qū)。數(shù)值模擬工況與模型試驗一致，具體見表1。

圖2 數(shù)值波浪水槽示意Fig. 2 Schematic diagram of a numerical wave flume

2.3 模型驗證

為驗證數(shù)值模擬的正確性，需將數(shù)值模擬結果與試驗結果進行對比。圖3表示在涵洞深度dS=0.177 m，涵洞高度S=0.096 m 的條件下，CFD結果與試驗結果的波面曲線對比，其中圖3(a)表示波高Hi=0.08 m,波周期Ti=1.57 s的波面對比，圖3(b)表示Hi=0.02 m,Ti=1.98 s的波面對比。從圖3中可以看出，圖3(a)的波浪透射試驗結果以及CFD模擬結果幾乎一致，圖3(b)的透射結果略微偏大主要是由于試驗消波區(qū)略短，存在部分反射波與透射波的疊加，故隨著時間的增長，透射值逐漸變大；反射結果兩者略有差異，但整體差異不大。上述結果表明雙液相VPM-THINC/QQ模型能較好重現(xiàn)實際波浪與涵洞式防波堤的相互作用，具有一定可靠性。

圖3 試驗及CFD波面對比(dS=0.177 m,S=0.096 m)Fig. 3 Comparisons of experimental results and CFD wave surface (dS=0.177 m,S=0.096 m)

為分析涵洞內(nèi)部水體輸移規(guī)律，試驗中，涵洞內(nèi)部水體用染色劑進行染色。該部分染色水體在一段時間內(nèi)，能定性的反映出波浪作用下的水體運動規(guī)律。圖4表示在波高Hi=0.06 m，周期Ti=2.4 s，涵洞深度dS=0.177 m，涵洞高度S=0.046 m條件下試驗與數(shù)值模擬的涵洞內(nèi)部水體輸移軌跡對比。圖4(a)、4(b)分別表示初始時刻以及波浪運動初期數(shù)值模擬條件下的雙液相水體分布情況；圖4(c)、4(d)分別表示初始時刻以及波浪運動初期試驗條件下的染色水體分布情況。從圖4中可以看出，當波浪傳到結構物時，涵洞內(nèi)部水體在波浪的作用下會向防波堤內(nèi)側輸移，并在涵洞口形成近乎對稱的兩個渦旋，逐漸把涵洞內(nèi)部水體帶出。對比試驗與數(shù)值結果發(fā)現(xiàn)，兩者現(xiàn)象幾乎一致，渦旋軌跡及大小也近乎相同，表明雙液相模型能夠較好地模擬出波浪作用下的水體輸移規(guī)律。

圖4 雙液相驗證(Hi=0.06 m,Ti=2.4 s ,dS=0.177 m,S=0.046 m)Fig. 4 Verifications of two-liquid phase (Hi=0.06 m,Ti=2.4 s ,dS=0.177 m,S=0.046 m)

模擬與試驗結果對比表明，文中所提出的雙液相VPM-THINC/QQ模型能夠真實反映波浪作用下涵洞式直立堤的水體交換特性。

3 結果與討論

由于已經(jīng)驗證文中數(shù)值模擬能夠較好重現(xiàn)波浪與涵洞式直立堤的相互作用，故下文結果均為數(shù)值所得。為探究涵洞式直立堤的水體交換特性，定量分析涵洞內(nèi)部的流體運動規(guī)律，定義涵洞內(nèi)14.25 m處的水平振蕩流量Qn如下：

(11)

3.1 涵洞內(nèi)水體交換的機理分析

為分析涵洞式直立堤的水體交換機理，通過分析涵洞附近的雙液相分布、渦量分布以及流線分布規(guī)律對其展開研究。其中雙液相為了揭示涵洞式直立堤的水體交換作用，渦量分布用來說明波浪與防波堤作用機理，流線分布則是渦量的一個補充，用來說明涵洞內(nèi)部水體的振蕩規(guī)律。圖5、6表示波浪與結構物作用穩(wěn)定后，不同周期的4個典型相位下(時間間隔為Ti/4)雙液相分布、渦量分布以及流線分布。

圖5 一個周期內(nèi)涵洞內(nèi)部水體交換特性(Ti=2.40 s,Hi=0.08 m,dS=0.177 m,S=0.096 m)Fig. 5 Water exchange characteristics inside the culvert in a period (Ti=2.40 s,Hi=0.08 m,dS=0.177 m,S=0.096 m)

圖5的試驗條件為波周期Ti=2.40 s，波高Hi=0.08 m，涵洞深度dS=0.177 m，涵洞高度S=0.096 m。由圖5可看出，自入射波從波谷到達波峰的半個周期內(nèi)(圖5(a))，港池內(nèi)側出現(xiàn)了未染色水體而港池外側出現(xiàn)部分染色水體，表明涵洞內(nèi)部以及港池內(nèi)側(透射波方向)水體逐漸被帶入港池外部(來波方向)，期間水體大量摻混，涵洞內(nèi)外水體發(fā)生交換。對比渦旋軌跡(圖5(b))能發(fā)現(xiàn)，當?shù)糖安ɡ颂幱诓ü葧r，涵洞入口(來波方向)附近會形成上下兩個渦，涵洞出口(透射波方向)水面以及水底附近也會形成渦旋。此后，在波浪從波谷到波峰的運動過程中，涵洞入口附近的兩個渦逐步被帶到堤前水面附近，并逐漸消散，同時帶走了大量涵洞內(nèi)部的水體；而防波堤外側水面及水底附近的渦逐漸運動到涵洞出口處，并進入涵洞。觀察流線圖能發(fā)現(xiàn)，涵洞內(nèi)部流線較為平順，僅在波浪轉向(波峰、波谷)時內(nèi)部流線發(fā)生改變，表明在波浪與其作用的過程中，整個水體交換在半個周期內(nèi)(波峰到波谷，或者相反)是穩(wěn)定而平順的，在整個周期內(nèi)則是呈現(xiàn)出振蕩的特性。當入射波從波峰到達波谷的半個周期內(nèi)，涵洞內(nèi)部以及港池外側水體逐漸被帶入港池內(nèi)部，對比渦旋軌跡能發(fā)現(xiàn)，此時渦旋運動規(guī)律與波谷到波峰時的規(guī)律相反。堤前水面以及水底形成的渦逐漸運動到涵洞內(nèi)部，而涵洞出口附近的渦則逐漸被帶入到水面以及水底附近。上述結果表明，涵洞內(nèi)部水體在波浪的驅動下，往復環(huán)流，形成一股振蕩的水流，進而完成水體交換。在此過程中，雖然單次波浪周期所完成的水體交換量不算太大，如圖5的雙液相可以發(fā)現(xiàn)港池外側逐漸出現(xiàn)港池內(nèi)測的水體，但在長時間的波浪作用下，水體交換范圍將逐漸擴大，進而可以較好地實現(xiàn)港池內(nèi)外水體交換的目的。

圖6試驗條件與圖5類似，只是波周期變?yōu)?.18 s，可以看出，其分布規(guī)律與長周期大體相似，但仍有一定差異，主要體現(xiàn)在水體交換強度弱于長周期波浪，同時，涵洞附近的渦脫軌跡也有一些不同。主要由于，在短周期波浪作用下，水體垂向能量部分不均勻，水面附近能量遠遠大于水底，故涵洞入口處僅在涵洞上部形成渦，而后逐漸被帶入到水體表面，如此來回振蕩；涵洞出口處與長波類似，形成上下兩個渦，但由于透過波能較小，其運動軌跡遠遠小于長波，僅在涵洞出口附近來回振蕩。波浪作用下涵洞內(nèi)部的流線分布規(guī)律幾乎一致。

圖6 一個周期內(nèi)涵洞內(nèi)部水體交換特性(Ti=1.18 s,Hi=0.08 m,dS=0.177 m,S=0.096 m)Fig. 6 Water exchange characteristics inside the culvert in a period (Ti=1.18 s,Hi=0.08 m,dS=0.177 m,S=0.096 m)

3.2 涵洞參數(shù)對水體交換的影響

圖7表示當波高Hi=0.02 m，涵洞深度dS=0.177 m時不同涵洞高度S下的涵洞內(nèi)部振蕩流Qn隨時間t的變化曲線，圖7(a)、7(b)分別表示波周期Ti=1.00 s、Ti=2.40 s時的Qn變化曲線。從圖7中可以看出，Qn隨t正弦變化；隨著S的增大，Qn在短周期下隨著高度的增長近乎同比例增大，在長周期下，增大倍數(shù)略小于高度增長率。

圖7 涵洞高度對涵洞內(nèi)部振蕩流Qn的影響(Hi=0.02 m,dS=0.177 m)Fig. 7 The effects of the culvert heights on the oscillatory flow Qn inside the culvert (Hi=0.02 m,dS=0.177 m)

圖8為當波高Hi=0.02 m，涵洞高度S=0.096 m時不同涵洞深度dS下的涵洞內(nèi)部振蕩流Qn隨時間t的變化曲線，圖8(a)、8(b)分別表示波周期Ti=1.00 s、Ti=2.40 s時的Qn變化曲線。從圖8中可以看出，隨著dS的增大，短周期下的Qn不斷變?。婚L周期下的Qn也在變小，但變化幅度很小。主要由于短周期波浪能量集中于水面附近，越靠近水面，涵洞所能透過的波能越多，進而Qn變大；而長周期下的波浪水質(zhì)點速度整個垂向分布較為均勻，改變涵洞位置對Qn的影響不大。

圖8 涵洞深度對涵洞內(nèi)部振蕩流Qn的影響(Hi=0.02 m,S=0.096 m)Fig. 8 The effects of the culvert depths on the oscillatory flow Qn inside the culvert (Hi=0.02 m,S=0.096 m)

上述結果表明，增大涵洞高度能夠明顯地提高涵洞內(nèi)部的流量，促進水體交換；減小涵洞深度對短周期波浪能夠一定程度上提高涵洞內(nèi)部的流量，但對長周期波浪作用很小。

3.3 波浪參數(shù)對水體交換的影響

由前文可知，在小波高情況下，涵洞內(nèi)部振蕩流隨時間呈現(xiàn)出正弦變化，可以將其表示為Qmaxcos(σt)，其中Qmax為振蕩流振幅，σ為波浪圓頻率，t為時間。故在一個周期內(nèi)涵洞內(nèi)部水體的總輸移量(周期積分)等于0。所以應該考慮的是半個周期的平均凈流出，而半個周期的平均凈流出流量Qave與振蕩流振幅存在線性關系：

(12)

圖9為在波高Hi=0.02 m，涵洞高度dS=0.177 m，涵洞深度S=0.096 m條件下，不同波周期Ti下涵洞內(nèi)部振蕩流頻譜。

圖9 不同波周期下涵洞內(nèi)部振蕩流頻譜(Hi=0.02 m,dS=0.177 m,S=0.096 m)Fig. 9 Spectrum of the oscillatory flow in different wave periods (Hi=0.02 m,dS=0.177 m,S=0.096 m)

可看出，隨著Ti的增大，振蕩流基頻振幅AQ不斷增大；除了基頻振幅外，涵洞內(nèi)部存在部分倍頻項，但倍頻振幅較小，與基頻振幅相比幾乎可忽略。同時，發(fā)現(xiàn)波周期Ti=2.40 s (kh=0.55,k為波數(shù),h為水深) 所對應的基頻振幅約為Ti=0.80 s (kh=2.55)的7倍。上述結果表明，波周期對水體交換有較大影響，長周期波浪的振蕩流量遠遠大于短周期波浪。

圖10 波高對涵洞內(nèi)部振蕩流的影響(Ti=0.80 s,dS=0.177 m,S=0.096 m)Fig. 10 The effects of wave heights on the oscillatory flow inside the culvert (Ti=0.80 s,dS=0.177 m,S=0.096 m)

圖11為50 s后不同波高下的雙液相水體分布范圍。可見隨著波浪波高的增大，水體交換范圍明顯擴大；當作用時間足夠時，水體交換能夠逐漸擴散到整個港區(qū)，進而提高水體質(zhì)量。上述結果表明，涵洞式直立堤在波浪作用下能夠有效地實現(xiàn)港池內(nèi)外的水體交換。

圖11 50 s后不同波高的水體交換對比(dS=0.177 m,S=0.096 m)Fig. 11 Comparisons of the water exchange at different wave heights after 50 s (dS=0.177 m,S=0.096 m)

4 結 語

在VPM-THINC/QQ模型的基礎上，提出一種能夠標記涵洞內(nèi)外水體變化的歐拉兩相VOF方法，并結合物理模型試驗，開展了波浪作用下涵洞式直立堤水體交換特性的研究，得到如下結論：

1) 基于VPM-THINC/QQ算法的雙液相VOF方法能夠較好反映真實流體的內(nèi)部流動，具有較高的精度。

2) 涵洞內(nèi)部水體交換機理如下：當防波堤堤前波浪從波谷向波峰運動時，涵洞入口處和涵洞出口的水面以及水底附近形成渦旋，帶動涵洞內(nèi)部水體向堤前運動；當?shù)糖安ɡ藦牟ǚ逑虿ü冗\動時，涵洞入口的水面附近和涵洞出口處形成渦旋，帶動涵洞內(nèi)部水體向堤后運動。故在波浪作用下，涵洞內(nèi)部會形成往復的振蕩流，且該振蕩流主導著港池內(nèi)外的水體交換。

3) 涵洞高度對涵洞內(nèi)部振蕩流有較大影響；涵洞深度對短周期波浪有較大影響，而對長周期波浪影響不大；隨著波浪非線性的增強或波周期的增大，涵洞內(nèi)外水體交換作用加強。同時，隨著非線性的增強，涵洞內(nèi)部存在一個類似于斯托克斯漂移的凈輸移。在波浪長時間的作用下，涵洞式直立堤能夠有效地促進港池內(nèi)外水體交換。