楊師宇 吳靜萍 汪 敏* 陳昌哲

(武漢理工大學(xué)船海與能源動(dòng)力工程學(xué)院1) 武漢 430063)(上海交通大學(xué)船舶海洋與建筑工程學(xué)院2) 上海 200240)

0 引 言

對(duì)于預(yù)報(bào)波浪與結(jié)構(gòu)物的相互作用，隨著計(jì)算機(jī)和數(shù)值模擬技術(shù)的發(fā)展，建立數(shù)值波浪水池(numerical wave tank，NWT)進(jìn)行數(shù)值模擬試驗(yàn)是一種重要的研究手段.

早期，數(shù)值波浪水池技術(shù)忽略流體黏性，流動(dòng)無旋，基于勢(shì)流理論，采用Green函數(shù)方法[1-2]或高階譜(higher-order spectral, HOS)方法[3]求解Laplace方程的初邊值問題，自由面捕捉技術(shù)從線性、弱非線性發(fā)展到強(qiáng)非線性、完全非線性，自由面邊界條件非線性問題一直是提高計(jì)算波浪運(yùn)動(dòng)問題精度的關(guān)鍵.忽略流體黏性的波浪運(yùn)動(dòng)計(jì)算方法，除了求解Laplace方程的邊值問題之外，還引入人工渦黏模型的非靜壓模型淺水方程方法[4-5]、Boussinesq 水波方程方法[6-7]和適用于強(qiáng)非線性波浪問題的層析波浪理論(green-naghdi理論)方法[8-9].

隨著CFD(computational fluid dynamics)的發(fā)展，考慮流體黏性，求解黏性流體運(yùn)動(dòng)的N-S(Navier-Stokes)微分方程，采用VOF(volume of fluid)等方法捕捉自由面，實(shí)現(xiàn)了強(qiáng)非線性的波浪流動(dòng)的模擬[10-11].由于微分方程常用差分方法、有限元方法或有限體積法等方法離散求解，波浪問題需做非定常計(jì)算，計(jì)算量大且耗時(shí).雖然粘-勢(shì)流耦合模型[12]有效提高了計(jì)算效率，但是還是相當(dāng)耗時(shí).雖然波浪與小尺度構(gòu)件相互作用黏性的影響比較重要，但是在處理波浪與大尺度結(jié)構(gòu)物相互作用的問題中，流體的黏性往往被忽略，所以勢(shì)流理論在大尺度的波固耦合問題中仍然在發(fā)展和使用.

基于勢(shì)流理論的數(shù)值計(jì)算方法，發(fā)展較早的是采用Green函數(shù)在流場(chǎng)邊界布置源匯分布的邊界元方法.其按照是否考慮流動(dòng)隨時(shí)間的變化，分為頻域和時(shí)域兩大類方法.相對(duì)于頻域方法對(duì)弱非線性、周期性問題有效的局限性而言，時(shí)域方法應(yīng)用范圍更加廣泛.按照求解的積分方程是否具有奇異性，又出現(xiàn)奇異邊界元方法和去奇異邊界元方法.奇異邊界元方法發(fā)展了高階邊界元法(higher-order boundary element model,HOBEM)[13]和泰勒展開邊界元法(taylor expansion boundary element method，TEBEM)[14]以降低處理奇異積分的難度.去奇異邊界元方法(desingularized boundary integral equation method，DBIEM)最早由Cao等[15]提出，隨后的30多年來，相繼被應(yīng)用于求解與波浪運(yùn)動(dòng)相關(guān)的問題[16-19].

去奇異邊界元方法將源(匯)分布移至流域邊界之外，從而解決了邊界積分方程的奇異問題.與奇異常值邊界元方法相比，去奇異邊界元法在形狀突變處的計(jì)算精度得到了提高[20]，編程與計(jì)算的效率也高[21].但是去奇異邊界元法的計(jì)算精度對(duì)去奇異距離的變化較為敏感，且去奇異距離在不同問題中的最佳取值不盡相同.除此之外，源點(diǎn)數(shù)量、時(shí)域計(jì)算時(shí)的時(shí)間步長(zhǎng)，在不同問題中的取值也對(duì)計(jì)算精度影響較大.

文中采用Cao等[22]的時(shí)域DBIEM構(gòu)建了一種二維數(shù)值波浪水池.自由面邊界條件考慮了弱非線性，采用混合歐拉-拉格朗日法追蹤瞬時(shí)自由面流體質(zhì)點(diǎn)，采用四階Adams-Bashforth-Moulton(ABM)預(yù)測(cè)-校正法對(duì)下一時(shí)間步的波面抬高與自由面速度勢(shì)進(jìn)行更新.同時(shí)，采用人工阻尼層(damping layer)避免水池尾端壁面的反射，采用斜坡函數(shù)(ramp function)避免造波入口處水流速度變化過快造成的計(jì)算不穩(wěn)定，采用非均勻布點(diǎn)的Chebyshev五點(diǎn)光順法[23]對(duì)自由面上的速度勢(shì)以及波面進(jìn)行光順，避免鋸齒狀不穩(wěn)定現(xiàn)象(saw-tooth instability).將文中構(gòu)建的數(shù)值波浪水池的數(shù)值計(jì)算波形與二階Stokes波解析解進(jìn)行對(duì)比，依次探討了源點(diǎn)數(shù)量、去奇異距離以及時(shí)間步長(zhǎng)對(duì)計(jì)算精度的影響，并綜合考慮計(jì)算精度與計(jì)算效率，給出了計(jì)算參數(shù)建議值.

1 數(shù)學(xué)模型

二維數(shù)值波浪水池示意圖見圖1，建立笛卡爾坐標(biāo)系oxz，坐標(biāo)原點(diǎn)o位于靜水面與造波入口面的交點(diǎn)處，ox水平向右，oz豎直向上，圖1中h為水池靜水深；V為流域；造波入口面、自由面、尾端壁面和水底壁面四個(gè)流域邊界分別為ΓU，ΓF，ΓD，ΓB；xu為人工阻尼層起始端點(diǎn)的x軸坐標(biāo)；Ldl為人工阻尼層的長(zhǎng)度.

圖1 二維數(shù)值波浪水池示意圖

假定流體為理想流體且其運(yùn)動(dòng)無旋，則可采用勢(shì)流理論來描述水池內(nèi)的水體運(yùn)動(dòng)，其邊界條件與初始條件為

2φ=0，inV

(1)

(2)

(3)

(4)

(5)

φ=0，ast≤0，onΓU&ΓF&ΓD&ΓB

(6)

η=0，ast≤0，onΓF

(7)

式中：φ(x,z,t)為速度勢(shì)；n=(ni,nk)為流域邊界外法線方向的單位向量；η(x,t)為自由面上的波面抬高；g為重力加速度；ρ為流體密度；pa為自由面上的大氣壓力，假定恒等于0；φU(x,z,t)為造波入口面的輸入速度勢(shì)；Rm(t)和vu(x)分別為斜坡函數(shù)與阻尼系數(shù).

(8)

(9)

式中：Tm為斜坡函數(shù)的作用時(shí)間，文中取值為2倍的入射波周期；ω為入射波角頻率；α為控制參數(shù)，用于控制阻尼層的阻尼強(qiáng)度，參照文獻(xiàn)[19]取1，阻尼層長(zhǎng)度Ldl至少為1倍波長(zhǎng)，參照文獻(xiàn)[17]取1倍入射波波長(zhǎng)；xu相應(yīng)取尾端壁面前端1倍波長(zhǎng)位置的x軸坐標(biāo).

2 數(shù)值方法

2.1 去奇異邊界元法

去奇異邊界元法將原本位于流域邊界上的源點(diǎn)移至流域外一定距離，即將積分表面Sd移到流域V外，以保證源點(diǎn)和配置點(diǎn)不重合，從而避免奇異性，示意圖見圖2.

圖2 源點(diǎn)與配置點(diǎn)示意圖

去奇異邊界元法分為直接法與間接法，直接法直接運(yùn)用格林第二定理.本文采用間接法，流域內(nèi)任意點(diǎn)的速度勢(shì)為積分表面上源點(diǎn)對(duì)該點(diǎn)引起的速度勢(shì)的線性疊加，表示為Rankine源形式.

(10)

式中：q=(xq,zq)為流域中的任意點(diǎn)(包括流域邊界上的配置點(diǎn)與流域內(nèi)的場(chǎng)點(diǎn))坐標(biāo)向量；p=(xp,zp)為流域外的源點(diǎn)坐標(biāo)向量；σ(p)為對(duì)應(yīng)源點(diǎn)的源強(qiáng)；G采用簡(jiǎn)單格林函數(shù)，對(duì)于該二維問題，G(p,q)可定義為

G(p,q)=ln|p-q|

(11)

對(duì)于式(8)中的源強(qiáng)，需運(yùn)用已知的邊界條件來求解.應(yīng)用相關(guān)邊界條件，用來求解未知源強(qiáng)度的去奇異化邊界積分方程為

(12)

(13)

配置點(diǎn)與源點(diǎn)一一對(duì)應(yīng)，兩個(gè)對(duì)應(yīng)點(diǎn)之間的距離稱為去奇異距離Ld，為

Ld=ld(Dm)β

(14)

式中：ld和β為常數(shù)；Dm為局部網(wǎng)格尺寸，一般取網(wǎng)格大小的平方根，以源點(diǎn)代替網(wǎng)格，網(wǎng)格大小取自由面處源點(diǎn)之間的距離乘以造波入口面處源點(diǎn)之間的距離.其中l(wèi)d和β對(duì)計(jì)算精度影響很大，后文中將對(duì)這兩個(gè)系數(shù)對(duì)計(jì)算精度的影響進(jìn)行研究.

2.2 混合歐拉-拉格朗日法

為了追蹤流體微粒在瞬態(tài)自由面上隨時(shí)間的變化，將歐拉形式下的自由面運(yùn)動(dòng)學(xué)邊界條件式(3)及動(dòng)力學(xué)邊界條件式(4)轉(zhuǎn)換為拉格朗日形式，其轉(zhuǎn)換方法為采用隨體導(dǎo)數(shù).

(15)

式中：v為瞬態(tài)自由面上流體微粒的運(yùn)動(dòng)速度.

將式(13)代入自由表面邊界條件式(3)與式(4)，自由面邊界條件可改寫為如下形式.

(16)

(17)

(18)

(19)

2.3 四階Adams-Bashforth-Moulton預(yù)測(cè)-校正法

(20)

(21)

式中：Δt為時(shí)間步長(zhǎng).

然后將四階Adams-Bashforth方法中得到的預(yù)測(cè)值代入隱式四階Adams-Moulton公式中進(jìn)行校正，最終得到η與φ的校正值.

(22)

(23)

四階Adams-Bashforth-Moulton預(yù)測(cè)-校正法還需要前四個(gè)時(shí)間步的波面抬高與自由面速度勢(shì)才能自啟動(dòng)，其中第一個(gè)時(shí)間步的波面抬高與自由面速度勢(shì)已由初始條件式(6)與式(7)給出，其他三個(gè)時(shí)間步的波面抬高與自由面速度勢(shì)的計(jì)算表達(dá)式為

ηt+Δt=Δt·ft

(24)

φt+Δt=Δt·gt

(25)

3 計(jì)算參數(shù)對(duì)數(shù)值計(jì)算結(jié)果的影響

建立二維數(shù)值波浪水池，水池靜水深h為0.6 m.目標(biāo)波形的參數(shù)分別為波長(zhǎng)L=1 m，波高H=0.03 m，周期T=0.801 1 s.

取數(shù)值波浪水池中x1=L與x2=4L兩個(gè)監(jiān)測(cè)點(diǎn)，通過將兩個(gè)監(jiān)測(cè)點(diǎn)的穩(wěn)定波形與二階Stokes波解析解進(jìn)行比較，來分析計(jì)算參數(shù)對(duì)數(shù)值計(jì)算結(jié)果的影響，穩(wěn)定波形取15T～18T時(shí)間段.討論的計(jì)算參數(shù)包括源點(diǎn)數(shù)量、去奇異距離和時(shí)間步長(zhǎng)，其中源點(diǎn)數(shù)量又包括自由面單位波長(zhǎng)距離的源點(diǎn)數(shù)量nx和造波入口面的源點(diǎn)數(shù)量nz，去奇異距離又包括常系數(shù)β和ld，時(shí)間步長(zhǎng)用Δt表示.計(jì)算參數(shù)的初始取值參照文獻(xiàn)[17]，取nx=60，nz=30，β=0.5，ld=0.5，Δt=T/100.

一般來說，判別水池中生成波形的精度，應(yīng)從波高、波長(zhǎng)和周期三個(gè)方面進(jìn)行綜合評(píng)估，然而經(jīng)過預(yù)計(jì)算，在本文的計(jì)算參數(shù)取值范圍內(nèi)，所得波形的平均波長(zhǎng)和平均周期均與解析解相近，誤差小于1%，因此后文僅討論波高在不同計(jì)算參數(shù)下的精度.為更為直觀地觀察波高與解析解之間的誤差，定義x1監(jiān)測(cè)點(diǎn)處穩(wěn)定波形的平均波高為H1，與解析解波高H之間的相對(duì)誤差δ1為

(26)

定義x2監(jiān)測(cè)點(diǎn)處穩(wěn)定波形的平均波高為H2，與解析解波高H之間的相對(duì)誤差δ2為

(27)

為計(jì)量波高沿傳播方向的變化，定義x1監(jiān)測(cè)點(diǎn)與x2監(jiān)測(cè)點(diǎn)之間穩(wěn)定波形的平均波高相對(duì)變化δ3為

(28)

由式(12)可知，去奇異距離同時(shí)受源點(diǎn)數(shù)量nx、nz、常系數(shù)β和常系數(shù)ld的影響，后文先討論源點(diǎn)數(shù)量對(duì)計(jì)算精度的影響，確定了合適的源點(diǎn)數(shù)量后再討論兩個(gè)常系數(shù)的影響.

3.1 源點(diǎn)數(shù)量的影響

改變自由面單位波長(zhǎng)距離的源點(diǎn)數(shù)量nx，分別取nx=20，40，60和80.圖3為不同nx時(shí)x1與x2監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線，表1為nx對(duì)波高的影響.由圖3和表1可知，當(dāng)源點(diǎn)數(shù)量nx在40～80之間時(shí)，x1與x2監(jiān)測(cè)點(diǎn)的波高相對(duì)誤差δ1與δ2均在1%左右，波高相對(duì)變化δ3均在1%以下，nx為20時(shí)，x1監(jiān)測(cè)點(diǎn)的波高相對(duì)誤差δ1小于1%，但波高沿傳播方向的損失較嚴(yán)重，波高相對(duì)變化δ3達(dá)到了11%.選取源點(diǎn)數(shù)量nx=40.

圖3 x1與x2監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線

表1 nx對(duì)波高的影響

取上文中建議的源點(diǎn)數(shù)量nx=40，改變?cè)觳ㄈ肟诿娴脑袋c(diǎn)數(shù)量.由于兩個(gè)邊界交點(diǎn)處的邊界條件無法確定，本文在邊界交點(diǎn)處均未布置源點(diǎn)，同時(shí)為了確保自由面附近的計(jì)算精度，在接近自由面與造波入口面交點(diǎn)的位置z=-H位置在造波入口面額外增加了一個(gè)源點(diǎn)，因此造波入口面源點(diǎn)的布點(diǎn)方案為在z軸坐標(biāo)(-h,-H]的區(qū)間內(nèi)均勻布點(diǎn)，數(shù)量分別取nz=2，3，4，6，8和10.由于x1與x2監(jiān)測(cè)點(diǎn)的波形幾乎相同，因此圖4為不同nz時(shí)x1監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線，表2為nz對(duì)波高的影響.由圖4和表2可知，z方向的源點(diǎn)數(shù)量nz并非越大波形精度就越高，當(dāng)源點(diǎn)數(shù)量nz在2～10時(shí)，數(shù)值水池生成波浪的波高相對(duì)誤差δ1與δ2均小于2%，波高相對(duì)變化δ3小于1%.選取源點(diǎn)數(shù)量nz=3.

圖4 x1監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線

表2 nz對(duì)波高的影響

鑒于上文中造波入口僅需布置少量的源點(diǎn)，就能得到足夠精度的結(jié)果，對(duì)造波入口處源點(diǎn)的z軸位置對(duì)波形精度的影響進(jìn)行補(bǔ)充研究.僅在入口布置單個(gè)源點(diǎn)，分別取z=-H，z=-(h+H)/2和z=-h，所得的x1監(jiān)測(cè)點(diǎn)波面時(shí)歷曲線見圖5.由圖5可知，自由面附近，即z=-H處的源點(diǎn)對(duì)計(jì)算結(jié)果的精度影響很大，其單個(gè)源點(diǎn)所生成的波形已與解析解相差較小，而與自由面距離較遠(yuǎn)的源點(diǎn)所生成的波形與解析解相差很大，上文中入口源點(diǎn)均包含了z=-H處的源點(diǎn)，這也是造波入口僅布置少量源點(diǎn)就能獲得足夠波形精度的主要原因.

圖5 x1監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線

3.2 去奇異距離的影響

取上文中建議的源點(diǎn)數(shù)量nx=40，源點(diǎn)數(shù)量nz=3，改變參數(shù)β，分別取β=0.4，0.5，0.6，0.7，0.8，0.9和1.0.由于x1與x2監(jiān)測(cè)點(diǎn)的波形幾乎相同，圖6為不同β時(shí)x1監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線，表3為β對(duì)波高的影響.由圖6和表3可知，當(dāng)參數(shù)β在0.4～1.0時(shí)，波高相對(duì)變化δ3均小于1%，且參數(shù)β與波高相對(duì)誤差δ1與δ2的關(guān)系較復(fù)雜，當(dāng)β=0.5和0.8時(shí)，波高相對(duì)誤差δ1與δ2最小，當(dāng)β=0.4，0.7，0.9和1.0時(shí)，波高相對(duì)誤差δ1與δ2大于2%.當(dāng)β<0.4時(shí)，源點(diǎn)與配置點(diǎn)過近，邊界積分方程接近奇異，計(jì)算結(jié)果誤差過大或計(jì)算無法進(jìn)行，因此并未給出β<0.4時(shí)的計(jì)算結(jié)果，綜上選取參數(shù)β=0.5.

圖6 x1監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線

表3 β對(duì)波高的影響

取上文中建議的源點(diǎn)數(shù)量nx=40，源點(diǎn)數(shù)量nz=3，參數(shù)β=0.5，改變參數(shù)ld，分別取ld=0.4，0.5，0.6，0.7，0.8，0.9和1.0.由于x1與x2監(jiān)測(cè)點(diǎn)的波形幾乎相同，圖7為不同ld時(shí)x1監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線，表4為ld對(duì)波高的影響.由圖7和表4可知，當(dāng)參數(shù)ld在0.4～1.0時(shí)，波高相對(duì)變化δ3均小于1%，且參數(shù)ld與波高相對(duì)誤差δ1與δ2之間的關(guān)系較復(fù)雜，當(dāng)ld=0.5和1.0時(shí)，波高相對(duì)誤差δ1與δ2最小，當(dāng)β=0.7，0.8和0.9時(shí)，波高相對(duì)誤差δ1與δ2大于2%.當(dāng)ld=0.3時(shí)，波高相對(duì)誤差δ1與δ2在4%左右，ld<0.3時(shí)，源點(diǎn)與配置點(diǎn)過近，邊界積分方程接近奇異，計(jì)算結(jié)果誤差過大或計(jì)算無法進(jìn)行，因此并未給出ld<0.4時(shí)的計(jì)算結(jié)果，綜上選取參數(shù)ld=0.5.

圖7 x1監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線

表4 ld對(duì)波高的影響

3.3 時(shí)間步長(zhǎng)的影響

取上文中建議的源點(diǎn)數(shù)量nx=40，源點(diǎn)數(shù)量nz=3，參數(shù)β=0.5，參數(shù)ld=0.5，改變時(shí)間步長(zhǎng)Δt，分別取Δt=T/50，T/75，T/100，T/125和T/150.由于x1與x2監(jiān)測(cè)點(diǎn)波形幾乎相同，圖8為不同Δt時(shí)x1監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線，表5為Δt對(duì)波高的影響.由圖8和表5可知，當(dāng)時(shí)間步長(zhǎng)Δt在T/150～T/50時(shí)，波高相對(duì)誤差δ1與δ2以及波高相對(duì)變化δ3均小于1%.時(shí)間步長(zhǎng)越大，計(jì)算效率越高，綜上選取時(shí)間步長(zhǎng)Δt=T/50.

圖8 x1監(jiān)測(cè)點(diǎn)的波面時(shí)歷曲線

表5 Δt對(duì)波高的影響

4 結(jié) 論

1) 構(gòu)建的數(shù)值波浪水池對(duì)二維弱非線性計(jì)算有足夠的精度，為后續(xù)研究水面式防波堤的水動(dòng)力性能問題奠定基礎(chǔ).

2) 造波入口面的源點(diǎn)數(shù)量nz取3.造波入口處的源點(diǎn)中，距離自由面較近的源點(diǎn)對(duì)計(jì)算精度影響很大，在造波入口設(shè)置源點(diǎn)時(shí)應(yīng)至少在自由面附近布置一個(gè)源點(diǎn)，否則將產(chǎn)生較大計(jì)算誤差.

3) 源點(diǎn)數(shù)量中，自由面單位波長(zhǎng)距離的源點(diǎn)數(shù)量nx取40，時(shí)間步長(zhǎng)Δt取T/50，有足夠的計(jì)算精度同時(shí)兼顧了計(jì)算效率，其對(duì)應(yīng)的CFL數(shù)為0.8，實(shí)際計(jì)算中為保持計(jì)算穩(wěn)定，應(yīng)使CFL數(shù)≤0.8.

4) 去奇異距離中，建議參數(shù)β取0.5，參數(shù)ld取0.5.參數(shù)β與ld對(duì)計(jì)算精度影響很大，且兩個(gè)參數(shù)與計(jì)算精度之間的關(guān)系比較復(fù)雜，在計(jì)算中應(yīng)慎重取值.當(dāng)β為0.5與0.8，ld為0.5與1.0時(shí)，計(jì)算精度最高