張慶海 李陽

(浙江大學數(shù)學科學學院, 杭州 310027)

不可壓Navier-Stokes方程是流體力學的基本控制方程, 其高精度數(shù)值模擬具有重要的科學意義.本綜述性文章回顧了求解Navier-Stokes方程的投影方法, 重點介紹了時空一致四階精度的GePUP方法.該方法用一個廣義投影算子對不可壓Navier-Stokes方程進行了重新表述, 使得投影流速的散度由一個熱方程控制, 保持了UPPE方法的優(yōu)點.與UPPE方法不同的是, GePUP方法的推導不依賴于Leray-Helmholtz投影算子的各種性質, 并且GePUP表述中的演化變量無需滿足散度為零的條件, 因此數(shù)值近似Leray-Helmholtz投影算子的誤差對精度和穩(wěn)定性的影響非常透明.在GePUP方法中, 時間積分和空間離散是完全解耦的, 因此對這兩個模塊都能以“黑匣子”的方式自由替換.時間積分模塊的靈活性實現(xiàn)了時間上的高階精度, 并使得GePUP算法能同時適用于低雷諾數(shù)流體和高雷諾數(shù)流體.空間離散模塊的靈活性使得GePUP算法能很好地適應不規(guī)則邊界.理論分析和數(shù)值測試結果都顯示, 相對于二階投影方法, GePUP方法無論在精度上還是效率上都具有巨大優(yōu)勢.

1 概 述

1900年, 著名數(shù)學家希爾伯特在國際數(shù)學聯(lián)盟IMU (International Mathematical Union)巴黎會議上提出了23個數(shù)學問題, 這個列表極大地影響了20世紀的數(shù)學發(fā)展.1998年, 菲爾茲獎得主Smale[1]應IMU的邀請, 給出了18個他認為在21世紀非常重要的數(shù)學問題.2000年5月, 克雷數(shù)學研究所發(fā)布了七個千禧年數(shù)學問題[2], 由于這些問題的重要性, 解決其中任何一個都能獲得100萬美元的獎金.在千禧年數(shù)學問題和Smale的18個問題中, 有一個問題同時存在, 它就是三維空間中Navier-Stokes方程全局解的存在性和正則性.

數(shù)學家們對Navier-Stokes方程的重視源于該方程極其廣泛的應用, 尤其是該方程的解對于理解湍流很可能具有根本性的重要意義.如果一個牛頓流體是性質均勻且各向同性的, 其不可壓Navier-Stokes方程 (INSE)的基本形式為

其中t是時間變量,x∈RD是空間變量( D =2,3 ),g是外力項,ν是動力黏性系數(shù),p(x,t) 是壓強, 而u(x,t) 是流速.假定流速場的初始條件u0∈C∞光滑無源, 并且其空間導數(shù)滿足一定的有界條件, 對于外力項的空間時間偏導數(shù)也做類似假定.在千禧年問題中, 空間域?是 R3或3-環(huán)面( R3/Z3), 因此?沒有邊界.但是在實際物理問題中, 往往需要指定(1)式中的邊界條件.除周期邊界條件外, 本文主要討論最普遍也是最重要的無滑移邊界條件:

在無滑移邊界條件下, 中高雷諾數(shù)不可壓流體往往含有多個時間尺度和多個空間尺度的流動結構.因此在進行數(shù)值模擬時, 數(shù)值方法也必須擁有足夠的精度來捕捉所有相關尺度的特征, 這對計算效率提出了極大的挑戰(zhàn).一階和二階方法的計算結果雖然可以收斂, 但是為了達到給定的計算精度,有限的計算資源可能會很快地消耗殆盡.另外, 至關重要的湍流特征往往取決于流速場的梯度或高階導數(shù).例如, 渦量是流速的一階導數(shù), 而轉捩點是切向流速的二階法向導數(shù)為零的位置.低階方法的流速場計算結果雖然可以收斂, 但是對于流速導數(shù)決定的重要物理量往往存在較大的誤差.綜上所述, 我們希望發(fā)展一個時空一致高階方法, 達到以下五個目標:

(Goal-1) 在時間上和空間上同時達到四次收斂階;

(Goal-2) 在每個時間步內只需求解常數(shù)個關于壓強或關于流速的線性方程組, 并且方程組求解的復雜度為O(N) , 這里N是控制體個數(shù);

(Goal-3) 時間積分方法的具體細節(jié)不出現(xiàn)在整體算法中;

(Goal-4) 時間積分與空間離散完全解耦, 也就是說用戶可以在一定范圍內自由選擇空間離散和時間積分模塊組合成一個整體;

(Goal-5) 同時適用于高雷諾數(shù)和低雷諾數(shù)流體.

眾所周知, 當動力黏性系數(shù)很大時, INSE的半離散系統(tǒng)具有剛性, 需要對(1a)式中的擴散項采用隱式時間積分, 否則時間步長的選取將受到數(shù)值穩(wěn)定性非?？量痰南拗?而當動力黏性系數(shù)很小時, 對流占優(yōu)使得INSE的半離散系統(tǒng)不具有剛性, 此時用顯式時間積分方法也可以使時間步長免受穩(wěn)定性的限制.如果不能滿足(Goal-3), 那么針對不同情況更換時間積分方法將是一件非常繁復的任務.反之, 如果能以“黑匣子”的方式調用高階時間積分模塊, 將會實現(xiàn)(Goal-4), 同時能為整體算法提供最大的靈活性, 從而輕松實現(xiàn)(Goal-5).因此, (Goal-5)可以看成是(Goal-3)和(Goal-4)的副產(chǎn)品.

在第2節(jié)預備知識和第3節(jié)經(jīng)典投影方法的基礎上, 本文重點介紹了一個滿足上述所有目標((Goal-1)—(Goal-5))的四階投影方法GePUP[3?5], 該方法達到核心目標((Goal-3)和(Goal-4))的關鍵是策略性地利用投影算子處理不可壓約束條件(1b)式.第4節(jié)推導INSE的另一種表述, 作為GePUP方法的控制方程, 并簡要總結了在有限體積框架下的數(shù)值算法.由于(Goal-4)帶來的靈活性, 用針對不規(guī)則區(qū)域的離散格式替換經(jīng)典有限體積離散格式就得到了一個時空一致四階精度, 適用于不規(guī)則區(qū)域的投影方法.第5節(jié)中的數(shù)值實驗結果驗證了GePUP算法的時空一致四階精度及其捕捉真實物理過程的能力, 其中5.4節(jié)系統(tǒng)地論述了GePUP算法相對二階方法[6,7]在精度和效率上的優(yōu)越性.

2 基礎知識

2.1 Helmholtz分解定理

與流速場不同, 動量方程(1a)中的壓強p既不需要初始條件也不需要邊界條件.這是因為Helmholtz分解定理指出, 有界區(qū)域??RD上任何充分光滑的向量場v?都可以被唯一地分解為一個無源場v和一個無旋場??之和:

該分解可以通過求解Neumann邊界條件下的Poisson方程實現(xiàn):

其中n表示區(qū)域邊界??的單位外法向量.

如果定義歐拉加速度

則可將動量方程(1a)寫成(3a)式的形式:

因此壓強可由歐拉加速度a?在相差一個常數(shù)的意義下由Helmholtz定理唯一確定, 并且在這個確定過程中無需任何p的邊界條件, 或者說p的邊界條件也是由歐拉加速度a?和a在邊界上的法向分量確定的.另外一個從流速場u得到壓強的方法是對動量方程(1a)求散度, 應用(1b)式后可以得到關于壓強的Poisson方程, 從而求解出壓強及其梯度.

2.2 Leray-Helmholtz投影算子

下面的定理在很多偏微分方程的教科書上都可以找到, 如文獻[8]:

定理1如果f和g足夠光滑, 具有Neumann邊界條件的Poisson方程

解存在的充要條件是

并且?在相差一個常數(shù)的意義下是唯一確定的.

容易驗證(4)式滿足可解條件(8)式且??唯一確定.因此對于任意給定的向量場v?都可以唯一確定(3a)式中的無源場v.由這個唯一性可以定義一個線性算子, 稱為Leray-Helmholtz投影算子, 從而把Helmholtz定理中的分解步驟整合其中,

由(3b)式可知, 在任意時刻P都滿足

將Leray-Helmholtz投影算子作用在歐拉加速度a?上可以得到a, 從而得到壓強梯度這個壓強梯度的提取方式與時間無關.

在周期邊界條件下, Leray-Helmholtz投影算子和Laplace算子 ? 是可交換的, 將作用于(1a)式, 再結合(10)式中最后一個性質可得

在應用Method-of-lines (MOL)方法數(shù)值求解INSE時, 由于(1b)式不是一個動力過程而是一個限制條件, 直接對(1)式進行離散得到的常微分方程組很難在求解過程中保證收斂階[9].相比來說, (11)式是關于流速場的一個動力過程, 可以直接應用MOL方法.這是(11)式相對于(1)式的一個顯著優(yōu)勢.

2.3 Laplace-Leray交換子

在無滑移條件下, Leray-Helmholtz投影算子和Laplace算子 ? 是不可交換的.為了度量這兩者的不可交換性并推導出不帶約束條件的關于流速場的動力過程, 定義Laplace-Leray交換子為

為了推導其表達式, 先定義梯度-散度交換子

滿足

其中(14a)式是由(13)式以及 ? 和?·的交換性得到的, 而(14b)式是由(9)式, ? 和?的交換性,(4a)式以及(13)式得到的:

應用于(15)式的第一和第三項, 結合(10)式中最后一個性質, 可得

因此Laplace-Leray交換子滿足

其中I表示恒等算子.由 (9)式知因此是某標量場的梯度, 若令v?為(1)式中的流速u, 得到的標量場稱為 Stokes壓強, 其梯度滿足

由(17)式和(18)式可得

結合(10)式中第二個性質以及(14a)式可得?ps=0 , 即Stokes壓強既無旋也無源, 是調和的.

定義標量pc滿足

應用Leray-Helmholtz投影算子于(1a)式, 再由(6)式, (17)式, (18)式及(20)式可得

因此INSE中的壓強梯度由兩部分構成, 其中?pc平衡外力項和非線性對流項的散度, 而?ps度量了 Laplace算子和Leray-Helmholtz投影算子的非交換性.在ν→0 和ν→+∞這兩種極限情況下, 壓強梯度相應地由?pc和ν?ps所主導, 因而INSE的性質在這兩種情況下會有很大的差別.

3 一些經(jīng)典方法

3.1 鞍點法

鞍點法[10]將 (1)式中的流速和壓強耦合起來同時求解.例如, 由Crank-Nicolson方法可得

由于矩陣A具有鞍點結構, 這個方法被稱為“鞍點法”.同時求解流速和壓強使得線性方程組中的未知量個數(shù)最大化, 因此這種方法對方程組的求解效率提出了很大的挑戰(zhàn)[11].另外, (23)式中矩陣A的推導依賴于空間離散和具體時間積分方法的細節(jié).對于較為復雜的時間積分方法, 如高階Additive Runge-Kutta (ARK)方法, 矩陣A的復雜性往往是難以接受的.因此, 這個方法難以實現(xiàn)(Goal-3)和(Goal-4), 從而也難以做到(Goal-1)和(Goal-2).

3.2 Chorin的一階投影方法

1968年, Chorin[12]提出了以下一階投影方法:

其中C(u?,un) 是對對流項的近似,P=I?GL?1D是對Leray-Helmholtz投影算子的一個近似.該方法首先在忽略壓強梯度項的情況下計算一個輔助流速場u?, 再將u?投影到散度為零的空間上來實現(xiàn)不可壓約束(1b)式.由于流速和壓強在計算域內解耦, 在每個時間步內, 都只需要求解一系列關于流速或壓強的邊值問題, 因此線性方程組的規(guī)模較鞍點法會小很多, 這種效率優(yōu)勢可能是投影方法得到廣泛應用的重要原因之一.

3.3 分步二階投影方法

在Chorin的基礎上, 學者們提出了許多二階投影方法[13,6,14?17], 出發(fā)點是先對INSE進行時間上的二階隱式離散:

這里 [(u·?)u]n+1/2表示對流項在時刻tn+1/2的顯式二階近似.如果再進行空間離散就得到了鞍點法的(22)式.與Chorin的投影方法類似, 可以引入輔助變量將壓強和流速解耦.不同的是, 達到二階精度需要引入兩個輔助變量.

定義2分步二階投影方法由以下三個步驟構成[17]:

(Proj-1) 選擇一個邊界條件B(u?)=0 以求解輔助變量u?:

(Proj-2) 利用投影算子實現(xiàn)不可壓約束(25b)式:

注意?n+1的邊界條件要與un+1的物理邊界條件以及u?的邊界條件B(u?)=0 相容.(Proj-3) 更新壓強pn+1/2:

一個具體的分步二階投影方法由以上步驟中的三個要素確定, 即壓強近似函數(shù)q、輔助變量u?的邊界條件B(u?) 以及壓強更新函數(shù)U.其中投影步驟(Proj-2)通過求解Neumann邊界條件下的Poisson方程實現(xiàn):

將(27a)式代入(26)式, 結合(25a)式知(Proj-3)步驟中的壓強更新公式(28)式可寫為

以下簡要回顧在定義2的框架下的兩個經(jīng)典的二階投影方法.

3.3.1 Bell, Colella & Glaz (BCG)

Bell, Colella和Glaz[6]指定

在文獻[17]中, 作者運用正則模式分析和數(shù)值實驗驗證了在L∞范數(shù)度量下, 該方法對于流速的求解是二階收斂的, 而壓強的計算結果則僅有一階收斂.造成壓強降階的原因是(31)式中對壓強梯度的更新公式與(30)式是不匹配的, 忽略(30)式中的高階項會導致壓強p的精度損失(通常會出現(xiàn)邊界層)以及壓強梯度在邊界上的不準確.的確,(29b)式和(31)式意味著

而這一般是不對的.但若修改(31)式中的壓強梯度更新公式使得其與(30)式一致, 即

則對流速和壓強的求解都能達到二階精度[17].

Martin, Colella和Graves提出了BCG方法的改進版, 稱為MCG方法, 見5.4節(jié).

3.3.2 無壓投影

在Kim和Moin[13]提出的無壓投影方法中,壓強近似變量q不出現(xiàn)在輔助變量u?的求解過程(Proj-1)中:

其中τ為??的單位切向量, 而q=0 這一選擇導致u?與un+1之間差別較大; 如果達到二階精度, 其邊界條件不能簡單地提為u?|??=0.由于n·un+1的正確性總能由投影步驟(Proj-2)保證, 因此從理論上說u?的法向分量n·u?的選擇具有很大的靈活性.但從數(shù)值計算角度看,u?的邊界條件選取方式會影響u?在靠近計算域邊界處的性質, 也會對壓強p在邊界處的性質產(chǎn)生影響.的確, (29a)式和(33c)式意味著

在文獻[17]中, Brown等用數(shù)值實驗展示了除非選擇u?的邊界條件使得其在靠近邊界處光滑, 否則該方法求解壓強時不能達到二階精度.

分步二階投影方法雖然取得了巨大的成功, 但要直接將其推廣到更高階精度是非常困難的, 主要表現(xiàn)在定義2中三個不確定因素的耦合: 若改變時間積分方法, 則三個因素的選取方式全都需要重新推導, 這無疑是極其繁瑣的.其次, 輔助變量u?不具有物理意義, 其邊界條件必須通過具體時間積分方法的細節(jié)來進行推導.另外, 壓強和流速盡管在計算域內部是解耦的, 但它們在邊界上仍然通過輔助變量u?耦合在一起.綜上所述, 用分步二階投影方法很難達到第1節(jié)中的五個目標.也正是因為這些原因, 之前很多學者把投影方法稱為時間分步法, 并認為它是不可能達到四階精度的.

3.4 Gauge方法

在二階投影方法中, 對輔助變量u?提恰當?shù)倪吔鐥l件通常比較困難.Weinan和Liu[15]提出了Gauge方法, 該格式的出發(fā)點是用規(guī)范變量χ來替代壓強p, 他們還引入輔助變量

并策 略性地選擇χ和m使得由此導出的流速u和壓強p滿足INSE.考慮

由(34)式, (35a)式和(1a)式可得

(34)式和(35)式的優(yōu)勢在于可以用規(guī)范變量的自由度來對m和χ提明確的邊界條件, 例如

或對控制方程(34)式, (35)式及邊界條件(37)式離散后得到以下定義.

定義3Gauge方法包括以下步驟:

(GM-1) 求解規(guī)范變量mn+1:

注意,mn+1的邊界條件要與(mn+1??χn+1)|??=0相容,

其中第二個等式是由外插公式(χn+1≈2χn?χn?1)近似后得到的.

(GM-2) 求解流速un+1:

其中χn+1滿足

(GM-3) 如果需要的話, 壓強p可由(36)式得到:

該方法對流速u和壓強p都能達到二階收斂,并且(40)式中使用的外插公式是必要的, 若僅使用滯后值χn則只有一階精度[15].由邊界條件(40)式可知, 該方法與前面提到的方法具有同樣的問題, 即具體算法依賴于時間積分方法的細節(jié).

3.5 PPE方法

Pressure Poisson Equation (PPE)方法[18?24]的核心理念是用一個壓力Poisson方程替換流速的不可壓條件:

其中(44c)式是對動量方程(1a)取散度后結合不可壓約束(1b)式得到的, 而(44d)式是由(44a)式的法向分量及無滑移邊界條件得到的.若將?·u=0作為(44)式的邊界條件, 則(44)式等價于INSE[19].

相對于前面四小節(jié)中的方法, PPE方法最大的優(yōu)勢在于可將時間積分模塊完全看成“黑匣子”:流速u是(44)式中唯一的演化變量, 而壓強p可視為u的隱函數(shù).在任一時刻, 壓強都可以按照(44c)式從流速中提取出來, 并且p和u的邊界條件不再與時間積分方法的細節(jié)相耦合.當應用一個時間積分方法的時候, 只要給出了(44a)式右端項在該方法指定時刻的“樣本”值, 就可以根據(jù)該方法的“契約”得到流速在時間上的高階解.在這個過程中, 無需知道該時間積分方法的任何細節(jié).

PPE方法(44)式也存在一些不足之處.首先,其推導過程假設了流速的不可壓條件總是成立的;由于(44c)式的推導過程依賴于(1b)式, 要分析近似(1b)式產(chǎn)生的誤差對INSE解的影響是很困難的.更重要的是, (44a)式和(44c)式意味著即流速散度的耗散是退化的, 也就是說在應用該方法時對流速散度是沒有控制的.數(shù)值實驗結果顯示, 直接對PPE方法的控制方程進行MOL離散得到的四階方法是不穩(wěn)定的.

3.6 UPPE方法

利用2.3節(jié)中的Laplace-Leray交換子估計,Liu等[22]提出了以下UPPE (Unconstrained PPE)方法:

該方法相對PPE方法的優(yōu)勢在于對(45a)式取散度后結合(45c)式可得因此流速的散度由一個熱方程所控制, 根據(jù)熱方程的極大值原理, UPPE格式(45)式對流速散度?·u的控制要比原始PPE格式(44)式可靠得多.在有限元離散的框架下, 數(shù)值實驗結果[23]印證了UPPE方法的穩(wěn)定性.

4 GePUP: 時空一致四階精度的廣義投影方法

在把(45)式中的UPPE表述應用于基于MOL的有限差分和有限體積方法的時候我們遇到了困難.首先, 直接對控制方程(45)進行MOL離散再求解常微分方程組得到的算法是不穩(wěn)定的.雖然在時間積分的過程中對流速的中間變量進行投影有時可以得到穩(wěn)定的收斂結果, 但是由于控制方程(45)中并沒有出現(xiàn)任何投影算子, 這個對流速中間變量進行投影的額外步驟是直接違背原方程的.另外, 即使可以把控制方程(45)中的一些u用進行替換, 也很難找到一個離散的四階投影算子P滿足Leray-Helmholtz投影性質 (10)式.換言之,離散四階投影P難以滿足以下任何條件:

其中〈u〉和〈?〉是定義在控制體上的平均值,D和G分別是離散散度算子和離散梯度算子.由于計算域的非周期性和四階精度要求, 這一困難在交錯網(wǎng)格上仍然存在.

UPPE表述(45)式的推導過程依賴于Leray-Helmholtz投影的性質(10)式.在這種情況下, 如果仍然把不可壓流速場作為動力系統(tǒng)的演化變量,那么用一個不滿足(47)式的離散投影算子來近似Leray-Helmholtz投影算子產(chǎn)生的離散誤差將對數(shù)值穩(wěn)定性造成怎樣的影響呢? 這個問題在UPPE表述下是難以回答的, 我們認為這也是導致UPPE表述的MOL離散格式不穩(wěn)定的根本原因.

為了在理論上刻畫不滿足(47)式的離散投影算子, 4.1節(jié)定義了一個廣義投影算子(generic projection).在4.2節(jié)中, 我們首先推導Laplace算子和廣義投影算子的交換子, 然后在4.3節(jié)中將它作用在INSE上得到GePUP表述.這三小節(jié)理論推導背后的主要思想是, 既然滿足(47)式的離散投影算子難以找到, 投影以后的流速場也不一定能完全滿足不可壓條件, 那么索性將散度不為零的流速場作為動力系統(tǒng)的主變量.之后四小節(jié)內容討論在GePUP表述基礎上的離散算法.其中4.4節(jié)介紹經(jīng)典的有限體積離散, 4.5節(jié)列舉了一些四階離散算子的形式, 4.6節(jié)簡述了如何把4.4節(jié)和4.5節(jié)中對于規(guī)則區(qū)域的空間離散推廣到不規(guī)則區(qū)域上,4.7節(jié)引入時間積分, 總結了GePUP-IMEX算法的主要步驟.

4.1 廣義投影算子

廣義投影算子是向量空間上的一個線性算子,其輸入向量與輸出向量之差為一個梯度場:

其中?·v=0 不一定成立.對于非周期邊界條件的有界區(qū)域, 還要求廣義投影算子P保持輸入流場的法向邊界條件, 即

無穿透邊界上的Leray-Helmholtz投影算子顯然是一個廣義投影算子.

對一個給定的廣義投影算子P, 由(48)式和(49)式知以下Poisson方程

滿足定理1中的可解條件, 從而可唯一確定?ψ.但是通常

例如ψ??可以是2.3節(jié)中提及的Stokes壓強.

因此, (48)式中的Pv?=v???ψ并不是對v?的唯一分解, 而是由給定的P直接決定了梯度場?ψ.由于許多向量場無法寫成梯度場的形式, 因此(48)式和(49)式定義的廣義投影算子是一種特殊的線性算子, 其特殊性在于總是為輸入的v?選擇一個輸出Pv?, 使得兩者之差為一個梯度場.

為什么可以用廣義投影算子來刻畫那些不滿足(47)式的離散投影呢? 這是因為離散投影通常是Leray-Helmholtz投影的一個高階(p≥4 )近似, 從而有

和??相比, 高階小量O(hp) 往往可以忽略.

4.2 Laplace算子和廣義投影算子的交換子[?,P]

首先定義梯度-散度-投影交換子為

由(48)式、 ? 和?的交換性以及(13)式可得

因此

對(52)式中第三項和最后一項應用P, 有

整理后得到

由(53)式, (54)式和 (51)式可得 Laplace算子和廣義投影算子的交換子為

若P滿足(10)式, 則 [??·,P]≡0 , 此時(55)式簡化為(17)式.

4.3 INSE的GePUP表述

由4.1節(jié)中的討論可知, 當廣義投影算子P作用在一個散度為零的流場u上時, 可以看成是對u加一個梯度擾動項使其散度不一定為零.在文獻[5]中, 我們對動量方程(1a)應用廣義投影算子P并結合 (55)式和(51)式推導了 INSE的GePUP表述.本文的推導過程在形式上略有不同.

由(6)式, (53)式, (5)式以及(15)式可得

其中梯度場?q定義為

由(6)式及廣義投影算子的定義可知 (58)式的右端確實是個梯度場.標量q可由如下Neumann邊值問題確定:

其中(59a)式是對(58)式取散度后由?·和 ? 的交換性得到的, 而(59b)式是 (58)式的法向分量.在無穿透邊界條件下, 由(5)式和(49)式可得

另外, 無滑移條件(2)式意味著對流項在邊界上為零:

再結合(5)式知可用g+ν?u替換 (59b)式中的a?.

定義投影流速

并假設其散度滿足

可以得到無滑移邊界條件下INSE的GePUP表述:

其中(64a)式是由 (57)式和假定(63a)式得到的,(64e)式是由(59a)式和假定(63b)式得到的, 邊界條件(64f)式來自于(59b)式, (60)式以及(61)式.由于廣義投影算子對不可壓流場u的擾動是增加了一個梯度項, 而任何梯度都在Leray-Helmholtz投影算子的零空間中, 因此(64c)式顯然成立.

在(64)式這個三變量系統(tǒng)中, 投影流速w是唯一演化變量, 無源流速u由w唯一確定, 而壓強梯度?q又由u唯一確定.因此(64)式可以看成是關于w的動力過程, 對其進行空間MOL離散之后得到一個常微分方程組, 而求解這個常微分方程組的時間積分方法可作為“黑匣子”使用, 這樣就可以實現(xiàn)第1節(jié)中的(Goal-3).

上述GePUP表述保持了UPPE表述的核心優(yōu)勢, 即投影流速的散度由一個熱方程所控制.確實, 對(64a)式取散度并結合(64e)式可得

若在??上指定邊界條件?·w=0 , 則由熱方程的極大值原理可知?·w在?上隨時間呈指數(shù)衰減①在數(shù)值算法中, 將 ? ·w=0 和 ? ·u=0 附加到(64)式中對流項以及擴散項的計算中.在不規(guī)則邊界附近, 這些條件耦合了不同流速分量之間的多維高階多項式擬合..

GePUP表述相對UPPE表述的優(yōu)勢有以下幾點:

(I) 推導過程不依賴于任何Leray-Helmholtz投影算子P的性質(10)式;

(II) 演化變量w沒有不可壓限制條件;

(III) 投影算子P出現(xiàn)在(64)式的右端, 因此用離散投影算子來近似P時對數(shù)值精度和算法穩(wěn)定性的影響非常清晰.

對一個滿足不可壓條件的w的初始值來說, 如果INSE的解存在且唯一, 熱方程(65)將使得假定(63)式成立.但如果INSE的解本身將爆破的話, 其散度必先爆破.因此只有當INSE在理論上確實有解的時候我們的計算才有意義, 換言之, 我們無法通過重新表述INSE的方式回避其解存在性的千禧年問題.這是假定(63)式背后的主要思想.

4.4 矩形區(qū)域上的空間離散

用正方形網(wǎng)格劃分計算域?, 并以多元索引i∈ZD代表第i個控制體Ci, 其區(qū)域為

控制體i在維度d上的高側界面表示為

其中xO∈RD是固定的坐標原點,h是均勻網(wǎng)格間距, 1∈ZD是分量全為1的多元索引,ed∈ZD表示第d個分量為1, 其余分量全為0的多元索引.

如圖1所示,?在控制體Ci上的平均值定義為

圖1 面平均值和控制體平均值.面平均值用帶 〈·〉 且下標為分數(shù)的符號表示, 而控制體平均值對應下標為整數(shù).因此 和 表 示 面 平 均值, 而 〈 ?〉i 是控制體平均值[5]Fig.1.Notation of face-averaged and cell-averaged values.A symbol with angled brackets denotes either a cell-averaged value if the subscript is an integer multi-index, or a face-averaged value if the subscript is a fractional multiindex.Hence and are face-averaged values while 〈 ?〉i is a cell-averaged value.Horizontal and vertical hatches represent the averaging processes over a vertical cell face and a horizontal cell face,respectively.Light gray area represents averaging over a cell[5].

而?在控制體邊界面上的平均值定義為

這兩種平均值之間的關系可通過微積分基本定理得到[25]:

由于控制體中心的點值和控制體平均值之間相差O(h2) , 二階有限體積方法不需要對它們進行區(qū)分, 而在四階方法中則必須考慮這一差異.此外,作用于平均值上的有限體積算子與作用于點值上的有限差分算子是不同的.

4.5 四階離散算子

離散梯度、離散散度以及離散Laplace算子分別定義如下:

離散散度算子D還可作用在張量平均值上

其中兩標量乘積的面平均值可近似為

而表示橫向離散梯度算子

命題4(71)式—(75)式中的離散算子都具有四階精度:

證明見文獻[4].

如圖2所示, 我們在計算域邊界之外再設置兩層虛擬單元.如果是周期性邊界條件則直接把計算域內對應控制體的值復制到虛擬單元中, 否則由給定邊界條件外插得到虛擬單元的平均值.例如, 對于Dirichlet邊界條件ψ=0 , 用以下五階公式來填充虛擬單元i+ed及i+2ed:

圖2 計算域邊界上的兩層虛擬單元 ( d=1).粗實線代表域邊界, 細實線表示內部控制體, 而虛線表示虛擬單元[5]Fig.2.An example ( d =1 ) of ghost filling for the domain boundary.The thick line represents a physical boundary,the thin solid lines interior cells, and the dotted lines ghost cells[5].

在確定了所有虛擬單元上的平均值之后, 再用格式(71)式—(74)式計算離散算子.這樣做的好處是將邊界條件和連續(xù)算子的離散形式完全解耦.

用(71)式—(73)式可定義離散投影算子

由可知.在周期邊界條件下,P是對Leray-Helmholtz投影算子P的四階近似, 并且可以證明其譜半徑和2-范數(shù)都是1[4].在無滑移邊界下,P不是對稱矩陣, 但是數(shù)值實驗驗證了P的譜半徑不大于1[5].總之,P可以看成是有限維線性空間上的廣義投影算子, 作為對GePUP表述(64)式中Leray-Helmholtz投影算子的近似.

4.6 不規(guī)則區(qū)域上的空間離散

如果計算域?是不規(guī)則的, 我們仍然使用正方形控制體離散一個包含了?的矩形區(qū)域, 并定義控制體流相Vi:=Ci∩?.當Vi=Ci,?時, 分別稱Ci為純流相和無流相控制體, 這兩種以外的情況稱為邊界控制體.記‖V‖為控制體流相V的體積, 定義?在V上的平均值為

對于遠離邊界的純流相控制體, 空間算子L在其上的平均值〈L?〉i可以用4.5節(jié)中的離散算子近似到四階精度.在圖3(a)的例子中,L=?.對于邊界控制體Ci或其附近的純流相控制體, 我們也希望用一些控制體平均值的線性組合來近似〈L?〉Ci.如圖3(b)所示, 主要思路是用選定控制體的平均值以及邊界條件擬合一個多維完全多項式.其核心難點是如何選擇Vi附近的一組控制體流相使得多項式擬合的方程組有唯一解,這個問題最終歸結為在ZD中找到個點滿足一些適定條件.我們提出了一個基于對稱群作用的回溯算法[27], 見圖4.

圖3 在不規(guī)則計算域上用正交網(wǎng)格離散Laplace算子 ? [26]Fig.3.Discretizing the Laplacian on irregular domains with rectangular grids[26].

圖4 從給定起始點開始選擇一組格點以擬合多維n次完全多項式.用擬合的完全多項式可以離散空間算子, 并且能很方便地滿足邊界條件[26]Fig.4.For a given starting point, we choose a set of nodes to fit a high-order, multi-variate, complete polynomial, which facilitates the discretization of spatial operators and the enforcement of boundary conditions[26].

4.7 時間積分

對(64)式在控制體上積分, 并用離散算子近似連續(xù)算子的積分, 可得常微分方程組:

其中外力項〈g〉已知, 未知量是控制體平均值向量〈w〉, 而〈u〉和〈q〉都可以通過〈w〉得到.(81)式的初始條件為〈w〉(t0)=〈u〉(t0) , 邊界條件可由(64)式得到.

由于P是對Leray-Helmholtz投影算子的四階近似, 結合命題4可知常微分方程組(81)是對 GePUP表述(64)式的四階離散.

在求解半離散系統(tǒng)(81)式時, 時間積分方法可以完全作為黑匣子使用.在文獻[5]中我們把GePUP和IMEX方法[28]耦合起來得到了一組求解INSE的半隱式時間推進算法, 稱為GePUPIMEX.這里以Kennedy與Carpenter[29]提出的ERK-ESDIRK方法作為一個例子來闡述GePUPIMEX算法.

對于

ERK-ESDIRK方法在每個時間步的求解步驟如下:

其中上標 (s) 代表中間階段,t(s)=tn+cs?t是該階段對應時刻;ns是總階段數(shù);A,b,c是Butcher表中的系數(shù).A[E],b和c表示一個ERK方法, 而A[I],b和c表示一個ESDIRK方法.在每個中間階段,ψ(s)通過求解線性方程組(83b)得到, 其中右端項是結合前面步驟值計算得到的, 故(83)式確實屬于半隱式格式.另外(83c)式是由得到的.有關該方法的更多細節(jié)請參考文獻[29].

對GePUP半離散系統(tǒng)(81)定義

并套用ERK-ESDIRK方法(83)式的契約, 我們得到如下的GePUP-IMEX算法:

在初始時刻t=t0, 有〈w〉0=〈u〉0≈〈u(t0)〉.

GePUP-IMEX算法(85)式的具體步驟直接來自于將ERK-ESDIRK方法應用到GePUP的半離散系統(tǒng)(81)式上.其中(85b)式的第一步對應于(81a)式和(81b)式, 而(85b)式的第二步對應于

(81c)式, 即用離散投影P將〈w〉(s)映成〈u〉(s), 以

便為下一階段計算X[E](〈u〉(j),t(j)) 做好準備.在(85c)式中, 算子X[E]中的非線性對流項導致〈w〉?的散度要遠大于中間階段〈w〉(s)的散度, 因此我們將該時間步的最終解〈w〉n+1設置為〈u〉n+1.

根據(jù)(85b)式和(81)式, 在每個中間階段, 需要求解三個線性方程組.

(HLS-1) Neumann邊界條件下的Poisson方

程(81b): 用來從〈u〉中提取〈q〉;

(HLS-2)無滑移邊界條件下的Helmholtz方程(85b): 用來在時間上推進〈w〉;

(HLS-3) Neumann邊界條件下的Poisson方程(81c): 用來對〈w〉?進行投影.

求解這些方程組的具體算法見文獻[5].在(85c)式中, 需要求解(HLS-1)和(HLS-3), 因此在每個時間步內, GePUP-IMEX算法(ns=6 )總共需要求解17個線性方程組.我們使用的方法是標準的多重網(wǎng)格V循環(huán), 其中光滑迭代為加權Jacobi方法(ω=2/3 ), 網(wǎng)格轉移算子為單射算子和常值限制算子[30].由于?L屬于本性正型[31,32],具有M-矩陣的良好性質, 例如正定性等.(HLS-2)中的算子的特征值的實部全為正,因此經(jīng)典幾何多重網(wǎng)格方法對三個線性方程組的求解都相當高效.

為了估算穩(wěn)定的庫朗數(shù)范圍, 我們將對流項線性化, 再用與處理對流擴散方程[25]同樣的方法可估算穩(wěn)定的庫朗數(shù)范圍為

綜上所述, GePUP-IMEX算法是將IMEX時間積分方法直接應用到GePUP半離散系統(tǒng)上得到的.這個簡單的思路對于其他時間積分方法如顯式Runge-Kutta方法也是適用的, 具體步驟往往會更簡單.例如, 我們將經(jīng)典四階Runge-Kutta方法應用到GePUP半離散系統(tǒng)上得到的算法稱為GePUP-ERK算法[5].

5 數(shù)值實驗

5.1 不規(guī)則區(qū)域上求解INSE

在不規(guī)則區(qū)域(傾斜正方形)

上給定流速和壓強的精確解

可以根據(jù)動量方程(1a)推導出本節(jié)測試中使用的外力項g.

在計算域(87)上, 用(88)式設置流速場的初始條件和無滑移邊界條件, 用GePUP-ERK算法將INSE的解從t0=0 推進到te=0.1.表1列出了最終時刻結果的誤差和收斂階, 顯然流速u和壓強p都能達到四階收斂.

表1 用GePUP-ERK算法在不規(guī)則計算域(87)上求解以(88)式和(89)式為解析解的INSE得到的誤差和收斂階.Re = 1000, t 0=0 , t e=0.1 , C r=0.2 [5]Table 1.Errors and convergence rates of GePUP-ERK for solving the INSE on the irregular domain (87) with Eq.(88)and Eq.(89) as the analytic solutions.Re = 1000, t 0=0 , t e=0.1 , C r=0.2 [5].

5.2 三維頂蓋驅動方腔流

三維頂蓋驅動方腔內的不可壓縮流動幾乎表現(xiàn)出流體力學中所有常見的現(xiàn)象: 漩渦、二次流、混沌粒子運動、不穩(wěn)定性、過渡、湍流和如展向渦的復雜三維流動結構[33], 因此其數(shù)值模擬具有重要的科學意義.另外由于幾何設定非常簡單, 頂蓋驅動方腔流是檢驗INSE數(shù)值求解方法最常使用的基準測試之一.

在規(guī)則計算區(qū)域 [0,1]×[0,1]×[0,2] 上, 用GePUP-IMEX算法(85)式對三維頂蓋驅動方腔流進行數(shù)值模擬.選擇的具體參數(shù)為: 流速場u=(u,v,w) 的 邊界條件為在x=1 處,v=1 , 其余邊界上u=0 ; 而在初始時刻,u在計算域內全為0,t0=0 及te=12.雷諾數(shù)Re= 1000, 庫朗數(shù)Cr=0.5 , 均勻網(wǎng)格間距h=1/128.下面的數(shù)值結果可以在文獻[5]中找到.

圖5給出了在時刻t=2,4,6,8,10,12 , 對稱平面z=1 內流速場數(shù)值解的流線.數(shù)值模擬剛開始時, 蓋子的移動會沿腔蓋產(chǎn)生很大的剪切應力, 并在t=2 之前形成一個漩渦: 位于漩渦中心控制體的流速要比附近所有控制體流速小50倍左右.隨著模擬繼續(xù)進行, 逆時針旋轉漩渦附近的高速流體逐漸穿透最初靜止的流體, 在t=4 之后不久, 又有另外一個順時針旋轉的漩渦開始形成, 并在t=6變得顯著.然后二次渦開始增強并移動到計算域左上角.最終在t=12 時刻, 主渦、二次渦和流速場達到穩(wěn)定狀態(tài).這些特征與文獻[34]中的實驗結果完全符合.三維與二維頂蓋驅動方腔流動的一個主要區(qū)別是沿翼展方向壁端渦旋的演化和對稱面附近的Taylor-G?rtler型渦旋, 圖6清晰地顯示了時刻t=8,10,12的這些突出特征.

圖5 三維頂蓋驅動方腔流測試中均勻網(wǎng)格上( h =1/128 )對稱平面 z =1 內的流線快照.不同顏色代表不同的 ‖ u‖2 [5]Fig.5.Snapshots of streamlines for the 3D lid-driven cavity test inside the symmetry plane z =1 on a uniform grid with h =1/128.Different colors represent different values of ‖ u‖2 [5].

圖6 三維頂蓋驅動方腔流測試中均勻網(wǎng)格上( h =1/128 )渦量的x-分量等值面.紅色, 橙色, 藍色和青色對應的渦量x-分量值分別為–0.50, –0.25, 0.25和0.50[5]Fig.6.Isosurfaces of the x component of the vorticity vector for the 3D lid-driven cavity test on a uniform grid with h =1/128.The values of the vorticity component for the red, orange, blue, cyan surfaces are –0.50, –0.25, 0.25, and 0.50, respectively[5].

本文還對GePUP-IMEX算法結果和 Guermond等[34]的實驗和計算結果進行了定量比較:圖7和圖8分別給出了主漩渦的中心位置和在平面z=1,0.5 內的流速剖面圖.可以看到, GePUPIMEX算法的計算結果和Guermond等[34]的實驗和計算結果都非常符合.此外, GePUP-IMEX算法的某些計算結果, 例如流速的v分量, 要比Guermond等[34]的計算結果更加符合實驗數(shù)據(jù).

圖7 時刻 t =1,2,4,6,8,10,12 對 稱平面 z =0 上主旋渦的中心位置比較.“ ? ”表示本文的數(shù)值結果(由最小化‖u‖2 得到), 而“+”和“ □ ”分別是Guermond等[34]得到的實驗和計算結果[5]Fig.7.A comparison of the center locations of the primary eddy in the symmetry plane z =0 at time instances t=1,2,4,6,8,10,12.The circles represent numerical results of this work, which are determined by minimizing the speed ‖ u‖2.The experimental and computational results of Guermond et.al.[34] are represented by the crosses and the squares, respectively[5].

圖8 三維頂蓋驅動方腔流在時刻 t =12 的流速剖面圖比較( h =1/128 ).實線表示本文的數(shù)值結果, 即 分別作為 的函數(shù).“ × ”和“ · ”分別表示Guermond等[34]得到的實驗和計算結果[5]Fig.8.A comparison of velocity profiles for the three dimensional lid-driven cavity test ( h =1/128 ) at t =12.Solid lines represent computational results of this work, i.e. as a function of and as a function of The computational and experimental results of Guermond et.al.[34] are represented by the solid dots and crosses, respectively[5].

5.1 節(jié)驗證了基于GePUP表述的INSE投影方法能夠達到時空一致四階精度.在這個基礎上,本節(jié)的結果進一步驗證了GePUP方法精確捕捉實際物理過程的能力.

5.3 二維頂蓋驅動方腔流

本節(jié)在傾斜正方形((87)式)上模擬二維頂蓋驅動方腔流, 并將得到的數(shù)值結果與對應規(guī)則區(qū)域 [0,1]2上的數(shù)值結果進行比較.在由點P=確定的域邊界線段上, 流速場u的邊界條件為

其中

其余計算域邊界上u=0 , 而在初始時刻,u在計算域內全為0,t0=0 ,te=60 , C r=0.1.雷諾數(shù)Re= 1000, 均勻網(wǎng)格間距h=1/128.

圖9和圖10分別給出了在時刻t=4,8,32 時的流速場和渦量快照.在規(guī)則區(qū)域和不規(guī)則區(qū)域上, 這兩個物理量的特征完全一致.

圖9 二維頂蓋驅動方腔流測試中均勻網(wǎng)格上 ( h=1/128) 的流速場Fig.9.Numerical results of the velocity field in the two dimensional lid-driven cavity test on the unit box and the rotated box with h=1/128.

圖10 二維頂蓋驅動方腔流測試中均勻網(wǎng)格上 ( h=1/128) 的渦量快照Fig.10.Numerical results of the vorticity in the two dimensional lid-driven cavity test on the unit box and the rotated box with h=1/128.

5.4 GePUP與一個經(jīng)典二階方法的性能比較

盡管在許多應用中[35?37], 高階方法都被認為要優(yōu)于二階方法, 但有學者也報告了某些四階方法的效率事實上是不盡如人意的[38?40], 原因是這些“四階方法”[38?40]僅能在空間上達到四階精度, 在時間上的精度仍為二階.由于求解誤差主要是由二階時間誤差主導的, 這些方法在求解非定常流動時精度和效率都相對較低.因此即使一個方法在空間上具有高精度甚至譜精度, 如果時間上只有二階精度, 我們仍然認為該方法是一個二階方法.另一方面, 這表明要想獲得一個真正的四階方法并充分發(fā)揮其優(yōu)勢, 必須保證時間上的四階精度.現(xiàn)有文獻中對四階方法與二階方法的比較缺乏系統(tǒng)的論述,本節(jié)回顧文獻[5]中給出的一個性能比較公式.

在比較具有不同收斂階的INSE數(shù)值方法的效率時, 作以下兩個假定:

(ECA-1) 求解線性方程組的時間遠大于顯式計算離散算子所耗費的時間;

(ECA-2) 求解單個線性方程組消耗的CPU時間與其未知量個數(shù)成線性關系.基于這兩個假定, 有:

命題5定義一個ξ階方法相對于一個η階方法的性能加速比為

其中ξ,η∈R+,ξ≥η,Tξ(?) 和Tη(?) 分別是ξ階 方法和η階方法達到給定精度?需要的CPU時間.則對于任一給定問題有[5]:

其中nT是在同一網(wǎng)格下求解該問題時ξ階方法與η階方法所需CPU時間的比值, 而Eξ和Eη分別是該固定網(wǎng)格下相應的誤差,Eξ≤Eη.

證明設hξ(?) 和hη(?) 分別是該ξ階方法和η階方法為達到給定精度?所取的網(wǎng)格間距, 則

其中f=Θ(g) 是 指存在常數(shù)cl,cu>0 使得clg≤f≤cug.

根據(jù)假定(ECA-1)和(ECA-2), 可得

其中cξ和cη是兩個正常數(shù),的冪次中的“+1”是由于初始值需要被推進Θ(1/h) 個時間步才能得到最終解.由性能比Sξ|η(?) 的定義(92)式, (94)式和(95)式可得

其中cξ|η對于該給定問題是個常數(shù).在相同網(wǎng)格上,Eξ和Eη是運行ξ階方法和η階方法得到的精度, 而nT是這兩個方法消耗CPU時間的比值.假設η階方法的精度由Eη提升到Eξ需要對網(wǎng)格加密nη→ξ次, 且相鄰兩網(wǎng)格間距比率為r, 則

其中(97a)式中rη的冪次η是η階方法的收斂階,(97b)式是由(95)式得到的.在(96)式中取?為Eξ,再結合(97b)式可得

將(98)式代入 (96)式即可得到(92)式.

利用(93)式, 可以對四階投影方法GePUPIMEX和Martin等[7]提出的經(jīng)典二階投影方法MCG①MCG是在C++軟件包 Chombo[41]的基礎上實現(xiàn)的, 其源代碼見文獻[42].的效率作比較.

表2給出了在不同的測試中為了達到同樣的L2流速計算精度?, GePUP-IMEX相對于MCG的性能加速比S4|2.對于較低的相對精度, 如?=10?4, GePUP-IMEX相對于MCG在效率上的優(yōu)勢還不太明顯, 而要達到適中的相對精度, 如?=10?6,10?8, GePUP-IMEX的性能優(yōu)勢是相當明顯的.盡管對于其他問題,S4|2的具體值可能與表2中的結果不同, 但是S4|2(?) 是 1 /?的冪函數(shù)這一增長模式是始終保持不變的.此外, 即使數(shù)值精度不是首要考慮的因素, GePUP-IMEX也允許使用更粗糙的網(wǎng)格來獲得與MCG相同的精度, 因此在實際計算時GePUP-IMEX依然具有巨大的效率優(yōu)勢.

由于渦量?×u是由流速場的一階導數(shù)得到的, GePUP-IMEX和MCG對渦量的計算在L∞范數(shù)下分別具有三階和一階精度, 表3給出了在不同測試中為了達到同樣的L∞渦量計算精度?,GePUP-IMEX相對于MCG的性能加速比S3|1.即使對于很低的相對精度, 如?=10?2, GePUPIMEX相對于MCG在效率上的優(yōu)勢也是非常明顯的.

以上展示的測試均是在作者的個人電腦(Intel?i7-3930, 138 GFlop/s)上進行的, 截至2020年11月,世界上最快的超級計算機是日本的 Fugaku(7630848 核, peak performance 537212 TFlop/s)[43],因此由個人電腦切換到超級計算機上進行計算所帶來的性能提升最大為 4×106.而注意到表2和表3中許多結果是大于 4×106的, 這說明為了在一些測試問題中達到給定精度要求, 在作者的個人電腦上運行四階GePUP-IMEX算法比在國際上最快的超級計算機上運行二階MCG方法都要快!

表2 不同的測試中為了達到同樣的 L2 流速計算精度 ? , GePUP-IMEX相對于MCG[7]的性能加速比 S 4|2 [5]Table 2.The speedup S 4|2 of GePUP-IMEX over MCG[7] for achieving the same L 2 accuracy ? of the velocity[5].

在數(shù)值模擬高雷諾數(shù)不可壓流體的物理過程時, 往往要用到依賴于流速的一階和二階導數(shù)的物理量, 例如渦量和曲率等.表3中的結果說明了四階GePUP-IMEX投影方法對精確捕捉高雷諾數(shù)流體渦量場的演化過程具有重要的科學意義.

表3 不同的測試中為了達到同樣的 L ∞ 渦量計算精度 ? , GePUP-IMEX相對于MCG[7]的性能加速比 S 3|1 [5]Table 3.The speedup S 3|1 of GePUP-IMEX over MCG[7] for achieving the same L ∞ accuracy ? of the vorticity[5].

隨著流體力學的飛速發(fā)展, 我們越來越需要研究具有多長度尺度和多時間尺度的復雜物理問題,這對數(shù)值模擬在精度上和效率上都提出了越來越高的要求.四階和更高階精度的算法極大地提高了計算精度和計算效率, 若進一步將這些算法與并行計算和自適應網(wǎng)格結合, 我們將能解決更多的實際物理問題.

6 結 論

本文首先回顧了一些求解不可壓Navier-Stokes方程的經(jīng)典數(shù)值方法, 然后基于廣義投影算子和UPPE表述推導了GePUP表述和相應的數(shù)值方法.GePUP方法最突出的優(yōu)勢在于時間積分方法和空間離散格式是完全解耦的, 因此這兩個算法模塊都能夠作為“黑匣子”進行處理.時間積分方法的靈活自由選擇便于實現(xiàn)INSE解在時間上的高階精度, 并使得GePUP格式同時適用于高低雷諾數(shù)流體.空間離散的靈活自由選擇使得算法可以很好地適應不規(guī)則邊界.數(shù)值測試結果印證了該方法的四階精度及其捕捉實際物理過程的能力, 再結合一個簡單的效率比較公式就可以看出四階GePUP方法相對于現(xiàn)有二階投影方法在效率上和精度上都具有巨大優(yōu)勢.

目前我們已經(jīng)將GePUP方法推廣到靜止不規(guī)則區(qū)域上, 我們計劃將其與高階界面追蹤方法[44,45]耦合, 針對動邊界不可壓流體提出一個時空一致四階精度的投影算法.