流體-流體相互作用模型的無條件穩(wěn)定格式?

李偉，黃鵬展

(新疆大學數(shù)學與系統(tǒng)科學學院，新疆烏魯木齊 830017)

0 引言

在計算流體力學中,兩種互不相溶流體的多區(qū)域多物理場耦合數(shù)值模擬是許多科學和工業(yè)應用中的一個重要內(nèi)容, 例如環(huán)境工程中的大氣海洋相互作用問題、血流異質性的模擬等等. 這些實際問題可以由兩個不可壓縮牛頓流體的Navier-Stokes 方程通過摩擦型界面條件耦合來模擬, 這里要求界面滿足剛性耦合條件, 即沒有流體穿過交界面.

我們考慮? ?R2是一個有界的區(qū)域, 該區(qū)域有兩個子域?1和?2且兩個區(qū)域的界面為I. 設計算時間為t ∈(0,T], 求速度 ui:(0,T]×?i→R2和壓力 pi:(0,T]×?i→R,i=1,2,滿足如下方程 (t ∈(0,T])[1]

其中|·|表示歐幾里得范數(shù), ui,0∈H1(?i)2是初始速度. 這里的νi>0 表示運動粘性, κ>0 是摩擦系數(shù), 外力項fi:[0,T]→ H1(?i)2. 向量 ni是 ??i上的單位法線, τ 是 I 上的任意切向量, 使得 τ ·ni=0. 注意到該流體-流體相互作用模型在I 上采用了非線性界面條件.

由于模型(1)在實際應用中的重要性, 許多研究工作致力于構造高效的解耦方法[2?7]. 最近, 基于解耦時間步長方法, Connors 等[1]提出了兩種解耦時間步長方法,第一種是隱式/顯式格式,它是最簡單和最自然的解耦方法. 隨后, Zhang 等[8]證明了該格式是條件穩(wěn)定的. 另一種是幾何平均格式, 該格式是無條件穩(wěn)定[1]. 但是這些方法都不是以顯式格式處理界面項.我們知道顯格式雖然計算簡便,但是會產(chǎn)生穩(wěn)定化條件.故構造無條件穩(wěn)定的顯格式就顯得尤為重要. 因此, 本文將使用全顯式的方法來解耦非線性界面條件.

本文主要構建并研究了一種無條件穩(wěn)定的有限元方法來求解非線性流體-流體相互作用模型. 新方法具有以下優(yōu)點: (1) 與隱式/顯式格式和幾何平均格式不同, 非線性界面條件采用顯格式, 計算簡便易于實現(xiàn); (2) 新方法是無條件穩(wěn)定的. 首先介紹了所研究問題的一些符號、函數(shù)空間和一些重要結果. 其次提出了具有顯式格式的界面解耦方法, 描述了求解過程, 并推導了算法的無條件穩(wěn)定性. 最后設計了一些數(shù)值實驗來驗證該方法的穩(wěn)定性和有效性.

1 預備知識

在本節(jié)中, 將給出一些必要的定義和等式. 對于L2(?i)空間, 我們分別用‖·‖0和(·,·)來表示其范數(shù)及其內(nèi)積. 接著, 我們用H1(?i)表示Sobolev 空間W2,1(?i). 對于非線性流體-流體相互作用模型(1), 引入以下壓力和速度函數(shù)空間: 對于i=1,2, 有

接下來, 我們在 Xi×Xi和 Xi×Mi上定義兩個連續(xù)的雙線性形式 a(·,·)和 d(·,·)如下:

以及在 Xi×Xi×Xi上定義一個非線性項 b(·,·,·)[9]

在以上函數(shù)空間的定義下, 我們給出問題(1)的變分形式如下, 求ui:(0,T]→Xi和pi:(0,T]→Mi, 使得對于任意的 (vi,qi)∈Xi×Mi,i,j=1,2,i/=j 滿足

2 流體-流體相互作用模型的全離散有限元算法

2.1 修正的標量輔助變量技巧

標量輔助變量方法最早基于能量不變二次化方法提出, 應用于梯度流問題[10]. 通過引入標量輔助變量, 可以對原偏微分方程進行修改. 在時間離散系統(tǒng)中, 通過動量方程與輔助變量方程的求和可以消除修改后動量方程中的非線性項. 本文基于標量輔助變量方法思想, 構造一個修正的標量輔助變量. 這時, 非線性界面耦合條件可以在時間離散系統(tǒng)中抵消.

這時, 修正的標量輔助變量Q 對于t 的導數(shù)寫作

注意到, 對(4)中界面項求和等于零,

在連續(xù)情形下(1+t)Q=1, 故(1)式結合修改的界面條件(3)的解與原問題(1)的解相同.

對應的(1)與(3)組成系統(tǒng)的變分形式如下, 求ui:(0,T]→Xi和pi:(0,T]→Mi, 對于?(vi,qi)∈Xi×Mi,i,j=1,2, i/=j 滿足

2.2 基于修正標量輔助變量的全離散有限元方法

首先, 對于N>0, 令{tn}Nn=0是對[0,T]均勻劃分得到的時間剖分, 其中時間步長?t= NT并有 tn=n?t. 對于 i=1, 2, 令 πih是子區(qū)域?i上的三角形剖分并定義為πh=π1h∪π2h, 其h 是πh上三角形單元的最大外接圓直徑. 此外, 在 πih上速度和壓力的離散空間分別記為Xih?Xi和Mih?Mi. 其有限元離散子空間定義如下:

其中: 對于?Ki∈πih,P1(Ki)表示定義在Ki上的最高次數(shù)為1 次的多項式,?b ∈H01(Ki)為在單元Ki上重心取值為1 的函數(shù),并滿足0 ≤?b ≤1.眾所周知,該有限元空間配對Mih和 Xih滿足離散的 Ladyzenskaja-Babu?ska-Brezzi條件

其中: βi>0 為只依賴于區(qū)域?i的正常數(shù). 更進一步, 記(uni,h,pni,h)為模型(1)在t=tn的有限元近似解. 最后,記fin=fi(tn).

基于修正標量輔助變量Q, 時間離散采用向后Euler 方法, 得到全離散的有限元方法:

給定(uni,h,pni,h)∈ Xih×Mih和 Qn∈ R, 對于 1 ≤ n ≤ N ?1, 求 (uni,+h1,pni,+h1)∈ Xih×Mih和 Qn+1∈ R 使得對于?(vi,h,qi,h)∈Xih×Mih, i,j=1,2, i/=j,滿足

注意到,在算法(7)～(8)中,(un+1i,h, pn+1i,h)和(un+1j,h, pn+1j,h)是解耦計算的. 但它們?nèi)耘cQn+1耦合在一起,因此構建的離散系統(tǒng)需要解耦才能達到較高的計算效率.

現(xiàn)在, 開始描述如何有效地應用全離散方法(7)求解模型(1). 記Sn+1=(1+tn+1)Qn+1. 令

將式(9)代入式(7), 并整理得:

基于(10)～(11)得到的?un+1i,h和ˉun+1i,h, i=1,2, Sn+1可以從標量方程(8)計算得出, 其計算方程如下:

最后, un+1i,h和pn+1i,h(i=1,2)可以由式(9)得出.

2.3 穩(wěn)定性分析

在本節(jié)中, 研究所提出的全離散有限元算法(7)～(8)的穩(wěn)定性.

則有如下不等式成立

證明 令(7)式中(vn+1i,h,qn+1i,h)=2?t(un+1i,h,pn+1i,h)得

這里我們使用了等式2(a?b,a)=|a|2?|b|2+|a?b|2, ?a, b ∈Rd和非線性項的反對稱形式性質.

標量方程(8)乘2?tQn+1可得

將(14)式和(15)式結合得

此外, 對于(16)的右端項, 有

結合(16)式和(17)式得

最后, 對 (18)式從 n=0, 1, 2,···N ?1 求和, 即完成了證明.

3 數(shù)值模擬

在這一節(jié)中, 通過一些精確解數(shù)值實驗來測試式(7)～(8)的穩(wěn)定性和收斂性. 除此之外, 通過Aggul 等[11]提出的一個實際問題(海岸山或懸崖問題)來測試算法長時間計算的有效性.

3.1 無條件穩(wěn)定性驗證

在本節(jié)中, 通過一個初值問題測試當前算法(7)～(8) 的穩(wěn)定性. 假設? = ?1∪?2且 ?1= [0,1]×[0,1],?2=[0,1]×[?1,0]. 顯而易見, I=(0,1)×{0}. 在 I 上, n1=[0,?1]T, n2=[0,1]T. 在本算例中除界面外,其余邊界選擇為齊次的Dirichlet 邊界條件, 并且設置其外力項fi=0, i=1,2. 選取初始值如下以滿足無散度條件,

固定空間步長 h=1/32,參數(shù)選取為 ν1=2.0e?2, ν2=5.0e?2, κ=1.0. 記 E(ui)=?i(u2i,1+u2i,2)dx. 選擇最終時間為T =5. 圖1 展示了當前算法在不同時間步長?t=0.1, 0.05, 0.01 時能量衰減結果.由圖1 可知,所有的能量曲線都呈現(xiàn)單調的衰減. 此外, 我們發(fā)現(xiàn)該算法的數(shù)值結果并不會隨著時間步長變長而發(fā)散, 從而在數(shù)值上證實了當前方法的無條件穩(wěn)定性.

圖 1 能量耗散圖:?1 (左)和?2 (右)

3.2 收斂率驗證

本節(jié)考慮一個解析解問題來驗證算法的收斂率. 令?1=[0,1]×[0,1],?2=[0,1]×[?1,0]且I=(0,1)×{0}. 模型(1)的解析解如下所示:

右端項f1=(f1,1(t,x,y),f1,2(t,x,y))和f2=(f2,1(t,x,y),f2,2(t,x,y))的選擇必須滿足(u1,p1)和(u2,p2)是式(1)的解. 為了簡便起見,誤差簡記為:

固定參數(shù)ν1=0.05, ν2=0.5,κ=100, 選取5 個不同的空間步長h=1/4, 1/8, 1/16, 1/32 和1/64 驗證空間收斂率. 此外, 固定時間步長?t=h 驗證時間收斂率.

圖2 和圖3 分別展示了當前算法在?1和?2子域上所得的誤差及收斂率. 由圖2 和圖3 可以看到當前的算法表現(xiàn)良好, 并保持了速度、壓力的最優(yōu)收斂速度. 一階時間精度也得到了驗證.

圖 2 空間收斂率: 速度場(左)和壓力場(右)

圖 3 時間收斂率: 速度場(左)和壓力場(右)

3.3 海岸山或懸崖問題

考慮一個海岸山或懸崖實際問題用于驗證當前算法的長時間穩(wěn)定性. 該問題描述大氣層中以拋物線流入海岸山或懸崖上空的氣流與海水相遇時的現(xiàn)象. 計算區(qū)域與文獻[11]一致. 在此區(qū)域中, 在海岸山或懸崖邊界, 以及海洋底部施加了齊次的Dirichlet 邊界條件. 同時, 大氣中的流體流動是以拋物線入口流入, 對其它邊界不設置邊界條件.

固定ν1=0.005,ν2=0.05,κ=0.001,h=1/10 和τ =1/5,取不同的最終時間所得到的速度矢量圖和壓力輪廓圖如圖4～圖7 所示. 從圖中可以看出, 當前算法是穩(wěn)定的且并未出現(xiàn)不符合物理的振蕩. 此外,使用該方法所得出的模擬結果與Aggul 等在文獻[11]中所得到的結果保持一致.

圖 4 速度矢量圖: T =5(左)和T =10(右)

圖 5 速度矢量圖: T =20(左)和T =40(右)

圖 6 壓力輪廓圖: T =5(左)和T =10(右)

圖 7 壓力輪廓圖: T =20(左)和T =40(右)