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

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

0 引 言

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

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

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

基于勢流理論的數(shù)值計算方法，發(fā)展較早的是采用Green函數(shù)在流場邊界布置源匯分布的邊界元方法.其按照是否考慮流動隨時間的變化，分為頻域和時域兩大類方法.相對于頻域方法對弱非線性、周期性問題有效的局限性而言，時域方法應用范圍更加廣泛.按照求解的積分方程是否具有奇異性，又出現(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多年來，相繼被應用于求解與波浪運動相關(guān)的問題[16-19].

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

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

1 數(shù)學模型

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

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

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

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

(8)

(9)

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

2 數(shù)值方法

2.1 去奇異邊界元法

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

圖2 源點與配置點示意圖

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

(10)

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

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

(11)

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

(12)

(13)

配置點與源點一一對應，兩個對應點之間的距離稱為去奇異距離Ld，為

Ld=ld(Dm)β

(14)

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

2.2 混合歐拉-拉格朗日法

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

(15)

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

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

(16)

(17)

(18)

(19)

2.3 四階Adams-Bashforth-Moulton預測-校正法

(20)

(21)

式中：Δt為時間步長.

然后將四階Adams-Bashforth方法中得到的預測值代入隱式四階Adams-Moulton公式中進行校正，最終得到η與φ的校正值.

(22)

(23)

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

ηt+Δt=Δt·ft

(24)

φt+Δt=Δt·gt

(25)

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

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

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

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

(26)

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

(27)

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

(28)

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

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

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

圖3 x1與x2監(jiān)測點的波面時歷曲線

表1 nx對波高的影響

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

圖4 x1監(jiān)測點的波面時歷曲線

表2 nz對波高的影響

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

圖5 x1監(jiān)測點的波面時歷曲線

3.2 去奇異距離的影響

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

圖6 x1監(jiān)測點的波面時歷曲線

表3 β對波高的影響

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

圖7 x1監(jiān)測點的波面時歷曲線

表4 ld對波高的影響

3.3 時間步長的影響

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

圖8 x1監(jiān)測點的波面時歷曲線

表5 Δt對波高的影響

4 結(jié) 論

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

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

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

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