基于非線性壓力損失邊界條件的半圓型開孔潛堤水動(dòng)力特性解析研究*

呂祚睿，劉 勇，李華軍

(中國(guó)海洋大學(xué)海岸與海洋工程研究所，山東 青島 266100)

潛堤淹沒于水下，不影響海岸景觀，具有抵御波浪、防止岸線侵蝕等作用。半圓型潛堤的堤身為半圓型拱圈，具有穩(wěn)定性好、適用軟弱基礎(chǔ)，構(gòu)件受力性能好等優(yōu)點(diǎn)。借鑒Jarlan[1]提出的開孔沉箱結(jié)構(gòu)，對(duì)半圓型潛堤進(jìn)行表面開孔，可以加強(qiáng)其對(duì)波浪能量的耗散作用，有效降低結(jié)構(gòu)承受的波浪力。

日本于1990年代初在宮崎港建成了半圓型開孔防波堤的試驗(yàn)段，Tanimoto和Takahashi[2]建議利用修正后的合田公式計(jì)算半圓型防波堤的波浪力。1997年，我國(guó)在天津港設(shè)計(jì)建成了半圓型防波堤，采用合田公式修正后的經(jīng)驗(yàn)公式計(jì)算結(jié)構(gòu)波浪力；隨后，在長(zhǎng)江口深水航道治理工程中修建了半圓型導(dǎo)堤，并基于5階斯托克斯波的數(shù)學(xué)模型分析和物理模型試驗(yàn)，提出了適用于半圓型潛堤的結(jié)構(gòu)波浪力計(jì)算公式[3]。一些學(xué)者通過物理模型試驗(yàn)和數(shù)值模擬對(duì)半圓型潛堤的水動(dòng)力特性進(jìn)行了深入研究。Dhinakaran等[4]試驗(yàn)研究了半圓型開孔潛堤在隨機(jī)波作用下的反射系數(shù)、波浪爬高及受力特性，建議結(jié)構(gòu)的最優(yōu)開孔率為11%。Knoll等[5]利用鋼筋網(wǎng)制作半圓型開孔潛堤模型，試驗(yàn)研究了結(jié)構(gòu)在不規(guī)則波作用下的受力特性。Forbes 和 Schwartz[6]以及Forbes[7]在不同水深和弗勞德數(shù)下研究了均勻流通過半圓型潛堤時(shí)的流動(dòng)特性。Kasem和Sasaki[8]通過數(shù)值求解粘性流體運(yùn)動(dòng)方程，研究了波浪對(duì)半圓型潛堤的作用。Yuan和Tao[9]將數(shù)值分析結(jié)果與模型試驗(yàn)結(jié)果進(jìn)行對(duì)比分析，總結(jié)了不同淹沒狀態(tài)下半圓型潛堤的受力特性，提出結(jié)構(gòu)波浪力的簡(jiǎn)化計(jì)算公式。

除上述模型試驗(yàn)和數(shù)值分析外，Liu和Li[10]利用多極子展開法建立了波浪對(duì)半圓型開孔潛堤作用的精確解析解，利用線性開孔邊界條件[11]描述波浪通過開孔板的能量損失(壓力損失)。多極子展開法由Ursell[12]首先應(yīng)用于求解水波問題，具有收斂性好、計(jì)算速度快的優(yōu)點(diǎn)，可用于流體對(duì)水平圓柱、圓球、橢圓球等結(jié)構(gòu)作用的解析求解。Ursell[12]通過在流場(chǎng)中布置奇點(diǎn)，建立了極坐標(biāo)系下拉普拉斯方程的奇異解，并在此基礎(chǔ)上推導(dǎo)了不同邊界條件下的速度勢(shì)表達(dá)式。Thorne[13]進(jìn)一步采用更為簡(jiǎn)潔的方法得到了Ursell的結(jié)果。Martin[14]給出多極子方法中滿足收斂條件的結(jié)構(gòu)物面方程。Linton 和McIver[15]對(duì)多極子展開法進(jìn)行了比較系統(tǒng)全面的總結(jié)。

Liu和Li[10]在解析研究中采用線性開孔邊界條件，求解過程中無法直接考慮波高變化對(duì)開孔板處能量損失的影響。本文在Liu和Li[10]的研究基礎(chǔ)上，進(jìn)一步采用非線性壓力損失邊界條件[16-18]，利用多極子展開法建立波浪對(duì)半圓型開孔潛堤作用的新解析解，深入分析半圓型潛堤的水動(dòng)力特性。本文所建立的新解析解及相應(yīng)的分析結(jié)果可為工程設(shè)計(jì)提供重要參考。

1 理論分析

1.1 控制方程和邊界條件

圖1為波浪對(duì)半圓型開孔潛堤作用的示意圖。半圓型潛堤的半徑為a，水深為h；波浪從左向右傳播，波高為H，周期為T。以潛堤中垂線與靜水面的交點(diǎn)為坐標(biāo)原點(diǎn)，建立直角坐標(biāo)系(x,o,z)，以(0,-h)為原點(diǎn)建立極坐標(biāo)系(r,θ)，定義極軸垂直向下時(shí)θ=0。極坐標(biāo)系與直角坐標(biāo)系之間滿足：rcosθ=-(z+h)，rsinθ=x。沿潛堤的半圓弧面將整個(gè)流域劃分為兩個(gè)部分，區(qū)域①為潛堤外部流場(chǎng)，區(qū)域②為潛堤內(nèi)部流場(chǎng)。

假定該流場(chǎng)內(nèi)為無粘、無旋、不可壓縮的理想流體，考慮線性規(guī)則波對(duì)結(jié)構(gòu)的作用，速度勢(shì)可以表示為：

Φ(x,z,t)=Re[φ(x,z)e-iωt]。

(1)

速度勢(shì)滿足拉普拉斯方程：

(2)

速度勢(shì)還滿足自由水面條件、水底條件和遠(yuǎn)場(chǎng)輻射條件：

(3)

(4)

(5)

其中k為波數(shù)，滿足K=ω2/g=ktanhkh(g為重力加速度)；φI為入射波速度勢(shì)，可以表示為：

(6)

入射波速度勢(shì)在極坐標(biāo)系下可以表示為：

(7)

在潛堤的開孔半圓弧處，內(nèi)外流場(chǎng)的法向速度連續(xù)：

(8)

其中，下標(biāo)1和2分別表示外部和內(nèi)部流場(chǎng)內(nèi)的變量。此外，流體運(yùn)動(dòng)滿足非線性(二次)壓力損失邊界條件[16-18]：

(9)

其中：ε為半圓弧的開孔率(開孔面積與總面積的比值)；μ為射流系數(shù)；C為開孔薄板的阻塞系數(shù)(Blockage coefficient)[19-20]。Molin[16]指出：通常情況下，射流系數(shù)μ在開孔率ε<50%時(shí)約為0.3～0.4。Suh等[21]的研究表明，C可以表示為

(10)

其中b為開孔薄板的厚度。

式(2)～(9)形成一個(gè)完整的邊值問題，下面利用多極子展開方法進(jìn)行求解。

1.2 速度勢(shì)求解

利用多極子展開法，流場(chǎng)中的復(fù)速度勢(shì)可以寫為對(duì)稱解和反對(duì)稱解之和?；贚iu和 Li[10]的研究，區(qū)域①內(nèi)的滿足方程(2)～(5)的速度勢(shì)φ1可以表示為：

(11)

其中

(12)

(13)

(14)

(15)

利用分離變量法，在區(qū)域②內(nèi)滿足方程(2)和(3)的速度勢(shì)φ2可以寫為：

(16)

需要利用開孔半圓弧處的邊界條件(8)和(9)來確定速度勢(shì)中待定系數(shù)。根據(jù)Linton 和 McIver[15]，利用加權(quán)積分和三角函數(shù)的正交性來求解待定系數(shù)的值。由于原速度勢(shì)表達(dá)式為無窮級(jí)數(shù)，需要先取截?cái)鄶?shù)M，使表達(dá)式變?yōu)橛邢揄?xiàng)。將式(11)和(16)代入式(8)后，取正整數(shù)p=1,2,3，…，M，在式(8)兩端分別乘以sin(2p-1)θ或cos2pθ，可以得到2M個(gè)方程；之后，在(0.5π,1.5π)區(qū)間上對(duì)θ積分，消去θ及部分正交項(xiàng)，將方程簡(jiǎn)化。同理，將式(11)和(16)代入式(9)后，在式(9)兩端分別乘以sin(2p-1)θ或cos2pθ，并在(0.5π,1.5π)區(qū)間上對(duì)θ積分，可以再得到2M個(gè)方程。為使方程個(gè)數(shù)與待定系數(shù)相同，對(duì)式(9)在(0.5π,1.5π)區(qū)間直接對(duì)θ積分，得到p=0時(shí)cos2pθ的加權(quán)方程。最后共得到以下4M+1個(gè)方程：

(17)

p=1,2,3，…，M。

(18)

p=1,2，3，…，M。

(19)

p=0,1,2，…，M。

(20)

其中：

(21)

1.3 水動(dòng)力參數(shù)

(22)

(23)

(24)

(25)

根據(jù)線性伯努利方程，流場(chǎng)內(nèi)任意點(diǎn)的動(dòng)水壓強(qiáng)可以寫為：

p(x,z)=iρωφ(x,z)。

(26)

將動(dòng)水壓強(qiáng)沿結(jié)構(gòu)物表面積分，可以得到波浪力。水平向波浪力為：

(27)

垂向波浪力Fz包括兩部分，即潛堤半圓弧面所受的力Fz1和底板所受的力Fz2：

Fz=Fz1+Fz2；

(28)

(29)

(30)

2 解析解驗(yàn)證

根據(jù)解析研究結(jié)果，確定截?cái)鄶?shù)M后，可通過式(17)～(20)構(gòu)成的非線性方程組求解速度勢(shì)表達(dá)式中的展開系數(shù)，進(jìn)而計(jì)算半圓型開孔潛堤在波浪作用下的各項(xiàng)水動(dòng)力參數(shù)。隨截?cái)鄶?shù)M的增加，需要確認(rèn)本文解析解是否能夠收斂。表1和2列出了無因次波數(shù)kh不同時(shí)，半徑分別為a=0.6h和a=0.8h的半圓型開孔潛堤在不同截?cái)鄶?shù)M下的反射系數(shù)CR。在算例中迭代精度取10-4，波高H=0.1h，開孔率ε=20%，開孔薄板的厚度b=0.03h，射流系數(shù)μ=0.35，阻塞系數(shù)C通過式(10)進(jìn)行計(jì)算。可以看到，當(dāng)M=6時(shí)CR在10-4的精度范圍內(nèi)達(dá)到收斂要求。因此，在后文計(jì)算中均取截?cái)鄶?shù)M=6。

表1 不同截?cái)鄶?shù)M下潛堤反射系數(shù)CR計(jì)算結(jié)果(a=0.6h)Table 1 Reflection coefficient CR with different truncation numberM (a=0.6h)

表2 不同截?cái)鄶?shù)M下潛堤反射系數(shù)CR計(jì)算結(jié)果(a=0.8h)Table 2 Reflection coefficient CR with different truncation numberM (a=0.8h)

將本文解析解的水動(dòng)力參數(shù)計(jì)算結(jié)果與Liu和Li[20]所建立的分區(qū)邊界元方法的數(shù)值計(jì)算結(jié)果進(jìn)行了對(duì)比。圖2給出了半圓型開孔潛堤的反射系數(shù)CR和透射系數(shù)CT在以上兩種方法中的計(jì)算結(jié)果，計(jì)算條件為：半徑a=0.7h，波高H=0.1h，開孔率ε=20%，板厚b=0.03h，射流系數(shù)μ=0.35。分區(qū)邊界元方法中將流場(chǎng)邊界離散為2 300個(gè)邊界單元(直線段)，計(jì)算結(jié)果收斂。圖2中曲線代表本文解析解的計(jì)算結(jié)果，點(diǎn)代表分區(qū)邊界元方法的計(jì)算結(jié)果?？梢钥闯觯簝煞N不同方法的計(jì)算結(jié)果一致，表明本文解析解的求解過程是正確的。

圖2 本文解析解與分區(qū)邊界元解[20]的對(duì)比Fig.2 Comparison between the present analytic solution and the solution of multi-domain BEM[20]

3 算例討論

本節(jié)通過算例分析，研究半徑a和開孔率ε這兩個(gè)主要參數(shù)對(duì)半圓型開孔潛堤水動(dòng)力特性的影響，計(jì)算條件為：板厚b=0.03h，射流系數(shù)μ=0.35，波高H=0.05h。

3.1 半徑a的影響

圖3和4分別給出了半圓型開孔潛堤的相對(duì)半徑a/h對(duì)反射系數(shù)CR和透射系數(shù)CT的影響，計(jì)算中取開孔率ε=20%。從圖3可以看出：隨著無因次波數(shù)kh的增大，潛堤的反射系數(shù)CR先增大后減小，在kh=0.5左右達(dá)到最大值；隨著a/h的增大，CR的最大值不斷增大。從圖4可以看出：隨著波數(shù)的增大，潛堤透射系數(shù)CT先減小后增大；隨著a/h的增大，透射系數(shù)CT顯著減小，并且這一趨勢(shì)在波數(shù)較大時(shí)更為明顯。這說明在水深恒定時(shí)，增大潛堤半徑可以顯著減小波浪透射，提高潛堤的掩護(hù)效果。

圖3 潛堤相對(duì)半徑a/h對(duì)CR的影響(ε=20%)Fig.3 Effect of the relative radius a/h onCRof breakwater

圖4 潛堤相對(duì)半徑a/h對(duì)CT的影響(ε=20%)Fig.4 Effect of the relative radius a/h onCT of breakwater

圖5和6給出半圓型開孔潛堤的相對(duì)半徑a/h對(duì)結(jié)構(gòu)所受波浪力的影響，計(jì)算條件與圖3和4相同。圖中，CFx=Fx/ρgHa和CFz=Fz/ρgHa分別表示水平方向和垂向的無因次波浪力。由圖5可以看出：隨著波數(shù)的增大，潛堤所受水平波浪力先增大后減小，在kh=0.7左右達(dá)到最大值；隨著a/h的增大，潛堤所受水平波浪力的最大值小幅增加。由圖6可知：隨著波數(shù)的增加，潛堤所受垂向波浪力顯著減?。浑S著a/h的增大，潛堤所受垂向波浪力不斷減小。對(duì)比圖5和6可以看出，對(duì)于較長(zhǎng)周期波(波數(shù)較小)，潛堤的垂向波浪力明顯大于水平波浪力；隨著波數(shù)的增大，二者在數(shù)值上的差距不斷縮小。

圖5 潛堤相對(duì)半徑a/h對(duì)CFx的影響(ε=20%)Fig.5 Effect of the relative radius a/h onCFxof breakwater

圖6 潛堤相對(duì)半徑a/h對(duì)CFz的影響(ε=20%)Fig.6 Effect of the relative radius a/h onCFz of breakwater

3.2 開孔率ε的影響

圖7和8分別表示半圓型潛堤開孔率ε對(duì)潛堤反射系數(shù)CR和透射系數(shù)CT的影響，計(jì)算條件為a=0.9h。從圖7可看出：隨著ε的增大，潛堤的反射系數(shù)CR不斷遞減；當(dāng)開孔率ε從 5% 增加到20%時(shí)，CR的最大值從0.5降低到0.2左右。這說明在潛堤表面開孔可以顯著降低堤前反射波高。從圖8可以看到：隨著ε的增大，潛堤的透射系數(shù)CT先減小后增大，存在一個(gè)透射系數(shù)最低的最優(yōu)開孔率。這是因?yàn)?當(dāng)開孔率很小時(shí)，總開孔面積很小，潛堤對(duì)波浪能量的耗散作用較弱；如果開孔率非常大，潛堤的半圓弧面對(duì)流場(chǎng)的影響非常小，能量耗散也非常小。在實(shí)際工程中，合理選擇開孔率，可以有效提高潛堤的掩護(hù)效果。根據(jù)圖8的計(jì)算結(jié)果，當(dāng)kh<2時(shí)，推薦開孔率為10%～15%。

圖7 潛堤開孔率ε對(duì)CR的影響(a=0.9h)Fig.7 Effect of porosityεonCRof breakwater

圖8 潛堤開孔率ε對(duì)CT的影響(a=0.9h)Fig.8 Effect of porosityεonCTof breakwater

圖9和10給出半圓型潛堤開孔率ε對(duì)結(jié)構(gòu)所受波浪力的影響，圖中CFx和CFz的定義與圖5和6相同，計(jì)算條件與圖7和8相同。根據(jù)圖9，隨著開孔率ε的增大，水平波浪力CFx顯著減小，當(dāng)開孔率ε從 5% 增加到20% 時(shí)，CFx的最大值從0.6減小到0.2。這說明在半圓型潛堤表面開孔可以顯著降低潛堤所受水平向波浪力，提高潛堤的穩(wěn)定性。對(duì)比圖10和9可以看出：與水平波浪力不同，垂向波浪力CFz受開孔率變化的影響較小。

圖9 潛堤開孔率ε對(duì)CFx的影響(a=0.9h)Fig.9 Effect of porosityεonCFx of breakwater

圖10 潛堤開孔率ε對(duì)CFz的影響(a=0.9h)Fig.10 Effect of porosityεonCFzof breakwater

為進(jìn)一步研究開孔率ε對(duì)半圓型開孔潛堤水動(dòng)力特性的影響規(guī)律，圖11直接給出潛堤水動(dòng)力參數(shù)隨開孔率ε的變化曲線，圖中計(jì)算條件為：kh=1.5，a/h=0.9。可以看到，當(dāng)ε=10%左右，能量損耗系數(shù)CL達(dá)到最大值，透射系數(shù)CT達(dá)到最小值(掩護(hù)效果最好)，這與圖8中的分析結(jié)果相符，也與Dhinakaran等[4]的試驗(yàn)研究結(jié)果基本一致。

圖11 潛堤開孔率ε對(duì)水動(dòng)力參數(shù)的影響(kh=1.5，a/h=0.9)Fig.11 Effect of porosity ε on the hydrodynamic coefficients of breakwater

4 結(jié)語(yǔ)

本文基于勢(shì)流理論，利用多極子展開法和分離變量法，建立了波浪對(duì)半圓型開孔潛堤作用的解析解。解析解中引入了非線性壓力損失邊界條件，從而可以直接考慮波高變化對(duì)開孔板處能量損耗的影響。本文解析解為級(jí)數(shù)解，具有很好的收斂性。本文解析解的計(jì)算結(jié)果與分區(qū)邊界元方法的計(jì)算結(jié)果一致。算例分析表明：增大半圓型開孔潛堤的半徑會(huì)增大潛堤的反射系數(shù)，降低潛堤的透射系數(shù)，在工程設(shè)計(jì)中應(yīng)根據(jù)波浪條件，合理選擇較大的潛堤半徑；增大半圓型開孔潛堤的開孔率可以顯著降低潛堤的反射系數(shù)和水平波浪力，但是對(duì)潛堤垂向波浪力影響較小；為了達(dá)到良好的消浪、掩護(hù)效果，當(dāng)kh<2時(shí)，半圓型開孔潛堤的推薦開孔率為10%～15%。