基于去奇異邊界元法的二維數值波浪水池計算參數影響分析

2021-12-24 01:43:12楊師宇吳靜萍陳昌哲

武漢理工大學學報(交通科學與工程版) 2021年5期

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

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

0 引言

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

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

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

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

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

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

1 數學模型

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

圖1 二維數值波浪水池示意圖

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

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)分別為斜坡函數與阻尼系數.

(8)

(9)

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

2 數值方法

2.1 去奇異邊界元法

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

圖2 源點與配置點示意圖

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

(10)

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

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

(11)

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

(12)

(13)

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

Ld=ld(Dm)β

(14)

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

2.2 混合歐拉-拉格朗日法

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

(15)

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

將式(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 計算參數對數值計算結果的影響

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

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

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

(26)

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

(27)

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

(28)

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

3.1 源點數量的影響

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

圖3 x1與x2監測點的波面時歷曲線

表1 nx對波高的影響

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

圖4 x1監測點的波面時歷曲線

表2 nz對波高的影響

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

圖5 x1監測點的波面時歷曲線

3.2 去奇異距離的影響

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

圖6 x1監測點的波面時歷曲線

表3 β對波高的影響

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

圖7 x1監測點的波面時歷曲線

表4 ld對波高的影響

3.3 時間步長的影響

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