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

1、基本方程組的數(shù)值求解第1頁，共40頁，2022年，5月20日，12點29分，星期二時間導數(shù)項 在控制容積 上積分，得（45）這里假定 、 在整個控制容積中是均勻的，網(wǎng)格不隨時間而變在時間導數(shù)項的差分式子中，包含了兩個時間層（在類時間坐標上是兩個截面）的參數(shù)值與此相對應，其它坐標方向的對流項、擴散項以及源項中的參數(shù)，在時間層上可以有不同的取法，因而可以得到顯式、隱式和克蘭克-尼科爾森（Crank-Nicholson）格式第2頁，共40頁，2022年，5月20日，12點29分，星期二1）顯式差分格式（46）2）隱式差分格式（47）3）克蘭克-尼科爾森格式 若對流項、擴散項和源項中的參數(shù)都取兩個時間

2、層上的算術(shù)平均值，可得差分方程（48） 克蘭克-尼科爾森格式的差分方程也需要與其它結(jié)點的差分方程聯(lián)立求解。當時間步長較小時，其精度比隱式格式高。第3頁，共40頁，2022年，5月20日，12點29分，星期二4、源項的差分 設(shè)源項在控制容積內(nèi)均勻，并等于中心點的值，則有（49）微分方程離散化后得到的是一個擬線性代數(shù)方程，但源項往往是參數(shù)的非線性函數(shù)為能夠用解線性代數(shù)方程的方法求解差分方程，需要將源項線性化源項線性化有兩種辦法: 一是用前一時刻的參數(shù)代入源項，算得的值作為常數(shù)參加下一時刻（或下一循環(huán)）的計算；二是假定源項與參數(shù) 的函數(shù)關(guān)系可以近似用式子（50） 表示 第4頁，共40頁，2022年，

3、5月20日，12點29分，星期二 代入差分方程得 （51） 其中Sc和Sp可能也是 的函數(shù)，因此式（50）也是擬線性關(guān)系。第二種辦法優(yōu)于第一種辦法為了使擬線性代數(shù)方程組收斂，得到有意義的解，Sc和Sp必須滿足 , （52）第5頁，共40頁，2022年，5月20日，12點29分，星期二 在滿足條件（52）的基礎(chǔ)上， 和 的確定有一定的任意性。一個方便的辦法是用切線使其線性化，即令 （53） 這時 ，按照切線法得到的 和 往往很復雜，而且不一定滿足條件（52）比較多的是用經(jīng)驗辦法來確定第6頁，共40頁，2022年，5月20日，12點29分，星期二五、代表方程組的求解1、概述采用有限容積法在已生成的

4、網(wǎng)格上將所求解變量的控制方程離散后，就形成了各求解變量的代數(shù)方程組求解離散化所得的代數(shù)方程組是流動與燃燒過程數(shù)值計算的最后一個環(huán)節(jié)求解的方法有兩大類：直接解法與迭代法迭代法的研究主要是研究收斂性以及如何加快收斂速度求解代數(shù)方程的解法包括Gauss消元法、三對角陣（TDMA）算法、五對角陣（PDMA）算法等現(xiàn)有文獻中求解由QUICK等高階格式形成的代數(shù)方程組時多采用TDMA方法求解，TDMA算法是Gauss消元法的一種特例（每一行上僅三個非零元素）第7頁，共40頁，2022年，5月20日，12點29分，星期二2. Gauss消元法3. TDMA算法 對圖1中S-N方向的網(wǎng)格線進行計算時，差分方程

5、需重新整理為（56）方程右端各項是鄰近網(wǎng)格線上的結(jié)點值（或源項），取前一循環(huán)得到的值，故為已知值每個節(jié)點的代數(shù)方程中最多只包含三個節(jié)點的未知值，可以認為其它節(jié)點上未知值的系數(shù)均為零如果把上述有限差分離散方程組寫成矩陣的形式，其系數(shù)陣是一個三對角陣 僅對角元素及其上下鄰位上的元素不為零，而其它元素均為零把Gauss消元法應用于這種情形，便構(gòu)成了稱為三對角陣算法的有效求解方法，簡記為TDMA（Tridiagonal Matrix Algorithm）第8頁，共40頁，2022年，5月20日，12點29分，星期二（57） 第9頁，共40頁，2022年，5月20日，12點29分，星期二對于這種系數(shù)矩陣

6、為三對角線矩陣的方程組，TDMA求解的具體步驟是： （1）第一個方程的各項除以B1，得 ， ， 。第二個方程減去化簡后的第一個方程乘以A2，并將所得新方程的各項除以（ ），最后得 對第三個以后的每個方程都作同樣的處理，得到方程組系數(shù)的通式為：（58）第10頁，共40頁，2022年，5月20日，12點29分，星期二（2）自下而上的解方程 解的通式為（59）若有兩個或三個坐標是橢圓型的，則需對各橢圓型方向都進行逐線迭代掃描對兩個或三個方向各作了一次逐線計算，稱為進行了一次雙重或三重掃描，也叫一個迭代循環(huán)重復多重掃描，直至兩次相繼迭代得出的值差別不大為止。把最后結(jié)果作為n時間層的值對每條網(wǎng)格線進行計

7、算時，邊界條件很快傳入流場內(nèi)部，因此收斂速率比逐點迭代快得多。隱式差分方程的逐線迭代解法對時間步長也無限制，因此在實際中得到廣泛的應用第11頁，共40頁，2022年，5月20日，12點29分，星期二六、壓力和速度之間的耦合1、引言若知道流場中壓力的分布，則動量方程的求解與其它方程完全一樣，不會產(chǎn)生困難但是壓力值一般都是預先不知道的在可壓流中，壓力與密度間的關(guān)系由狀態(tài)方程確定，需通過連續(xù)方程和狀態(tài)方程確定在不可壓流中，流場中的壓力分布對速度場有很強的影響，壓力梯度以源項形式出現(xiàn)在動量方程中，但壓力卻沒有獨立的控制方程在計算流體力學的發(fā)展過程中，提出了多種不同的方法，來解決在以速度、壓力為求解變量

8、的原始變量法中的這個問題大致可分為：渦量-流函數(shù)方法、壓力修正方法、壓力與速度之間的迭代算法等幾種。以壓力修正方法中的SIMPLE系列方法應用最為廣泛第12頁，共40頁，2022年，5月20日，12點29分，星期二2壓力修正方法的基本思想在求解不可壓縮流體的流場問題時，如果我們把從動量方程與連續(xù)性方程離散得到的代數(shù)方程組聯(lián)立起來直接求解，就可以得到各速度分量及相應的壓力值這樣的直接解法要占用大量的計算機內(nèi)存，對于目前大多數(shù)工程應用場合還不適用如果采用分離式的迭代求解方法（segregated method），即先求解u速度場，再求解v速度場，則對壓力場因其無獨立的方程而無法對其求解（或無法改進

9、其原先的假定值）另一方面，上述分離式求解過程中只利用了u, v動量方程的離散形式而未用到連續(xù)性方程要解決的問題: 如何使得用質(zhì)量守恒方程假定的壓力場能不斷地隨迭代過程的進行而得到改進，這就是所謂的壓力修正算法壓力修正算法源于1972年由Patankar與Spalding提出的SIMPLE算法16，在20世紀80年代初期，又相繼提出了SIMPLER與SIMPLEC等方法，并由此形成了SIMPLE系列算法第13頁，共40頁，2022年，5月20日，12點29分，星期二3SIMPLE算法以二維直角坐標中的對流換熱求解為例，其控制方程為（60）（61）（62）在交錯網(wǎng)格（圖3）上動量離散方程為：（63

10、）（64）質(zhì)量守恒方程為（65）第14頁，共40頁，2022年，5月20日，12點29分，星期二采用SIMPLE算法來求解時計算步驟如下：（1）假定一個速度場，記為u(0)，v(0)，由此計算ae，an，anb及 b；（2）假定一個壓力場，記為 ；（3）求解動量離散方程（63）、（64），得 ；（4）計算壓力的修正值 ，要求與 相對應的 、 能滿足連續(xù)性方程第15頁，共40頁，2022年，5月20日，12點29分，星期二（5）計算速度的修正值 。要求 仍滿足線性化了的動量方程，即 與式（63）相減，得（69） 此式表明，要據(jù) 值確定 ，需要解一個代數(shù)方程組。為能利用 值顯式地求解 ，此處略去式

11、（69）右端第一項， ， 于是得（70）（6）將 、 作為本迭代層次之解，開始下一層次的計算，重復15步直到流場收斂，即所解得之流速能同時滿足連續(xù)性方程及動量方程為止第16頁，共40頁，2022年，5月20日，12點29分，星期二SIMPLE算法中引入了以下假定或簡化處理：（1）速度場的假定（u（0），v（0）與壓力場的假定（ ）是相互獨立地進行的，u(0)，v(0)與 間沒有任何聯(lián)系；（2）在導出速度修正值計算式（70）時，未計及鄰點速度修正值的影響；（3）動量離散方程中的b(控制容積P的剩余質(zhì)量流量)在速度修正前后保持不變；（4）由式（63）、（64）解出的 滿足動量守恒但未必滿足質(zhì)量守恒

12、，而由式（69）解出的 決定 時，保證 、 滿足質(zhì)量守恒，但動量守恒則未必滿足。注意，雖然在導出式（69）過程中曾要求 、 滿足線性化了的動量方程，但由于以后略去了式（69）中的 ，因而所得的 ， 就未必可使 、 滿足上述動量守恒方程。以上四條假設(shè)或近似處理，是SIMPLE算法提出之后所出現(xiàn)的一些改進方案的著眼之處。第17頁，共40頁，2022年，5月20日，12點29分，星期二4、SIMPLER算法 SIMPLER（1980）算法主要用以改進SIMPLE方法中的第一項近似處理方法。一旦速度場給定，壓力場就可以從動量離散方程中予以求解，而不再任意假定，其主要計算步驟如下：（1）假定一個速度場，

13、記為u(0)，v(0)，由此計算動量離散方程系數(shù)ae，an，anb，b及 ， ：（71） 注意，引入 ， 后，動量離散方程便可寫為（72）第18頁，共40頁，2022年，5月20日，12點29分，星期二（2）據(jù) 計算相應的壓力場 。為此，將式（72）代入質(zhì)量守恒方程的離散形式（65），得（73） 其中，anb及b的計算式形式上與SIMPLE算法中的方程的一樣，只要將 代替 即可。（3）求解動量離散方程，得替 。（4）求解壓力修正方程，得 。（5）用 修正速度得 （但不修正壓力）。（6）用 ， 開始下一層次的迭代，重復15步直到收斂。為有利于非線性迭代的收斂，SIMPLER算法中流速應予以亞松弛

14、，但 及p則不作亞松弛。第19頁，共40頁，2022年，5月20日，12點29分，星期二5、SIMPLEC算法6SIMPLEX算法7、SIMPLET算法第20頁，共40頁，2022年，5月20日，12點29分，星期二8SIMPLE系列算法小結(jié)SIMPLE，SIMPLEC，SIMPLER，SIMPLET及SIMPLEX的主要計算步驟可以歸納為：（1）由假定的或上一次計算得到的速度u(0)，v(0)確定動量離散方程的系數(shù)及常數(shù)項；（2）確定壓力場：對于SIMPLE，SIMPLEC，SIMPLET及SIMPLEX采用假定的或用上一次的計算值 ；對于SIMPLER則由已知的速度場通過求解壓力Poiss

15、on方程而獲得；（3）求解動量離散方程以獲得 ；（4）求解壓力修正值 方程，要求與 相應的 能 使 、 滿足連續(xù)性方程。其中除SIMPLET外，其它算法只考慮壓力修正值對 的影響，而SIMPLET中還考慮了溫度修正對 的影響，該影響最后體現(xiàn)在方程的源項之中；第21頁，共40頁，2022年，5月20日，12點29分，星期二（5）確定速度修正值：SIMPLE，SIMPLER，SIMPLET： ，鄰點速度影響略去不計；SIMPLEC： ，鄰點速度修正值的影響已基本考慮在內(nèi)；SIMPLEX： ，de等是通過求解一個Poisson方程得出的，相當于考慮了鄰點的影響。（6）采用下列速度、壓力值作為初值開始

16、下一層次的迭代計算（設(shè)速度的亞松弛已組合在u，v的求解過程中）：SIMPLE，SIMPLET：SIMPLEC，SIMPLEX：SIMLER： 以上計算過程可以大致分為兩步：預估步（1）（3）及校正步（4）（6）。因而可以說上述五種算法都是兩步算法（一步預估，一步校正）。第22頁，共40頁，2022年，5月20日，12點29分，星期二9加速迭代收斂的一些方法1）選擇合適的松馳因子之值對于正交的網(wǎng)格以及雖非正交但網(wǎng)格傾斜不是嚴重的情況，建議au與ap的取值應滿足以下關(guān)系：（76） 式中常數(shù)c取為123或1.124，同時au之值應盡可能地取得大，一般可取0.70.8。上述取值原則對交錯網(wǎng)格、同位網(wǎng)格

17、、正交網(wǎng)格與非正交網(wǎng)格原則上都適用，但對網(wǎng)格嚴重傾斜的情況，ap之值應小于按式（76）得出之值。第23頁，共40頁，2022年，5月20日，12點29分，星期二2）其它加速迭代收斂的方法迭代循環(huán)方式對收斂速率的影響描述燃燒過程的是一組互相耦合的方程組，在用迭代法求解時，迭代循環(huán)方式對收斂速率有很大的影響可將流體力學方程（動量和修正壓力方程）和其它標量方程（能量、湍流量方程等）一起迭代另一種選擇是將流體力學方程首先迭代若干次，然后把其它標量方程（能量、湍流量方程等）一起迭代；流體力學方程首先迭代若干次，然后把其它標量方程再加進去一起循環(huán)兩者的收斂速率是不一樣的。但目前還給不出最佳循環(huán)方式的一般原

18、則網(wǎng)格均勻度對收斂速率的影響網(wǎng)格均勻度對收斂速率有很大的影響計算時間與結(jié)點數(shù)成正比但增加結(jié)點數(shù)，若能使差分網(wǎng)格更均勻，由于收斂速率加快，有時仍有可能減少總的計算時間第24頁，共40頁，2022年，5月20日，12點29分，星期二七、邊界條件的數(shù)值處理方法1、引言 初始條件和邊界條件: 對拋物型問題，只需給出一端的邊界條件，如時間坐標上給出初始的參數(shù)分布，邊界層流動中給定進口截面上的參數(shù)分布等等對橢圓型問題，則需給出求解域兩端的邊界條件，如果三個空間坐標都是橢圓型坐標，則求解域的四周各點都需給出邊界條件第25頁，共40頁，2022年，5月20日，12點29分，星期二2、進口邊界 模型計算湍流時，

19、一般無法獲得由實驗測得的進口處k及 k:（1） （2） : (1) (2) 進口截面上k與 取值的這種不確定性對計算結(jié)果的影響（亦即計算結(jié)果對進口條件的敏感性）應當在數(shù)值計算中進行考察。第26頁，共40頁，2022年，5月20日，12點29分，星期二3、出口邊界描寫有回流的流動與燃燒的控制方程是關(guān)于空間坐標的橢圓型方程，從數(shù)學上要求在每個坐標方向上都應當有兩個邊界條件。出口邊界應該置在什么位置及怎樣確定這些人為設(shè)定的計算邊界條件，就成了計算流體力學及燃燒過程模擬中的一個重要研究課題。出口邊界應設(shè)置在什么位置上:Patankar曾認為，出口邊界不能設(shè)置在有回流的截面上，否則得到的解是沒有意義的數(shù)

20、值實踐證明，對出口邊界位置選擇的這種苛刻的要求是不必要的文獻27提出對這種出口邊界條件定性的要求是：（1）應當對出口邊界附近，尤其是離開出口邊界較遠處的流場沒有影響；（2）出口邊界應當是“透明的”, 又稱為“無反射”（3）無論出口邊界設(shè)置在何處，應對其它區(qū)域中的數(shù)值解沒有影 響第27頁，共40頁，2022年，5月20日，12點29分，星期二圖11求解流場時出口邊界位置的確定第28頁，共40頁，2022年，5月20日，12點29分，星期二目前，文獻中已提出了多種關(guān)于設(shè)置出口邊界條件的方法，歸納起來大約有以下幾種情形：1）局部單向化假設(shè)只有當出口截面無回流時才能近似地成立2）充分發(fā)展的假設(shè)假設(shè)n為

21、出口截面的法線方向，這一條件的數(shù)學表達式為（84）3）法向速度局部質(zhì)量守恒、切向速度齊次Neumann條件以圖11（b）中出口邊界相鄰的控制容積為例，取出后放大而示于圖12中。在本例中速度u是與出口截面平行的流速，而v為與出口截面垂直的流速。對所示控制容積作質(zhì)量守恒運算，得（85a）而按齊次Neumann條件有（85b）第29頁，共40頁，2022年，5月20日，12點29分，星期二圖12 局部質(zhì)量守恒的圖示 圖13 說明切向速度按局部質(zhì)量守恒的計算圖示第30頁，共40頁，2022年，5月20日，12點29分，星期二4）出口截面上的法向速度與切向速度都按局部質(zhì)量守恒確定在這一方法中，對圖12所

22、示情形， 計算式同式（85a），但ui,M的更新則通過對圖13所示半個控制容積的質(zhì)量守恒來獲得。以0.5(vi,M+vi,M-1)作為該半個控制容積南側(cè)流速，則可以寫出（86a） 將邊界上的切向流速取本迭代層次之值而其余速度取上一迭代層次的值，得 （86b）由圖11（b）可見，在計算區(qū)域右邊界， 之值為零，由此可推得 ，這樣依次計算可以得出出口邊界切向速度的更新值。第31頁，共40頁，2022年，5月20日，12點29分，星期二5）對流邊界條件用以下公式來更新出口截面法向速度及其它標量： （87a） 其中，c為對流速度或相速度，一般取為通道平均流速。上式離散后，可得（87b） 這里，vm為出口

23、截面平均法向流速。相當于假定在出口邊界上無擴散作用, 為一個純對流問題。第32頁，共40頁，2022年，5月20日，12點29分，星期二6）內(nèi)點插值方法 在文獻36中，提出采用由內(nèi)點的速度來決定出口流速的方法：（88a）其中，為與出口邊界相鄰的控制容積的寬度，x為出口邊界上游第四個截面位置起算的坐標，當網(wǎng)格均分時（如圖14所示），有（88b）第33頁，共40頁，2022年，5月20日，12點29分，星期二圖14 插值公式（88a）的圖示第34頁，共40頁，2022年，5月20日，12點29分，星期二7）各種高階導數(shù)為零的方法（89） 其中，x為出口方向，上角標j及q表示求導的階數(shù)，要求 ，并建

24、議取j=3，q=2。這一類方法實際上相當于提出了由內(nèi)點來插值的一個方式。導數(shù)的階數(shù)越高，涉及的內(nèi)點數(shù)越多，該階導數(shù)取零對解的影響也越微弱。第35頁，共40頁，2022年，5月20日，12點29分，星期二 出口邊界上邊界條件的設(shè)定有多種看來相差很大的方法。在從事工程流動與傳熱問題的數(shù)值計算中實用的處理出口邊界條件的方法：出口邊界上無回流時，可采用局部單向化或法向一階導數(shù)為零的方法；當出口邊界有回流時，可以采用法向流速按質(zhì)量守恒、切向流速按齊次Neumann條件的方法；無論采用哪一種方法，出口截面上的法向速度分布必須滿足總體質(zhì)量守恒的條件。第36頁，共40頁，2022年，5月20日，12點29分，

25、星期二4固體邊界 1）對于溫度一類的變量 用附加源項法來處理 2）對于速度固體壁面是受限流動中最常見的邊界，原則上應采用無滑移條件，即流體的速度等于壁面的速度（第一類邊界條件），層流時無論是對于動量方程中的流速還是連續(xù)性方程中的流速都是這樣（后者表現(xiàn)在SIMPLE算法壓力修正方程中的剩余質(zhì)量流量的計算上）當采用標準 模型時，平行于壁面的流速及垂直于壁面的流速還存在不同處理方法。對于平行于壁面的流速（例如u），采用無滑移條件，同時壁面上的動力粘性系數(shù)由壁面函數(shù)法的公式來計算；但對垂直于該壁面的流速（v）則應區(qū)分動量方程還是連續(xù)方程。在連續(xù)性方程中仍然采用無滑移條件（靜止壁面上取v=0），而對v動量方程則采用法向一階導數(shù)為零的條件： 第37頁，共40頁，2022年，5月20日，12點29分，星期二 3）對于k采用低Reynolds數(shù) 模型時，可直接令壁面上k=