差分方法基礎(chǔ)

1、應(yīng)用計算流體力學(xué)講稿 羅東明 1 第二講有限差分法基本原理 一般的流體控制方程都是非線性的偏微分方程。在絕大多數(shù)情況下，這些偏微分 方程無法得到精確解；而 CFD 就是通過采用各種計算方法得到這些偏微分方程的數(shù)值 解，或稱近似解。當(dāng)然這些近似解應(yīng)該滿足一定的精度。目前，主要采用的 CFD 方法 是有限差分法和有限體積法。本講主要介紹有限差分法，它也是下一講中的有限體積 法的基礎(chǔ)1。 有限差分法求解流動控制方程的基本過程是：首先將求解區(qū)域劃分為差分網(wǎng)格， 用有限個網(wǎng)格點代替連續(xù)的求解域，將待求解的流動變量(如密度、速度等)存儲在 各網(wǎng)格點上，并將偏微分方程中的微分項用相應(yīng)的差商代替，從而將偏微分

2、方程轉(zhuǎn)化 為代數(shù)形式的差分方程，得到含有離散點上的有限個未知變量的差分方程組。求出該 差分方程組的解，也就得到了網(wǎng)格點上流動變量的數(shù)值解。 2.1差分和逼近誤差 由丁通常數(shù)字計算機只能執(zhí)行算術(shù)運算和邏輯運算， 因此就需要一種用算術(shù)運算 來處理函數(shù)微分運算的數(shù)值方法。 而有限差分法就是用離散網(wǎng)格點上的函數(shù)值來近似 導(dǎo)數(shù)的一種方法。 設(shè)有 x 的解析函數(shù)y=f(x),從微分學(xué)知道函數(shù)y對 x 的導(dǎo)數(shù)為 蟲=飄.=貯 f(x“x)-f(x) (2-1) dx x0 x 0 x dy、dx分別是函數(shù)及自變量的微分，dy/dx是函數(shù)對自變量的導(dǎo)數(shù)，乂稱微商。相 應(yīng)地，上式中的Ax、Ay分別稱為自變量及函

3、數(shù)的差分， Ay/Ax為函數(shù)對自變量的差 商。在導(dǎo)數(shù)的定義中 Ax 是以任意方式逼近丁零的，因而 Ax 是可正可負的。在差分方 法中，Ax總是取某一小的正數(shù)。這樣一來，與微分對應(yīng)的差分可以有三種形式： 向前差分 Ly = f (x ：x) - f (x) 向后差分 * = f (x) - f (x -，x) 1 . 1 .、 中心差分 y = f(x x)-f(x x) 2 2 上面談的是一階導(dǎo)數(shù)，對應(yīng)的稱為一階差分。對一階差分再作一階差分，就得到 二階差分，記為 A2y。以前向差分為例，有應(yīng)用計算流體力學(xué)講稿 羅東明 2 2 . 、 L y = Gy) =.：f (x wx) - f (x)

4、 1 =Af (x+Ax) - Af (x) (2-2) -f (x 2x) - f (x LX) L (x wx) - f (x) 1 依次類推，任何階差分都可以由低一階再作一階差分得到。 函數(shù)的差分與自變量的差分之比，即為函數(shù)對自變量的差商。如一階向前差商為 .：y f (x 匚,；x) f (x) 已x LX 一階后差商為 :y f(x)-f(xr；x) LX LX 一階中心差商為 , 1 , 1 .y f(x 2 “)-伯-2 x) x x 或 ：y f (x ：x) - f (x ：x) 二階差商多取中心格式，即 圖2.1差商與導(dǎo)數(shù)的關(guān)系 差商與導(dǎo)數(shù)的關(guān)系可見圖 2.1。由導(dǎo)數(shù)（微商

5、）和差商的定義知道，當(dāng)自變量的 差分（增量）趨近丁零時，就可以由差商得到導(dǎo)數(shù)。因此在數(shù)值計算中常用差商近似 代替導(dǎo)數(shù)。差商與 2 . 2y _ f(x x) 2f (x) f(x ：x) Gx 應(yīng)用計算流體力學(xué)講稿 羅東明 3 導(dǎo)數(shù)之間的誤差表明差商逼近導(dǎo)數(shù)的程度，稱為逼近誤差。由函數(shù) 的 Taylor 展開，可以得到逼近誤差相對丁自變量差分（增量）的量級，稱為用差商 代替導(dǎo)數(shù)的精度，簡稱為差商的精度。應(yīng)用計算流體力學(xué)講稿 羅東明 4 展開: f (x wx) = f (x) wx f (x) ( x) f (x) 2! 將上式代入一階向前差商表達式中，有 f(x*x)-f(x)=f(x).g

6、.：x l(.：x)2 0() 2! 3! 這里符號 0()表示與括號中的量有相同的量級。上式表明一階向前差商的逼近誤差與 自變量的增量為同一量級。把 0(Axn)中 Ax 的指數(shù)作為精度的階數(shù)。這里 n=1,故一 階向前差商具有一階精度。由丁 Ax是個小量，因此階數(shù)越大精度越高。采用同樣的 辦法可知一階向后差商也具有一階精度。 對丁一階中心差商，將函數(shù)f(x+Ax)與f(x-Ax)在 x 的Ax鄰域作 Taylor 展開并 代入一階中心差商的表達式中，有 心+叫一冷一取匚山豚碩2) (2-3) 2* 可見一階中心差商具有二階精度。同樣，二階中心差商的精度也為二階。 2.2差分方程、截斷誤差和

7、相容性 從上節(jié)所述可知，差分相應(yīng)丁微分，差商相應(yīng)丁導(dǎo)數(shù)。只不過差分和差商是用有 限形式表示的，而微分和導(dǎo)數(shù)則是以極限形式表示的。如果將微分方程中的導(dǎo)數(shù)用相 圖2.2網(wǎng)格劃分 應(yīng)的差商近似代替，就可以得到有限形式的差分方程?，F(xiàn)以對流方程( 2-4)為例, 列出相對應(yīng)的差分方程。現(xiàn)以一階向前差商為例來分析其精度。將函數(shù) f(x十Ax)在 x 的 Ax 鄰域作 Taylor 峙 f (x) O( x)4) x =f (x) O(. ：應(yīng)用計算流體力學(xué)講稿 羅東明 5 X = x0 Lx, i = 0,1,2, tn = n, n = 0,1,2,. 劃分出矩形網(wǎng)格，如圖 2.2 所示。這里 Ax 和

8、 At 取常數(shù)。直線 t=tn稱為第 n 層，網(wǎng)格 交義點稱為結(jié)點。 網(wǎng)格點劃定后，就可針對某一結(jié)點，例如圖 2.2 中的結(jié)點(xi,tn),用差商近似代 替導(dǎo)數(shù)。現(xiàn)用()in表示括號內(nèi)函數(shù)在(xi,tn )點的值，則對流方程在該點為 (2-5) 如果時間導(dǎo)數(shù)用一階向前差商近似代替: -n n 1 n :u Ui -Ui ft i t 空間導(dǎo)數(shù)用一階中心差商近似代替： Z n n n 絲)七Uf Ui 頃x 2Ax 則對流方程在(xi,tn)點對應(yīng)的差分方程為 n T n n n U+0( Ui* UI =0 (2-6) ：t 2 x 按照前面關(guān)丁逼近誤差的分析知道，用時間向前差商代替時間導(dǎo)數(shù)

9、的誤差為 O(At),用空間中心差商代替空間導(dǎo)數(shù)時的誤差為 O(Ax)2),因而對流方程與對應(yīng)的 差分方程之間也存在一個誤差，這一誤差可由 Taylor 展開確定，即 U(X,tn +At)-U(xi -tn) + U(x +Ax,tn) - U(xi - Ax,tn) a t 2 x =蝕 1+偵空】 +ae+【空n(Ax)2+ (2-7) I 2 32Llcx .J 3 檔 x3 尸 ) .J / =T7 a O( t,( x)2) ：t :x i 這種用差分方程近似代替微分方程所引起的誤差，稱為截斷誤差。這里誤差量級相當(dāng).:u .:t 用差商近似代替導(dǎo)數(shù)時，首先要選定 上用平行丁坐標軸的

10、兩族直線： 3 n a = 0 ；:x Ax和At，稱為步長。然后在x -1坐標平面 (2-4) 應(yīng)用計算流體力學(xué)講稿 羅東明 6 丁 At的一次式、Ax的二次式 一個與時間相關(guān)的物理問題，應(yīng)用微分方程表示時，還必須給定初始條件，從而 形成一個完整的初值問題。對流方程的初值問題為 +a=0 i ct ex u(x,0) =U(x) (2-8) 這里U(x)為某已知函數(shù)。同樣，差分方程也必須有初始條件: in -1 n n n Ui 7 : uiuiJ _ 0 .t 0 2 x Ui =u(x) (2-9) 初始條件是一種定解條件，差分方程和其定解條件一起，稱為相應(yīng)微分方程定解問題 的差分格式。

11、將式(2-9)中第n+1時間層的量放在等號左邊，將其余時間層的量放 在等號右邊，有 n：；1 n -t . n n . i - ui a (ui 1 5 言) 2x =姑) (2-10) 稱其為 FTCSB 式(時間前差、空間中差) 0 ui 若時間和空間都用向前差分，貝 U 得 n n 山 1 一 ui ： L = 0 n 1 n 山 7 :t 璀 0 ui =u(xi) (2-11) 同樣，將第n+1時間層的量放在等號左邊，將其余時間層的量放在等號右邊，有 n -1 n & / n n、 % =q a了(山蟲一山 ) x 0 -z 、 ui = u(x) (2-12) 該格式稱為

12、FTFS式。若時 I 可米用向前差分、空間米用向后差分，則得到 FTBS 格式: n -1 n =t,n n . ui =7 -a (ui -uj) x u0 2(xi) (2-13) 觀察這三種差分格式，可以看出若知道第 n時間層的u,則可以由一個差分式子 直接算出第n+1時間層的u,稱這類格式為顯式格式。 差分方程的相容性：如果當(dāng) Ax、AtT。時，此差分方程的截斷誤差的某種范數(shù) 也趨近丁零，則表明從截斷誤差的角度來看，此差分方程是能用來逼近微分方程的， 通常稱這樣的差分方程和相應(yīng)的微分方程相容。 應(yīng)用計算流體力學(xué)講稿 羅東明 7 2.3收斂性和穩(wěn)定性 當(dāng)步長趨丁零時，要求差分格式的解趨丁

13、微分方程的解，稱這種是否趨丁微分方 程定解問題的解的情況為差分格式的收斂性。 在有限差分法的具體運算中，計算誤差總是不可避免的，如舍入誤差，以及這種 誤差的傳播、積累。如果這一誤差對以后的影響越來越小，或是這個誤差保持在某個 限度內(nèi)，那么就稱這個差分格式在給定的條件下穩(wěn)定。根據(jù)理論分析可以知道，上面 介紹的幾種差分格式是條件穩(wěn)定的。 2.4差分格式介紹 2.4.1迎風(fēng)格式 前面已經(jīng)指出，微分問題 u ;:u c a =0 ；:t ；x u(x,0) = u(x) 的 FTBSB 式，在a0和a竺芻的條件下穩(wěn)定，而 FTFS 格式在a0時的波 形傳播方向沿 x 軸正向，上游的量經(jīng)過一段時間要傳播

14、到下游， n+1時刻i站的量要 受到上游站 n 時刻量的影響，故只可以迎著風(fēng)向取空間差分，而不可以順著風(fēng)向取空 間差分。這種格式是迎著風(fēng)向往上游作差分所得到的，稱為迎風(fēng)格式。上述 FTBS 格 式和 FTFSa 式都必須在迎風(fēng)時有條件穩(wěn)定。 2.4.2 隱式格式 前面介紹的顯式格式往往是有條件穩(wěn)定的，甚至完全不穩(wěn)定。如 FTCS 是完全不 穩(wěn)定的，F(xiàn)TBS式是條件穩(wěn)定的。對丁 FTB 潞式，在a0和a里弓1的條件下穩(wěn)定， x x 即要求x,當(dāng)要求空間步長 8 很小時，時間步長也必須取的很小，才能保證格 a 式穩(wěn)定，而 At 取得小，計算工作量就大大增加，經(jīng)濟上也不合算。而本節(jié)將要介紹的 隱式格

15、式常常是無條件穩(wěn)定的，因此在許多情況下受到重視并被廣泛應(yīng)用。 隱式格式相當(dāng)丁從(x,tn+At)點出發(fā)，用時間的向后差分把第n+1時間層的量與 已知時間層的量聯(lián)系起來。現(xiàn)以對流方程為例，從(x”tn+At)點出發(fā)取 BTCSH 分可得應(yīng)用計算流體力學(xué)講稿 羅東明 8 n 1 n n -1 n -1 Ui - u, U, .1 - U, 1 c - - :- - = 0 t 2.x 或改寫為 a.、：t n 1 n 1 a t n 1 n Ui-Ui - Ui 1 = -Ui 2 ：x 2 x 由丁該方程含有三個第n十1時間層上的函數(shù)值，即一個方程含有三個未知量，必須解 聯(lián)立方程才能得到第n+1

16、時間層上的未知量，故稱該格式為隱式格式。 可以證明，用丁對流方程的隱式格式是完全穩(wěn)定的。由丁完全穩(wěn)定，時間步長可 以取得大些，從這一點來說，工作量減少了。 一時間步長內(nèi)工作量有所增加。 2.5耗散與色散 n n Ui Ui =0 x 利用 Taylor 展開,得到: 一 一 . .2 -U :u 1 U a=-a(1 - r)2 ft 汶 2 改 (3一22一1)皂 6 :x 上式就是差分方程(2-16)實際所模擬的微分方程，與原對流方程相比，多了二次導(dǎo) 數(shù)項和三次導(dǎo)數(shù)項。一般說：截斷誤差中含偶次導(dǎo)數(shù)項時，將引起耗散。在流體力學(xué) 方程中，二次導(dǎo)數(shù)項是與粘性項相關(guān)的。不同的是，在這里該項是差分方程數(shù)值離散 的結(jié)果，因此純粹的數(shù)值引發(fā)，沒有物理意義。因此，在 CFM 法中，類似的項被稱 為數(shù)值耗散。相似的，該項的系數(shù)，即；aAx(1-r),其作用類似丁物理粘性，因而被 稱為人工粘性。雖然數(shù)值耗散損害了計算精度，但卻改善了計算的穩(wěn)定性。因此，彳艮 多不穩(wěn)定的計算格式或顯式、或隱含地添加了人工粘性后，就變得穩(wěn)定了。同樣，在 截斷誤差中還含奇次導(dǎo)數(shù)項，該項將引起數(shù)值色散。 最后，總結(jié)一下有限差分法的優(yōu)缺