波浪模式的變分多尺度方法模擬

劉小剛，任睿超，王震，惠小健，章培軍

(1.西京學院 理學院，陜西 西安 710123；2.西北大學 數學學院，陜西 西安 710127)

1 引言

波浪模式是流體力學中的一個研究熱點[1-5]，水波自由晃動和波浪傳播是典型的波浪模式問題，它的性質和規律很大程度上反映了自由面問題的性質和規律，所以對其數值方法進行深入研究具有重大意義。用標準的有限元方法模擬自由面問題時，當對流項或反應項占優時，數值結果會伴隨著劇烈的數值偽振蕩。這種無任何物理意義的數值偽振蕩可以使數值結果嚴重失真[6]。變分多尺度方法通過對“細”尺度解的近似自動獲得穩定化結構的穩定化因子，并把“細”尺度上的近似解代入原偏微分方程對應的Galerkin 變分弱形式以修正“粗”尺度上的解，使得“粗”尺度上的解包含“細”尺度上解的特性，從而提高了數值求解過程的數值穩定性[7-8]，并基于變分多尺度方法，研究波浪模式問題。

2 波浪模式的模型描述

波浪模式的模型方程是一個非定常黏性不可壓Navier-Stokes 方程，即

式中，v=(u,v)T為流體的速度向量；p為壓力；f=(f1,f2)T為體積力；ρ為密度；μ是運動黏度系數；t為時間；σ為應力張量；Ω為二維求解區域；Γ1為區域 Ω的自由面邊界；Γ2為區域 Ω的固壁邊界；Γ=Γ1∪Γ2為求解區域的邊界；n是邊界Γ的單位外法向量。公式中ρ=1 000 kg/m3，μ=1.0×10－2m2/s，重力加速度f2=9.8 m/s2，f1=0.0 m/s2。其中，速度v=(u,v)T、應力張量 σ和法向量n為矢量，其他變量為標量。

3 數值算法

空間離散采用正規網格剖分T={Ωe|,e=1,2,···,nel}，Ω表示計算區域，nel為網格數，h=max{hΩe}，hΩe是網格Ωe的直徑。試驗函數空間取為

相應的檢驗函數空間取為

式中，PK(Ωe)是網格 Ωe上的K階完全多項式。壓力空間的試驗函數和檢驗函數空間取為

式（1）和式（2）的變分弱形式為

式中，q是壓力的試驗函數和檢驗函數；w、v分別為速度的檢驗函數和試驗函數。

式（2）是一個非線性方程組[9]，對其進行線性化可得

式中，v0為Picard 循環迭代中的前一步數值解。

下面根據變分多尺度的思想建立求解方程（9）的“細”尺度模型。首先把檢驗函數和試驗函數分解到“粗” “細”兩種尺度上，即

取如下的空間B={b|b∈(HK/PK)(Ωe),b=0,on?Ωe;e=1,2,···,nel}，并要求w′(x),v′(x)的每個分量均屬于B，即用高階“泡”近似“細”尺度上的解[10]。

把式（12）代入式（8）可得

式（12）中的細尺度v′(x,t)可用在時間域上間斷、空間域上連續的高階分段多項式表述。特別地，細尺度v′(x,t)可用時間域上為常數的分段函數表示。因此令于是把式（11）式（12）代入式（9），并分解[10]，可得

式（14）和式（15）就是動量方程在“粗” “細”兩種尺度上滿足的控制方程。

下面通過Petrov-Glerkin 方法求解“細”尺度方程（15）。

為簡單起見，“細”尺度上的試驗函數和檢驗函數均取一個基函數[10]，分別為則將其代入式（15），由 β的任意性可得“細”尺度上的近似解為

式（16）至式（18）即為基于原偏微分方程殘差的“細”尺度上的近似解。

再建立總體變分多尺度方程。把式（16）至式（18）代入式（14）中，由散度定理，可得

將式（16）代入式（13），并運用散度定理，有

方程（19）和（20）就是要求解的整體變分多尺度方程。

這里采用修正迭代技術確定自由面的位置。自由面上具有如下3 個邊界條件：

式中，u,v分別為x軸（、y軸）上的速度；nx,ny為垂直于自由面的單位法線，γx,γy是垂直于自由表面力；S為表面張力系數；(ρ1,ρ2)是曲面的主曲率半徑。定義h=h(x,t)為t時刻x處的自由面高度函數，則微分形式下的自由面運動條件為

式中，U=(u,v)。選取一階Euler 格式

作為自由面運動條件式（24）的更新格式。

4 數值模擬

4.1 水波自由晃蕩

計算中采用的時間步長為Δt=0.01，這里計算了一個周期T內自由面的變化趨勢，計算區域為1 m×3 m，網格步長為0.1 m，上邊界為自由面邊界，其他邊界為固壁邊界，初始波形為η=cos(kx)，從初始波形到恢復初始波形的晃動為1 個周期T。圖1 至圖5 分別給出了初始時刻和T時刻的速度及壓力。

前面在初始時刻指定的初始速度為0，而圖1a 和圖1b 卻給出了非零的速度，這是由于在保持初始自由面位置不變的條件下迭代求解了15 次控制方程，使之達到穩態。

從圖2 可以看出，在水體重力的作用下，微幅波左側水面下降，右側水面在左側水體的推動和擠壓下上升，并在時刻達到水平狀態。由于動力慣性的作用，左側水面繼續下降而右側水面繼續上升，經過約周期的運動后呈現圖3 所示的狀態。在這一時刻水的速度方向發生改變，開始了新一輪半個周期的運動，如圖4 和圖5 所示，其分析結果同前所述。由圖1至圖5 可以看出，水的運動速度呈現周期性變化的規律，壓力隨著水深的增加而逐漸增大。

4.2 波浪的傳播

波浪的傳播是海浪模式的一個主要研究內容。圖6 至圖8 給出了波浪傳播的數值模擬結果，其中每個圖分別給出了波浪的水平速度、垂直速度和壓力。在初始時刻，波浪位于計算區域的中間（圖6），隨后波浪向左傳播。由圖7 可見，在波浪向左傳播的過程中，波峰向左移動，波峰兩邊的垂直速度方向相反。圖8 顯示了波峰傳播到左邊界的模式。這個算例可以推廣到研究海浪與海岸相互作用實際問題中。

5 結論

圖1 初始時刻波浪的水平速度（u）、垂直速度（v）和壓力Fig.1 The horizontal velocity (u),the vertical velocity (v),and the pressure of the wave at the initial moment

圖2 T/4時刻波浪的水平速度（u）、垂直速度（v）和壓力Fig.2 The horizontal velocity (u),the vertical velocity (v),and the pressure of the wave at the T moment

圖3 T/2時刻波浪的水平速度（u）、垂直速度（v）和壓力Fig.3 The horizontal velocity (u),the vertical velocity (v),and the pressure of the wave at the T/2 moment

圖4 3T/4時刻波浪的水平速度（u）、垂直速度（v）和壓力Fig.4 The horizontal velocity (u),the vertical velocity (v),and the pressure of the wave at the 3T/4 moment

圖5 T時刻波浪的水平速度（ u ）、垂直速度（v）和壓力Fig.5 The horizontal velocity (u),the vertical velocity (v),and the pressure of the wave at the T moment

圖6 初始時刻波浪的水平速度（u）、垂直速度（v）和壓力Fig.6 The horizontal velocity (u),the vertical velocity (v),and the pressure of the wave at the initial moment

圖7 t=0.8 s 時刻波浪的水平速度（u）、垂直速度（v）和壓力Fig.7 The horizontal velocity (u),the vertical velocity (v),and the pressure of the wave at t=0.8 s moment

圖8 t=1.6 s 時刻波浪的水平速度（u）、垂直速度（v）和壓力Fig.8 The horizontal velocity (u),the vertical velocity (v),and the pressure of the wave at t=1.6 s moment

變分多尺度方法可以有效地處理具有多尺度效應的數學物理問題。本文在變分多尺度的理論框架內將變分多尺度方法與自由面技術相合，提出了模擬波浪模式的一種新的數值方法，其中自由表面位置的確定由自由面運動條件控制。作為算例，應用此技術數值求解了水波自由晃動問題和波浪傳播問題，得到了如下結論：

（1）通過Petrov-Galerkin 方法，用“泡”函數近似“細”尺度上的解，變分多尺度方法可以消除由對流項占優和速度-壓力失耦引起的數值偽振蕩，并得到滿意的數值結果。

（2）對水波自由晃動問題，結合自由面運動條件的變分多尺度方法能很好地模擬自由面的周期性變化趨勢，而且速度、壓力的分布也呈現周期性變化的規律。

（3）變分多尺度方法能夠有效地模擬波浪傳播問題，可以進一步推廣到研究海浪與海岸相互作用的實際工程問題。