速度邊界造波法與推板造波法的差異化分析

程友良，許 強，焦慎俐

（ 1.華北電力大學動力工程系，河北 保定 071003；2.華北電力大學 河北省低碳高效發電技術重點實驗室，河北 保定 071003； 3.華北電力大學 保定市低碳高效發電技術重點實驗室，河北 保定 071003）

數值波浪水池是用于模擬自由表面波、水動力和浮體運動的數值工具的總稱[1]。在海洋工程領域，物理模型和數值模擬是兩種主要的實驗方法。物理模型可以在實驗室中準確的再現波的傳播和變形過程，但這種方法成本較高，而且受到實驗設備設置的限制。隨著計算機技術的飛速發展，數值模擬已成為研究波浪及其對海洋結構物作用的有效工具。數值模擬的成本更低，且能夠得到更多難以通過物理實驗測量的細節[2]。BOO S Y[3]將線性規則波和不規則波的仿真結果與理論值進行了比較，線性波仿真結果與理論值吻合良好，但對非線性波的模擬要進一步研究。RYU S等[4]為研究非線性波浪開發了一種具有完全非線性自由表面邊界條件的數值波浪水池。此后，數值波浪水池的優化一直伴隨著數值波浪水池的發展。CHEN Y P 等[5]提出了一種基于變換的N-S方程求解的全非線性數值波浪水池，數值計算結果與線性或非線性解析解吻合較好，驗證了該數值波浪水池的準確性。ZHANG X T等[6]研究了在完全非線性數值波浪水池中的波傳播問題，解決了每個時間步長的潛在流動邊界值問題。隨著數值波浪水池可用性以及功能性快速增長，其精確性越來越重要。BUNNIK T等[7]對基于勢流理論和VOF（Volume of Fluid）理論的數值波浪水池準確性進行評估。LIU Z等[8]建立了基于Fluent的不可壓縮粘性流體兩相VOF數值波浪水池，以模擬波浪的產生和傳播，且波浪運動誘發了氣相的運動。LAL A等[9]使用ANSYS CFX模擬了襟翼類型造波機，并對各種湍流模型進行了比較，模擬采用層流模型、k-ε模型和SST（Shear-Stress Transport）切應力輸運模型，結果顯示這3種方法沒有顯著差異。GOMES M N等[10]采用了基于VOF模型的二維數值波浪水池，模擬了兩種不同形式的規則重力波，表明VOF是一種能夠正確模擬水與空氣之間相互作用的模型。無網格方法也適用于數值波浪水池。SENTURK U等[11]研究了在自由表面上受到非線性邊界條件的二維數值波浪水池數值模擬，采用局部無網格RBF（Radial Basis Function）方法。這項研究表明了無網格方法的靈活性，特別是對于移動邊界問題，其中求解過程中的每個時間步僅需要更新節點位置而不是重新整理整個域。隨著計算機能力的進步，波浪模擬精確性的提高，好的消波方式能進一步提高數值波浪水池的性能。KIM M W等[12]評估了用于漸進水波的各種人工阻尼方案。將五種類型的阻尼方案應用于計算域的末端，并比較了它們的阻尼能力，提出了在一定范圍的阻尼區域上與適當的阻尼形狀相關聯的最優數值阻尼方案。但反射通常是無法避免的，PERIC R等[13]提出了一種通過數值流動模擬來量化這種影響，使數值波浪水池的波反射最小化。3年后，PERIC R等[14]進一步提出了一種解析預測反射系數的理論，該理論可用于在進行仿真之前優化選擇源項參數，與規則自由表面波的有限體積流動模擬結果進行了對比驗證，結果表明，該方法具有令人滿意的實用精度。隨著造波與消波方法的優化發展，造波的其他方面也引起人們重視。GALERA-CALERO L等[15]對采用VOF的k-ε、k-ω和LES（Large Eddy Simulation）模型分別進行了分析和比較，在物理波浪水池上進行了實測，并與理論進行了比較，驗證了各種湍流模型的正確性。UDDIN M N等[16]提出了一種二維數值波浪水池來計算實驗波浪水槽下壁的靜壓變化，得出對于小型二維實驗波浪水槽，非定常無粘流體法可以較準確地預測表面壓力分布。

數值波浪水池模型有多種，但國內對不同種類造波方法的討論極少，本文對速度邊界造波法和推板造波法的波幅歷時曲線、歷程曲線、壓力分布圖和u、v速度分布圖進行了分析，得出的結論對海洋工程具有重要意義。

1 波浪要素

1.1 求解波浪要素

設周期T=1.94 s，靜水深d=1 m，波高H=0.1 m。由彌散關系對波長進行迭代求解，利用MATLAB的求解程序求得波長L=5 m。

1.2 驗證波浪要素

1.2.1 二階stokes波的理論驗證

根據梅沃特圖[17]采用兩個無因次獨立參數。

該波浪要素符合理論要求。

1.2.2 模擬條件驗證

（1）二階stokes波只適用于中、深水深的規則波。對于淺水波，；對于中、深水波，。該波浪要素d/L=1/5＞1/20，符合條件。

（2）波函數方程中的二次項必須要小于線性項，從而確保收斂，沒有波浪破碎。要滿足以上條件必須要滿足式（4）和式（5）[18]。

式中，A為振幅；k為波數；Ur為努塞爾數。

該波浪要素Ur=0.005＜1/3，A=0.05＜0.26，符合條件。

（3）波浪破碎極限

當波浪破碎時，沒有任何波浪理論是有效的。最有效的波動特征值在水深較大時為波陡H/L，在水深較小時為相對波高H/d。在深水中波浪的破碎極限如下。

式中，H0為深水波波高；L0為深水波波長。

該波浪要素未超過破碎極限，符合Stokes二階波。

2 速度邊界造波法

2.1 二階Stokes波

流體運動時，必須遵循的另一規律為動量守恒定理。在研究液體波浪運動時，一般認為流體是無粘性的理想流體，此時可以忽略流體的粘性效應，只考慮壓強而不必考慮切向應力。

Stokes二階波的速度勢函數見式（7）。

Stokes二階波的波面方程見式（8）。

Stokes波的水質點運動速度見式（9）和式（10）。

式中，u、v分別為x、y方向的水質點速度；為波頻。

2.2 控制方程

連續性方程見式（11）。

Navier-Stokes方程（N-S方程）見式（12）至式（14）。

式中，ρ為液體密度（kg/m3）；g為重力加速（m/s）；t為時間（s）；p為壓強（Pa）；μ為動力粘度（Pa·s）。

2.3 VOF模型

2.3.1 VOF機理

VOF模型中界面跟蹤是通過求解相連續方程實現的，界面的位置是通過求出體積分量中急劇變化的點來確定的。VOF模型可以對兩種或更多種不相溶的流體進行建模。典型的應用包括射流破裂的預測，液體中大氣泡的運動，潰壩后液體的運動以及任何液—氣界面的穩定或瞬態跟蹤。

VOF模型中，不同的流體組份共用著一套動量方程，通過引進相體積分數這一變量，實現對每一個計算單元相界面的追蹤。在每個控制體積內，所有相體積分數總和為1。因此，在任何確定單元中的變量和屬性要么純粹代表其中一個相，要么代表各相混合，這取決于體積分數值。故在單元中，若第q相流體體積分數為αq，那么可能存在以下三種情況：

（1）αq=0：單元里不存在第q相流體。

（2）αq=1：單元里充滿了第q相流體。

（3）0＜αq＜1：單元里包含了第q相流體和一相或者其他多相流體的界面。

基于αq的值，為計算域中每個計算單元分配合適的屬性與變量。

2.3.2 VOF方程

跟蹤相之間的界面是通過求解一相或者多相的容積比率的連續方程來完成的。對第q相，有

式中，S為質量源項；V為速度矢量。在默認情況下該方程右端源項為0，但當你給每一項制定常數或用戶定義的質量源，則右端不為0。初相的體積分數的計算基于公式（16）。

2.3.3 VOF使用限制

（1）VOF只能應用于壓力求解器，在密度基求解器中無法使用。

（2）每一計算單元必須充滿一種或多種流體，VOF模型不允許區域中不存在任何流體的情況。

（3）僅有一相可定義為可壓縮理想氣體，但對于使用UDF定義的可壓縮液體則無限制。

（4）流向周期運動（指定質量流率或指定壓力降）無法與VOF一起使用。

（5）二階隱式時間步格式無法與VOF顯示格式一起使用。

（6）當利用并行計算進行粒子追蹤時，在共享內存選項被激活的情況下，DPM模型無法與VOF模型一起使用[19]。

2.4 消波

2.4.1 阻尼消波

在消波區的動量方程添加阻尼源項。

式中，μ（x）為消波系數，根據公式（19）確定。

式中，a為常數，根據經驗選??；x0為消波區起始位置。

2.4.2 拉伸網格

在消波區增大網格尺寸，通過計算過程中的數值耗散來達到消波的目的。

2.4.3 混合消波

阻尼消波與拉伸網格方法相結合的混合網格消波方法。本文中所用消波方法都為此法。

2.5 模型建立

數值波浪水池如圖1所示，長40 m，高2 m，靜水深為1 m，20～40 m處為消波區。

圖1 速度邊界造波水池模型

2.5.1 網格劃分

網格劃分為表1的3種方式，網格劃分方法參考了相關文獻[20-21]。Δx為x方向單個單元格長度；加密處Δy為靜水面±H區域內y方向單個單元格長度，x方向不變；非加密區Δy為加密區之外區域y方向單個單元格長度，x方向不變。此外，消波區網格在x方向數量取150，最小為Δx，并以1.1的網格尺寸比率增加。

表1 網格劃分方式

網格的全局示意圖如圖2所示。為了降低壓力傳遞過程中的數值耗散，需要對波面附近的網格進行加密，加密區（圖3）。

圖2 全局網格

圖3 局部網格

2.5.2 模擬設置

該模型采用表2所示的設置。

表2 模型設置

在消波區添加編譯好的阻尼源項，和拉伸后的網格一起達到消波的目的。文獻[18]中提到使用層流模型與常用的模擬結果沒有顯著差異，故選用層流模型。在距進口5 m、10 m和15 m的波面處設置浪高儀，監測波幅的變化。

2.5.3 獨立性驗證

（1）網格獨立性驗證

在進行數值計算過程中，為了降低網格密度對計算結果的影響，通常要對多套疏密程度不同的網格系統進行計算，并比較不同網格系統下的計算結果。

從圖4可得出，3種網格模擬的波幅隨時間變化的趨勢相同，隨著模擬時間的增加，模擬結果與理論波幅的誤差逐漸增大。這是因為在消波區的混合消波方法雖然起到很好的消波效果，但并不能完全消除反射。故隨著模擬時間的增加，波浪的反射會越來越大，影響模擬的正確性。因此為了確保結果的準確性，需要控制模擬波浪的時間。經過多次模擬，該數值波浪水池模型下的模擬時間控制在30～50 s。由圖5可得波幅波峰和波谷的變化趨勢。

圖4 波幅歷時曲線

圖5（a）中3種網格下模擬的波幅與理論波幅在波峰相一致。在波谷：Mesh3與理論波幅近乎一致；Mesh2與Mesh3在10～16 s相一致，但在16～20 s會與Mesh3拉開距離。相比其他兩種網格，Mesh1下波谷處一直與理論結果相差較大。

圖5（b）中3種網格下模擬的波幅與理論波幅較符合。波峰：Mesh1與Mesh2下的波峰的波幅與理論波幅的波峰相一致；Mesh3下波峰的波幅比理論波幅的波峰要高。波谷：Mesh3下的波谷的波幅與理論波幅幾乎一致，Mesh2與Mesh1下波谷的波幅比理論波幅的波谷高，且Mesh1下最高。

圖5 波峰與波谷的局部對比圖

故Mesh1的結果與理論波幅相差較大，Mesh2與Mesh3較為接近。但Mesh2在波峰表現更好，Mesh3在波谷表現更佳。

為進一步驗證，引入波幅的相對誤差wc。

式中，y為模擬波幅下的縱坐標；y理論為理論波幅下的縱坐標。

由圖6可以看出：1）在波幅穩定后的35 s內，相對誤差符合wc1＞wc2＞wc3；2）在35 s后，Mesh1的相對誤差最大，Mesh2和Mesh3的相對誤差相似；3）由圖中的最大相對誤差值小于0.6%可驗證該模型造波情況良好。

圖6 相對誤差

根據以上結論和計算資源考慮，選取Mesh2的網格劃分方式。

（2）網格時間獨立性驗證

對于瞬態計算，其結果不僅依賴于計算網格，還依賴于時間步長。在進行瞬態計算過程中，時間步長的選擇固然是按照所需的時間分辨率來進行，但依照CFL條件，瞬態計算對于時間步長也存在明確的要求。

時間獨立性驗證方法與網格獨立性方法類似，取不同大小的時間步長進行瞬態計算，考察相同時間點上的計算結果，分析其誤差，最終確定最大時間步長用于計算[19]。

如圖7中局部放大圖所示，在波峰處，時間步長t=0.004 s下波幅較理論波幅下有較大的增大，t=0.01 s與t=0.02 s下波幅較吻合，且與理論波幅相比增大較小。在波谷處，3種時間步長下波幅與理論波長的差距：0.02＞0.01＞0.004，但t=0.01 s與t=0.004 s差距接近，都接近理論波長的波谷。故將對比情況與計算資源結合考慮，最終選擇Mesh2網格，時間步長為0.01 s，計算時間為50 s。

圖7 x=5 m處波幅歷時曲線

2.5.4 模擬結果

為了驗證模擬結果的正確性，需要將模擬結果進行多方面論證。

在計算前監測距進口5 m、10 m和15 m處波幅的歷時曲線（圖8），并與理論歷時曲線進行對比。波幅隨位置的變化并未發生波長方向的偏移，波幅方向與理論波幅吻合良好。故從歷時曲線的結果可以驗證模擬結果波幅的正確性。

圖8 5 m、10 m和15 m處波幅歷時曲線

消波區加入阻尼項的阻尼系數過小會導致消波不完全，影響波浪模擬的準確性，而通過加長消波區域的長度，則會增大計算機資源消耗。但阻尼系數過大，消波效果仿佛在20 m處放置了一堵墻，因此造成的反射同樣會影響波浪模擬的正確性。如圖9和圖10所示，波浪在10 s左右到達20 m處的消波區，且在20～35 m之間消波效果良好,在35～40 m趨于平穩。整個區域的波浪情況在20 s時基本穩定，20～50 s消波區域保持穩定，且在整個模擬時間段內消波區末端水平面一直保持水平狀態。故驗證了波浪模擬的準確性和消波區良好的消波效果。

圖9 5～25 s波幅歷程曲線

圖10 30～50 s波幅歷程曲線

如圖11所示，壓力分布圖隨時間的變化過程。在波浪傳遞過程中，壓力隨著波浪水池深度的增加而增大，且壓力數值與理論計算值相符。圖中虛線為該時刻理論波幅的狀態，通過與壓力分布圖對比得出，靜壓等值線與理論波幅相符。且靜壓強必然是沿著作用面的內法線法向，即可驗證流動狀態的波形與理論波形相符合。

圖11 10～50 s壓力分布圖

一個周期內x方向水質點的運動速度u（圖12）：波峰處最大；在波幅與靜水面的交界處為零；在交界處運動到波谷處速度增至最大，方向為負；在波谷運動到下一個交界處速度再次降為0；最終回到波峰時速度再次回到正向最大，以此循環往復。

圖12 x方向10～50 s 速度分布圖

一個周期內y方向水質點的運動速度v（圖13）：波峰處為零；在波幅與靜水面的交界處達到正向最大；在交界處運動到波谷處降為0；在波谷運動到下一個交界處速度增至最大，方向為負；最終回到波峰時速度再次降為0，以此循環往復。

圖13 y方向10～50 s速度分布圖

以上模擬結果的分析皆符合理論上一個周期內水質點運動的速度情況。

經過以上論證，該模型的精確性良好。

3 推板造波

如圖14所示，數值波浪水池與速度邊界模型相同，但進口、出口邊界條件及消波區的設置不同。

圖14 推板造波水池模型

3.1 推板造波機理

推板造波方法模仿的是實驗中的推板式造波機，通過動邊界實現?？刂品匠倘缦拢?/p>

式中，xp為推板位移；S為推板運動的沖程。

3.2 模擬設置

在模擬過程中發現推板造波方法的一些問題，對設置進行了一些調整。

（1）進口邊界采用動網格layering方法。為預防計算中負體積的計算，需要根據分裂因子α及合并因子對初始高度h0進行設置，即可在邊界運動中控制相鄰邊界的網格的分裂及合并[22]。

（2）阻尼項的適用范圍從4L變為L，網格拉伸區域仍為4L。因為在模擬過程中發現適用邊界速度法的設置會發生庫朗數過大的錯誤。庫朗數表達式如式（25）所示[19]：

在相同網格數量時減小時間步長，不改變阻尼項適用長度，在25 s左右發現計算不會收斂，經過調整松弛因子發現波幅會突然增大。最終將阻尼項適用長度改為L且網格拉伸區長度不變，得到的結果符合預期且收斂。

（3）時間步長的改變。原時間步長在推板造波方法上會發生庫朗數過大，故降低時間步長，并對其進行時間獨立性驗證。

如圖15所示，在波谷時y0.0025＞y＞y0.005＞y0.01（y為理論波幅），且y0.0025和y0.005與理論波幅相差較小，y0.01與理論波幅相差最大；在波峰時，y＞y0.0025＞y0.005＞y0.01，且y0.0025與y0.005幾乎一致，y0.01與理論波幅的相差最大；y0.01波幅雖然也較好，但在30 s左右會發生不收斂的情況。故選擇0.005 s為時間步長。

圖15 x=5 m波幅歷時曲線

3.3 模擬結果

同樣，為了驗證模擬結果的正確性，需要將模擬結果進行多方面論證。

同速度邊界造波法一樣，模擬結果波幅的歷時曲線與波幅的理論歷時曲線進行對比（圖16）。波幅隨位置的變化并未發生波長方向的偏移，波幅方向與理論波幅吻合良好。故從歷時曲線的結果可以驗證模擬結果波幅的正確性。

圖16 x=5 m、10 m和15 m處幅面歷時曲線

如圖17和18所示，波浪在10 s左右到達20 m處的消波區，且可以發現只有網格拉伸方法區域的消波效果并不理想，只靠拉伸網格方法達到消波的作用消波區的長度要數以倍計，故起到主要消波作用的是添加阻尼項的消波段。37.5～40 m處在20～50 s期間阻尼項消波區一直保持水平狀態，消波效果較好。

圖17 5～25 s波幅歷程曲線

圖18 30～50 s波幅歷程曲線

如圖19所示靜壓云圖隨時間（10 s、20 s、30 s、40 s、50 s）的變化過程。在波浪運動過程中，壓力隨波浪水池深度的增加而增大，且壓力的分層現象與壓力數值皆與理論相符。

圖19 10～50 s壓力分布圖

一個周期內x方向水質點的運動速度u如圖20所示。波面以下與速度邊界造波法相同，但波面以上氣相的u在波浪運動到消波區后部分區域速度會有較大升高，隨計算時間增加，影響到其他部分。但未影響波面以下。

圖20 x方向5～25 s速度分布圖

一個周期內y方向水質點的運動速度v如圖21所示，波面以下與速度邊界造波法，波面以上與速度u的發生的變化相同，氣相不穩定，隨計算時間增加影響會增大。

圖21 y方向5～25 s速度分布圖

以上模擬結果的分析皆符合理論上一個周期內水質點運動的速度情況。

4 兩種造波法方法的比較

基于兩種造波方法模擬結果的分析，可以得到以下結果：

（1）如圖20和圖21所示，因為出口邊界設置的不同，氣相速度分布會有明顯的差別：推板造波方法在波面以下的速度分布圖與速度邊界造波方法的速度分布圖幾乎相同，但推板造法方法在波面上方氣體部分的速度分布圖更加無序，速度量級也較高，但對波面下影響不大。

（2）如圖22所示，三個位置波幅的歷時曲線表明：推板造波方法從造波開始到波浪形成這個階段比速度邊界造波方法更加平穩，且圖23所示速度邊界方法的波幅50 s時在接近消波區在波長方向會有輕微的偏移。因此推板造波方法的波幅相比速度邊界方法更加貼合理論波幅，但邊界速度方法的波幅與理論波幅的相對誤差極小。

圖22 x=5 m、10 m和15 m處波幅歷時曲線的對比

圖23 10～50 s波幅歷程曲線的對比

（3）相同網格下：推板造波對時間步長更加敏感，需要更小的時間步長才能達到收斂要求；推板造波采用動網格方法，對網格的要求比速度邊界方法更高。

（4）適當網格拉伸對波幅的影響極小，故混合消波方式有更好的發展空間。通過兩種消波方式區域長度的不同組合，不僅可以改善收斂，還可以在計算資源受限和水池長度較長時得到應用。

5 結 論

從推板造波法結果可證明，純拉伸網格方式消波效果不佳，消波主力仍是阻尼消波。但兩種方法的組合可以在大幅降低網格數量的同時保持較好的計算精度，且促進了推板造波的收斂。兩種造波方法的波幅與理論波幅的吻合度都很好，需要根據情況選擇，具體如下。

（1）計算資源與時間充足的情況下，推板造波法比速度邊界造波法更好。但推板造波法對網格要求更高，需要更小的時間步長或更精密的網格劃分。

（2）推板造波法出口設置為wall，需要考慮長時間計算波面上方氣相速度變化帶來的影響，在考慮氣相速度變化的數值模擬中采用開放邊界的數值波浪水池更合適。

（3）在計算資源受限或需要加長數值水槽長度時，可將主要模擬區域之后的網格逐漸拉伸，在消波區末端加入阻尼項?？杉骖櫽嬎阗Y源和精度兩方面。