2.3　偏微分方程

在CFD和海洋數(shù)值模擬中，與其他許多包含流體流場計算的研究一樣，我們主要使用偏微分方程求解，這些方程通常是非線性的（也就是會出現(xiàn)包含因變量和（或）因變量的導數(shù)項），它們不服從解析解，因此需要使用數(shù)值方法進行求解。對于由Navier-Stokes方程控制的流體流問題而言，非線性平流項是產(chǎn)生困難的主要根源，然而，這些項包含一階導數(shù)，不是方程中的最高階項，最高階項通常是包含二階導數(shù)的擴散項，因此是這些擴散項決定了方程的性質(zhì)，也決定了對其進行求解的合適方法。這可以用一個簡單的單一因變量φ（x，t）（其中t可以是時間也可以是空間）的線性二階偏微分方程進行說明

該方程的性質(zhì)只與包含最高階導數(shù)的項有關(guān)，即前三項，其余項無關(guān)緊要（即使較低階的項是非線性的也一樣）。在數(shù)學上，x-t空間中存在兩個基本方向，用特征方程的解對其進行定義為

這兩個方向稱為特征方向，由下面的式子得到

如果兩個特征是不同實數(shù)［（B2-4AC）＞0］，則方程是雙曲線；如果是兩個相同實數(shù)（B2-4AC=0），則方程是拋物線；如果是復數(shù)［（B2-4AC）＜0］，則方程是橢圓。

在物理學中，特征定義了一條線，信息在物理區(qū)域中（x-t）沿著該線進行傳播。波傳播問題（超音速流是一個典型的實例）是由具有實特征的雙曲線方程問題。穩(wěn)定狀態(tài)問題（Stommel問題和Munk問題是經(jīng)典的例子）通常由具有虛特征的橢圓方程描述。特征的實數(shù)值或復數(shù)值性質(zhì)對問題的性質(zhì)非常重要，同時對用來進行求解的方法也很重要。

對于含有一個耦合偏微分方程組和所包含因變量不止一個（假設(shè)為n）的真實的流問題來說，其控制方程可以轉(zhuǎn)化為一個一階偏微分方程組如下

其中Ф是一個包含因變量的矢量，矩陣A和B是因變量x，y和時間t的函數(shù)，P為另一個矢量。該方程的性質(zhì)與矩陣A和B的特征值有關(guān)，如果A的特征值是不同的實數(shù)，則問題是關(guān)于x和t的雙曲線，如果是復數(shù)，則是關(guān)于x和t的橢圓，而B則控制著y-t空間中解的性質(zhì)。

式（2.3.4）中只含自變量x和y的子集定義了二維空間中的穩(wěn)定狀態(tài)問題，它的特征方向由下式（Hoffman和Chiang，1993）的解來決定

該方程特征曲面法向量具有n個解的n階代數(shù)方程，特征曲面就是特征方向my/mx。如果所有特征都是實數(shù)，則系統(tǒng)為雙曲線；如果都是復數(shù)，則為橢圓。

2.3.1　一致性、收斂性以及穩(wěn)定性

該問題的適定性（Hadamard定義）也是很重要的，對于與CFD相關(guān)的數(shù)學問題來說，適定性就是不僅要求解存在并是唯一的，還必須持續(xù)依賴于輔助數(shù)據(jù)。換句話說，控制方程并不能單獨決定唯一解，而是和輔助條件（初始條件和邊界條件）一起決定唯一解。一般數(shù)值問題要求滿足一致性、收斂性和穩(wěn)定性幾個條件，也就是對控制微分方程的離散逼近必須在一種意義上是一致的，即令離散間隔任意趨于0時可以對原始微分方程進行還原，這一點保證所用的形式逼近所考慮的微分方程，而不是其他方程。可以通過離散過程求逆來對一致性進行檢查，具體方法是用鄰近格網(wǎng)點替換泰勒級數(shù)展開式，并要求時間步長Δt和格網(wǎng)大小Δx，Δy以任意方式趨于0，結(jié)果應該是原始偏微分方程，通常使用泰勒級數(shù)展開式求控制微分方程的有限差分近似。因此，通常情況下一致性不構(gòu)成問題，只有在極少的形式中會遇到一致性問題，一個罕見個例是拋物線方程的Dufort-Frankel方案。

計算方案也必須在一種意義上滿足收斂性，即離散間隔減小時，數(shù)值解應該收斂于真實解。通常情況下，對于復雜系統(tǒng)來說收斂性很難證明，雖然通過使用不同格網(wǎng)大小（只受限于增長的累計舍入誤差）進行重復計算經(jīng)驗性地檢查收斂性的原理比較簡單，但由于資源需求大，一般很少從這方面進行驗證。因此，很難保證收斂性。然而，Lax提出的定理對此有很大幫助，Lax等價定理表明，如果一個方案滿足一致性和穩(wěn)定性，那么它也滿足收斂性。對一個滿足一致性的方案來說，只需要對該方案的穩(wěn)定性進行檢查，推理如下：如果一個方案是一致的，那么可以使用泰勒級數(shù)展開式對離散過程求逆從而得到控制偏微分方程，其次，如果該方案也是穩(wěn)定的，那么格網(wǎng)的逐步求精應該能使有限差分解至少在理論上收斂到精確解，因為截斷誤差單調(diào)遞減到0。實際上，舍入誤差會產(chǎn)生干預，但是，收斂性的檢查依然依賴于該理論，而很少使用經(jīng)驗性手段。大多數(shù)時候，模型分辨率的選擇與可用的計算機資源以及穩(wěn)定性需求有關(guān)。

計算算法也必須滿足穩(wěn)定并不易受到舍入誤差的積累和增大帶來的數(shù)值問題的影響。在求解過程中，解的任何自發(fā)性擾動必須在任一點處無增大趨勢，否則，將會導致計算誤差增大，無法得到真實解，從而使最終結(jié)果失去意義。這通常對允許的時間步長強加了一個限制并降低了方案的效率。一個離散方案的穩(wěn)定性非常重要，最常用的方案是Von Neumann方案，然而它的適用范圍僅限于線性系統(tǒng)和域中的內(nèi)點，非線性問題必然要調(diào)用局部離散化方法。

Von Neumann的方法依賴于一個事實：由于方程具有線性性質(zhì)，解的誤差即數(shù)值解和真實解的差異和數(shù)值解本身一樣，也受到相同方程的控制。現(xiàn)在如果認為誤差由具有不同波長nΔx的分量構(gòu)成，并且如果在所有分量在時間上有界的條件下進行研究的話，將得到該方案的一個穩(wěn)定狀態(tài)。L.F.Richardson在第一次世界大戰(zhàn)期間最早嘗試了數(shù)值天氣預測（Richardson，1922），他使用了非穩(wěn)定方案（以及控制方程的不合適形式），因此雖然他的預測完全是開拓性的，卻沒有多大用處。1928年Courant，F(xiàn)riedrichs和Lewy推導出了顯式方案穩(wěn)定性的時間步長約束，20世紀30年代Carl-Gustaf Rossby的工作讓大氣預測渦度守恒方程的應用成為可能（取代原始方程慣性——重力波問題引起的困擾），該方程將棘手的重力波問題去除掉。然而，成功的數(shù)值天氣預測直到40年代晚期電子計算機問世后才成為可能（Richardson估算過，別說預測，僅僅跟上全球范圍內(nèi)最新的天氣情況就需要64,000個真人“計算機”），Charney、Fjortoft和von Neumann在普林斯頓使用渦度守恒方程和第一臺電子計算機（ENIAC）首次非常成功地對天氣進行了預測（雖然從某種程度上說只是對于重力勢的預測）。

在控制變量u（x）的時間演變的一維問題中，Von Neumann技術(shù)包括=ξnexpi（jkΔx），k=1，2，…的替換，其中ξ（=ξn+1ξ-n）是一個復數(shù)，表示第k個空間模式的振幅，稱為放大因子，如果ξ保持有界，即對于波數(shù)k的所有可能取值|ξ|≤1，計算誤差或者減小，或者不再增大，此時解是穩(wěn)定的。這一標準的檢查非常簡單，盡管其代數(shù)方法比較麻煩。需要注意的是，ξ通常是復數(shù)。對于因變量大于1的方程系來說（比如CFD中的方程），Von Neuman方法得到一個放大矩陣ξ，而不是一個標量放大因子，然后需要求取該矩陣的特征值λm，其中m為變量個數(shù)，穩(wěn)定條件為：對于所有m和波數(shù)k的所有可能取值，|λm|≤1。

可以用一維線性平流方程對Von Neuman方法進行說明，該方程的精確解是u=f（x-ct），u=f（x）表示初始條件，隨著時間的增長，初始形式在x的正方向上以平流速度c進行不變的傳播。利用時間步長Δt和格網(wǎng)大小Δx，如果cΔt=Δx或者Courant數(shù)C=cΔt/Δx，很容易得到穿過一個格網(wǎng)的真實解等于1。很顯然如果對數(shù)值解進行求取，那么將解從一個格網(wǎng)點到另一格網(wǎng)點進行傳播的速度不可能超過真實解本身的傳播速度，也就是說，在解變?yōu)榉俏锢淼闹埃珻不可能大于1，通俗地說，所考慮的點的計算依賴域必須包含依賴的物理域（圖2.3.1），實際上，對于隨機舍入誤差來說變得不穩(wěn)定的解會導致計算結(jié)果的失敗。

要在數(shù)學上對其進行證明，則可以將解看成不同波數(shù)的傅里葉級數(shù)的重疊并確保所有級數(shù)都不隨時間而增大。由于方程的線性性質(zhì)，該誤差也滿足相同的方程，這也意味著數(shù)值解中的誤差保持有界，可以用平流方程對此進行說明。考慮波數(shù)k的單一級數(shù)為

將其帶入式（2.5.3）得到振蕩方程為

它的解為U（t）=Re［U（0）exp（-ikct）］。將式（2.3.6）和式（2.3.7）結(jié)合很容易看出，解是傳播保持不變的初始形式。在此，考慮式（2.3.6）的有限差分等價方程為

使用式（2.3.7）的中央（跳步）時間差分為

如果用U（n+1）=ξU（n）定義一個放大因子ξ，則數(shù)值解的穩(wěn)定性要求|ξ|n=|U（n）|/| U（0）|有界，或者說當Δt/t較小時，|ξ|≤o（Δt/t）或|ξ|≤1+o（Δt/t），這是穩(wěn)定性的Von Neumann必要條件。然而，對于隨機的舍入誤差和其他計算誤差而言，充要條件是|ξ|≤1，因為ξ通常是復數(shù)，這意味著在復數(shù)平面（Reξ-Imξ）中，對于方案的穩(wěn)定性，ξ的所有可能取值都必須位于以原點為中心的半徑1的圓內(nèi)。式（2.3.9）得到跳步中心的空間差分方案為

如果對于kΔx的所有容許值有| Csin（kΔx）|≤1，則滿足穩(wěn)定性的充要條件|ξ|≤1，因為kΔx可假設(shè)值為0～2π，則sin（kΔx）的最大可能取值為1，因此方案的穩(wěn)定性要求C≤1，方案是條件穩(wěn)定的。

雖然Von Neumann方法在概念上比較簡單，但實際上，它只有在內(nèi)部域并且不把邊界條件的影響考慮在內(nèi)時有效。矩陣方法確實依賴于將時間步長為n+1時的解用步長為n時的解進行表示的形式為

考慮了空間域中所有格網(wǎng)點處的解。此時，U是一個長度為N的向量，N為格網(wǎng)點個數(shù)，A是N×N矩陣，A的形式可以很容易地用邊界條件下的差分方程得到，很顯然U（n）=A（n）U（0），因此A是放大矩陣，該方案的穩(wěn)定性取決于該矩陣的性質(zhì)，解有界的充要條件是該矩陣的值大小不超過1，即|λn≤1|，n=1，…，N，N個特征值由特征方程|A-λI|的根得到，一般而言，就像上面對Neumann的分析那樣，它們會通過較小的O（Δt/t）超過單位一。而且，如果放大的本征模在初始條件下不存在，則解是不穩(wěn)定的，但事實上，求解過程中的舍入誤差必不可少地要求特征值都不大于1。這個過程可以用擴散方程［見下文式（2.5.21）］進行說明，為了簡單，使用兩級FTCS（時間上為前向，空間上為中央）的方案，即

（N-1）×（N-1）階放大矩陣A變?yōu)?/p>

其中D=vΔt/Δx2，使用的邊界條件=0。該矩陣是三對角矩陣，它的元素為a（=D），b（=1-2D）和c（=D），它的特征值用下面式子得到（Haltiner和Williams，1980）

如果D≤1/2，則所有特征值有界，因此該方案是條件穩(wěn)定的。

Nuemann方法和矩陣方法都是適合于線性方程的最好方法。第三種方法適合于非線性方程，它也是能量方法，對因變量（相當于能量）的平方和的有界性進行檢查。如果能量有界，那么在每個點處的解也有界，并且具有穩(wěn)定性。這可以用式（2.5.3）的FTUS（時間上向前，上游差分）進行說明，因為對它來說，與其他方案相比這種方法更為簡單（Mesiniger和Arakawa, 1976；Haltiner和Williams, 1980）。

由此得到

使用循環(huán)邊界條件得到利用它和式（2.3.16）中要求（1-C）C＞0的Schwarz不等式∑ab＜（∑a2）1/2（∑b2）1/2，得到

因此，只要滿足該不等式，那么方案就是穩(wěn)定的，當然，這要求CFL條件C＜1。

2.3.2　橢圓系統(tǒng)、雙曲線系統(tǒng)及拋物線系統(tǒng)

橢圓方程問題包含物理域邊界條件，在該鄰域中對方程進行求解。這些邊界條件可以是狄利克雷邊界條件，其中規(guī)定因變量的值位于邊界上，也可以是Neumann條件，對其導數(shù)（通常與邊界垂直）進行描述，在Nuemann條件情形中，可以確定解是積分的一個未知常數(shù)。邊界條件還可以是包含因變量與其導數(shù)線性結(jié)合的混合（或Robin）條件。所需邊界條件的數(shù)量取決于方程的階數(shù)。調(diào)和方程只包含二階導數(shù)，它只需要因變量的一個條件。例如，在橢圓拉普拉斯方程描述下的穩(wěn)定狀態(tài)流體流問題中，拉普拉斯方程為

或者等式右側(cè)被f（x，y）替換的泊松方程，它足以對φ或者其邊界上的導數(shù)進行描述。此外，還必須為以下的雙調(diào)和方程描述兩個條件

需要在邊界上對φ及其導數(shù)進行描述。

雙曲線方程多少包含較少耗散或無耗散的傳播。物理域中某一給定點處的擾動與兩個特征一起，從該點以一個有限的速度向該域的“下游”進行傳播，因此這種擾動的影響只有在被這兩個特征劃定界限的楔形域中才會產(chǎn)生。同樣地，物理域中的一個點只會受到以兩個特征為邊界的楔形域中向“上游”傳播的擾動的影響（圖2.3.1，考慮x正方向上的超音速流，對此的理解便很清晰）。由于這種方向性，可以在受穩(wěn)定性標準控制的一致向下游的適當增大條件下，用已知的上游條件和控制方程的逐步積分進行求解，此處的穩(wěn)定性標準包含媒介中擾動（c）的傳播速度。涉及時間和空間的雙曲線問題是初始值問題，需要在已知的初始條件前向時間下用積分逐步對它進行求解，尋找所要求的滿足物理域中邊界條件的解，求解過程中的時間步長受到微擾動的傳播速度（在超音速流中為聲速）的影響，從物理學上來看，就是對于選定的格網(wǎng)大小來說，數(shù)值解的前向時間傳播速度不能小于擾動本身在該格網(wǎng)中的傳播速度：時間步長Δt以Δx/c為界，Δx為物理空間中的格網(wǎng)大小，相當于柯朗數(shù)C=cΔt/Δx必須小于1。更確切地說，計算域必須包含進行時間步長點的物理域，這就是著名的Courant-Friedrichs-Lewy（CFL）條件，它用包含偏微分方程的時間步進式問題對解的效率進行約束，在顯式時間步進方案中，每到下一時間步長t+Δ t時使用一個物理點（一個x）對解進行確定。

可以使用隱式時間步進方法來避開這一限制，但是這一方法需要在下一個時間步長時求取所有物理點處（所有x）的解。因此，盡管時間步長不受穩(wěn)定性的約束，僅受到過渡解的精度要求，但是解通常要求對一個大而稀疏的矩陣求逆，通常很少使用蠻力求逆手段比如傳統(tǒng)的高斯消元法，而是使用迭代方法從初始估計開始，最終收斂到正確解。使用哪一種方法主要取決于簡化程度和計算效率的要求。通常在涉及對方程的二階導數(shù)項進行二階差分的問題中（例如，涉及垂直方向上混合的問題），所得矩陣為三對角矩陣形式或者其他同樣簡單的矩陣形式，因為這種形式可以使用快速算法。否則，能否采用快速矩陣解算器或者能否在所用的迭代方案下使解快速地收斂就成為一個問題，可用的計算機體系結(jié)構(gòu)也是一個問題。因為不同過程之間具有相互聯(lián)系的需求，大型并行計算機的體系結(jié)構(gòu)本質(zhì)上不太適合對涉及全局依賴的問題進行求解，只包含局部依賴的顯式方法使用這種體系結(jié)構(gòu)則會更加有效。

拋物線問題是一種典型問題，它涉及空間域中的性質(zhì)隨時間增長而產(chǎn)生擴散的問題，這種擴散的經(jīng)典實例有熱量擴散、動量的黏性擴散和渦度的黏性擴散等。解的趨勢是最終消除由物理過程的擴散性質(zhì)引起的任何初始不連續(xù)性（而在雙曲線問題中，低階的非線性平流項能產(chǎn)生不連續(xù)性，即使初始狀態(tài)下沒有任何不連續(xù)性存在）。時間依賴問題的穩(wěn)定狀態(tài)限制涉及的空間方向不止一個，并受一個雙曲線偏微分方程控制，將其降低為一個橢圓方程。

因此，拋物線問題和雙曲線問題都是“時間依賴”問題，而橢圓問題是一個“穩(wěn)定狀態(tài)”問題。數(shù)值上，拋物線問題的求解方法與雙曲線問題的求解方法沒有多大差異，它們都是時間依賴初始值問題（初始值問題s），這是最重要的。而橢圓問題采用的求解方法則有很大不同，因為它是穩(wěn)定狀態(tài)邊界值問題（邊界值問題s），它要求整個物理域上的解在邊界處服從一些條件，通常從一個初始估計開始進行迭代求解，最終收斂到真實解。在橢圓問題中，物理域中任何點處的一個擾動也會同時在該域中的其他所有點產(chǎn)生，這與時間依賴的拋物線問題和雙曲線問題不同，這兩個問題中可以每次對一個物理點進行求解。在數(shù)值上，這意味著將橢圓偏微分方程轉(zhuǎn)化為一個代數(shù)方程組，因此也“轉(zhuǎn)化”一個大的矩陣。實際上，迭代方法被用來進行求解，在此，解的收斂速度再次成為一個關(guān)鍵問題。

當進行數(shù)值求解時，初始值問題與邊界值問題之間的差異通常變得模糊起來，因為它們都可以使用迭代矩陣解算器進行求解。實際上，通過對邊界值問題添加含有因變量的時間導數(shù)的項以及從一個估計的初始狀態(tài)開始進行積分直到一個穩(wěn)定的漸進狀態(tài)，可以將一個邊界值問題轉(zhuǎn)化為一個初始值問題，這稱為偽瞬態(tài)方法。事實上，邊界值問題的迭代松弛法可以看成是在時間上進行迭代的方法，換句話說，就是一個初始值問題。