李杰權(quán)

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

0 引 言

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

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

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

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

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

(1)

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

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

(2)

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

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

(3)

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

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

(4)

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

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

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

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

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

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

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

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

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

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

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

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

(8)

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

方程組(8)對光滑流動是有效的，可由偏微分方程的知識建立時空的關(guān)聯(lián)性。注意這里的通量函數(shù)F(u,?u,…)與式(7)略有不同，它涉及雷諾輸運定理等有關(guān)變換，這里不做詳細討論。以后歐拉控制體邊界上的通量用大寫字母表示，否則用小寫字母。對于實際流動來說，這一轉(zhuǎn)化是非常粗糙的，微妙之處在于Cauchy通量C(?Ω;t)的數(shù)學(xué)性質(zhì)，即正則性?？蓧嚎s流場蘊涵豐富的現(xiàn)象，如沖擊波、物質(zhì)界面、湍流等效應(yīng)，流場(或稱方程的解)的正則性非常弱，人們對它的認識很少，故而該定理的應(yīng)用并不顯然。歷史上有很多研究[13-16]，直到最近，這方面的研究仍是如火如荼[17]，遺憾的是所取得的結(jié)果與實際的期望仍然相差甚遠。正如最近的專著[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對應(yīng)于這里的f(u)·n，A對應(yīng)于一個給定的應(yīng)力張量。也就是說：人們還不知道如何從應(yīng)力張量中計算出通量密度函數(shù)。

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

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

(9)

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

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

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

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

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

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

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

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

(10)

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

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

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

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

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

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

(iii) 積分平衡律成立:

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

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

2 有限體積方法及其時空耦合性

2.1 半離散有限體積方法

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

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

將之帶入(13)即得:

(19)

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

2.2 全離散有限體積方法

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

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

(21)

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

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

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

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

(24)

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

2.3 有限體積方法的相容性

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

對半離散格式(18)來說，在每個固定時刻給出的通量誤差至多為:

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

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

(26)

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

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

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

(27)

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

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

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

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

?tu+?xf(u)=0

(28)

其2p+1點守恒性差分格式為:

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

(30)

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

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

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

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

2.4 再訪Godunov方法

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

?tu+?xF(u)=0

(31)

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

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

由于

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

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

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

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

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

(38)

(i) 相對于界面的法向黎曼問題

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

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

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

(ii) 頂點處二維黎曼問題

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

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

(x,y)∈2,t>0

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

(40)

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

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

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

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

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

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

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

其中Φ(x)=0是定向的n-1維光滑曲面，經(jīng)過適當?shù)淖鴺俗儞Q，可記該曲面為x=0 (記空間坐標為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)對應(yīng)的相容性要求。具體地，式(44)有:

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

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

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

2.6 幾個時空耦合算法的例子

下面給出幾個時空耦合算法的例子。

2.6.1 線性輸運方程

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

從而由(46)可得:

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

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

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

t∈(tn,tn+1)

(51)

直接計算可得:

以及

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

2.6.2 線性對流擴散方程

考慮下列方程:

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

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

(55)

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

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

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

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

u(x,0)=u0(x)

(56)

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

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

(59)

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

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

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

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

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

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

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

(63)

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

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

(64)

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

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

3.1 線性GRP解法器

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

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

其中A、C、D是常數(shù)矩陣，A可以對角化，Λ=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)

特別當uL=uR=u*時有:

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

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

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

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

?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)

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

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

對于非線性系統(tǒng)來說，當流場是光滑的或者波的強度不大時，使用近似解法器進行通量的逼近就可以取得不錯的效果，即聲波近似或線性化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)

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

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

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

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

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

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

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

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

3.3 真正非線性GRP解法器

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

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

(80)

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

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

對于多維情形，仍然可以把切向的變化和間斷面的變化看作擾動，考慮有下列形式問題:

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

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

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

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

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

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

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

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

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

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

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

所以對于間斷解來說，通量的截斷誤差為:

誤差階q=0。不過，對于光滑解來說，式(84)中的差賦予了一個“階數(shù)”，

從而誤差階為q=1。

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

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

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

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

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

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

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

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

3.7 高階時空耦合算法

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

最近，文獻[29]發(fā)展了兩步四階時間推進方法，成功地融合了Runge-Kutta型方法簡單性優(yōu)點和基于Lar-Wendroff型解法器的時空耦合性質(zhì)。文獻[4]詳述了該類方法的特性： (i) 計算效率。同等條件下其計算效率是Runge-Kutta方法的一半(二維情形為例)[57]。(ii) 緊致性。由于時間推進步驟減少一半，模版就會減少一半；時空耦合的HWENO型重構(gòu)[22, 58]可以使計算格式更加緊致。(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)“強剛性”等多尺度問題的求解。

3.8 投影過程時空耦合性的一點注釋

從而有:

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

4 數(shù)值算例

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

4.1 大壓力比(大密度比)問題

這個問題又稱為LeBlanc問題，它是經(jīng)典的激波管問題的極端情形，從中可以看出數(shù)值方法對強稀疏波的分辨能力以及它對激波產(chǎn)生的影響?？刂品匠虨橐痪S歐拉方程(31)，初始數(shù)據(jù)具有如下形式：

(90)

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

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

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

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

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

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

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

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

(a) m=50

4.2 激波和熵波相互作用的問題

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

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

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

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

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

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

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

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

圖11 激波撞擊氣泡不同時刻的陰影圖，每個子圖的左邊是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é)的時空關(guān)聯(lián)模型刻畫了物理量空間變化在時間方向上的演化的傳播，如各種波的傳播等。當進行數(shù)值模擬時，這種性質(zhì)理應(yīng)得到繼承，相應(yīng)地所用算法應(yīng)該具備時空耦合特性。實踐過程中這個看似簡單課題理應(yīng)得到關(guān)注。遺憾的是，相對于時空解耦方法，時空耦合算法需要更好理解模型的物理性質(zhì)或數(shù)學(xué)理論，人們自然“避難就易”。一個直白的問題是——即使物理問題的數(shù)學(xué)模型十分完美，當相應(yīng)的算法缺乏相應(yīng)完美性質(zhì)時，其數(shù)值結(jié)果怎么令人信服呢？

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

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

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

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

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

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

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

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

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

后記。作為一個數(shù)學(xué)工作者，很忐忑地接受邀請在力學(xué)的專業(yè)期刊《空氣動力學(xué)學(xué)報》奉獻此類稿件，畢竟隔行如隔山。但想到數(shù)學(xué)、力學(xué)本就是“一家”，向力學(xué)家“學(xué)習(xí)”本就是“數(shù)學(xué)”的一大研究動機，心里稍微坦然點。

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