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

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

1 引言

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

2 波浪模式的模型描述

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

式中，v=(u,v)T為流體的速度向量；p為壓力；f=(f1,f2)T為體積力；ρ為密度；μ是運(yùn)動(dòng)黏度系數(shù)；t為時(shí)間；σ為應(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}，Ω表示計(jì)算區(qū)域，nel為網(wǎng)格數(shù)，h=max{hΩe}，hΩe是網(wǎng)格Ωe的直徑。試驗(yàn)函數(shù)空間取為

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

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

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

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

式（2）是一個(gè)非線性方程組[9]，對(duì)其進(jìn)行線性化可得

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

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

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

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

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

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

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

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

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

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

將式（16）代入式（13），并運(yùn)用散度定理，有

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

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

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

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

作為自由面運(yùn)動(dòng)條件式（24）的更新格式。

4 數(shù)值模擬

4.1 水波自由晃蕩

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

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

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

4.2 波浪的傳播

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

5 結(jié)論

圖1 初始時(shí)刻波浪的水平速度（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時(shí)刻波浪的水平速度（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時(shí)刻波浪的水平速度（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時(shí)刻波浪的水平速度（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時(shí)刻波浪的水平速度（ u ）、垂直速度（v）和壓力Fig.5 The horizontal velocity (u),the vertical velocity (v),and the pressure of the wave at the T moment

圖6 初始時(shí)刻波浪的水平速度（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 時(shí)刻波浪的水平速度（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 時(shí)刻波浪的水平速度（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ùn)動(dòng)條件控制。作為算例，應(yīng)用此技術(shù)數(shù)值求解了水波自由晃動(dòng)問題和波浪傳播問題，得到了如下結(jié)論：

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

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

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