李杰權(quán)

(1. 北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所 計(jì)算物理重點(diǎn)實(shí)驗(yàn)室, 北京 100088; 2. 北京大學(xué) 應(yīng)用物理與技術(shù)中心, 北京 100871)

0 引 言

1959年，華羅庚先生在《大哉數(shù)學(xué)之為用》(數(shù)與量)[1]一文中關(guān)于“時(shí)空”有以下論述：四維空間聽(tīng)來(lái)好象有些神秘，其實(shí)早已有之，即以“宇宙”二字來(lái)說(shuō)，“往古來(lái)今謂之宙，四方上下謂之宇”(《淮南子·齊俗訓(xùn)》)，就是宇是東西、南北、上下三維擴(kuò)展的空間，而宙是一維的時(shí)間。

接著他引用相對(duì)論和生活的常識(shí)敘述了時(shí)空既獨(dú)立又統(tǒng)一的性質(zhì)：愛(ài)因斯坦不再把“宇”、“宙”分開(kāi)來(lái)看，也就是時(shí)間也在進(jìn)行著。每一瞬間三維空間中的物質(zhì)都占有它一定的位置。他根據(jù)麥克斯韋——洛倫茲的光速不變假定，并繼承牛頓的相對(duì)性原理，提出了狹義相對(duì)論。狹義相對(duì)論中的洛倫茲變換把時(shí)空聯(lián)系在一起，但并不消滅時(shí)空特點(diǎn)。如向東走三里再向西走三里就回到原處，但時(shí)間則不然，共用了走六里的時(shí)間。時(shí)間是一去不復(fù)返地流逝著。

華先生意在說(shuō)明數(shù)學(xué)在描述“宇宙之大”中的作用。本文的引用，旨在強(qiáng)調(diào)時(shí)空之間的關(guān)系深刻地植根于計(jì)算流體力學(xué)的模型理解和算法構(gòu)造。離開(kāi)了時(shí)空演化，一切歸于“穩(wěn)態(tài)”。正因?yàn)闀r(shí)空演化，流體的世界才色彩斑瀾、波瀾壯闊。

定性甚至定量描述流體世界的時(shí)空演化并不容易，數(shù)值模擬同樣是貌似簡(jiǎn)單，實(shí)則不易。下面舉一個(gè)眾所周知的例子:

?tu(x,t)+a?xu(x,t)=0

(1)

其中a是一個(gè)常數(shù)，x表示空間坐標(biāo)，t是時(shí)間變量。此模型表示了變量u(x,t)的空間變化(擾動(dòng)、波)與時(shí)間變化(傳播)的關(guān)系，描述了一個(gè)時(shí)空關(guān)聯(lián)的性質(zhì)。

u(x,t+Δt)=u(x-aΔt,t)

(2)

式(2)事實(shí)上比式(1)更加根本和直接，并且不牽涉函數(shù)u(x,t)本身更多性質(zhì),比如正則性。即使u(x,t)是一個(gè)方波或脈沖, 式(2)仍然有效。用平移算子e-aΔt?x來(lái)表示式(2)，則有：

u(x,t+Δt)=u(x-aΔt,t)=e-aΔt?xu(x,t)

(3)

從這里出發(fā)，如果想精確表達(dá)式(1)的輸運(yùn)性質(zhì)，則需要進(jìn)行展開(kāi)(如泰勒展開(kāi)):

u(x,t+Δt)=u(x-aΔt,t)=e-aΔt?xu(x,t)

(4)

為了設(shè)計(jì)實(shí)際工程中的計(jì)算方法(簡(jiǎn)稱算法)，需進(jìn)行必要操作：一是截?cái)嗵幚?，涉及算法與控制方程的相容性，以及精度的概念；二是對(duì)空間變化?xu(x,t)的處理(如有限差分)，不同的處理得到不同的算法。這些都與函數(shù)u(x,t)的正則性和模型(1)的內(nèi)在性質(zhì)息息相關(guān)：正則性越好，精度越高; 模型(1)的時(shí)空關(guān)聯(lián)性需要通過(guò)式(2)精確描述。這些處理充滿技巧和無(wú)限可能，怎樣“令人信服”地設(shè)計(jì)時(shí)空耦合算法且精準(zhǔn)描述相應(yīng)模型的內(nèi)在時(shí)空關(guān)聯(lián)性，是一個(gè)需要深入思考的本質(zhì)問(wèn)題。

近年來(lái)，計(jì)算流體力學(xué)蓬勃發(fā)展，各種算法和相應(yīng)的應(yīng)用軟件目不暇接。有些算法具有堅(jiān)實(shí)的理論基礎(chǔ)，如基于變分原理的有限元方法等。而流體力學(xué)中關(guān)于可壓縮流動(dòng)的研究，由于流場(chǎng)結(jié)構(gòu)非常復(fù)雜(激波、物質(zhì)界面、渦團(tuán)和湍流等),人們對(duì)流動(dòng)(解的)性質(zhì)缺乏很好理解，模型的適定性和相關(guān)數(shù)值模擬中算法的探討相當(dāng)工程化，形式化的表達(dá)往往導(dǎo)致似是而非的結(jié)果，常常是模型具有很好的性質(zhì)，而相應(yīng)數(shù)值模擬的離散模型失去了該保持的性質(zhì)，比如模型的時(shí)空關(guān)聯(lián)性和算法的時(shí)空耦合性等。不同的出發(fā)點(diǎn)導(dǎo)出的數(shù)值算法往往截然不同。本文試圖在這方面做一點(diǎn)嘗試: 基于連續(xù)介質(zhì)假定探討流體力學(xué)模型時(shí)空變化的內(nèi)在關(guān)聯(lián)性，闡述相應(yīng)算法應(yīng)保持的時(shí)空耦合性質(zhì)。

從流體力學(xué)微團(tuán)法的直接建模開(kāi)始思考模型的時(shí)空關(guān)聯(lián)性。通過(guò)對(duì)通量的研究，發(fā)現(xiàn)經(jīng)典的瞬時(shí)Cauchy通量無(wú)法反映時(shí)空關(guān)聯(lián)性，提出了某一時(shí)間段內(nèi)時(shí)間區(qū)間通量，在任意特定時(shí)間段內(nèi)反映流體的動(dòng)力學(xué)過(guò)程[2-3]，進(jìn)而提出有限體積格式的基本原理：在連續(xù)介質(zhì)假定性下，任一時(shí)間段內(nèi)流體通過(guò)控制體邊界的時(shí)間區(qū)間通量關(guān)于邊界的擾動(dòng)是Lipschitz連續(xù)的。這一Lipschitz連續(xù)性恰好體現(xiàn)流場(chǎng)時(shí)空關(guān)聯(lián)性的本質(zhì)特征，實(shí)現(xiàn)了數(shù)學(xué)上流體力學(xué)方程組弱解和物理上積分型平衡律之間的統(tǒng)一，后者其實(shí)就是全離散有限體積格式的物理起源。通過(guò)對(duì)流體力學(xué)方程組時(shí)空關(guān)聯(lián)性質(zhì)的研究，實(shí)現(xiàn)有限體積框架下的時(shí)空耦合。在算法的實(shí)現(xiàn)過(guò)程中，物理量(質(zhì)量、動(dòng)量、能量等)以及它的變化率這一對(duì)量起著重要的作用。一個(gè)量的變化率通過(guò)它的空間變化關(guān)系來(lái)反映。變化率不是一類常規(guī)的數(shù)學(xué)演算，而是反映了內(nèi)在的物理規(guī)律。例如，加速度是速度的變化率，它等于控制體表面受到的力，即應(yīng)力，這是牛頓第二定律。在流體力學(xué)計(jì)算方法中使用這樣的量是自然的也是必須的，這反映了時(shí)空變化的聯(lián)系。在文獻(xiàn)[4]中，筆者倡導(dǎo)的時(shí)空耦合高精度算法正是這一思想的具體表現(xiàn)。本文再次描述實(shí)現(xiàn)時(shí)空耦合算法的一種基本流程，并通過(guò)算例展示時(shí)空耦合算法的數(shù)值表現(xiàn)。

從偏微分方程理論上的Cauchy-Kowalevski方法到數(shù)值解中的Lax-Wendroff方法，都是某種意義下的時(shí)空耦合方法?，F(xiàn)代計(jì)算流體力學(xué)的算法時(shí)空觀常有體現(xiàn)，如迎風(fēng)格式、廣義黎曼解法器[5-7]、守恒元/解元(CE/SE)方法[8]、氣體動(dòng)理學(xué)解法器[9-10]等。最近幾年，筆者致力于思考高精度時(shí)空耦合算法的必要性與理論基礎(chǔ)，但遠(yuǎn)遠(yuǎn)不夠深入。隨著計(jì)算技術(shù)的提高和工程需求的精細(xì)化，這方面的研究應(yīng)該得到重視。應(yīng)當(dāng)建立相應(yīng)的時(shí)空觀，這對(duì)算法的設(shè)計(jì)與實(shí)現(xiàn)很有幫助。故作此文，拋磚引玉。

1 流體力學(xué)方程組的時(shí)空關(guān)聯(lián)性

眾所周知，流體力學(xué)建模過(guò)程常用微團(tuán)法，即在某一時(shí)間段(t,t+δt)，研究控制體Ω(t)=(x,x+δx)內(nèi)流體質(zhì)團(tuán)的運(yùn)動(dòng)。時(shí)間間隔δt非常重要，給出了微團(tuán)運(yùn)動(dòng)的動(dòng)力學(xué)過(guò)程響應(yīng)時(shí)間。為了描述方便，基于連續(xù)介質(zhì)假定，流體力學(xué)方程組寫成如下形式：

給定一個(gè)控制體Ω(t),流體運(yùn)動(dòng)滿足下列的積分平衡律(略去外力場(chǎng)等效應(yīng))：

這里u是密度函數(shù)，n是?Ω(t)的單位外法向。相應(yīng)地，稱:

為控制體Ω(t)上的質(zhì)量，右邊:

為Cauchy通量，f為通量密度函數(shù)，簡(jiǎn)稱為通量函數(shù)，它是u以及空間變化量?u等的函數(shù)f=f(u,?u,…)。方程(5)反映了質(zhì)量時(shí)間變化率與通過(guò)控制體邊界?Ω(t)的通量之間的瞬時(shí)關(guān)系。

Cauchy在特定連續(xù)介質(zhì)假定下，論述了Cauchy通量的性質(zhì)[11-12]，這是他對(duì)連續(xù)介質(zhì)力學(xué)最重要的貢獻(xiàn)[11]。要深入理解這一關(guān)系并不容易，這與Gauss-Green定理密切相關(guān)?；谶@一定理，可以將式(5)寫成散度型偏微分方程(PDE)形式：

?tu+?·F(u,?u,…)=0.

(8)

這里，略去雷諾輸運(yùn)定理等細(xì)致討論，并設(shè)Ω不隨流場(chǎng)運(yùn)動(dòng)，即歐拉型控制體。相應(yīng)地，式(1)中的控制體為拉格朗日型控制體。

方程組(8)對(duì)光滑流動(dòng)是有效的，可由偏微分方程的知識(shí)建立時(shí)空的關(guān)聯(lián)性。注意這里的通量函數(shù)F(u,?u,…)與式(7)略有不同，它涉及雷諾輸運(yùn)定理等有關(guān)變換，這里不做詳細(xì)討論。以后歐拉控制體邊界上的通量用大寫字母表示，否則用小寫字母。對(duì)于實(shí)際流動(dòng)來(lái)說(shuō)，這一轉(zhuǎn)化是非常粗糙的，微妙之處在于Cauchy通量C(?Ω;t)的數(shù)學(xué)性質(zhì)，即正則性?？蓧嚎s流場(chǎng)蘊(yùn)涵豐富的現(xiàn)象，如沖擊波、物質(zhì)界面、湍流等效應(yīng)，流場(chǎng)(或稱方程的解)的正則性非常弱，人們對(duì)它的認(rèn)識(shí)很少，故而該定理的應(yīng)用并不顯然。歷史上有很多研究[13-16]，直到最近，這方面的研究仍是如火如荼[17]，遺憾的是所取得的結(jié)果與實(shí)際的期望仍然相差甚遠(yuǎn)。正如最近的專著[18]的1.3節(jié)有這樣一段論述,“Regarding the right-hand side of (1.3) one needs to keep in mind the following comment concerning the identification of the boundary flux: the drawback of this functional analytic, demonstration is that it does not provide any clues on how the qDmay be computed from A” ，其中qD對(duì)應(yīng)于這里的f(u)·n，A對(duì)應(yīng)于一個(gè)給定的應(yīng)力張量。也就是說(shuō)：人們還不知道如何從應(yīng)力張量中計(jì)算出通量密度函數(shù)。

仔細(xì)審視可以看出Cauchy通量C(?Ω;t)是表示一個(gè)特定時(shí)刻t流動(dòng)的瞬時(shí)行為，方程(5)只是表示一個(gè)非常弱的“時(shí)-空”關(guān)聯(lián)性質(zhì)。現(xiàn)在常用的一個(gè)觀點(diǎn)是下面的“弱解”概念，它來(lái)自于偏微分方程(8)，利用了分布理論。

定義1.1設(shè)Ω是n中的一個(gè)有界區(qū)域，n=1,2,3, 0≤t1≤t2。 對(duì)任意試驗(yàn)函數(shù)如果函數(shù)u(x,t)滿足:

(9)

則稱u(x,t)是方程(8)的弱解。

這個(gè)表達(dá)式比式(5)前進(jìn)了一大步：如果解u(x,t)是光滑的，式(8)和式(9)是等價(jià)的。近年來(lái)偏微分方程理論得到極大的發(fā)展，人們對(duì)于式(8)的認(rèn)識(shí)比直接對(duì)式(5)的認(rèn)識(shí)要深刻得多，且偏微分方程(8)的時(shí)空關(guān)聯(lián)性一目了然：流場(chǎng)的空間攝動(dòng)可以通過(guò)式(8)反映到時(shí)間的變化中, 線性波的傳播式(1)正是這種時(shí)空關(guān)聯(lián)性的特例。另外，當(dāng)一個(gè)間斷面在無(wú)限薄的假定下，可以導(dǎo)出Rankine-Hugoniot間斷關(guān)系，即間斷面(線)兩側(cè)解的“跡”的關(guān)系。

定義1.1并沒(méi)有對(duì)函數(shù)u(x,t)做過(guò)多假設(shè)，只要式(9)成立即可。一旦流場(chǎng)出現(xiàn)復(fù)雜間斷或其它流場(chǎng)結(jié)構(gòu)時(shí)，這個(gè)定義并不能給出太多的信息。我們回避不了的事實(shí)是：在可壓縮流動(dòng)中，流場(chǎng)即使初始性質(zhì)非常“好”，但在不久的將來(lái)由于復(fù)雜的非線性相互作用，可能變得非?！霸愀狻薄＿@可以從最簡(jiǎn)單的Burgers方程看到[19-20]。因此，人們不得不面對(duì)復(fù)雜的情形，深入思考為什么應(yīng)用Gauss-Green定理理解Cauchy通量非常困難。

任何流動(dòng)都有一個(gè)動(dòng)力學(xué)過(guò)程，這也許是個(gè)常識(shí), 但仍需要嚴(yán)格論證。最近的研究表明[3-4]，流動(dòng)的時(shí)空關(guān)系既相互獨(dú)立又彼此關(guān)聯(lián)。下面的原理深刻地反映這一事實(shí)。

有限體積方法基本原理。設(shè)u(x,t)是定義1.1下方程(8)的弱解，Ω是n中任一緊集。假定:

(i)u(x,t)是局部有界；

(ii) 質(zhì)量m(t)是時(shí)間t的連續(xù)函數(shù)。

那么有下列結(jié)論：

(10)

(ii) 對(duì)任意δ>0，時(shí)間段(t,t+δt)內(nèi)流過(guò)Ω邊界?Ω的時(shí)間區(qū)間通量:

關(guān)于Ω的邊界?Ω擾動(dòng)是Lipschitz連續(xù)的。

由此，我們給出方程(8)的弱解另一種表述，它恰恰是物理上常見(jiàn)的積分型平衡律。

定義1.2設(shè)Ω是n中任一緊集,δt>0。如果它滿足下列條件：

(i)u(x,t)局部有界，且質(zhì)量m(t)是時(shí)間t的連續(xù)函數(shù)；

(ii) 整體通量(4)有意義且關(guān)于Ω的邊界?Ω擾動(dòng)是連續(xù)的；

(iii) 積分平衡律成立:

其中n是?Ω的單位外法向, 則稱函數(shù)u(x,t)是偏微分方程(8)的弱解。

事實(shí)上，已經(jīng)證明定義1.1和定義1.2是等價(jià)的[3-4]，并且給出了通量(11)的正則性質(zhì)。這個(gè)定義說(shuō)明后面討論的有限體積方法就是計(jì)算流體力學(xué)的基本格式，和物理最原始的建模一致。流體力學(xué)方程組(8)可理解為：如果流場(chǎng)光滑，用偏微分方程組(8)刻畫；但當(dāng)流場(chǎng)失去正則性時(shí)，解需要滿足積分平衡律式(12)。流體力學(xué)有限體積格式的構(gòu)造過(guò)程就是基于式(12)的直接建模過(guò)程：在偏微分方程(8)誘導(dǎo)出的時(shí)空關(guān)聯(lián)解輔助下，構(gòu)造出和式(12)相容的離散模型，精準(zhǔn)反映真實(shí)的物理過(guò)程。

2 有限體積方法及其時(shí)空耦合性

2.1 半離散有限體積方法

在Ωj上，對(duì)式(8)的空間散度項(xiàng)應(yīng)用Gauss-Green公式,可得半離散有限體積方程:

它直接對(duì)應(yīng)于流體力學(xué)方程組(5)。在初始時(shí)刻t=0, 給定u(x,t)的分布:

將之帶入(13)即得:

(19)

半離散有限體積方法屬于線方法，從某種意義上來(lái)說(shuō)是時(shí)空解耦方法。在光滑流場(chǎng)中，流體力學(xué)方程組的偏微分方程關(guān)系可以直接反映時(shí)空關(guān)聯(lián)性。但對(duì)間斷解來(lái)說(shuō)，式(19)中的局部誤差的累積會(huì)給數(shù)值模擬帶來(lái)巨大傷害。

2.2 全離散有限體積方法

本文將把式(12)應(yīng)用到時(shí)空控制體Ωj×[tn,tn+1],tn+1=tn+Δtn,Δtn>0為時(shí)間步長(zhǎng)，得到有限體積格式的全離散形式:

也可以把式(20)看作是式(13)在時(shí)間段(tn,tn+1)上的積分。與半離散情形下求解瞬時(shí)通量(16)不同，這里需要利用偏微分方程(8)的時(shí)空關(guān)聯(lián)性質(zhì)逼近[tn,tn+1]上時(shí)間區(qū)間通量:

(21)

將之代入(20)，得到有限體積公式:

這里記逼近解在控制體Ωj上的平均值為:

由有限體積方法基本原理，可以得到:

=O(Δtn)q+1|Γjl|

(24)

其中q>0。值得注意的是，式(24)中的誤差和式(19)中的誤差非常不同,式(19)中的誤差是用解的局部跳躍(變差)來(lái)測(cè)量，而式(24)中的誤差直接用時(shí)間步長(zhǎng)(等價(jià)于網(wǎng)格尺度)來(lái)測(cè)量。當(dāng)流場(chǎng)相對(duì)光滑時(shí)，這兩者是等價(jià)的,但當(dāng)流場(chǎng)中有劇烈間斷、脈動(dòng)時(shí)，兩者差別是明顯的。即使網(wǎng)格加密，式(19)也得不到收斂解(在標(biāo)量方程的研究中，網(wǎng)格加密是可以得到收斂解的，原因在于全局變差有界性質(zhì)。到目前為止還沒(méi)有例證可證明該性質(zhì)在流體力學(xué)方程組成立[7])。后文算例4.1可以看出數(shù)值通量的重要性。

2.3 有限體積方法的相容性

所謂相容性，描述的是離散格式與背景方程之間的關(guān)系。傳統(tǒng)上常常使用Taylor展開(kāi)的方式研究數(shù)值格式與偏微分方程(如式(8))之間的相容性，這樣的運(yùn)算只對(duì)光滑流場(chǎng)成立。當(dāng)流場(chǎng)比較復(fù)雜時(shí)，目前大部分的數(shù)值分析某種意義上只是啟發(fā)性的，比如以標(biāo)量方程為模型建立相關(guān)的理論[21]。

對(duì)半離散格式(18)來(lái)說(shuō)，在每個(gè)固定時(shí)刻給出的通量誤差至多為:

=O((Δu)jl)|Γjl|

=O((Δu)j)|Γj|dj

(26)

其中(Δu)j表示在控制體附近解的局部跳躍，|Γj|是Ωj邊的最大長(zhǎng)度，dj=diam(Ωj)。并且在每個(gè)時(shí)間層，隨著Runge-Kutta步的增加，誤差也會(huì)累積。在式(26)中，dj來(lái)源于不同方向上通量差的獲利，在3.5節(jié)可以進(jìn)一步看到。

對(duì)于全離散方法(25)，我們討論每個(gè)時(shí)間層數(shù)值通量的逼近，并期望有下列的估計(jì):

=O(Δtn)q+2|Γj|.

(27)

其中q>0。比較式(26)和式(27)，它們有本質(zhì)的差別，即(Δu)j與Δtn的差別。對(duì)光滑解來(lái)說(shuō)它們是等價(jià)的。這樣我們給出下列的定義。

定義2.1 (有限體積格式的相容性)。設(shè)Ωj是任一控制體，j=1,…,N, 如果對(duì)于某個(gè)q>0,式(27)成立，則有限體積格式(25)相容于平衡律(12)，并具有q階相容性。

如上所述，相容性概念常常相對(duì)與偏微分方程所言。如果式(8)是雙曲守恒律，Lax和Wendroff 給出了有限差分方法的相容性，常稱為L(zhǎng)ax相容性[23]。由于Lax相容性概念影響深遠(yuǎn)，有必要進(jìn)行回顧。

定義2.2 (Lax相容性[23])?？紤]雙曲守恒律方程組:

?tu+?xf(u)=0

(28)

其2p+1點(diǎn)守恒性差分格式為:

g(u,…,u)=f(u)

(30)

則稱差分方法(29)與雙曲守恒律(28)相容。

基于這個(gè)定義，建立了影響深遠(yuǎn)的Lax-Wendroff收斂定理，直至今日，該定理仍被廣泛應(yīng)用。隨著高階精度格式的發(fā)展以及工程應(yīng)用的需求，這個(gè)定理常常被誤用，典型表現(xiàn)為：

(ii) 在此基礎(chǔ)上建立的Lax-Wendroff定理只適用于一致網(wǎng)格或結(jié)構(gòu)網(wǎng)格，對(duì)非結(jié)構(gòu)網(wǎng)格并不適用[24-25]。數(shù)值驗(yàn)證過(guò)程中使用的網(wǎng)格加密方法測(cè)試收斂階，需要倍加小心。

實(shí)際上，定義2.1第一次給出高精度有限體積數(shù)值格式與積分平衡律之間的相容性[3]。有限體積方法與網(wǎng)格的結(jié)構(gòu)無(wú)關(guān)，因此適用于無(wú)結(jié)構(gòu)網(wǎng)格。特別地，收斂階測(cè)試時(shí)一般針對(duì)光滑流進(jìn)行，數(shù)據(jù)重構(gòu)部分能夠達(dá)到應(yīng)有的精度，通量逼近階事實(shí)上就是收斂階。但是，當(dāng)流場(chǎng)含有間斷或其他復(fù)雜結(jié)構(gòu)時(shí)，數(shù)據(jù)重構(gòu)仍存在諸多爭(zhēng)議，通量的逼近效果往往被忽略。通過(guò)以上分析可以看到，如果通量相容性(27)不成立，逼近解不可能收斂。也就是說(shuō)，式(27)是有限體積格式相容性的一個(gè)必要條件。

2.4 再訪Godunov方法

談到流體力學(xué)中的有限體積方法，需要回顧Godunov方法[26]。經(jīng)過(guò)半個(gè)多世紀(jì)的發(fā)展和檢驗(yàn)，它已成為現(xiàn)代計(jì)算流體力學(xué)基石?？紤]一維可壓縮歐拉方程組:

?tu+?xF(u)=0

(31)

即積分平衡律(12)的形式。對(duì)應(yīng)于式(22)有限體積格式為:

?tu+?xF(u)=0,x∈,t>tn,

由于

可知式(35)和式(12)完全一致。因此，從積分平衡律的逼近來(lái)說(shuō)，基于分片常數(shù)的初值逼近，式(35)與式(12)完全相容或通量有無(wú)窮階相容性。從整體逼近式(25)來(lái)說(shuō)，所有的誤差來(lái)自于投影過(guò)程。因此習(xí)慣上稱Godunov格式為一階格式。

如果用逼近的黎曼解法解，界面上的值一般可以寫為:

這個(gè)誤差在強(qiáng)間斷的情形下對(duì)計(jì)算結(jié)果有巨大的傷害。如式(19)所述，這個(gè)誤差并不隨著網(wǎng)格加密而減小。

在多維情形下，例如二維雙曲守恒律情形:

?tu+?xF(u)+?yG(u)=0

(38)

(i) 相對(duì)于界面的法向黎曼問(wèn)題

由于流體力學(xué)方程組的伽利略不變性，通過(guò)旋轉(zhuǎn)變換，總可以設(shè)x-方向?yàn)棣l的外法向，Γjl兩側(cè)單元的值記為uL和uR，法向黎曼問(wèn)題定義為:

?tu+?xF(u)=0, (x,y)∈2,t>0

它的求解和上面的一維黎曼問(wèn)題(34)沒(méi)有本質(zhì)差異。顯然，此解只反映了沿Γjl法向流場(chǎng)的變化情況。

(ii) 頂點(diǎn)處二維黎曼問(wèn)題

這是真正的多維問(wèn)題，可以寫成如下形式[27]:

?tu+?xF(u)+?yG(u)=0,

(x,y)∈2,t>0

u(x,y,0)=u, (x,y)∈Θ

(40)

這里Θ，=1,…,K,是從原點(diǎn)出發(fā)的角形區(qū)域,見(jiàn)圖1。由于有限傳播性質(zhì)，從中心點(diǎn)發(fā)出的復(fù)雜波結(jié)構(gòu)只會(huì)影響有限的區(qū)域。僅從通量的構(gòu)造來(lái)說(shuō)，頂點(diǎn)處解對(duì)界面通量的貢獻(xiàn)占比很小,可以適當(dāng)限制Courant數(shù)，減少頂點(diǎn)周邊流場(chǎng)對(duì)數(shù)值通量的影響，從而可以忽略此問(wèn)題的求解。在移動(dòng)網(wǎng)格方法中，特別是拉格朗日方法中，需要用頂點(diǎn)解確定頂點(diǎn)運(yùn)動(dòng)速度，于是頂點(diǎn)處二維黎曼問(wèn)題(40)的解非常重要。但由于解的結(jié)構(gòu)過(guò)于復(fù)雜，倒不如通過(guò)頂點(diǎn)近似黎曼解法器來(lái)處理更為方便可靠[28]。

圖1 二維黎曼問(wèn)題的初始數(shù)據(jù)分布

黎曼問(wèn)題的(逼近)解被廣泛用到雙曲守恒律，并延展到更一般的流體力學(xué)方程組半離散數(shù)值方法中。與下面廣義黎曼問(wèn)題的解進(jìn)行比較可以發(fā)現(xiàn)，這個(gè)解不能充分反映出流體的動(dòng)力學(xué)過(guò)程；即使從微觀角度，歐拉方程反映了粒子無(wú)窮碰撞的結(jié)果。再者，在每個(gè)控制單元上及初始時(shí)間層，流場(chǎng)處于常狀態(tài)，空間的變化依賴相鄰單元之間的間斷來(lái)實(shí)現(xiàn)。因此，Godunov格式的解不能充分反映瞬時(shí)行為。對(duì)長(zhǎng)時(shí)間、多物理以及多尺度現(xiàn)象，數(shù)值模擬結(jié)果常常出現(xiàn)似是而非的現(xiàn)象。

到目前為止黎曼問(wèn)題只限于對(duì)雙曲守恒律進(jìn)行研究, 相關(guān)的拓展可以從下面的廣義黎曼問(wèn)題研究中看出。

2.5 高階數(shù)值通量與廣義黎曼問(wèn)題

一般地，方程組(8)的廣義黎曼問(wèn)題可以寫成如下形式:

?tu+?·F(u,?u,…)=0，x∈n,t>0

其中Φ(x)=0是定向的n-1維光滑曲面，經(jīng)過(guò)適當(dāng)?shù)淖鴺?biāo)變換，可記該曲面為x=0 (記空間坐標(biāo)為x=(x,y,z))，即式(41)可化為如下形式:

根據(jù)有限體積方法基本原理，需要直接逼近:

這里Aε={(0,y,z);y∈(y0-ε,y0+ε),z∈(z0-ε,z0+ε)}是面x=0上的一部分。在一維情況下，(43)即為:

根據(jù)已有流體力學(xué)方程組相關(guān)的偏微分方程理論與應(yīng)用研究，可以有效地逼近或求解式(42)，滿足與式(27)對(duì)應(yīng)的相容性要求。具體地，式(44)有:

這里需要被積函數(shù)F(u(0;t))滿足關(guān)于時(shí)間t的正則性。對(duì)于式(42)中的初始數(shù)據(jù)，這一點(diǎn)可以得到保障，從而可以對(duì)式(44)進(jìn)行相容性逼近。

定義2.3 (有限體積方法的時(shí)空耦合性)?？紤]時(shí)空關(guān)聯(lián)模型(8)的有限體積方法(25)。如果方程(8)的解可以用來(lái)逼近數(shù)值通量，并使式(25)在定義2.1的意義下相容, 數(shù)據(jù)重構(gòu)也充分利用式(8)的時(shí)空關(guān)聯(lián)信息，則算法(25)是時(shí)空耦合的。

注意文獻(xiàn)[22]給出的Hermite型數(shù)據(jù)重構(gòu)方法是時(shí)空耦合的，該方法已用在兩步四階方法[4,29]中。文獻(xiàn)中數(shù)據(jù)重構(gòu)的時(shí)空耦合性研究還不充分，上述定義中“數(shù)據(jù)重構(gòu)也充分利用式(8)的時(shí)空關(guān)聯(lián)信息”這句話比較模糊，需要做更深入的研究。

2.6 幾個(gè)時(shí)空耦合算法的例子

下面給出幾個(gè)時(shí)空耦合算法的例子。

2.6.1 線性輸運(yùn)方程

對(duì)于線性方程(1), 其有限體積格式可以寫成:

從而由(46)可得:

這就是迎風(fēng)格式。數(shù)值通量的構(gòu)造完全使用了解的時(shí)空關(guān)聯(lián)信息，因此迎風(fēng)格式(49)是唯一的一階時(shí)空耦合格式。眾所周知，迎風(fēng)格式在所有一階穩(wěn)定格式中具有“最優(yōu)”性質(zhì)。

如果初始數(shù)據(jù)是分片多項(xiàng)式，特別是分片線性函數(shù)時(shí):

則式(1)的解式(2)可寫為:

t∈(tn,tn+1)

(51)

直接計(jì)算可得:

以及

更一般地，我們可以利用文獻(xiàn)[29，22]的方式來(lái)構(gòu)造高階時(shí)空耦合方法，或者Lax-Wendroff型的單步高階方法，只是推廣到實(shí)際工程問(wèn)題時(shí)，真正非線性和多維性會(huì)給單步方法帶來(lái)實(shí)質(zhì)性的困難。

2.6.2 線性對(duì)流擴(kuò)散方程

考慮下列方程:

這里物理黏性系數(shù)μ>0。把它寫成散度形式有:

?tu+?x(au-μ?xu)=0,μ>0

(55)

在時(shí)空控制體上，把它表示為式(32)的形式，進(jìn)行通量逼近。文獻(xiàn)[8]中的守恒元/解元方法展示其時(shí)空耦合的屬性。由于其構(gòu)造需要一點(diǎn)篇幅展示，有興趣的讀者可以查閱文獻(xiàn)[8]及其后來(lái)的一系列文章。

相對(duì)來(lái)說(shuō), Navier-Stokes方程組及其相關(guān)時(shí)空關(guān)聯(lián)模型的時(shí)空耦合算法研究較少. 盡管形式上可以進(jìn)行簡(jiǎn)單的時(shí)空對(duì)換處理，但對(duì)耗散項(xiàng)等高階空間導(dǎo)數(shù)項(xiàng)的處理仍需要嚴(yán)格數(shù)學(xué)論證。當(dāng)然，GKS方法間接提供了Navier-Stokes方程的數(shù)值算法，具有時(shí)空耦合性[9,30]。

2.6.3 線性松弛系統(tǒng)

考慮下列簡(jiǎn)單的松弛系統(tǒng):

u(x,0)=u0(x)

(56)

其中τ>0是松弛參數(shù)，a為常數(shù)，v也是常數(shù)。這個(gè)方程的有限體積形式為：

為了構(gòu)造時(shí)空耦合算法，可以用式(56)的解表達(dá)式:

(59)

很顯然，格式(59)～式(61)是時(shí)空耦合的，并具有時(shí)空二階精度。

3 基于廣義黎曼解法器 (GRP solver)的時(shí)空耦合算法

數(shù)值通量的構(gòu)造是執(zhí)行有限體積格式的兩個(gè)核心步驟之一。對(duì)于非線性問(wèn)題我們一般無(wú)法像上述線性方程那樣,得到解的表達(dá)式。下面將要利用廣義黎曼問(wèn)題(41)或(42)解的局部正則性質(zhì)，發(fā)展廣義黎曼解法器。其主要思想可以概括為：

為了方便敘述，這一節(jié)統(tǒng)一把界面放在x=0, 把廣義黎曼問(wèn)題的求解點(diǎn)平移到坐標(biāo)原點(diǎn)(0,0), 初始數(shù)據(jù)表示成式(42)的形式。記:

廣義黎曼問(wèn)題(GRP)解法器：給定控制方程以及形如式(42)的分片光滑函數(shù)作為初始數(shù)據(jù)，求解式(42) 并得到形如式(62)的極限值。

一旦有了極限值(62)，受益于解u(x,t)在界面上關(guān)于時(shí)間t的連續(xù)性，我們有:

u(0,δt)=u*+(?tu)*δt+O(δt2)

(63)

由于u*只表示了一種瞬時(shí)行為，其解由相應(yīng)的黎曼解給定:

u*=R(0;uL(0),uR(0))

(64)

接下來(lái)的關(guān)鍵任務(wù)是求值(?tu)*。“非常不嚴(yán)格”地由控制方程(42)直接得到:

然后求得極限值。這樣解u(x,t)的變化率可以通過(guò)空間的變化反映，就是經(jīng)典的Lax-Wendroff方法[23]，或Cauchy-Kowalevski方法[31]。通過(guò)這種途徑把模型的時(shí)空關(guān)聯(lián)性質(zhì)嵌入到數(shù)值格式中，所得的算法具有時(shí)空耦合性質(zhì)。

3.1 線性GRP解法器

對(duì)于線性系統(tǒng):

?tu+A?xu=Cu+D,t>0

其中A、C、D是常數(shù)矩陣，A可以對(duì)角化，Λ=R-1AR, 其中Λ是由A的特征值λi組成的矩陣，Λ=diag{λi}。記w=R-1u, 有:

?tw+Λ?xw=R-1Cu+R-1D,t>0

(I+R-1CuL+I-R-1CuR)+R-1D

(68)

(RI+R-1CuL+RI-R-1CuR)+D

(69)

特別當(dāng)uL=uR=u*時(shí)有:

從中可以看出如何應(yīng)用Lax-Wendroff解法器和間斷追蹤計(jì)算(?tu)*。

多維情形是類似的。例如考慮二維線性方程組：

?tu+A?xu+B?yu=Cu+D,t>0

這里A、B、C、D是常矩陣，且A、B可對(duì)角化。切向變化量B?yu可看作相對(duì)法向的一個(gè)擾動(dòng), 將線性方程組寫成如下形式：

?tu+A?xu=-B?yu+Cu+D,t>0

與上面類似方法，我們可以得到:

(?tu)*=-[A+(?xu)L+A-(?xu)R]-

RI+R-1B(?yu)L+RI-R-1B(?yu)R]+

(RI+R-1CuL+RI-R-1CuR)+D

(73)

特別強(qiáng)調(diào)，通過(guò)GRP解法器把切向變化自然地蘊(yùn)含到數(shù)值通量中，這是法向(一維)黎曼解法器做不到的[38]。換句話說(shuō)，如果只使用黎曼問(wèn)題解法器，導(dǎo)致格式失去多維性，該格式應(yīng)用到多維問(wèn)題模擬時(shí)自然會(huì)有數(shù)值缺陷，不得不用其他方式進(jìn)行補(bǔ)償。再者，GRP解法器是保證數(shù)值格式真正多維性的有效手段，在合理的數(shù)據(jù)投影基礎(chǔ)上，算法具有渦量保持等性質(zhì)。

3.2 聲波近似 ——線性化GRP解法器

對(duì)于非線性系統(tǒng)來(lái)說(shuō)，當(dāng)流場(chǎng)是光滑的或者波的強(qiáng)度不大時(shí)，使用近似解法器進(jìn)行通量的逼近就可以取得不錯(cuò)的效果，即聲波近似或線性化GRP解法器。Toro等的ADER解法器就是GRP的聲波近似版本[35]。

先看下面的一維雙曲平衡律:

?tu+?xF(u)=H(x,u),x∈,t>0

?tv+A(u*)?xv=H(x,u*),u=u*+v

(75)

這是一個(gè)線性系統(tǒng)。類似于式(66)，從中可以計(jì)算(?tu)=(?tv)*,

這里A±(u*)=R(u*)Λ±(u*)R-1(u*),Λ±(u*)是由Λ(u*)的“±”部分組成Λ+=diag{max(λi),0)}，Λ-=diag{min(λi),0)}。從而

當(dāng)uL(0-0)≠uR(0+0)時(shí)，聲波近似策略如下：以局部黎曼解u*=R(0;uL(0),uR(0))為背景解, 線性化式(74)得到式(75)，從而可得線性化GRP解法器式(77),具體見(jiàn)文獻(xiàn)[32]。

對(duì)于多維系統(tǒng)(仍然以二維為例):

?tu+?xF(u)+?yG(u)=H(x,y,u),t>0

采用聲波近似的策略，線性化這個(gè)方程組得到:

?tv+A(u*)?xv=-B(u*)?yv+H(x,y,u*),t>0

從而得到類于式(73)的解法器, 這里B(u*)=?uG(u)。

3.3 真正非線性GRP解法器

在式(74)中，如果初始數(shù)據(jù)在原點(diǎn)(相對(duì)網(wǎng)格尺寸)有“很大”跳躍，就會(huì)出現(xiàn)強(qiáng)非線性波。線性化GRP解法器不足以分辨強(qiáng)非線性波，需要使用真正非線性GRP解法器，否則就會(huì)出現(xiàn)明顯的數(shù)值缺陷[7,34-35]。一般的界定標(biāo)準(zhǔn)為：

‖uL(0)-uR(0)‖?αΔx

(80)

參數(shù)α>0是一個(gè)重要的經(jīng)驗(yàn)參數(shù)。求解GRP解法器的基本思想是：解析強(qiáng)稀疏波，追蹤強(qiáng)間斷。特別需要指出，在強(qiáng)間斷的情形下，熱力學(xué)效應(yīng)起著重要作用，真正非線性的GRP解法器反映了熱力學(xué)相容的性質(zhì)[7]。

本文不具體給出真正非線性GRP解法器。對(duì)于可壓縮歐拉方程組，可參閱文獻(xiàn)[5，36]。后期發(fā)展的不依賴于坐標(biāo)系統(tǒng)的GRP解法器，詳見(jiàn)文獻(xiàn)[6]。對(duì)一般的雙曲方程組，可參見(jiàn)文獻(xiàn)[37]。

對(duì)于多維情形，仍然可以把切向的變化和間斷面的變化看作擾動(dòng)，考慮有下列形式問(wèn)題:

?tu+?xF(u)=-?yG(u)+H(x,y,u),t>0

從而利用3.1和3.2節(jié)中法向解法器求解式(72)和式(73)，這主要源自于一個(gè)關(guān)鍵的事實(shí)：切向擾動(dòng)在法向的傳播是連續(xù)的。由此切向的變化在通量中得到充分反映，得到的數(shù)值方法具有真正的多維性[38]。

正如前面所討論的那樣，除非涉及(依賴于網(wǎng)格移動(dòng))中心拉格朗日型數(shù)值方法，這里不關(guān)注頂點(diǎn)多維黎曼解法器。由于真正多維黎曼問(wèn)題的理論基礎(chǔ)很不成熟[27,39-41]，真正多維頂點(diǎn)黎曼解法器常常帶來(lái)不穩(wěn)定的數(shù)值結(jié)果[42]。在實(shí)踐中，人們更傾向使用健壯的逼近GRP解法器, 例如，Maire等利用新的框架并結(jié)合拉格朗日型GRP解法器[43]，發(fā)展了穩(wěn)定的中心型高階拉氏方法。

3.4 一般時(shí)空關(guān)聯(lián)模型的GRP解法器

對(duì)于一般的時(shí)空關(guān)聯(lián)模型，如多相流、多介質(zhì)模型[44]和Navier-Stokes方程等，廣義黎曼問(wèn)題解法器可以進(jìn)行類似的構(gòu)造，簡(jiǎn)述如下。

對(duì)于多相流、多介質(zhì)模型，廣義黎曼問(wèn)題解法器可以歸結(jié)為一般的非線性平衡律的范疇，模型的時(shí)空關(guān)聯(lián)性質(zhì)在GRP解法器中可以得到充分體現(xiàn)[45-47]。困難在于介質(zhì)組份以及混合引起的數(shù)值振蕩、物質(zhì)分?jǐn)?shù)為負(fù)等非物理解，涉及物理建模本身的缺陷以及有限體積格式的投影(平均)過(guò)程。由于熱力學(xué)關(guān)系的高度非線性，投影過(guò)程不能反映精確的物理過(guò)程, 例如內(nèi)能的平均過(guò)程導(dǎo)致物質(zhì)界面附近的壓力振蕩. 因此，在投影步實(shí)現(xiàn)時(shí)空耦合，充分利用解的信息也許可以提高相關(guān)數(shù)值質(zhì)量。

對(duì)于Navier-Stokes方程組等宏觀模型，基本的想法是類似的。對(duì)于對(duì)流項(xiàng)(歐拉部分)，使用上述標(biāo)準(zhǔn)的廣義黎曼解法器。需要指出的是：初始數(shù)據(jù)式(41)需要進(jìn)行適當(dāng)?shù)钠ヅ?。也就是說(shuō)該初始數(shù)據(jù)至少是二階以上分片多項(xiàng)式。另外可使用的策略為，對(duì)于對(duì)流占優(yōu)問(wèn)題，采用松弛策略，把相應(yīng)模型雙曲化[48-49]。這樣的策略是基于能量原理——在對(duì)流占優(yōu)情形下，能量集中于初始能量傳播影響區(qū)域內(nèi)(雙曲性質(zhì))。

從Boltzmann方程直接出發(fā)的動(dòng)理學(xué)方法是設(shè)計(jì)流體力學(xué)數(shù)值方法的一條有效途徑。一般來(lái)說(shuō)，很難寫出Boltzmann方程的解，但在特定的近似假定下，如BGK模型[50]，可以類似于式(58)那樣隱式得出解的表達(dá)式，并將之用于數(shù)值通量的構(gòu)造，得出的數(shù)值方法如氣體動(dòng)理學(xué)格式(GKS)[9]、統(tǒng)一氣體動(dòng)理學(xué)格式UGKS[10]等。特別地，GKS可以用來(lái)模擬Navier-Stokes方程描述的流動(dòng)，可以作為Navier-Stokes的解法器來(lái)使用。按照這樣的定義，在通量的構(gòu)造部分，GKS和UGKS解法器顯然是時(shí)空耦合的。

3.5 高精度方法中黎曼解法器和GRP解法器的比較

在高精度數(shù)值方法中，黎曼解法器和GRP解法器都可用來(lái)構(gòu)造數(shù)值通量。前者在高階線方法中常用，如Runge-Kutta型高階方法，后者用在兩步四階方法中。下面仍以一維雙曲守恒律方程組(31)來(lái)比較兩種解法器的差別。

由解關(guān)于時(shí)間t的正則性得到:

所以對(duì)于間斷解來(lái)說(shuō)，通量的截?cái)嗾`差為:

誤差階q=0。不過(guò)，對(duì)于光滑解來(lái)說(shuō)，式(84)中的差賦予了一個(gè)“階數(shù)”，

從而誤差階為q=1。

與上面討論類似，對(duì)于間斷解，GRP解法器有一階精度q=1，而對(duì)于光滑解有二階精度q=2。

3.6 時(shí)空耦合數(shù)值邊界條件

邊界條件的數(shù)值處理是計(jì)算流體力學(xué)中的一個(gè)核心問(wèn)題，大部分的數(shù)值處理基于對(duì)虛擬區(qū)域的延拓(外插技巧)或特征傳播理論[51]。最近逆Lax-Wendroff方法[52]用來(lái)數(shù)值處理邊界條件，這種方法蘊(yùn)含了時(shí)空耦合性。我們認(rèn)識(shí)到，流體力學(xué)方程組的數(shù)值邊界條件等價(jià)于單邊廣義黎曼問(wèn)題的數(shù)值求解[53]。事實(shí)上，先前的工作已經(jīng)認(rèn)識(shí)到單邊黎曼問(wèn)題在數(shù)值邊界條件處理上的重要性[54-55]。這些工作與直接的逆Lax-Wendroff方法相比，有以下特點(diǎn)：

(ⅰ) 在有限體積框架下，計(jì)算區(qū)域的邊界自然為控制體的邊界，所以數(shù)值邊界條件處理等價(jià)于邊界上的通量數(shù)值逼近，它歸結(jié)為單邊黎曼解法器的構(gòu)造。

(ⅱ) 單邊黎曼解法器中的(?tu)*其實(shí)反映了邊界條件上的受力情況，即牛頓第二定律的數(shù)學(xué)表現(xiàn)，這是時(shí)空耦合算法的體現(xiàn)。在文獻(xiàn)[53]中，這一事實(shí)是單邊黎曼解法器的基石。

(ⅲ) 在技術(shù)上，如此處理數(shù)值邊界條件可以盡量少使用虛擬網(wǎng)格。在時(shí)空二階格式中無(wú)需使用虛擬網(wǎng)格； 相比較其它更高階方法，至多使用一半的虛擬網(wǎng)格。

(ⅳ) 無(wú)需使用更高階導(dǎo)數(shù)，避免了人為的復(fù)雜性和虛假物理性質(zhì)等[52]。

單邊黎曼解法器的使用將極大簡(jiǎn)化邊界條件的數(shù)值處理，特別是解決無(wú)結(jié)構(gòu)網(wǎng)格下虛擬區(qū)域的插值問(wèn)題[56]。

3.7 高階時(shí)空耦合算法

高階精度數(shù)值方法對(duì)小尺度問(wèn)題的數(shù)值模擬非常重要，如湍流等?；谟邢摅w積格式的框架式(25)，高階方法需要在數(shù)值通量和數(shù)據(jù)重構(gòu)兩部分具有高階逼近性質(zhì)。在通量的逼近部分，可以采用Lax-Wendroff方法，但流體力學(xué)方程組(8)的非線性性質(zhì)導(dǎo)致高階展開(kāi)過(guò)于復(fù)雜，缺乏實(shí)際工程價(jià)值。更糟糕的是，解的間斷性質(zhì)意味著高階展開(kāi)的運(yùn)算沒(méi)有意義, 因此需要尋求其他方式。在過(guò)去的幾十年中，各種空間重構(gòu)技術(shù)以及Runge-Kutta型時(shí)間推進(jìn)方法取得了巨大的進(jìn)步，可以容易地查閱到相關(guān)文獻(xiàn)。

最近，文獻(xiàn)[29]發(fā)展了兩步四階時(shí)間推進(jìn)方法，成功地融合了Runge-Kutta型方法簡(jiǎn)單性優(yōu)點(diǎn)和基于Lar-Wendroff型解法器的時(shí)空耦合性質(zhì)。文獻(xiàn)[4]詳述了該類方法的特性： (i) 計(jì)算效率。同等條件下其計(jì)算效率是Runge-Kutta方法的一半(二維情形為例)[57]。(ii) 緊致性。由于時(shí)間推進(jìn)步驟減少一半，模版就會(huì)減少一半；時(shí)空耦合的HWENO型重構(gòu)[22, 58]可以使計(jì)算格式更加緊致。(iii) 穩(wěn)定性。已經(jīng)證明其穩(wěn)定區(qū)域比Runge-Kutta大[59]。(iv) 兼容性。該方法很容易和其它方法相結(jié)合，如GKS解法器[60]、多矩方法[61]以及CRP重構(gòu)技術(shù)[62]等。

目前這類方法還局限在“顯式”框架內(nèi)，相關(guān)的研究處于起步階段。為了應(yīng)用的需要，亟須發(fā)展“隱式”兩步四階方法，以適應(yīng)“強(qiáng)剛性”等多尺度問(wèn)題的求解。

3.8 投影過(guò)程時(shí)空耦合性的一點(diǎn)注釋

從而有:

在最近的兩步四階方法中，我們也使用了這一策略[22]，得到的數(shù)據(jù)更加緊致，從而格式也具有緊致性。

4 數(shù)值算例

在可壓縮流體的算法設(shè)計(jì)及數(shù)值模擬中，有大量基準(zhǔn)算例。隨著算法進(jìn)步和工程需求的提高,基準(zhǔn)算例應(yīng)該與時(shí)俱進(jìn)，以面對(duì)相當(dāng)極端的環(huán)境。在文獻(xiàn)[65]中，一系列基準(zhǔn)算例值得用新的數(shù)值算法進(jìn)行測(cè)試。下面幾個(gè)算例，可以看出時(shí)空耦合性的重要性。

4.1 大壓力比(大密度比)問(wèn)題

這個(gè)問(wèn)題又稱為L(zhǎng)eBlanc問(wèn)題，它是經(jīng)典的激波管問(wèn)題的極端情形，從中可以看出數(shù)值方法對(duì)強(qiáng)稀疏波的分辨能力以及它對(duì)激波產(chǎn)生的影響。控制方程為一維歐拉方程(31)，初始數(shù)據(jù)具有如下形式：

(90)

多方指數(shù)取1.4，計(jì)算的終止時(shí)間t=0.12。文獻(xiàn)[66]利用這個(gè)例子檢測(cè)了多個(gè)數(shù)值格式的性能。這里我們同樣用這個(gè)算例針對(duì)性討論本文涉及的一些觀點(diǎn)，數(shù)值結(jié)果來(lái)源于文獻(xiàn)[7,29]。

首先在圖2中展現(xiàn)使用不同解法器的一階格式?？梢钥闯鲆浑AGodunov方法隨著網(wǎng)格加密，可以收斂到精確解；使用逼近解法器，收斂相對(duì)較慢，但仍然具有收斂性。圖3使用了空間二階重構(gòu)，而時(shí)間離散使用一階向前歐拉方法，解法器分別使用黎曼解法器和逼近的黎曼解法器(HLLC,Roe), 數(shù)值表現(xiàn)很差，不具有收斂性。正如在3.5節(jié)所述,當(dāng)空間具有高階精度，使用一階黎曼解法器構(gòu)造數(shù)值通量等，算法不具有相容性；類似的情形反映在圖4, 即使時(shí)間離散使用二階Runge-Kutta方法。

(a) 200網(wǎng)格點(diǎn)

(a) 1000網(wǎng)格點(diǎn)

(a) 1000網(wǎng)格點(diǎn)

圖5考察了在這種極端情況下使用線性化方法的數(shù)值表現(xiàn)，黎曼解使用聲波近似解法器。其實(shí)文獻(xiàn)[34-35]已經(jīng)指出線性化解法器的數(shù)值缺陷。圖6使用了非線性GRP解法器，數(shù)值結(jié)果立即改善，說(shuō)明了強(qiáng)非線性出現(xiàn)后真正非線性GRP解法器的重要性。圖7展示的結(jié)果是基于GRP解法器的兩步四階方法[29]，從中看到用100網(wǎng)格點(diǎn)可以很好分辨間斷，300網(wǎng)格點(diǎn)可以得到完美效果。這是許多數(shù)值方法很難達(dá)到的，從而說(shuō)明了時(shí)空耦合算法的必要性。

(a) 1000網(wǎng)格點(diǎn)

(a) 1000網(wǎng)格點(diǎn)

(a) m=50

4.2 激波和熵波相互作用的問(wèn)題

這個(gè)算例是Shu-Osher算例的推廣，考慮了更極端的情形，又稱為Titarev-Toro問(wèn)題[67]?？刂品匠桃廊粸闅W拉方程，初始數(shù)據(jù)為:

這里使用了GKS解法器[9]得到圖8的數(shù)值結(jié)果，詳細(xì)說(shuō)明見(jiàn)文獻(xiàn)[60]。

圖8 Titarev-Toro問(wèn)題. 這里使用1000網(wǎng)格點(diǎn)數(shù)，空間重構(gòu)采用了不同的重構(gòu)策略，計(jì)算終止時(shí)間t=5

4.3 激波和懸浮剛體的相互作用

圖9 激波和懸浮剛體的相互作用后壓力的分布情況.

4.4 多介質(zhì)激波和氣泡的相互作用

這個(gè)例子展示了激波和氣泡相互作用的情況(圖10)，該問(wèn)題已經(jīng)成為多相流領(lǐng)域一個(gè)標(biāo)準(zhǔn)算例[68]。下面的結(jié)果取自文獻(xiàn)[45]，分別用能量分裂的Godunov方法和GRP方法模擬所得結(jié)果(圖11)。顯然，GRP很好地提高了流場(chǎng)的分辨率。

圖10 激波撞擊氦氣泡的裝置示意圖. 初始時(shí)刻激波的馬赫數(shù)Ms=1.22,計(jì)算網(wǎng)格為2500×890, 其他參數(shù)詳見(jiàn)文獻(xiàn)[45]的算例D

圖11 激波撞擊氣泡不同時(shí)刻的陰影圖，每個(gè)子圖的左邊是Godunov格式的結(jié)果，右邊是GRP的結(jié)果. (a) t=32 μs; (b) t=62 μs; (c) t=72 μs; (d) t=102 μs; (e) t=427 μs; (f) t=674 μs

4.5 各向同性可壓縮湍流的模擬

圖12 各向同性可壓縮湍流的模擬[69]. 其中湍流馬赫數(shù)Mat=1.2, 初始泰勒微尺度雷諾數(shù)Reλ=72, 速度梯度張量第二不變量的等值面Q=25.

5 討論與展望

流體力學(xué)的時(shí)空關(guān)聯(lián)模型刻畫了物理量空間變化在時(shí)間方向上的演化的傳播，如各種波的傳播等。當(dāng)進(jìn)行數(shù)值模擬時(shí)，這種性質(zhì)理應(yīng)得到繼承，相應(yīng)地所用算法應(yīng)該具備時(shí)空耦合特性。實(shí)踐過(guò)程中這個(gè)看似簡(jiǎn)單課題理應(yīng)得到關(guān)注。遺憾的是，相對(duì)于時(shí)空解耦方法，時(shí)空耦合算法需要更好理解模型的物理性質(zhì)或數(shù)學(xué)理論，人們自然“避難就易”。一個(gè)直白的問(wèn)題是——即使物理問(wèn)題的數(shù)學(xué)模型十分完美，當(dāng)相應(yīng)的算法缺乏相應(yīng)完美性質(zhì)時(shí)，其數(shù)值結(jié)果怎么令人信服呢？

本文首先嚴(yán)格論述有限體積方法。正如引言所解釋的原因：(i) 無(wú)論從物理定律的形成還是數(shù)值算法的構(gòu)造，有限體積方法是解流體力學(xué)問(wèn)題最自然的框架，本文總結(jié)了這個(gè)框架形成的數(shù)學(xué)理論和內(nèi)在涵義；(ii) 流體現(xiàn)象的復(fù)雜性(如間斷等)客觀地阻礙了有限體積方法嚴(yán)格數(shù)學(xué)理論的形成，對(duì)于這一框架的認(rèn)識(shí)常出現(xiàn)諸多似是而非的爭(zhēng)論; (iii) 對(duì)于許多實(shí)際工程問(wèn)題，基于無(wú)結(jié)構(gòu)網(wǎng)格的算法設(shè)計(jì)是一個(gè)自然要求，在此之上許多方法相互沖突[24]，有必要從最基本的地方出發(fā)建立一個(gè)根本準(zhǔn)則。幸運(yùn)地是，許多工程軟件，如Fluent, 自然選擇在有限體積框架下生成；許多工程人員當(dāng)遇到難以跨越的困難時(shí)，往往借用有限體積框架度過(guò)難關(guān)。從論證過(guò)程可以看到時(shí)空耦合是算法的根本屬性。

有限體積格式的步驟很簡(jiǎn)單：投影過(guò)程和數(shù)值通量的構(gòu)造。目前CFD算法的大部分研究者將注意力集中于前者，基本上在空間逼近論的范疇內(nèi)進(jìn)行探索。盡管黎曼不變量等物理量以及其他技術(shù)用于重構(gòu)過(guò)程，但是在隨機(jī)選取方法中重要的時(shí)空耦合性質(zhì)相當(dāng)缺乏。 黎曼問(wèn)題的解被用在隨機(jī)選取中，在某種意義下，它可以看成是一種時(shí)空耦合投影方法。本文側(cè)重于后者的討論，給出了有限體積格式與積分平衡律(12)之間的相容性準(zhǔn)則。特別對(duì)于數(shù)值通量介紹了黎曼解法器和GRP解法器，并在3.5節(jié)中進(jìn)行了對(duì)比。簡(jiǎn)單地總結(jié)如下，具體內(nèi)容可以到文獻(xiàn)[2]中查閱。

(i)關(guān)于Riemann解法器。對(duì)于雙曲守恒律(歐拉方程)，當(dāng)初始數(shù)據(jù)是分片常數(shù)時(shí)，由Riemann解法器產(chǎn)生的數(shù)值通量具有無(wú)窮階相容性，即Godunov格式就是積分平衡律；當(dāng)初始數(shù)據(jù)是非常數(shù)的分片光滑函數(shù)時(shí)，Riemann解法器產(chǎn)生的數(shù)值通量對(duì)光滑解是一階相容的，但對(duì)間斷解是不相容的(見(jiàn)3.5節(jié))。也就是說(shuō)對(duì)于高階空間重構(gòu)，如果不能有匹配的數(shù)值通量，產(chǎn)生的數(shù)值格式是不相容的。值得注意的是，如果控制方程包含外力項(xiàng)時(shí)，Godunov格式永遠(yuǎn)是不相容的。

(ii)關(guān)于GRP解法器。GRP解法器給出時(shí)空耦合的數(shù)值通量。對(duì)于光滑流場(chǎng)，GRP解法器得到二階相容的數(shù)值通量；但對(duì)于間斷解只有一階時(shí)間精度。GRP解法器是保證有限體積格式的逼近解收斂的一個(gè)必要條件。

這里之所以再次強(qiáng)調(diào)這點(diǎn)并細(xì)致給出第算例4.1節(jié)的, 原因在于為了看清它與傳統(tǒng)理解上的差異。事實(shí)上，算法時(shí)空耦合性也體現(xiàn)在其他方面，比如最近研究氣體動(dòng)理學(xué)的漸近性質(zhì)時(shí)，提出了統(tǒng)一保持性質(zhì) (unified preserving property, UP) 的概念[71]。對(duì)于一個(gè)動(dòng)理學(xué)方法，不僅應(yīng)該具有歐拉層面的漸近保持性質(zhì)(Asymptotic preserving property, AP)，還應(yīng)該根據(jù)需要具有Navier-Stokes或更深層面的UP性質(zhì)，這個(gè)過(guò)程實(shí)際上是通過(guò)時(shí)空耦合的思想實(shí)現(xiàn)的。該文分別用時(shí)空耦合DUGKS方法以及時(shí)空解耦CLR格式，對(duì)二維不可壓縮Taylor渦進(jìn)行數(shù)值模擬，見(jiàn)圖13，具體見(jiàn)文獻(xiàn)[71]。特別指出的是，這一問(wèn)題本質(zhì)上是多尺度問(wèn)題。在多尺度問(wèn)題及其相關(guān)的研究中，時(shí)空耦合策略是一條有效途徑。

圖13 關(guān)于二維不可壓縮Taylor渦的DUGKS(時(shí)空耦合)和CLR(時(shí)空解耦)模擬的比較

由于篇幅限制，本文只給出了有限體積方法基本原理的相容性準(zhǔn)則以及通過(guò)GRP解法器實(shí)現(xiàn)相容性的技術(shù)細(xì)節(jié)，對(duì)很多根本性質(zhì)還沒(méi)有涉及，如時(shí)空相容格式的熵穩(wěn)定性等。對(duì)于可壓縮流體力學(xué)來(lái)說(shuō)，熵穩(wěn)定性對(duì)應(yīng)于熱力學(xué)相容性[7]，這部分留在將來(lái)的文章中探討。

客觀地說(shuō)，時(shí)空耦合算法的研究非常不成熟，目前只在相對(duì)比較“純粹”的模型上進(jìn)行分析和驗(yàn)證，但是這一思想應(yīng)該加以推廣與應(yīng)用，比如如何將時(shí)空耦合算法的思想應(yīng)用于一般的時(shí)空關(guān)聯(lián)湍流模型[72]、粒子建模和時(shí)空耦合算法等。就算法本身來(lái)說(shuō)，時(shí)空耦合的隱式格式研究還沒(méi)有開(kāi)展，這一領(lǐng)域的研究應(yīng)該是下一階段的重點(diǎn)。

在工程應(yīng)用中時(shí)空耦合的程序顯得更少，一方面是受限算法模塊的歷史延續(xù)性制約，另一方面是工程領(lǐng)域內(nèi)的“理性”力學(xué)越來(lái)越少。加強(qiáng)工程算法的科學(xué)性是需要目前亟待解決的觀念問(wèn)題。

后記。作為一個(gè)數(shù)學(xué)工作者，很忐忑地接受邀請(qǐng)?jiān)诹W(xué)的專業(yè)期刊《空氣動(dòng)力學(xué)學(xué)報(bào)》奉獻(xiàn)此類稿件，畢竟隔行如隔山。但想到數(shù)學(xué)、力學(xué)本就是“一家”，向力學(xué)家“學(xué)習(xí)”本就是“數(shù)學(xué)”的一大研究動(dòng)機(jī)，心里稍微坦然點(diǎn)。

致謝:此文的想法是在與下列合作者合作過(guò)程中逐步形成的，在此一并致謝，他們是Matania Ben-Artzi、齊進(jìn)、汪玥、杜知方、雷昕、徐昆、郭照立、李啟兵，等。特別感謝曹貴瑜博士提供的關(guān)于各向同性可壓縮湍流模擬算例4.5, 感謝徐昆教授、郭照立教授和汪玥副研究員提出許多根本性改進(jìn)意見(jiàn)，感謝陳海波給予了文字上的潤(rùn)色。