基本方程組的數(shù)值求解

基本方程組的數(shù)值求解第一頁，共四十頁，編輯于2023年，星期日時間導(dǎo)數(shù)項在控制容積上積分，得

（45）這里假定、在整個控制容積中是均勻的，網(wǎng)格不隨時間而變在時間導(dǎo)數(shù)項的差分式子中，包含了兩個時間層（在類時間坐標(biāo)上是兩個截面）的參數(shù)值與此相對應(yīng)，其它坐標(biāo)方向的對流項、擴(kuò)散項以及源項中的參數(shù)，在時間層上可以有不同的取法，因而可以得到顯式、隱式和克蘭克-尼科爾森（Crank-Nicholson）格式第二頁，共四十頁，編輯于2023年，星期日1）顯式差分格式

（46）2）隱式差分格式

（47）3）克蘭克-尼科爾森格式若對流項、擴(kuò)散項和源項中的參數(shù)都取兩個時間層上的算術(shù)平均值，可得差分方程

（48）克蘭克-尼科爾森格式的差分方程也需要與其它結(jié)點的差分方程聯(lián)立求解。當(dāng)時間步長較小時，其精度比隱式格式高。第三頁，共四十頁，編輯于2023年，星期日4、源項的差分

設(shè)源項在控制容積內(nèi)均勻，并等于中心點的值，則有

（49）微分方程離散化后得到的是一個擬線性代數(shù)方程，但源項往往是參數(shù)的非線性函數(shù)為能夠用解線性代數(shù)方程的方法求解差分方程，需要將源項線性化源項線性化有兩種辦法:一是用前一時刻的參數(shù)代入源項，算得的值作為常數(shù)參加下一時刻（或下一循環(huán)）的計算；二是假定源項與參數(shù)的函數(shù)關(guān)系可以近似用式子

（50）表示

第四頁，共四十頁，編輯于2023年，星期日代入差分方程得

（51）其中Sc和Sp可能也是的函數(shù)，因此式（50）也是擬線性關(guān)系。第二種辦法優(yōu)于第一種辦法為了使擬線性代數(shù)方程組收斂，得到有意義的解，Sc和Sp必須滿足, （52）第五頁，共四十頁，編輯于2023年，星期日在滿足條件（52）的基礎(chǔ)上，和的確定有一定的任意性。一個方便的辦法是用切線使其線性化，即令

（53）這時，按照切線法得到的和往往很復(fù)雜，而且不一定滿足條件（52）比較多的是用經(jīng)驗辦法來確定第六頁，共四十頁，編輯于2023年，星期日五、代表方程組的求解1、概述采用有限容積法在已生成的網(wǎng)格上將所求解變量的控制方程離散后，就形成了各求解變量的代數(shù)方程組求解離散化所得的代數(shù)方程組是流動與燃燒過程數(shù)值計算的最后一個環(huán)節(jié)求解的方法有兩大類：直接解法與迭代法迭代法的研究主要是研究收斂性以及如何加快收斂速度求解代數(shù)方程的解法包括Gauss消元法、三對角陣（TDMA）算法、五對角陣（PDMA）算法等現(xiàn)有文獻(xiàn)中求解由QUICK等高階格式形成的代數(shù)方程組時多采用TDMA方法求解，TDMA算法是Gauss消元法的一種特例（每一行上僅三個非零元素）第七頁，共四十頁，編輯于2023年，星期日2.Gauss消元法3.TDMA算法對圖1中S-N方向的網(wǎng)格線進(jìn)行計算時，差分方程需重新整理為

（56）方程右端各項是鄰近網(wǎng)格線上的結(jié)點值（或源項），取前一循環(huán)得到的值，故為已知值每個節(jié)點的代數(shù)方程中最多只包含三個節(jié)點的未知值，可以認(rèn)為其它節(jié)點上未知值的系數(shù)均為零如果把上述有限差分離散方程組寫成矩陣的形式，其系數(shù)陣是一個三對角陣—僅對角元素及其上下鄰位上的元素不為零，而其它元素均為零把Gauss消元法應(yīng)用于這種情形，便構(gòu)成了稱為三對角陣算法的有效求解方法，簡記為TDMA（TridiagonalMatrixAlgorithm）第八頁，共四十頁，編輯于2023年，星期日（57）第九頁，共四十頁，編輯于2023年，星期日對于這種系數(shù)矩陣為三對角線矩陣的方程組，TDMA求解的具體步驟是：（1）第一個方程的各項除以B1，得，，。第二個方程減去化簡后的第一個方程乘以A2，并將所得新方程的各項除以（），最后得對第三個以后的每個方程都作同樣的處理，得到方程組系數(shù)的通式為：

（58）第十頁，共四十頁，編輯于2023年，星期日（2）自下而上的解方程……解的通式為

（59）若有兩個或三個坐標(biāo)是橢圓型的，則需對各橢圓型方向都進(jìn)行逐線迭代掃描對兩個或三個方向各作了一次逐線計算，稱為進(jìn)行了一次雙重或三重掃描，也叫一個迭代循環(huán)重復(fù)多重掃描，直至兩次相繼迭代得出的值差別不大為止。把最后結(jié)果作為n時間層的值對每條網(wǎng)格線進(jìn)行計算時，邊界條件很快傳入流場內(nèi)部，因此收斂速率比逐點迭代快得多。隱式差分方程的逐線迭代解法對時間步長也無限制，因此在實際中得到廣泛的應(yīng)用第十一頁，共四十頁，編輯于2023年，星期日六、壓力和速度之間的耦合1、引言若知道流場中壓力的分布，則動量方程的求解與其它方程完全一樣，不會產(chǎn)生困難但是壓力值一般都是預(yù)先不知道的在可壓流中，壓力與密度間的關(guān)系由狀態(tài)方程確定，需通過連續(xù)方程和狀態(tài)方程確定在不可壓流中，流場中的壓力分布對速度場有很強(qiáng)的影響，壓力梯度以源項形式出現(xiàn)在動量方程中，但壓力卻沒有獨立的控制方程在計算流體力學(xué)的發(fā)展過程中，提出了多種不同的方法，來解決在以速度、壓力為求解變量的原始變量法中的這個問題大致可分為：渦量-流函數(shù)方法、壓力修正方法、壓力與速度之間的迭代算法等幾種。以壓力修正方法中的SIMPLE系列方法應(yīng)用最為廣泛第十二頁，共四十頁，編輯于2023年，星期日2．壓力修正方法的基本思想在求解不可壓縮流體的流場問題時，如果我們把從動量方程與連續(xù)性方程離散得到的代數(shù)方程組聯(lián)立起來直接求解，就可以得到各速度分量及相應(yīng)的壓力值這樣的直接解法要占用大量的計算機(jī)內(nèi)存，對于目前大多數(shù)工程應(yīng)用場合還不適用如果采用分離式的迭代求解方法（segregatedmethod），即先求解u速度場，再求解v速度場，則對壓力場因其無獨立的方程而無法對其求解（或無法改進(jìn)其原先的假定值）另一方面，上述分離式求解過程中只利用了u,v動量方程的離散形式而未用到連續(xù)性方程要解決的問題:如何使得用質(zhì)量守恒方程假定的壓力場能不斷地隨迭代過程的進(jìn)行而得到改進(jìn)，這就是所謂的壓力修正算法壓力修正算法源于1972年由Patankar與Spalding提出的SIMPLE算法[16]，在20世紀(jì)80年代初期，又相繼提出了SIMPLER與SIMPLEC等方法，并由此形成了SIMPLE系列算法第十三頁，共四十頁，編輯于2023年，星期日3．SIMPLE算法以二維直角坐標(biāo)中的對流換熱求解為例，其控制方程為

（60）

（61）

（62）在交錯網(wǎng)格（圖3）上動量離散方程為：

（63）

（64）質(zhì)量守恒方程為

（65）第十四頁，共四十頁，編輯于2023年，星期日采用SIMPLE算法來求解時計算步驟如下：（1）假定一個速度場，記為u(0)，v(0)，由此計算ae，an，anb及b；（2）假定一個壓力場，記為；（3）求解動量離散方程（63）、（64），得；（4）計算壓力的修正值，要求與相對應(yīng)的、能滿足連續(xù)性方程第十五頁，共四十頁，編輯于2023年，星期日（5）計算速度的修正值。要求仍滿足線性化了的動量方程，即與式（63）相減，得

（69）此式表明，要據(jù)值確定，需要解一個代數(shù)方程組。為能利用值顯式地求解，此處略去式（69）右端第一項，，于是得

（70）（6）將、作為本迭代層次之解，開始下一層次的計算，重復(fù)1~5步直到流場收斂，即所解得之流速能同時滿足連續(xù)性方程及動量方程為止第十六頁，共四十頁，編輯于2023年，星期日SIMPLE算法中引入了以下假定或簡化處理：（1）速度場的假定（u（0），v（0））與壓力場的假定（）是相互獨立地進(jìn)行的，u(0)，v(0)與間沒有任何聯(lián)系；（2）在導(dǎo)出速度修正值計算式（70）時，未計及鄰點速度修正值的影響；（3）動量離散方程中的b(控制容積P的剩余質(zhì)量流量)在速度修正前后保持不變；（4）由式（63）、（64）解出的滿足動量守恒但未必滿足質(zhì)量守恒，而由式（69）解出的決定時，保證、滿足質(zhì)量守恒，但動量守恒則未必滿足。注意，雖然在導(dǎo)出式（69）過程中曾要求、滿足線性化了的動量方程，但由于以后略去了式（69）中的，因而所得的，就未必可使、滿足上述動量守恒方程。以上四條假設(shè)或近似處理，是SIMPLE算法提出之后所出現(xiàn)的一些改進(jìn)方案的著眼之處。第十七頁，共四十頁，編輯于2023年，星期日4、SIMPLER算法

SIMPLER（1980）算法主要用以改進(jìn)SIMPLE方法中的第一項近似處理方法。一旦速度場給定，壓力場就可以從動量離散方程中予以求解，而不再任意假定，其主要計算步驟如下：（1）假定一個速度場，記為u(0)，v(0)，由此計算動量離散方程系數(shù)ae，an，anb，b及，：

（71）注意，引入，后，動量離散方程便可寫為

（72）第十八頁，共四十頁，編輯于2023年，星期日（2）據(jù)計算相應(yīng)的壓力場。為此，將式（72）代入質(zhì)量守恒方程的離散形式（65），得

（73）其中，anb及b的計算式形式上與SIMPLE算法中的方程的一樣，只要將代替即可。（3）求解動量離散方程，得替。（4）求解壓力修正方程，得。（5）用修正速度得（但不修正壓力）。（6）用，開始下一層次的迭代，重復(fù)1~5步直到收斂。為有利于非線性迭代的收斂，SIMPLER算法中流速應(yīng)予以亞松弛，但及p則不作亞松弛。第十九頁，共四十頁，編輯于2023年，星期日5、SIMPLEC算法6．SIMPLEX算法7、SIMPLET算法第二十頁，共四十頁，編輯于2023年，星期日8．SIMPLE系列算法小結(jié)SIMPLE，SIMPLEC，SIMPLER，SIMPLET及SIMPLEX的主要計算步驟可以歸納為：（1）由假定的或上一次計算得到的速度u(0)，v(0)確定動量離散方程的系數(shù)及常數(shù)項；（2）確定壓力場：對于SIMPLE，SIMPLEC，SIMPLET及SIMPLEX采用假定的或用上一次的計算值；對于SIMPLER則由已知的速度場通過求解壓力Poisson方程而獲得；（3）求解動量離散方程以獲得；（4）求解壓力修正值方程，要求與相應(yīng)的能使、滿足連續(xù)性方程。其中除SIMPLET外，其它算法只考慮壓力修正值對的影響，而SIMPLET中還考慮了溫度修正對的影響，該影響最后體現(xiàn)在方程的源項之中；第二十一頁，共四十頁，編輯于2023年，星期日（5）確定速度修正值：SIMPLE，SIMPLER，SIMPLET：，鄰點速度影響略去不計；SIMPLEC：，鄰點速度修正值的影響已基本考慮在內(nèi)；SIMPLEX：，de等是通過求解一個Poisson方程得出的，相當(dāng)于考慮了鄰點的影響。（6）采用下列速度、壓力值作為初值開始下一層次的迭代計算（設(shè)速度的亞松弛已組合在u，v的求解過程中）：SIMPLE，SIMPLET：SIMPLEC，SIMPLEX：SIMLER：以上計算過程可以大致分為兩步：預(yù)估步（1）~（3）及校正步（4）~（6）。因而可以說上述五種算法都是兩步算法（一步預(yù)估，一步校正）。第二十二頁，共四十頁，編輯于2023年，星期日9．加速迭代收斂的一些方法1）選擇合適的松馳因子之值對于正交的網(wǎng)格以及雖非正交但網(wǎng)格傾斜不是嚴(yán)重的情況，建議au與ap的取值應(yīng)滿足以下關(guān)系：

（76）式中常數(shù)c取為1[23]或1.1[24]，同時au之值應(yīng)盡可能地取得大，一般可取0.7~0.8。上述取值原則對交錯網(wǎng)格、同位網(wǎng)格、正交網(wǎng)格與非正交網(wǎng)格原則上都適用，但對網(wǎng)格嚴(yán)重傾斜的情況，ap之值應(yīng)小于按式（76）得出之值。第二十三頁，共四十頁，編輯于2023年，星期日2）其它加速迭代收斂的方法①迭代循環(huán)方式對收斂速率的影響描述燃燒過程的是一組互相耦合的方程組，在用迭代法求解時，迭代循環(huán)方式對收斂速率有很大的影響可將流體力學(xué)方程（動量和修正壓力方程）和其它標(biāo)量方程（能量、湍流量方程等）一起迭代另一種選擇是將流體力學(xué)方程首先迭代若干次，然后把其它標(biāo)量方程（能量、湍流量方程等）一起迭代；流體力學(xué)方程首先迭代若干次，然后把其它標(biāo)量方程再加進(jìn)去一起循環(huán)兩者的收斂速率是不一樣的。但目前還給不出最佳循環(huán)方式的一般原則②網(wǎng)格均勻度對收斂速率的影響網(wǎng)格均勻度對收斂速率有很大的影響計算時間與結(jié)點數(shù)成正比但增加結(jié)點數(shù)，若能使差分網(wǎng)格更均勻，由于收斂速率加快，有時仍有可能減少總的計算時間第二十四頁，共四十頁，編輯于2023年，星期日七、邊界條件的數(shù)值處理方法1、引言

初始條件和邊界條件:

對拋物型問題，只需給出一端的邊界條件，如時間坐標(biāo)上給出初始的參數(shù)分布，邊界層流動中給定進(jìn)口截面上的參數(shù)分布等等對橢圓型問題，則需給出求解域兩端的邊界條件，如果三個空間坐標(biāo)都是橢圓型坐標(biāo)，則求解域的四周各點都需給出邊界條件第二十五頁，共四十頁，編輯于2023年，星期日2、進(jìn)口邊界模型計算湍流時，一般無法獲得由實驗測得的進(jìn)口處k及k:（1）

（2）

:(1)

(2)

進(jìn)口截面上k與取值的這種不確定性對計算結(jié)果的影響（亦即計算結(jié)果對進(jìn)口條件的敏感性）應(yīng)當(dāng)在數(shù)值計算中進(jìn)行考察。第二十六頁，共四十頁，編輯于2023年，星期日3、出口邊界描寫有回流的流動與燃燒的控制方程是關(guān)于空間坐標(biāo)的橢圓型方程，從數(shù)學(xué)上要求在每個坐標(biāo)方向上都應(yīng)當(dāng)有兩個邊界條件。出口邊界應(yīng)該置在什么位置及怎樣確定這些人為設(shè)定的計算邊界條件，就成了計算流體力學(xué)及燃燒過程模擬中的一個重要研究課題。出口邊界應(yīng)設(shè)置在什么位置上:Patankar曾認(rèn)為，出口邊界不能設(shè)置在有回流的截面上，否則得到的解是沒有意義的數(shù)值實踐證明，對出口邊界位置選擇的這種苛刻的要求是不必要的文獻(xiàn)[27]提出對這種出口邊界條件定性的要求是：（1）應(yīng)當(dāng)對出口邊界附近，尤其是離開出口邊界較遠(yuǎn)處的流場沒有影響；（2）出口邊界應(yīng)當(dāng)是“透明的”,又稱為“無反射”（3）無論出口邊界設(shè)置在何處，應(yīng)對其它區(qū)域中的數(shù)值解沒有影響第二十七頁，共四十頁，編輯于2023年，星期日圖11 求解流場時出口邊界位置的確定第二十八頁，共四十頁，編輯于2023年，星期日目前，文獻(xiàn)中已提出了多種關(guān)于設(shè)置出口邊界條件的方法，歸納起來大約有以下幾種情形：1）局部單向化假設(shè)只有當(dāng)出口截面無回流時才能近似地成立2）充分發(fā)展的假設(shè)假設(shè)n為出口截面的法線方向，這一條件的數(shù)學(xué)表達(dá)式為

（84）3）法向速度局部質(zhì)量守恒、切向速度齊次Neumann條件以圖11（b）中出口邊界相鄰的控制容積為例，取出后放大而示于圖12中。在本例中速度u是與出口截面平行的流速，而v為與出口截面垂直的流速。對所示控制容積作質(zhì)量守恒運算，得

（85a）而按齊次Neumann條件有

（85b）第二十九頁，共四十頁，編輯于2023年，星期日圖12局部質(zhì)量守恒的圖示圖13說明切向速度按局部質(zhì)量守恒的計算圖示第三十頁，共四十頁，編輯于2023年，星期日4）出口截面上的法向速度與切向速度都按局部質(zhì)量守恒確定在這一方法中，對圖12所示情形，計算式同式（85a），但ui,M的更新則通過對圖13所示半個控制容積的質(zhì)量守恒來獲得。以0.5(vi,M+vi,M-1)作為該半個控制容積南側(cè)流速，則可以寫出

（86a）

將邊界上的切向流速取本迭代層次之值而其余速度取上一迭代層次的值，得

（86b）由圖11（b）可見，在計算區(qū)域右邊界，之值為零，由此可推得，這樣依次計算可以得出出口邊界切向速度的更新值。第三十一頁，共四十頁，編輯于2023年，星期日5）對流邊界條件用以下公式來更新出口截面法向速度及其它標(biāo)量：

（87a）其中，c為對流速度或相速度，一般取為通道平均流速。上式離散后，可得

（87b）這里，vm為出口截面平均法向流速。相當(dāng)于假定在出口邊界上無擴(kuò)散作用,為一個純對流問題。第三十二頁，共四十頁，編輯于2023年，星期日6）內(nèi)點插值方法

在文獻(xiàn)[36]中，提出采用由內(nèi)點的速度來決定出口流速的方法：

（88a）其中，為與出口邊界相鄰的控制容積的寬度，x為出口邊界上游第四個截面位置起算的坐標(biāo)，當(dāng)網(wǎng)格均分時（如圖14所示），有

（88b）第三十三頁，共四十頁，編輯于2023年，星期日圖14插值公式（88a）的圖示第三十四頁，共四十頁，編輯于2023年，星期日7）各種高階導(dǎo)數(shù)為零的方法

（89）其中，x為出口方向，上角標(biāo)j及q表示求導(dǎo)的階數(shù)，要求，并建議取j=3，q=2。這一類方法實際上相當(dāng)于提出了由內(nèi)點來插值的一個方式。導(dǎo)數(shù)的階數(shù)越高，涉及的內(nèi)點數(shù)越多，該階導(dǎo)數(shù)取零對解的影響也越微弱。第三十五頁，共四十頁，編輯于2023年，星期日出口邊界上邊界條件的設(shè)定有多種看來相差很大的方法。在從事工程流動與傳熱問題的數(shù)值計算中實用的處理出口邊界條件的方法：出口邊界上無回流時，可采用局部單向化或法向一階導(dǎo)數(shù)為零的方法；當(dāng)出口邊界有回流時，可以采用法向流速按質(zhì)量守恒、切向流速按齊次Neumann條件的方法；無論采用哪一種方法，出口截面上的法向速度分布必須滿足總體質(zhì)量守恒的條件。第三十六頁，共四十頁，編輯于2023年，星期日4．固體邊界

1）對于溫度一類的變量用附加源項法來處理2）對于速度固體壁面是受限流動中最常見的邊界，原則上應(yīng)采用無滑移條件，即流體的速度等于壁面的速度（第一類邊界條件），層流時無論是對于動量方程中的流速還是連續(xù)性方程中的流速都是這樣（后者表現(xiàn)在SIMPLE算法壓力修正方程中的剩余質(zhì)量流量的計算上）當(dāng)采用標(biāo)準(zhǔn)模型時，平行于壁面的流速及垂直于壁面的流速還存在不同處理方法。對于平行于壁面的流速（例如u），采用無滑移條件，同時壁面上的動力粘性系數(shù)由壁面函數(shù)法的公式來計算；但對垂直于該壁面的流速（v）則應(yīng)區(qū)分動量方程還是連續(xù)方程。在連續(xù)性方程中仍然采用無滑移條件（