計(jì)算流體力學(xué)基礎(chǔ)簡(jiǎn)介

計(jì)算流體力學(xué)基礎(chǔ)簡(jiǎn)介本碩111班王鵬主要內(nèi)容：

流體力學(xué)的控制方程組適合CFD使用的控制方程組偏微分方程的離散化——有限差分計(jì)算流體力學(xué)的基本方法流體力學(xué)的控制方程組一.粘性流動(dòng)的N-S方程

1.連續(xù)性方程

非守恒形式守恒形式

2.動(dòng)量方程

非守恒形式

流體力學(xué)的控制方程組守恒形式3.能量方程非守恒形式

流體力學(xué)的控制方程組

守恒形式

流體力學(xué)的控制方程組一.無粘流Euler方程

1.連續(xù)性方程

非守恒形式守恒形式

2.動(dòng)量方程

非守恒形式

流體力學(xué)的控制方程組守恒形式3.能量方程非守恒形式守恒形式

流體力學(xué)的控制方程組對(duì)于牛頓流體，斯托克斯在1845年得到，在上面的方程組中：

根據(jù)傅里葉熱傳導(dǎo)定律，熱傳導(dǎo)產(chǎn)生的熱流與當(dāng)?shù)氐臏囟忍荻瘸烧?其中k為熱導(dǎo)率。適合CFD使用的控制方程組

守恒形式的連續(xù)性方程、動(dòng)量方程和能量方程可以使用同一個(gè)通用方程來表達(dá)，這有助于計(jì)算程序的簡(jiǎn)化和結(jié)構(gòu)組織，為算法設(shè)計(jì)和編程提供了方便。

守恒形式的通用方程

（1）

如果將上式中的U、F、G、H和J看成列向量，則上式就可代表整個(gè)守恒形式的控制方程。這些列向量為：

適合CFD使用的控制方程組

（1）式中的F、G和H稱為通量項(xiàng)(通量向量），J代表源項(xiàng)，U稱為解向量。對(duì)于一般的時(shí)間推進(jìn)方法，我們可將（1）式變換一下以方便求解

適合CFD使用的控制方程組

針對(duì)Euler方程也可進(jìn)行類似分析。偏微分方程的離散化——有限差分

偏微分方程的解析解是封閉形式的表達(dá)式，它給出了未知函數(shù)在區(qū)域內(nèi)的連續(xù)變化。偏微分方程的離散化是用另外一個(gè)類似的表達(dá)式來近似這些偏微分方程，但是這個(gè)近似表達(dá)式只在去域內(nèi)有限多個(gè)離散點(diǎn)或控制體上規(guī)定了取值。

利用泰勒級(jí)數(shù)展開可以推導(dǎo)出各階導(dǎo)數(shù)的有限差分的一般形式。（2）（3）利用（2）、（3）式可以得到以下一些差分公式：偏微分方程的離散化——有限差分一階導(dǎo)數(shù)的一階向前差分一階導(dǎo)數(shù)的二階中心差分一階導(dǎo)數(shù)的一階向后差分二階導(dǎo)數(shù)的二階中心差分一階混合導(dǎo)數(shù)的二階中心差分……..偏微分方程的離散化——有限差分利用多項(xiàng)式也可以實(shí)現(xiàn)有限差分。例如，如右圖所示的邊界上，u可以表示成多項(xiàng)式：

（4）（5）（6）（7）偏微分方程的離散化——有限差分解出b,再對(duì)（4）式求導(dǎo)，邊界上y=0,代入上式有，最終得到上式就是單側(cè)差分表達(dá)式，之所以被稱為單側(cè)，是因?yàn)樗挥玫搅诉吔缫粋?cè)的網(wǎng)格點(diǎn)上的信息。計(jì)算流體力學(xué)的基本方法

計(jì)算流體力學(xué)包括很多方法，其中較基礎(chǔ)的有拉克斯-溫德羅夫方法、麥考馬克方法、空間推進(jìn)方法、松弛法、交替方向隱式方法、壓力修正法等。1.拉克斯-溫德羅夫方法

以非定常二維無粘流為例，控制方程如下

（8）（9）（10）（11）計(jì)算流體力學(xué)的基本方法

（12）（13）（14）（15）

計(jì)算流體力學(xué)的基本方法

（17）

再將方程（8）對(duì)時(shí)間求導(dǎo)：

上式中出現(xiàn)二階混合導(dǎo)數(shù)，可由方程式（8）到（11）對(duì)相應(yīng)的空間自變量求導(dǎo)得到。例如（16）

將右邊所有的項(xiàng)還是由t時(shí)刻的二階中心差分表示，即計(jì)算流體力學(xué)的基本方法

計(jì)算流體力學(xué)的基本方法

二.麥考馬克（Macormack）方法

麥考馬克方法是從拉克斯-溫德羅夫方法中變化過來的一種方法，但更簡(jiǎn)單，在1969年一經(jīng)提出，就成為解決流動(dòng)問題最流行的顯示有限差分方法，一直流行了15年，現(xiàn)在它已被更先進(jìn)的方法取代了。

（18）計(jì)算流體力學(xué)的基本方法

預(yù)估步：在連續(xù)方程（8）中，用向前差分代替方程右邊的空間導(dǎo)數(shù)：

（19）同理可求出u、v和e的預(yù)估值，即計(jì)算流體力學(xué)的基本方法（20）

式（18）中密度的時(shí)間導(dǎo)數(shù)的平均值

計(jì)算流體力學(xué)的基本方法（21）

松弛法是一種適合于求解橢圓型偏微分方程的迭代法。例如，亞聲速無粘流動(dòng)由橢圓型偏微分方程控制松弛法經(jīng)常被用來求解亞聲速的低速流動(dòng)。松弛法可以是顯式的也可以是隱式的，這里介紹顯式松弛法。

我們將在右圖所示的網(wǎng)格上數(shù)值求解方程（21），利用二階導(dǎo)數(shù)的二階中心差分式代替上式的偏導(dǎo)數(shù)，得

上式中，上標(biāo)n和n+1表示迭代次數(shù)。（22）（23）

三.松弛法計(jì)算流體力學(xué)的基本方法（24）

在對(duì)所有內(nèi)部網(wǎng)格點(diǎn)都使用了方程（23）后，就完成了第一次迭代，n=1。接著進(jìn)行第二次迭代，n=2。這個(gè)迭代過程反復(fù)不斷，直至收斂到解。（為什么一定是朝解的方向收斂？）

再具體些，將方程（23）應(yīng)用到網(wǎng)格點(diǎn)21，進(jìn)行第n+1次迭代，有

（25）計(jì)算流體力學(xué)的基本方法

計(jì)算流體力學(xué)的基本方法（26）

四.交替方向隱式（ADI）方法講這個(gè)方法之前，有必要先介紹一下顯示方法和隱式方法。所有的有限差分方法都可以歸結(jié)為這兩種不同的通用方法——顯示方法或隱式方法。顯示方法。以一位熱傳導(dǎo)方程為例等式左邊用向前差分，等是右邊用二階中心差分，即整理此式后可寫成方程（26）是拋物線型偏微分方程，可推進(jìn)求解。此時(shí)，推進(jìn)變量是時(shí)間t，為更簡(jiǎn)明起見，考慮右圖所示的有限差分網(wǎng)格。

假設(shè)在第n個(gè)時(shí)間層上的每個(gè)網(wǎng)格點(diǎn)上溫度T是已知的，時(shí)間推進(jìn)就意味著第n+1個(gè)時(shí)間層每個(gè)網(wǎng)格點(diǎn)上的T值都要用第n個(gè)時(shí)間層上的已知量計(jì)算出來。（28）（27）計(jì)算流體力學(xué)的基本方法（29）

例如，在網(wǎng)格點(diǎn)2處，式（28）將寫成

隱示方法

我們可以再大膽些，把式（27）右邊的空間差分寫成第n個(gè)時(shí)間層的量與第n+1個(gè)時(shí)間層上的量的平均值，即

計(jì)算流體力學(xué)的基本方法（30）

考慮右圖所示的有限差分網(wǎng)格，以七點(diǎn)空間網(wǎng)格點(diǎn)作為具體例子。整理式（29），將未知量放在等式左邊，已知量放到等式右邊，得

于是，式（30）可寫成

（31）計(jì)算流體力學(xué)的基本方法

把方程（31）~（35）寫成矩陣形式：（32）（33）（34）（35）

這是一個(gè)三對(duì)角矩陣，可用托馬斯算法（國(guó)內(nèi)稱為追趕法）求解。解出之后即得到第n+1個(gè)時(shí)間層上的溫度T值。重復(fù)以上過程，即可得到最終的數(shù)值解。計(jì)算流體力學(xué)的基本方法（37）

現(xiàn)在，介紹交替方向隱式（AlternatingDirectionImplicit,ADI）方法。以二維熱傳導(dǎo)方程為例，

對(duì)上式運(yùn)用克蘭克-尼克爾森方法，有

（36）（38）計(jì)算流體力學(xué)的基本方法（39）

方程（38）可化為三對(duì)角形式

式中

計(jì)算流體力學(xué)的基本方法（40）

方程（40）可化為三對(duì)角形式

式中

（41）計(jì)算流體力學(xué)的基本方法（46）五.壓力修正法

壓力修正法本質(zhì)上也是一種迭代法，思路如下：

1）迭代開始時(shí)，先給定流場(chǎng)壓力的初始近似值p*。

2）用p*的值從動(dòng)量方程中求解出u,v,w.因?yàn)檫@些速度都與p*有關(guān)，所以用u*,v*,w*表示它們。

3）

u*,v*,w*不一定滿足連續(xù)方程。因此，要用連續(xù)方程構(gòu)造壓力的修正量p’，加到p*上，使速度滿足連續(xù)方程。設(shè)修正后的壓力為

相應(yīng)的速度修正量u’,v’,w’可以從p’得到，使得

壓力修正法的基本原理（42）（43）（44）（45）計(jì)算流體力學(xué)的基本方法（47）

此時(shí)我們要做到就是求出壓力修正量p’，為簡(jiǎn)單起見，考慮二維流動(dòng)，并忽略體積力。對(duì)不可壓流動(dòng)，x方向和y方向的動(dòng)量方程分別是式（43）和（44），這些方程都是非守恒型的。若寫成守恒形式就是4）用方程（46）左邊的p作為新的p*，回到步驟2）。重復(fù)這個(gè)過程，直到速度場(chǎng)滿足連續(xù)方程為止。壓力修正公式（48）計(jì)算流體力學(xué)的基本方法（49）

考慮如右圖所示的交錯(cuò)網(wǎng)格，壓力在實(shí)心圓點(diǎn)上計(jì)算，速度在空心圓點(diǎn)上計(jì)算。在右圖中的點(diǎn)（i+1/2,j)點(diǎn)對(duì)x方向動(dòng)量方程（47）進(jìn)行中心差分（右圖中陰影部分相當(dāng)于有限體積法中的控制體）。中心差分要用到陰影區(qū)上邊a點(diǎn)和下邊b點(diǎn)處的v的平均值，我們用相鄰兩點(diǎn)的線性插值來定義這些平均值，即

以（i+1/2,j)點(diǎn)為中心，方程（47）的差分方程為

即式中計(jì)算流體力學(xué)的基本方法（50）

y方向動(dòng)量方程（48）的中心差分方程也可以用相同的方式得到。中心差分要用到陰影區(qū)左邊的c點(diǎn)和右邊的d點(diǎn)處u的平均值，定義為

以（i,j+1/2)點(diǎn)為中心，對(duì)方程（48）進(jìn)行差分，有式中計(jì)算流體力學(xué)的基本方法（51）根據(jù)前面給出的步驟，在迭代開始時(shí)，給定p=p*。此時(shí)，式（49）和（50）分別為（52）從式（49）中減去式（51），得式中從式（50）中減去式（52），得（53）（54）計(jì)算流體力學(xué)的基本方法（55）（56）式中

式（53）和式（54），就是用壓力修正量p’和速度修正量u’,v’表示的動(dòng)量方程。

則，式（55）可表達(dá)為計(jì)算流體力學(xué)的基本方法（58）（57）

則，式（56）可表達(dá)為

回到二維連續(xù)方程

在網(wǎng)格點(diǎn)（i,j)處寫出相應(yīng)的差分方程，即（59）

去掉式（57）（58）的上標(biāo)，然后代入式（59），得到

整理上式，得到計(jì)算流體力學(xué)的基本方法（60）式中方程（60）就是壓力修正公式，它具有橢圓性質(zhì)。因此用前面講過的松弛法