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

程友良，許 強(qiáng)，焦慎俐

（ 1.華北電力大學(xué)動(dòng)力工程系，河北 保定 071003；2.華北電力大學(xué) 河北省低碳高效發(fā)電技術(shù)重點(diǎn)實(shí)驗(yàn)室，河北 保定 071003； 3.華北電力大學(xué) 保定市低碳高效發(fā)電技術(shù)重點(diǎn)實(shí)驗(yàn)室，河北 保定 071003）

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

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

1 波浪要素

1.1 求解波浪要素

設(shè)周期T=1.94 s，靜水深d=1 m，波高H=0.1 m。由彌散關(guān)系對(duì)波長(zhǎng)進(jìn)行迭代求解，利用MATLAB的求解程序求得波長(zhǎng)L=5 m。

1.2 驗(yàn)證波浪要素

1.2.1 二階stokes波的理論驗(yàn)證

根據(jù)梅沃特圖[17]采用兩個(gè)無因次獨(dú)立參數(shù)。

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

1.2.2 模擬條件驗(yàn)證

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

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

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

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

（3）波浪破碎極限

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

式中，H0為深水波波高；L0為深水波波長(zhǎng)。

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

2 速度邊界造波法

2.1 二階Stokes波

流體運(yùn)動(dòng)時(shí)，必須遵循的另一規(guī)律為動(dòng)量守恒定理。在研究液體波浪運(yùn)動(dòng)時(shí)，一般認(rèn)為流體是無粘性的理想流體，此時(shí)可以忽略流體的粘性效應(yīng)，只考慮壓強(qiáng)而不必考慮切向應(yīng)力。

Stokes二階波的速度勢(shì)函數(shù)見式（7）。

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

Stokes波的水質(zhì)點(diǎn)運(yùn)動(dòng)速度見式（9）和式（10）。

式中，u、v分別為x、y方向的水質(zhì)點(diǎn)速度；為波頻。

2.2 控制方程

連續(xù)性方程見式（11）。

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

式中，ρ為液體密度（kg/m3）；g為重力加速（m/s）；t為時(shí)間（s）；p為壓強(qiáng)（Pa）；μ為動(dòng)力粘度（Pa·s）。

2.3 VOF模型

2.3.1 VOF機(jī)理

VOF模型中界面跟蹤是通過求解相連續(xù)方程實(shí)現(xiàn)的，界面的位置是通過求出體積分量中急劇變化的點(diǎn)來確定的。VOF模型可以對(duì)兩種或更多種不相溶的流體進(jìn)行建模。典型的應(yīng)用包括射流破裂的預(yù)測(cè)，液體中大氣泡的運(yùn)動(dòng)，潰壩后液體的運(yùn)動(dòng)以及任何液—?dú)饨缑娴姆€(wěn)定或瞬態(tài)跟蹤。

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

（1）αq=0：?jiǎn)卧锊淮嬖诘趒相流體。

（2）αq=1：?jiǎn)卧锍錆M了第q相流體。

（3）0＜αq＜1：?jiǎn)卧锇说趒相流體和一相或者其他多相流體的界面。

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

2.3.2 VOF方程

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

式中，S為質(zhì)量源項(xiàng)；V為速度矢量。在默認(rèn)情況下該方程右端源項(xiàng)為0，但當(dāng)你給每一項(xiàng)制定常數(shù)或用戶定義的質(zhì)量源，則右端不為0。初相的體積分?jǐn)?shù)的計(jì)算基于公式（16）。

2.3.3 VOF使用限制

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

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

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

（4）流向周期運(yùn)動(dòng)（指定質(zhì)量流率或指定壓力降）無法與VOF一起使用。

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

（6）當(dāng)利用并行計(jì)算進(jìn)行粒子追蹤時(shí)，在共享內(nèi)存選項(xiàng)被激活的情況下，DPM模型無法與VOF模型一起使用[19]。

2.4 消波

2.4.1 阻尼消波

在消波區(qū)的動(dòng)量方程添加阻尼源項(xiàng)。

式中，μ（x）為消波系數(shù)，根據(jù)公式（19）確定。

式中，a為常數(shù)，根據(jù)經(jīng)驗(yàn)選取；x0為消波區(qū)起始位置。

2.4.2 拉伸網(wǎng)格

在消波區(qū)增大網(wǎng)格尺寸，通過計(jì)算過程中的數(shù)值耗散來達(dá)到消波的目的。

2.4.3 混合消波

阻尼消波與拉伸網(wǎng)格方法相結(jié)合的混合網(wǎng)格消波方法。本文中所用消波方法都為此法。

2.5 模型建立

數(shù)值波浪水池如圖1所示，長(zhǎng)40 m，高2 m，靜水深為1 m，20～40 m處為消波區(qū)。

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

2.5.1 網(wǎng)格劃分

網(wǎng)格劃分為表1的3種方式，網(wǎng)格劃分方法參考了相關(guān)文獻(xiàn)[20-21]。Δx為x方向單個(gè)單元格長(zhǎng)度；加密處Δy為靜水面±H區(qū)域內(nèi)y方向單個(gè)單元格長(zhǎng)度，x方向不變；非加密區(qū)Δy為加密區(qū)之外區(qū)域y方向單個(gè)單元格長(zhǎng)度，x方向不變。此外，消波區(qū)網(wǎng)格在x方向數(shù)量取150，最小為Δx，并以1.1的網(wǎng)格尺寸比率增加。

表1 網(wǎng)格劃分方式

網(wǎng)格的全局示意圖如圖2所示。為了降低壓力傳遞過程中的數(shù)值耗散，需要對(duì)波面附近的網(wǎng)格進(jìn)行加密，加密區(qū)（圖3）。

圖2 全局網(wǎng)格

圖3 局部網(wǎng)格

2.5.2 模擬設(shè)置

該模型采用表2所示的設(shè)置。

表2 模型設(shè)置

在消波區(qū)添加編譯好的阻尼源項(xiàng)，和拉伸后的網(wǎng)格一起達(dá)到消波的目的。文獻(xiàn)[18]中提到使用層流模型與常用的模擬結(jié)果沒有顯著差異，故選用層流模型。在距進(jìn)口5 m、10 m和15 m的波面處設(shè)置浪高儀，監(jiān)測(cè)波幅的變化。

2.5.3 獨(dú)立性驗(yàn)證

（1）網(wǎng)格獨(dú)立性驗(yàn)證

在進(jìn)行數(shù)值計(jì)算過程中，為了降低網(wǎng)格密度對(duì)計(jì)算結(jié)果的影響，通常要對(duì)多套疏密程度不同的網(wǎng)格系統(tǒng)進(jìn)行計(jì)算，并比較不同網(wǎng)格系統(tǒng)下的計(jì)算結(jié)果。

從圖4可得出，3種網(wǎng)格模擬的波幅隨時(shí)間變化的趨勢(shì)相同，隨著模擬時(shí)間的增加，模擬結(jié)果與理論波幅的誤差逐漸增大。這是因?yàn)樵谙▍^(qū)的混合消波方法雖然起到很好的消波效果，但并不能完全消除反射。故隨著模擬時(shí)間的增加，波浪的反射會(huì)越來越大，影響模擬的正確性。因此為了確保結(jié)果的準(zhǔn)確性，需要控制模擬波浪的時(shí)間。經(jīng)過多次模擬，該數(shù)值波浪水池模型下的模擬時(shí)間控制在30～50 s。由圖5可得波幅波峰和波谷的變化趨勢(shì)。

圖4 波幅歷時(shí)曲線

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

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

圖5 波峰與波谷的局部對(duì)比圖

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

為進(jìn)一步驗(yàn)證，引入波幅的相對(duì)誤差wc。

式中，y為模擬波幅下的縱坐標(biāo)；y理論為理論波幅下的縱坐標(biāo)。

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

圖6 相對(duì)誤差

根據(jù)以上結(jié)論和計(jì)算資源考慮，選取Mesh2的網(wǎng)格劃分方式。

（2）網(wǎng)格時(shí)間獨(dú)立性驗(yàn)證

對(duì)于瞬態(tài)計(jì)算，其結(jié)果不僅依賴于計(jì)算網(wǎng)格，還依賴于時(shí)間步長(zhǎng)。在進(jìn)行瞬態(tài)計(jì)算過程中，時(shí)間步長(zhǎng)的選擇固然是按照所需的時(shí)間分辨率來進(jìn)行，但依照CFL條件，瞬態(tài)計(jì)算對(duì)于時(shí)間步長(zhǎng)也存在明確的要求。

時(shí)間獨(dú)立性驗(yàn)證方法與網(wǎng)格獨(dú)立性方法類似，取不同大小的時(shí)間步長(zhǎng)進(jìn)行瞬態(tài)計(jì)算，考察相同時(shí)間點(diǎn)上的計(jì)算結(jié)果，分析其誤差，最終確定最大時(shí)間步長(zhǎng)用于計(jì)算[19]。

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

圖7 x=5 m處波幅歷時(shí)曲線

2.5.4 模擬結(jié)果

為了驗(yàn)證模擬結(jié)果的正確性，需要將模擬結(jié)果進(jìn)行多方面論證。

在計(jì)算前監(jiān)測(cè)距進(jìn)口5 m、10 m和15 m處波幅的歷時(shí)曲線（圖8），并與理論歷時(shí)曲線進(jìn)行對(duì)比。波幅隨位置的變化并未發(fā)生波長(zhǎng)方向的偏移，波幅方向與理論波幅吻合良好。故從歷時(shí)曲線的結(jié)果可以驗(yàn)證模擬結(jié)果波幅的正確性。

圖8 5 m、10 m和15 m處波幅歷時(shí)曲線

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

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

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

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

圖11 10～50 s壓力分布圖

一個(gè)周期內(nèi)x方向水質(zhì)點(diǎn)的運(yùn)動(dòng)速度u（圖12）：波峰處最大；在波幅與靜水面的交界處為零；在交界處運(yùn)動(dòng)到波谷處速度增至最大，方向?yàn)樨?fù)；在波谷運(yùn)動(dòng)到下一個(gè)交界處速度再次降為0；最終回到波峰時(shí)速度再次回到正向最大，以此循環(huán)往復(fù)。

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

一個(gè)周期內(nèi)y方向水質(zhì)點(diǎn)的運(yùn)動(dòng)速度v（圖13）：波峰處為零；在波幅與靜水面的交界處達(dá)到正向最大；在交界處運(yùn)動(dòng)到波谷處降為0；在波谷運(yùn)動(dòng)到下一個(gè)交界處速度增至最大，方向?yàn)樨?fù)；最終回到波峰時(shí)速度再次降為0，以此循環(huán)往復(fù)。

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

以上模擬結(jié)果的分析皆符合理論上一個(gè)周期內(nèi)水質(zhì)點(diǎn)運(yùn)動(dòng)的速度情況。

經(jīng)過以上論證，該模型的精確性良好。

3 推板造波

如圖14所示，數(shù)值波浪水池與速度邊界模型相同，但進(jìn)口、出口邊界條件及消波區(qū)的設(shè)置不同。

圖14 推板造波水池模型

3.1 推板造波機(jī)理

推板造波方法模仿的是實(shí)驗(yàn)中的推板式造波機(jī)，通過動(dòng)邊界實(shí)現(xiàn)。控制方程如下：

式中，xp為推板位移；S為推板運(yùn)動(dòng)的沖程。

3.2 模擬設(shè)置

在模擬過程中發(fā)現(xiàn)推板造波方法的一些問題，對(duì)設(shè)置進(jìn)行了一些調(diào)整。

（1）進(jìn)口邊界采用動(dòng)網(wǎng)格layering方法。為預(yù)防計(jì)算中負(fù)體積的計(jì)算，需要根據(jù)分裂因子α及合并因子對(duì)初始高度h0進(jìn)行設(shè)置，即可在邊界運(yùn)動(dòng)中控制相鄰邊界的網(wǎng)格的分裂及合并[22]。

（2）阻尼項(xiàng)的適用范圍從4L變?yōu)長(zhǎng)，網(wǎng)格拉伸區(qū)域仍為4L。因?yàn)樵谀M過程中發(fā)現(xiàn)適用邊界速度法的設(shè)置會(huì)發(fā)生庫(kù)朗數(shù)過大的錯(cuò)誤。庫(kù)朗數(shù)表達(dá)式如式（25）所示[19]：

在相同網(wǎng)格數(shù)量時(shí)減小時(shí)間步長(zhǎng)，不改變阻尼項(xiàng)適用長(zhǎng)度，在25 s左右發(fā)現(xiàn)計(jì)算不會(huì)收斂，經(jīng)過調(diào)整松弛因子發(fā)現(xiàn)波幅會(huì)突然增大。最終將阻尼項(xiàng)適用長(zhǎng)度改為L(zhǎng)且網(wǎng)格拉伸區(qū)長(zhǎng)度不變，得到的結(jié)果符合預(yù)期且收斂。

（3）時(shí)間步長(zhǎng)的改變。原時(shí)間步長(zhǎng)在推板造波方法上會(huì)發(fā)生庫(kù)朗數(shù)過大，故降低時(shí)間步長(zhǎng)，并對(duì)其進(jìn)行時(shí)間獨(dú)立性驗(yàn)證。

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

圖15 x=5 m波幅歷時(shí)曲線

3.3 模擬結(jié)果

同樣，為了驗(yàn)證模擬結(jié)果的正確性，需要將模擬結(jié)果進(jìn)行多方面論證。

同速度邊界造波法一樣，模擬結(jié)果波幅的歷時(shí)曲線與波幅的理論歷時(shí)曲線進(jìn)行對(duì)比（圖16）。波幅隨位置的變化并未發(fā)生波長(zhǎng)方向的偏移，波幅方向與理論波幅吻合良好。故從歷時(shí)曲線的結(jié)果可以驗(yàn)證模擬結(jié)果波幅的正確性。

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

如圖17和18所示，波浪在10 s左右到達(dá)20 m處的消波區(qū)，且可以發(fā)現(xiàn)只有網(wǎng)格拉伸方法區(qū)域的消波效果并不理想，只靠拉伸網(wǎng)格方法達(dá)到消波的作用消波區(qū)的長(zhǎng)度要數(shù)以倍計(jì)，故起到主要消波作用的是添加阻尼項(xiàng)的消波段。37.5～40 m處在20～50 s期間阻尼項(xiàng)消波區(qū)一直保持水平狀態(tài)，消波效果較好。

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

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

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

圖19 10～50 s壓力分布圖

一個(gè)周期內(nèi)x方向水質(zhì)點(diǎn)的運(yùn)動(dòng)速度u如圖20所示。波面以下與速度邊界造波法相同，但波面以上氣相的u在波浪運(yùn)動(dòng)到消波區(qū)后部分區(qū)域速度會(huì)有較大升高，隨計(jì)算時(shí)間增加，影響到其他部分。但未影響波面以下。

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

一個(gè)周期內(nèi)y方向水質(zhì)點(diǎn)的運(yùn)動(dòng)速度v如圖21所示，波面以下與速度邊界造波法，波面以上與速度u的發(fā)生的變化相同，氣相不穩(wěn)定，隨計(jì)算時(shí)間增加影響會(huì)增大。

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

以上模擬結(jié)果的分析皆符合理論上一個(gè)周期內(nèi)水質(zhì)點(diǎn)運(yùn)動(dòng)的速度情況。

4 兩種造波法方法的比較

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

（1）如圖20和圖21所示，因?yàn)槌隹谶吔缭O(shè)置的不同，氣相速度分布會(huì)有明顯的差別：推板造波方法在波面以下的速度分布圖與速度邊界造波方法的速度分布圖幾乎相同，但推板造法方法在波面上方氣體部分的速度分布圖更加無序，速度量級(jí)也較高，但對(duì)波面下影響不大。

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

圖22 x=5 m、10 m和15 m處波幅歷時(shí)曲線的對(duì)比

圖23 10～50 s波幅歷程曲線的對(duì)比

（3）相同網(wǎng)格下：推板造波對(duì)時(shí)間步長(zhǎng)更加敏感，需要更小的時(shí)間步長(zhǎng)才能達(dá)到收斂要求；推板造波采用動(dòng)網(wǎng)格方法，對(duì)網(wǎng)格的要求比速度邊界方法更高。

（4）適當(dāng)網(wǎng)格拉伸對(duì)波幅的影響極小，故混合消波方式有更好的發(fā)展空間。通過兩種消波方式區(qū)域長(zhǎng)度的不同組合，不僅可以改善收斂，還可以在計(jì)算資源受限和水池長(zhǎng)度較長(zhǎng)時(shí)得到應(yīng)用。

5 結(jié) 論

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

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

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

（3）在計(jì)算資源受限或需要加長(zhǎng)數(shù)值水槽長(zhǎng)度時(shí)，可將主要模擬區(qū)域之后的網(wǎng)格逐漸拉伸，在消波區(qū)末端加入阻尼項(xiàng)。可兼顧計(jì)算資源和精度兩方面。