限制性航道中船行波傳播特性的數(shù)值研究

王孟飛，鄧 斌,2,3，蔣昌波,2,3，伍志元,2,3

(1.長(zhǎng)沙理工大學(xué) 水利工程學(xué)院，湖南 長(zhǎng)沙 410114； 2.洞庭湖水環(huán)境治理與生態(tài)修復(fù)湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410114； 3.水沙科學(xué)與水災(zāi)害防治湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410114)

1 研究背景

限制性航道是指因水域狹窄、斷面系數(shù)小而對(duì)船舶航行有明顯限制作用的航道。限制性航道中船行波具有波高大、流速急等特點(diǎn)[1]，易破壞航道環(huán)境[2]和航道邊坡[3]或?qū)ζ渌嫠ㄖ镌斐筛蓴_或損害[4]。目前，我國(guó)內(nèi)河通航標(biāo)準(zhǔn)中關(guān)于限制性航道尺度的確定，僅給出雙航線下的最小取值(限制性航道斷面系數(shù)大于6)，航道尺寸的確定僅依據(jù)航跡帶寬度附加一定安全距離，沒(méi)有充分考慮航道尺寸與船舶航行條件的影響。為防止船行波對(duì)航道及周圍船舶產(chǎn)生影響，往往僅規(guī)定航道內(nèi)船舶的限制航速[5]。近年來(lái)，隨著“一帶一路”的深入推進(jìn)，我國(guó)內(nèi)河航運(yùn)的任務(wù)日益繁重，船舶大型化、快速化的發(fā)展特點(diǎn)，使得現(xiàn)有航運(yùn)條件日益緊張，加大對(duì)限制性航道的開(kāi)發(fā)與恢復(fù)通航越來(lái)越受關(guān)注。因此，深入開(kāi)展限制性航道通航關(guān)鍵技術(shù)研究極為必要，特別是限制性水域條件中船行波的相關(guān)研究。

目前，國(guó)內(nèi)外很多學(xué)者通過(guò)現(xiàn)場(chǎng)調(diào)查[3,6]、物理實(shí)驗(yàn)[7-8]及數(shù)值模擬[9-10]等方法對(duì)船行波的傳播特性[11, 12]、船舶的穩(wěn)定性[13]以及船行波對(duì)周圍護(hù)岸[3]、泥沙[2]及船舶本身[4]的影響進(jìn)行了大量的研究。船行波有首波系、尾波系、散波和橫波等不同分類形式，通常是隨弗汝徳數(shù)Fr的不同呈現(xiàn)不同的運(yùn)動(dòng)形式[11]。部分學(xué)者針對(duì)開(kāi)闊水域條件建立了可描述船行波波高與弗汝德數(shù)、離船舷距離、船型等因素的相關(guān)關(guān)系的數(shù)學(xué)模型[14-15]。但限制性水域條件下船行波的傳播特性不同于開(kāi)闊水域[16]。因此，越來(lái)越多的學(xué)者研究限制性水域條件內(nèi)船行波的相關(guān)特性。如：Zhu等[17]解釋了淺水條件下單船體與雙船體船的首波系與尾波系的不同傳播規(guī)律。Ji等[1]通過(guò)求解Navier-Stokes方程和采用不同的湍流方法分析了不同船型和航速下船行波波高的演化規(guī)律，并討論了船舶尾流的水動(dòng)力特性。駱婉珍等[18]基于CFD技術(shù)應(yīng)用重疊網(wǎng)格的方法模擬分析狹窄、淺水航道內(nèi)船行波波形的疊加、反射等現(xiàn)象，但該類方法計(jì)算量大，難以求解大船速下船行波的傳播演化[2]。

近年來(lái)，隨著Boussinesq方程不斷的改進(jìn)，如消除了淺水的限制又考慮了波浪的破碎，減小了數(shù)值計(jì)算量[19]，越來(lái)越多的學(xué)者運(yùn)用Boussinesq方程模擬波浪的傳播變形等。如：Dam等[16]通過(guò)建立二維Boussinesq模型，模擬細(xì)長(zhǎng)船行波及其傳播過(guò)程，并得到不同船速下波高的變化規(guī)律。Shi等[20]基于完全非線性的Boussinesq模型FUNWAVE-TVD模擬了船行波引起的紊動(dòng)破碎，并針對(duì)兩種不同的破碎準(zhǔn)則(黏性破碎準(zhǔn)則和波高破碎準(zhǔn)則)進(jìn)行對(duì)比分析，結(jié)果表明兩種破碎標(biāo)準(zhǔn)均能準(zhǔn)確的模擬船行波在不同F(xiàn)r下的傳播過(guò)程。卓明泉等[21]基于COULWAVE開(kāi)源程序代碼實(shí)現(xiàn)淺水航道船行波的數(shù)值模擬，再現(xiàn)了淺水條件下船首孤立波現(xiàn)象。然而，上述研究大都基于Boussinesq模型下船行波傳播特性的研究，較少對(duì)船舶吃水以及跨臨界航速范圍內(nèi)船行波的特性進(jìn)行詳細(xì)描述。

綜上所述，盡管現(xiàn)有研究對(duì)限制性水域內(nèi)船行波的傳播規(guī)律有了一定的認(rèn)識(shí)，但基于限制性水域內(nèi)船舶航行誘導(dǎo)的船首孤立波波高以及橫向流速的研究相對(duì)較少。因此，為彌補(bǔ)現(xiàn)有研究的不足，本研究運(yùn)用完全非線性的Boussinesq模型FUNWAVE-TVD[22-23]，以壓力源項(xiàng)近似表示船舶引起自由面的波動(dòng)，采用Shi等[20]提出的破碎準(zhǔn)則，重點(diǎn)討論限制性水域內(nèi)船行波的相關(guān)特性。

2 數(shù)值模型建立

基于完全非線性Boussinesq方程建立了限制性航道內(nèi)船舶航行的數(shù)值模型，具體控制方程如下。

質(zhì)量守恒方程：

(1)

(2)

動(dòng)量守恒方程為：

(3)

渦流黏度如下所示：

v=Bδ2(h+η)ηt

(4)

式中：δ為混合長(zhǎng)度系數(shù)，本研究取δ=1；B為避免數(shù)值振蕩的常數(shù)[19]，本研究取B=1。

船舶模型的實(shí)現(xiàn)由壓力源確定，即在完全非線性Boussinesq模型中加入壓力源項(xiàng)。其中，船舶中心點(diǎn)為(x*,y*)的壓力分布可由公式(5)～(7)給出：

(5)

(6)

(7)

數(shù)值模型基于三階Runge-Kutta方法的自適應(yīng)時(shí)間步長(zhǎng)和MUSCL-TVD數(shù)值解法；干濕邊界網(wǎng)格最小深度為0.001 m，CFL數(shù)設(shè)置為0.25，通量計(jì)算采用HLL(Harten Lax & van Leer)構(gòu)造方法，數(shù)值消波技術(shù)采用耗散型海綿層消波方法。

3 數(shù)值模型工況設(shè)置與驗(yàn)證

3.1 數(shù)值模型工況設(shè)置

參考我國(guó)《內(nèi)河通航標(biāo)準(zhǔn)》(GB 50139-2014)規(guī)范的要求，模型計(jì)算區(qū)域采用Ⅵ級(jí)單線航道尺寸，船模采用Ⅵ級(jí)航道標(biāo)準(zhǔn)駁船(船長(zhǎng)L=32 m，寬B=7.0 m，滿載設(shè)計(jì)吃水p為1.0 m，船舶形狀系數(shù)α和β均為0.5)，并在計(jì)算區(qū)域左側(cè)設(shè)置消波區(qū)域。考慮到船行波的傳播在航道內(nèi)通常呈對(duì)稱分布，計(jì)算在航道一側(cè)布置15個(gè)測(cè)點(diǎn)，數(shù)值模型計(jì)算平面布置圖如圖1所示。

圖1 數(shù)值模型計(jì)算平面布置圖(單位：m)

考慮Ⅵ級(jí)航道內(nèi)標(biāo)準(zhǔn)駁船的空載吃水p≈0.4～0.5 m，滿載吃水p≈0.9～1.0 m，以及標(biāo)準(zhǔn)航行船速，共設(shè)置5種計(jì)算工況(A1～A5)，擬定的模型計(jì)算工況如表1所示。

表1 模型計(jì)算工況設(shè)置

3.2 數(shù)值模型驗(yàn)證

為驗(yàn)證數(shù)值模型計(jì)算限制性航道內(nèi)船行波的可靠性，分別從數(shù)值實(shí)例和網(wǎng)格收斂性驗(yàn)證兩個(gè)方面進(jìn)行詳細(xì)說(shuō)明。二維驗(yàn)證算例取自參考文獻(xiàn)[7]，實(shí)驗(yàn)水槽寬D為3.5 m，水深h為0.114 m；船舶長(zhǎng)度L為1.6 m，寬度B為0.4 m，吃水p為0.1 m，船舶的形狀系數(shù)α和β均為0.5；初始時(shí)刻船舶位于x=6 m、y=1.75 m的位置，其他參數(shù)設(shè)置均與參考文獻(xiàn)[7]保持一致。

圖1給出了3種工況下測(cè)點(diǎn)(x=33.125 m，y=0.72 m)波高歷時(shí)曲線的對(duì)比，結(jié)果顯示，數(shù)值模型計(jì)算波高與文獻(xiàn)[7]中船行波的變化趨勢(shì)相吻合。采用skill數(shù)驗(yàn)證模型精度要求：

(8)

圖2為船行波波高計(jì)算值的實(shí)驗(yàn)驗(yàn)證結(jié)果。計(jì)算結(jié)果表明各工況下skill數(shù)均大于0.90，均方根誤差均小于0.038，說(shuō)明模型能較好地模擬船行波的傳播。此外，圖2(a)～2(c)黑色方框內(nèi)為船首孤立波，這是典型淺水航道內(nèi)船行波的重要特征。

3.3 網(wǎng)格收斂性驗(yàn)證

為探究網(wǎng)格大小對(duì)計(jì)算結(jié)果的影響，確保所有航速下CFL數(shù)均小于0.25。選取滿載吃水工況A4、Fr=1.2下x=80 m、y=1 m處4種不同大小網(wǎng)格的船行波波高歷時(shí)曲線進(jìn)行收斂性驗(yàn)證，驗(yàn)證結(jié)果見(jiàn)圖3。從圖3可以看出，各波高歷時(shí)曲線基本吻合，其中，以dx×dy=0.2 m2網(wǎng)格大小為參考，網(wǎng)格dx×dy=0.19、0.21、0.22 m2的波高歷時(shí)曲線的均方根誤差分別為0.004，0.004，0.008 m；同一位置處不同網(wǎng)格下船行波最大波高相對(duì)誤差小于3%?？梢?jiàn)，模型在風(fēng)格大小dx×dy=0.2 m2下已達(dá)到計(jì)算精度要求。因此，模型計(jì)算采用dx×dy=0.2 m2的網(wǎng)格尺度。

圖2 不同F(xiàn)r數(shù)船行波波高模擬值的實(shí)驗(yàn)驗(yàn)證結(jié)果

圖3 4種不同大小網(wǎng)格的收斂性驗(yàn)證

4 結(jié)果分析與討論

4.1 波面的演化及最大波高

由于船行波以Fr=1為界呈現(xiàn)不同的傳播特性[7]，其中，在跨臨界航速(Fr=1)下船行波中橫波和散波疊加易形成較大波高。因此，選取工況A3、Fr=1分析船舶運(yùn)動(dòng)過(guò)程中船行波波面的歷時(shí)演化過(guò)程，結(jié)果見(jiàn)圖4。由圖4可看出，在船舶航行的起始階段(t=5 s)，船舶前方形成彎曲的孤立波波列面。隨著船舶的運(yùn)動(dòng)，彎曲的波列面開(kāi)始變直，船舶前方孤立波數(shù)量逐漸增多。與此同時(shí)，船舶后方的尾波呈現(xiàn)出開(kāi)爾文尾流模式，船尾波傳播受側(cè)壁反射、疊加等影響，在限制性航道兩側(cè)及航道中間各點(diǎn)處(如：x=100 m，t=15 s)易形成較大波高，從而對(duì)近岸航道邊坡等設(shè)施產(chǎn)生不利影響。

為深入分析船首孤立波的傳播特征，以A3、Fr=0.75下船行波波高歷時(shí)曲線為例，定義船行波中船首孤立波和船尾波的最大波高，如圖5(a)所示。由于限制性航道內(nèi)船行波的傳播易形成復(fù)雜的橫向流動(dòng)，對(duì)航道兩岸設(shè)施及船舶航行造成嚴(yán)重的影響。因此，為分析近岸橫向流動(dòng)的分布規(guī)律，定義最大近岸橫向流速如圖5(b)所示。

圖6為工況A3、Fr=0.75下距離近岸1 m處不同位置(x=80, 100, 120, 140, 160 m)的船行波波高歷時(shí)曲線圖。從圖6可以看出，在淺水條件下，船首孤立波波高在80 m

圖5 最大波高及近岸橫向流速的定義

圖6 距離近岸1 m處不同位置(x=80, 100, 120, 140, 160 m)的船行波波高歷時(shí)曲線(A3、Fr=0.75)

圖7為工況A3、Fr=0.75下距離近岸不同位置(y=1, 2, 3 m)的船行波波高歷時(shí)曲線圖。可見(jiàn)船尾波波高在最大值出現(xiàn)之后易受航道寬度及邊壁影響，且越靠近岸壁船尾波波高越大。另外，隨著船首孤立波不斷向前傳播，船首孤立波逐漸形成新的次孤立波，造成航道水位的上升。

4.2 船行波最大波高分布

基于上述最大船首孤立波波高的定義，圖8(a)給出了不同F(xiàn)r數(shù)條件下最大船首孤立波波高隨相對(duì)吃水的分布規(guī)律。可以看出，隨著相對(duì)吃水的增加，最大船首孤立波波高呈線性規(guī)律增加，斜率大致為0.4～0.7。圖8(b)給出了不同吃水條件下最大船首孤立波波高隨Fr數(shù)的變化規(guī)律?？梢钥闯觯S著Fr的增加，最大船首孤立波波高呈冪函數(shù)形式增加，且冪次方約為3.5。

基于上述分析，通過(guò)多元非線性回歸方程得到了p、Fr與最大船首孤立波波高的經(jīng)驗(yàn)關(guān)系式：

(9)

式中：ηmax為最大船首孤立波波高，m；h為水深，m；Fr為弗汝德數(shù)；p為船舶吃水，m。

方程式(9)各系數(shù)標(biāo)準(zhǔn)差見(jiàn)表2，多元非線性回歸方程擬合優(yōu)度R2=0.967。為了更加直觀地展現(xiàn)多元回歸方程的擬合精度，圖8(c)給出了數(shù)值計(jì)算與多元回歸方程式(9)的結(jié)果對(duì)比，可見(jiàn)方程式(9)可以較好地預(yù)測(cè)不同吃水條件下最大船首孤立波的波高。

表2 船首孤立波回歸方程的系數(shù)

圖9為不同吃水條件下最大船尾波波高隨Fr的變化規(guī)律。由圖9可以看出，當(dāng)p≤0.7 m時(shí)，最大船尾波波高以Fr=0.9為界呈現(xiàn)不同的分布規(guī)律。當(dāng)Fr<0.9時(shí)，最大船尾波波高逐漸增加；相反，當(dāng)Fr>0.9時(shí)，最大船尾波波高逐漸減小，這是由于當(dāng)Fr>0.9時(shí)船尾波破碎導(dǎo)致波高減小，這與David等[24]的研究大體一致，但David等并未討論船舶相對(duì)吃水與最大船尾波波高之間的影響關(guān)系。然而，當(dāng)p>0.7 m時(shí)，最大船尾波波高隨Fr增加逐漸增加，與上述p≤0.7 m下最大船尾波波高呈現(xiàn)不同的分布規(guī)律，這是由于當(dāng)船舶吃水較大時(shí)，同一波高下最大船尾波波長(zhǎng)相對(duì)較長(zhǎng)，船尾波不易破碎。

4.3 船行波最大近岸橫向速度

影響限制性航道設(shè)計(jì)的主要因素除最大船行波波高外，還受航道內(nèi)橫向速度大小的影響。圖10為限制性航道內(nèi)不同吃水條件下最大近岸橫向流速隨Fr的變化規(guī)律。從圖10中可以看出，受船行波破碎影響，最大近岸橫向流速僅在p≤0.7 m、Fr≤0.9下整體呈線性關(guān)系，具體為隨Fr，p的增加分別呈線性減少和線性增加的趨勢(shì)。

圖7 距離近岸不同位置(y=1, 2, 3 m)的船行波波高歷時(shí)曲線(A3、Fr=0.75)

圖8 最大船首孤立波波高變化規(guī)律

圖9 不同吃水條件下最大船尾波波高隨Fr的變化規(guī)律 圖10 不同吃水條件下最大近岸橫向流速隨Fr的變化規(guī)律

5 結(jié) 論

基于完全非線性的Boussinesq方程建立了限制性航道內(nèi)船舶運(yùn)動(dòng)的數(shù)值模型，重點(diǎn)研究了船首孤立波、船尾波以及近岸橫向流速在不同p、Fr下的分布規(guī)律。與傳統(tǒng)開(kāi)闊水域及深水條件下船舶航行的數(shù)值模型相比，該模型能更合理的描述限制性航道中船舶航行過(guò)程中波高、流速等分布規(guī)律。具體結(jié)論如下：

(1) 基于我國(guó)《內(nèi)河通航標(biāo)準(zhǔn)》(GB 50139-2014)規(guī)范的要求，采用FUNWAVE-TVD模型建立了精確描述船行波傳播及演化的數(shù)值模型。研究表明該模型能很好地模擬限制性航道條件下船行波的傳播特征。

(2) 本文研究條件下，最大船首孤立波波高與船舶吃水之間呈線性關(guān)系，斜率約為0.4～0.7之間；最大船首孤立波波高與Fr之間呈冪函數(shù)關(guān)系，冪次方約為3.5；此外，基于上述規(guī)律與分析建立了限制性水域條件下最大船首孤立波波高與船舶吃水、弗汝德數(shù)的冪函數(shù)經(jīng)驗(yàn)公式。

(3) 船尾波波高及近岸橫向流速受波浪破碎影響，以p=0.7 m和Fr=0.9為界分別呈現(xiàn)不同的分布規(guī)律。其中，F(xiàn)r主要通過(guò)船行波波高影響船行波的破碎，而p主要通過(guò)船行波波長(zhǎng)影響船行波的破碎。同樣，受波浪破碎影響最大近岸橫向流速僅在p≤0.7 m、Fr≤0.9內(nèi)整體呈線性分布規(guī)律。

本研究基于完全非線性的Boussinesq方程、壓力源項(xiàng)以及黏性破碎準(zhǔn)則模擬了限制性航道內(nèi)船首波系以及船尾波系的傳播演化過(guò)程，為進(jìn)一步開(kāi)展限制性航道中船舶安全通航提供了理論基礎(chǔ)，但本研究?jī)H考慮一定船型、航道斷面尺寸、靜水條件下船行波的演變規(guī)律，更復(fù)雜的航道條件、船舶運(yùn)動(dòng)以及水動(dòng)力特性之間的內(nèi)在關(guān)系還有待進(jìn)一步研究。