賈進(jìn)生

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

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

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

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

1 數(shù)值模型

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

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

動量方程表示為：

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

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

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

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

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

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

2 動網(wǎng)格造波

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

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

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

其中，C表示波速，對于孤立波，有：

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

造波邊界的沖程為S，造波邊界運動時間為T，

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

3 數(shù)值計算驗證

基于上述數(shù)值模型，建立數(shù)值波浪水槽進(jìn)行驗證，水槽長40 m，高1.7 m，如圖1 所示，整個水槽均劃分為長方形網(wǎng)格進(jìn)行計算，為了保證數(shù)值計算精度，在水面附近對網(wǎ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)計算，在水槽中x=5 m，10 m和15 m 處的波面歷時曲線如圖2 所示，其中x是與造波邊界的距離。

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

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

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

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

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

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

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

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

4 結(jié)語

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