管道開裂射流場的數(shù)值模擬
本章利用已有的流場求解程序，通過有限體積法求解氣體動力學(xué)基本方程。文中對開裂后的高壓輸氣管道氣體逸出引起的可壓縮高速自由射流流場進(jìn)行了數(shù)值模擬，得到了作用于管壁的分布?xì)怏w壓力和管道內(nèi)外氣流場的特征參數(shù)。
本章數(shù)學(xué)符號獨(dú)立，與其他章節(jié)沒有繼承關(guān)系。
5.1湍流場的數(shù)值模擬方法
5.1.1直接數(shù)值模擬（DNS)
直接模擬方法即不引入任何湍流模型，通過高精度、低色散、低耗散的差分格式求解完整的非定常N-S方程，對湍流的瞬時運(yùn)動進(jìn)行直接的數(shù)值模擬。
這種方法的優(yōu)點(diǎn)是：方程本身是精確的，僅有的誤差只是由數(shù)值方法引入的誤差；數(shù)值模擬可以提供每一瞬間所有流動量在流場上的全部信息。
但是，湍流的直接數(shù)值模擬一直受到計算機(jī)速度和容量的限制。它的難點(diǎn)及存在的問題主要有：計算量十分巨大；流動平均量一般要比擾動量大幾個量級。求解時，由于差分格式存在截斷誤差，從而可能將真實(shí)擾動量淹沒掉。
理論上，直接數(shù)值模擬是射流流場數(shù)值模擬最精確的方法，但由于該方法計算量巨大，對計算機(jī)的計算能力有很高的要求，目前的計算機(jī)性能還不能夠滿足復(fù)雜流動的直接數(shù)值模擬的要求，利用直接數(shù)值模擬求解射流流場問題還沒有見到和實(shí)驗(yàn)結(jié)果進(jìn)行比較的文獻(xiàn)。
5.1.2大渦模擬（LES)
大渦模擬的基本思想是把包括脈動運(yùn)動在內(nèi)的湍流瞬時運(yùn)動通過某種濾波方法分解成大尺度運(yùn)動和小尺度運(yùn)動兩個部分。大尺度量要通過數(shù)值求解運(yùn)動微分方程直接計算出來；小尺度運(yùn)動對大尺度運(yùn)動的影響將在運(yùn)動方程中表現(xiàn)為類似雷諾應(yīng)力一樣的應(yīng)力項，稱之為亞格子雷諾應(yīng)力，它們將通過建立模型來模擬。
大渦模擬的最主要的困難不在于計算機(jī)的限制，而是大渦模擬方法本身。例如，現(xiàn)有的亞格子應(yīng)力模型還很不完善，特別是近壁區(qū)內(nèi)的模型以及邊界條件的提法等等都是急待解決的問題。使用該方法進(jìn)行射流流場的數(shù)值模擬也還在探索之中。
5.1.3湍流模型理論
湍流模型理論在射流流場的數(shù)值模擬中已經(jīng)取得了較好的應(yīng)用，目前人們已經(jīng)廣泛的應(yīng)用湍流模型理論對射流的流動進(jìn)行分析，并取得了很大的進(jìn)展。
所謂的湍流模型理論就是以雷諾平均運(yùn)動方程為基礎(chǔ)，依靠理論與經(jīng)驗(yàn)的結(jié)合，引進(jìn)一系列模型假設(shè)，建立一組描寫湍流平均量的封閉方程組的理論計算方法。目前，在射流流場的數(shù)值模擬中常用的湍流模型有B-L混合長模型，S-A一方程湍流模型，RNGK-ε方程湍流模型，雷諾應(yīng)力湍流模型等。
在湍流模型的選擇方面：
AndrewT.T.和ChistopherK.W．T.（1996)應(yīng)用修正的k-ε二方程湍流模型計算了軸對稱及方形、橢圓形噴嘴的自由射流，馬赫數(shù)為0.4～2，計算結(jié)果與實(shí)驗(yàn)值進(jìn)行了比較，吻合較好：
EghlimiA．和SahajwallaV.(1999)分別采用了k-ε湍流模型、RNGk-ε湍流模型、Realizek-ε湍流模型和雷諾應(yīng)力湍流模型對軸對稱的亞聲速自由射流進(jìn)行了數(shù)值模擬。通過計算值與實(shí)驗(yàn)值的比較，可以發(fā)現(xiàn)除了標(biāo)準(zhǔn)的k-ε湍流模型對軸對稱的亞聲速自由射流流場數(shù)值模擬的結(jié)果與實(shí)驗(yàn)值有較大偏差外，采用其他三種湍流模型計算結(jié)果與實(shí)驗(yàn)值都吻合得比較好；
NASA的ReddyD.R.,SteffenC.J.,ZamanK,B.M.Q.(1999)應(yīng)用S-A湍流模型對三維方形噴嘴自由射流進(jìn)行了數(shù)值模擬，計算結(jié)果表明，在超聲速流動下，計算值與實(shí)驗(yàn)值吻合較好，在亞聲速流動時偏差較大。
5.2 S-A湍流模型
根據(jù)流場數(shù)值計算領(lǐng)域已有的研究經(jīng)驗(yàn)，結(jié)合本文問題的性質(zhì)，在大量試算的基礎(chǔ)上，本文選擇了S-A湍流模型。這一節(jié)介紹該模型的原理。
S-A湍流模型是Spalart和Allmaras在1992年提出的一方程湍流模型，具體描述如下：
湍流運(yùn)動粘性系數(shù)
滿足輸運(yùn)方程
S-A湍流模型常用于大梯度、近壁的氣體流動的數(shù)值模擬，在渦輪機(jī)械、高速氣體流動等方面的計算中有著廣泛的應(yīng)用。
5.3N-S控制方程
非定常可壓縮的射流滿足如下的積分形式的N-S方程：
上式中，Ω是控制體，
Ω是控制體邊界面，W是求解變量，F(xiàn)是無粘通量，G是粘性通量。
式中P為密度，
表示速度，E為總能，
分別為沿x，y，z方向的單位矢量。
為粘性應(yīng)力項，
是熱通量。
將以上各項寫成分量形式：
其中μv是分子粘性系數(shù)，μt是湍流粘性系數(shù)，Prv是分子普朗特數(shù)，Prt是湍流普朗特數(shù)。
為了便于粘性通量的計算，采用原始變量
為變量，
定義為
5.4非結(jié)構(gòu)性網(wǎng)格
相對于結(jié)構(gòu)性網(wǎng)格而言，非結(jié)構(gòu)性網(wǎng)格能夠更有效的適應(yīng)形狀不規(guī)則、邊界彎曲的復(fù)雜領(lǐng)域，以及依據(jù)不同的分辨率要求，分配求解域內(nèi)的網(wǎng)格。特別是三角形網(wǎng)格（對于二維問題）和四邊形網(wǎng)格（對于三維問題）。
5.4.1非結(jié)構(gòu)性網(wǎng)格的特點(diǎn)
在結(jié)構(gòu)網(wǎng)格中，網(wǎng)格點(diǎn)總是分布在某種坐標(biāo)變換后的坐標(biāo)線上，在拓?fù)渖鲜且?guī)范的（個別變換的奇點(diǎn)除外），而非結(jié)構(gòu)性網(wǎng)格則不要求拓?fù)渖系囊?guī)范性。
非結(jié)構(gòu)性網(wǎng)格的網(wǎng)格點(diǎn)在空間的分布比較自由，不要求每個網(wǎng)格點(diǎn)具有相同數(shù)量的鄰點(diǎn)，每一個網(wǎng)格點(diǎn)參與構(gòu)成的單元數(shù)量也可以互不相同，網(wǎng)格點(diǎn)之間的連接不再具有方向性。在非結(jié)構(gòu)性網(wǎng)格的劃分中沒有網(wǎng)格線的概念，因而能夠更為有效地適應(yīng)形狀不規(guī)則、彎曲邊界的復(fù)雜領(lǐng)域，以及依據(jù)不同的分辨率要求，分配求解域內(nèi)的網(wǎng)格。非結(jié)構(gòu)性網(wǎng)格的劃分問題都直接在物理平面（或物理空間）中討論。
以三角形單元網(wǎng)格為例，互相靠近的不在同一直線上的三個網(wǎng)格點(diǎn)（內(nèi)點(diǎn)或邊界點(diǎn)），均可以相互連接而構(gòu)成一個計算單元，只要在這個三角形內(nèi)部以及其邊上，沒有其他的網(wǎng)格點(diǎn)。
5.4.2非結(jié)構(gòu)性網(wǎng)格的數(shù)值求解方法
非結(jié)構(gòu)性網(wǎng)格的數(shù)值求解方法通常必須采用有限體積法建立差分方程。
有限體積法是介于有限差分法和有限元法之間的一種計算方法，它針對一個有限體積的單元體，通過差分離散建立差分方程。由于非結(jié)構(gòu)性網(wǎng)格的網(wǎng)格劃分更接近于有限元的網(wǎng)格劃分思想，直接采用微分方程出發(fā)的差分方法很不方便，而必須采用從積分型方程出發(fā)。本文5.5節(jié)有詳細(xì)的介紹。
5.4.3非結(jié)構(gòu)性網(wǎng)格的生成
非結(jié)構(gòu)性網(wǎng)格的生成大致可分為兩大類：兩步法和一步法。
兩步法就是將非結(jié)構(gòu)性網(wǎng)格的生成分為兩個步驟完成，第一步是在求解域內(nèi)生成網(wǎng)格點(diǎn)：第二步是把邊界上的網(wǎng)格點(diǎn)和內(nèi)部已生成的網(wǎng)格點(diǎn)，連接成合適的三角形網(wǎng)格劃分；
一步法則把網(wǎng)格點(diǎn)的生成與三角形單元的形成統(tǒng)籌考慮，在生成網(wǎng)格點(diǎn)的同時考慮連接關(guān)系，而在考慮三角形的連接時，又同時考慮網(wǎng)格點(diǎn)的增刪，這類工作的代表是Delaunay方法。對于三角形網(wǎng)格單元，每個三角形越接近于正三角形時，網(wǎng)格的素質(zhì)越好。
本文中的非結(jié)構(gòu)性網(wǎng)格的生成是通過網(wǎng)格生成軟件包GAMBIT來實(shí)現(xiàn)的，采用的是Delaunay方法。
5.5有限體積法求解過程
采用有限體積法，控制體選取為網(wǎng)格單元，將物理量置于網(wǎng)格單元中心，如圖5-1所示。
對方程（5-10）在單元i上積分，有
(5-14）式中，nb(i）表示單元i的相鄰單元數(shù)，下標(biāo)“i,j'”表示單元i和單元j的交界面，
表示界面（i,j）處的無粘通量，
表示界面（i,j）處的粘性通量，
表示界面（i，j）的外法線單位矢量。
在三維的情況下，V
i為單元i的體積，A
i，j為界面（i,j）面積：在二維的情況下，V
i為單元i面積，A
i，j為界線（i,j）長度；Γ
i表示單元i的Jacobian矩陣Γ。
5.5.1無粘通量
的確定
確定無粘通量只愁的方法較多，大體上可以分成兩大類。
一類是中心型格式，這類格式需要加入人工粘性來抑制間斷附近的非物理振蕩，因此，對于這類格式必須較好的控制人工粘性系數(shù)，這帶有相當(dāng)大的人為性，并且難于找到普適的人工粘性系數(shù)，因此很不方便；
另一類是迎風(fēng)型格式，應(yīng)用的較為廣泛的有TVD,ENO和NND等格式。非結(jié)構(gòu)性網(wǎng)格上的TVD和ENO格式在求解Euler方程的時候，取得了較大的成功，但在粘性流體計算中需要借助于網(wǎng)格單元中與控制面相對的頂點(diǎn)值，因而對網(wǎng)格的依賴性較大，要求網(wǎng)格的均勻性較好。
然而，在粘流的計算中，出于分辨邊界層的需要，網(wǎng)格在壁面法線方向相當(dāng)密集。在現(xiàn)有的計算條件下，采用均勻網(wǎng)格是做不到的。NND格式和基于Roe矢通量差分分裂的迎風(fēng)格式，只借助于控制面兩側(cè)單元中心點(diǎn)值，因而對網(wǎng)格均勻性的要求有所減弱。本文采用基于Roe矢通量差分分裂的迎風(fēng)格式。計算無粘通量的時候，采用Roe矢通量差分分裂，有
ε為略大于零的數(shù)，一般取ε=0.05aM。
A1由下式計算：A1=|UM|，A2，3=|UM±aM|。
這樣，無粘通量
的確定就轉(zhuǎn)化成了
L和
R的確定。若取
便可方便的得出一階精度格式。顯然，一階精度格式對于粘流的計算來說是沒有應(yīng)用價值的。為了構(gòu)造高精度格式，必須尋求確定
L和
R的高精度方法。對于二階精度格式，
L和
R的表達(dá)式可由下式確定：
其中，（▽
）
i和（▽
）
j，分別為Q在單元i和j內(nèi)的梯度；k為界面（i,j）的中心，
表示單元j中心到k的距離矢量，
表示單元i中心到k的距離矢量。將（5-18)代入（5-16)，即可得到無粘通量
�？梢宰C明，這樣確定的
具有二階精度。
在（5-18）式，我們用到了原始變量的梯度（▽
）
i和（▽
）
j。在非結(jié)構(gòu)網(wǎng)格中，物理量的梯度可以由高斯積分公式求得。令ф代表任一變量，則
其中，Ω為積分路徑構(gòu)成的封閉體，
Ω為該封閉體的邊界面，
為邊界面上的
外法線單位矢量。原則上積分路線可以任意選擇，只要包圍計算點(diǎn)即可。
以二維情況為例，計算圖5-1中A點(diǎn)梯度可選C→B→D→C作為積分路徑，但為了保證積分的精度，最好使積分路徑構(gòu)成的封閉體形心與計算點(diǎn)重合，或至少比較接近。
的確定
由控制面上的速度、壓力及密度的梯度項組成，所以粘性通量
