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

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

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

1 引言

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

2 波浪模式的模型描述

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

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

3 數(shù)值算法

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

相應(yīng)的檢驗函數(shù)空間取為

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

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

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

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

式中，v0為Picard 循環(huán)迭代中的前一步數(shù)值解。

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

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

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

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

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

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

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

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

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

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

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

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

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

4 數(shù)值模擬

4.1 水波自由晃蕩

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

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

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

4.2 波浪的傳播

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

5 結(jié)論

圖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

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

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

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

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