基于動(dòng)網(wǎng)格方法的孤立波數(shù)值造波研究

賈進(jìn)生

(山西大同大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，山西大同 037009)

在海洋中，地震、海嘯引起的波浪通常是孤立波，其波長(zhǎng)較長(zhǎng)且波高較大，對(duì)海岸附近的人類生產(chǎn)生活和建筑物安全有著重大影響。因此，針對(duì)孤立波開展研究，有著重要的實(shí)際意義。

目前，對(duì)于波浪的常用研究手段包括物理模型實(shí)驗(yàn)和數(shù)值模擬。物理模型實(shí)驗(yàn)是在實(shí)驗(yàn)室中按照一定的比尺對(duì)實(shí)際工程進(jìn)行還原，用造波機(jī)產(chǎn)生目標(biāo)波浪，通過實(shí)驗(yàn)研究波浪的相關(guān)問題。然而，物理模型實(shí)驗(yàn)對(duì)實(shí)驗(yàn)條件要求較高，并且實(shí)驗(yàn)成本較貴。數(shù)值模擬的成本相對(duì)較低，同時(shí)較物理模型實(shí)驗(yàn)更為靈活，因此近年來已成為相關(guān)領(lǐng)域的研究熱點(diǎn)。

在實(shí)驗(yàn)室中，孤立波通常由造波機(jī)用推板造波[1-4]的方法產(chǎn)生，而數(shù)值模型中通常采用速度入口法或質(zhì)量源方法造波[2-4]。基于動(dòng)網(wǎng)格方法，采用類似推板造波的方式在數(shù)值模型中產(chǎn)生孤立波，并對(duì)網(wǎng)格尺寸對(duì)模擬結(jié)果的影響進(jìn)行研究。

1 數(shù)值模型

流體的基本控制方程包括質(zhì)量守恒方程和動(dòng)量守恒方程。在研究波浪問題時(shí)，通常將水視為不可壓縮的流體，即水的密度為常數(shù)，此時(shí)質(zhì)量守恒方程可以簡(jiǎn)化為：

其中，U是水體的速度矢量，該方程通常稱為連續(xù)性方程。

動(dòng)量方程表示為：

其中Pd是動(dòng)水壓力，g表示重力加速度，X是位置矢量。

在數(shù)值模型中，水體和氣體的交界面采用VOF方法進(jìn)行計(jì)算，其控制方程為：

其中，α 為體積分?jǐn)?shù)，代表水體占網(wǎng)格體積的百分?jǐn)?shù)，當(dāng)α=0 時(shí)代表氣體，α=1 時(shí)代表水體，當(dāng)網(wǎng)格在界面上時(shí)α介于0和1之間。

在空間上，對(duì)整個(gè)計(jì)算域用有限體積法進(jìn)行離散。為了求解上述控制方程，將方程中的每一項(xiàng)進(jìn)行離散，轉(zhuǎn)化為代數(shù)方程，從而采用PISO 算法耦合求解速度和壓力。用庫(kù)朗數(shù)來限制時(shí)間步長(zhǎng)，從而保證數(shù)值計(jì)算的精度，庫(kù)朗數(shù)的表達(dá)式為：

屈哨兵：這個(gè)問題問得好。這也是我最近一直在思考的一個(gè)問題，就是好教育的質(zhì)量觀。這個(gè)問題我此前還沒有做很多思考，最近稍微梳理了一下，大概有幾個(gè)基本觀點(diǎn)可以先說一下。

其中，Δt表示時(shí)間步長(zhǎng)，δx是網(wǎng)格的長(zhǎng)度，計(jì)算中的最大庫(kù)朗數(shù)設(shè)為0.5，時(shí)間步長(zhǎng)在計(jì)算中由于庫(kù)朗數(shù)的限制不斷變化。

2 動(dòng)網(wǎng)格造波

采用動(dòng)網(wǎng)格技術(shù)，可以在數(shù)值計(jì)算過程中，實(shí)現(xiàn)網(wǎng)格大小和位置的不斷變化。對(duì)于造波邊界處的網(wǎng)格，可模擬造波機(jī)的往復(fù)運(yùn)動(dòng)，從而用類似于推板造波的方式產(chǎn)生目標(biāo)波浪。

對(duì)于推板造波，在造波邊界上應(yīng)當(dāng)滿足條件：

其中，u代表水質(zhì)點(diǎn)的水平速度，在淺水條件下，由連續(xù)性方程可得水質(zhì)點(diǎn)水平速度的關(guān)系如下：

其中，C表示波速，對(duì)于孤立波，有：

其中，H為孤立4波的波高，d為水深，將式(6)、(7)代入(5)，積分可得造波邊界的位移公式：

造波邊界的沖程為S，造波邊界運(yùn)動(dòng)時(shí)間為T，

在數(shù)值模型中，當(dāng)時(shí)間t＜T時(shí)，在每一時(shí)間步開始時(shí)迭代求解式(9)，求得的X(t)即當(dāng)前時(shí)刻造波邊界的位置，當(dāng)造波結(jié)束后，造波邊界在X(t)=S處?kù)o止不動(dòng)。

3 數(shù)值計(jì)算驗(yàn)證

基于上述數(shù)值模型，建立數(shù)值波浪水槽進(jìn)行驗(yàn)證，水槽長(zhǎng)40 m，高1.7 m，如圖1 所示，整個(gè)水槽均劃分為長(zhǎng)方形網(wǎng)格進(jìn)行計(jì)算，為了保證數(shù)值計(jì)算精度，在水面附近對(duì)網(wǎng)格加密。水槽長(zhǎng)度方向（x方向）的網(wǎng)格尺寸為Δx=0.25 m，在水面附近的網(wǎng)格加密區(qū)內(nèi)，水槽高度方向（y方向）的網(wǎng)格尺寸為Δy=0.025 m。造波邊界位于水槽左側(cè)，模擬的孤立波波高為0.25 m，水深1.2 m。經(jīng)計(jì)算，在水槽中x=5 m，10 m和15 m 處的波面歷時(shí)曲線如圖2 所示，其中x是與造波邊界的距離。

圖1 數(shù)值水槽示意圖

由圖2 可以看出，在距離造波邊界5 m 處，模擬的孤立波波高與目標(biāo)波高基本一致，說明動(dòng)網(wǎng)格造波方法是有效的。但在孤立波繼續(xù)傳播的過程中，波高有所衰減，尤其是到達(dá)距離造波邊界15 m 處時(shí)，衰減較為明顯。

圖2 距造波邊界5 m，10 m和15 m處的波面歷時(shí)曲線

將x方向上的網(wǎng)格加密，網(wǎng)格尺寸減小為Δx=0.08 m，y方向的網(wǎng)格仍為Δy=0.025 m，再次進(jìn)行計(jì)算。此外，為了探討y方向網(wǎng)格尺寸的影響，取網(wǎng)格尺寸Δx=0.08 m，Δy=0.065 m作為對(duì)照。由圖2中的結(jié)果可以看出，在x=5 m和10 m處，由于距離造波邊界較近，其沿程衰減并不明顯，因此僅對(duì)比x=15 m處的計(jì)算結(jié)果。如圖3所示，紅色和藍(lán)色虛線分別為網(wǎng)格尺寸Δx=0.08 m，Δy=0.025 m和Δx=0.08 m，Δy=0.065 m時(shí)，x=15 m 處的波面歷時(shí)曲線，黑色線是該位置處的理論波面，用公式(7)計(jì)算。

圖3 距造波邊界5 m，10 m和15 m處的波面歷時(shí)曲線

可以看出，當(dāng)Δx=0.08 m，Δy=0.025 m 時(shí)，計(jì)算結(jié)果與理論波面幾乎重合。與圖2 對(duì)比可以看出，x方向的網(wǎng)格尺寸對(duì)孤立波的沿程衰減程度有一定影響，當(dāng)x方向的網(wǎng)格足夠密時(shí)，沿程衰減幾乎可以忽略不計(jì)。通過與Δx=0.08 m，Δy=0.065 m 時(shí)的計(jì)算結(jié)果對(duì)比可以看出，當(dāng)y方向的網(wǎng)格較疏時(shí)，不僅會(huì)造成孤立波的沿程衰減，還會(huì)造成孤立波相位的差異，這可能是由于此時(shí)水面處的速度計(jì)算存在誤差，因此造成了孤立波的相位差。

需要指出的是，模擬孤立波時(shí)，網(wǎng)格尺寸的選擇還與波浪參數(shù)有關(guān)。當(dāng)波高較小時(shí)，網(wǎng)格尺寸也需要相應(yīng)減小。例如當(dāng)水深0.55 m 時(shí)，模擬波高為0.04 m 的孤立波，應(yīng)當(dāng)選取尺寸較小的網(wǎng)格，Δx=0.05 m，Δy=0.008 m，計(jì)算結(jié)果如圖4 所示?？梢钥闯?，模擬的孤立波波高與目標(biāo)波高一致，并且在傳播至距離造波邊界20 m 處時(shí)仍沒有明顯的衰減。此外，評(píng)判孤立波的造波效果時(shí)通常會(huì)觀察是否存在尾波，即孤立波傳播過后，是否存在后續(xù)的波動(dòng)。在以上多個(gè)算例中均沒有明顯的尾波，說明孤立波的數(shù)值模擬方法是行之有效的。

圖4 距造波邊界10 m，15 m和20 m處的波面歷時(shí)曲線

4 結(jié)語

基于動(dòng)網(wǎng)格方法，用模擬推板造波的方式在數(shù)值模型中產(chǎn)生了孤立波。在數(shù)值計(jì)算中，沿水槽長(zhǎng)度和高度方向的網(wǎng)格較疏時(shí)，會(huì)造成孤立波波高的沿程衰減，而水槽高度方向的網(wǎng)格較疏還會(huì)造成孤立波的相位變化。整體而言，產(chǎn)生的孤立波與目標(biāo)波高一致，并且沒有明顯的尾波，說明文中的方法是有效的。