管道開裂射流場的數(shù)值模擬
本章利用已有的流場求解程序，通過有限體積法求解氣體動力學(xué)基本方程。文中對開裂后的高壓輸氣管道氣體逸出引起的可壓縮高速自由射流流場進(jìn)行了數(shù)值模擬，得到了作用于管壁的分布?xì)怏w壓力和管道內(nèi)外氣流場的特征參數(shù)。
本章數(shù)學(xué)符號獨立，與其他章節(jié)沒有繼承關(guān)系。
5.1湍流場的數(shù)值模擬方法
5.1.1直接數(shù)值模擬（DNS)
直接模擬方法即不引入任何湍流模型，通過高精度、低色散、低耗散的差分格式求解完整的非定常N-S方程，對湍流的瞬時運動進(jìn)行直接的數(shù)值模擬。
這種方法的優(yōu)點是：方程本身是精確的，僅有的誤差只是由數(shù)值方法引入的誤差；數(shù)值模擬可以提供每一瞬間所有流動量在流場上的全部信息。
但是，湍流的直接數(shù)值模擬一直受到計算機速度和容量的限制。它的難點及存在的問題主要有：計算量十分巨大；流動平均量一般要比擾動量大幾個量級。求解時，由于差分格式存在截斷誤差，從而可能將真實擾動量淹沒掉。
理論上，直接數(shù)值模擬是射流流場數(shù)值模擬最精確的方法，但由于該方法計算量巨大，對計算機的計算能力有很高的要求，目前的計算機性能還不能夠滿足復(fù)雜流動的直接數(shù)值模擬的要求，利用直接數(shù)值模擬求解射流流場問題還沒有見到和實驗結(jié)果進(jìn)行比較的文獻(xiàn)。
5.1.2大渦模擬（LES)
大渦模擬的基本思想是把包括脈動運動在內(nèi)的湍流瞬時運動通過某種濾波方法分解成大尺度運動和小尺度運動兩個部分。大尺度量要通過數(shù)值求解運動微分方程直接計算出來；小尺度運動對大尺度運動的影響將在運動方程中表現(xiàn)為類似雷諾應(yīng)力一樣的應(yīng)力項，稱之為亞格子雷諾應(yīng)力，它們將通過建立模型來模擬。
大渦模擬的最主要的困難不在于計算機的限制，而是大渦模擬方法本身。例如，現(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)展。
所謂的湍流模型理論就是以雷諾平均運動方程為基礎(chǔ)，依靠理論與經(jīng)驗的結(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é)果與實驗值進(jìn)行了比較，吻合較好：
EghlimiA．和SahajwallaV.(1999)分別采用了k-ε湍流模型、RNGk-ε湍流模型、Realizek-ε湍流模型和雷諾應(yīng)力湍流模型對軸對稱的亞聲速自由射流進(jìn)行了數(shù)值模擬。通過計算值與實驗值的比較，可以發(fā)現(xiàn)除了標(biāo)準(zhǔn)的k-ε湍流模型對軸對稱的亞聲速自由射流流場數(shù)值模擬的結(jié)果與實驗值有較大偏差外，采用其他三種湍流模型計算結(jié)果與實驗值都吻合得比較好；
NASA的ReddyD.R.,SteffenC.J.,ZamanK,B.M.Q.(1999)應(yīng)用S-A湍流模型對三維方形噴嘴自由射流進(jìn)行了數(shù)值模擬，計算結(jié)果表明，在超聲速流動下，計算值與實驗值吻合較好，在亞聲速流動時偏差較大。
5.2 S-A湍流模型
根據(jù)流場數(shù)值計算領(lǐng)域已有的研究經(jīng)驗，結(jié)合本文問題的性質(zhì)，在大量試算的基礎(chǔ)上，本文選擇了S-A湍流模型。這一節(jié)介紹該模型的原理。
S-A湍流模型是Spalart和Allmaras在1992年提出的一方程湍流模型，具體描述如下：
湍流運動粘性系數(shù)
滿足輸運方程
S-A湍流模型常用于大梯度、近壁的氣體流動的數(shù)值模擬，在渦輪機械、高速氣體流動等方面的計算中有著廣泛的應(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)格的特點
在結(jié)構(gòu)網(wǎng)格中，網(wǎng)格點總是分布在某種坐標(biāo)變換后的坐標(biāo)線上，在拓?fù)渖鲜且?guī)范的（個別變換的奇點除外），而非結(jié)構(gòu)性網(wǎng)格則不要求拓?fù)渖系囊?guī)范性。
非結(jié)構(gòu)性網(wǎng)格的網(wǎng)格點在空間的分布比較自由，不要求每個網(wǎng)格點具有相同數(shù)量的鄰點，每一個網(wǎng)格點參與構(gòu)成的單元數(shù)量也可以互不相同，網(wǎng)格點之間的連接不再具有方向性。在非結(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)格點（內(nèi)點或邊界點），均可以相互連接而構(gòu)成一個計算單元，只要在這個三角形內(nèi)部以及其邊上，沒有其他的網(wǎng)格點。
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)格點：第二步是把邊界上的網(wǎng)格點和內(nèi)部已生成的網(wǎng)格點，連接成合適的三角形網(wǎng)格劃分；
一步法則把網(wǎng)格點的生成與三角形單元的形成統(tǒng)籌考慮，在生成網(wǎng)格點的同時考慮連接關(guān)系，而在考慮三角形的連接時，又同時考慮網(wǎng)格點的增刪，這類工作的代表是Delaunay方法。對于三角形網(wǎng)格單元，每個三角形越接近于正三角形時，網(wǎng)格的素質(zhì)越好。
本文中的非結(jié)構(gòu)性網(wǎng)格的生成是通過網(wǎng)格生成軟件包GAMBIT來實現(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)格單元中與控制面相對的頂點值，因而對網(wǎng)格的依賴性較大，要求網(wǎng)格的均勻性較好。
然而，在粘流的計算中，出于分辨邊界層的需要，網(wǎng)格在壁面法線方向相當(dāng)密集。在現(xiàn)有的計算條件下，采用均勻網(wǎng)格是做不到的。NND格式和基于Roe矢通量差分分裂的迎風(fēng)格式，只借助于控制面兩側(cè)單元中心點值，因而對網(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)成的封閉體，
Ω為該封閉體的邊界面，
為邊界面上的
外法線單位矢量。原則上積分路線可以任意選擇，只要包圍計算點即可。
以二維情況為例，計算圖5-1中A點梯度可選C→B→D→C作為積分路徑，但為了保證積分的精度，最好使積分路徑構(gòu)成的封閉體形心與計算點重合，或至少比較接近。
但是，由（5-18）直接得到的格式在激波附近是不穩(wěn)定的，將會出現(xiàn)非物理振蕩。迎風(fēng)格式通過限制物理量的梯度值，來保證物理量在單元內(nèi)的分布具有單調(diào)性，從而達(dá)到抑制間斷附近非物理振蕩的目的。這時，（5-18）式變?yōu)?/DIV>
5.5.2粘性通量
的確定
粘性通量
由控制面上的速度、壓力及密度的梯度項組成，所以粘性通量
的計算主要是求出控制面上相關(guān)物理量的梯度值。首先由（5-19）式求出各單元內(nèi)物理量ф的梯度值后，界面（i，j）上物理量的梯度值可簡單地由相鄰單元內(nèi)物理量ф的梯度值的算術(shù)平均得出：
5.5.3離散方程的求解
在計算出無粘通量
和粘性通量
后，代入（5-14）式可得到半離散的差分方程：
方程（5-24）可以采用顯式或隱式的方法進(jìn)行求解。但是由于受時間步長的限制，采用顯式的R-K時間推進(jìn)方法需要迭代的步數(shù)較多，計算很不經(jīng)濟，因此本文采用了Gauss-Seidel迭代隱式求解的方法。
式（5-24）在時間方向采用Euler后向差分，有
對方程（5-25),采用Gauss-Seidel點迭代法隱式求解。
5.5.4邊界條件
在邊界處理方面，采用傳統(tǒng)的邊界處理方法。
1.計算域的入口采用入流條件，亞聲速時給定總壓、總溫和速度的方向，其他參數(shù)值由內(nèi)點外推；超聲速時給定邊界上的所有參數(shù)值。
2.計算域的出口采用出流條件，亞聲速的時候在出口截面上給定環(huán)境壓力，其他參數(shù)值由內(nèi)點外推，超聲速出流邊界不提任何條件；
3.對稱軸處采用軸對稱條件，對稱軸上法向速度為O，其他參數(shù)值一階偏導(dǎo)為O；
4.固壁采用無滑移固壁條件，u=v=w=O；
5.人工邊界用無反射邊界條件。
其中無反射邊界條件是根據(jù)特征線的方法來求出邊界上的參數(shù)值。根據(jù)一維完全氣體無粘非定常運動的方程組，在亞聲速的條件下由邊界上的入射波和出射波關(guān)系可以得到兩個Riemann不變量：
其中Vn為沿邊界法向的速度分量，c為當(dāng)?shù)芈曀�，下�?biāo)∞表示無窮遠(yuǎn)處的參數(shù)值，下標(biāo)i表示計算域內(nèi)緊鄰邊界的點的參數(shù)值。
無窮遠(yuǎn)處的參數(shù)值為已知，計算域內(nèi)的參數(shù)值可以通過內(nèi)點的計算得出，因此我們可以得到邊界上的沿邊界法向方向的速度分量和當(dāng)?shù)氐穆曀僦担?/DIV>
邊界上的沿邊界切向方向的速度分量和熵的值可以通過內(nèi)點外推得到。
5.6結(jié)果分析
圖5-2是第二章穩(wěn)態(tài)裂紋擴展計算得到的管道壁道面在某上瞬時的變形狀態(tài)。
圖5-3是以圖5-2為依據(jù)建立的流場計算網(wǎng)格。設(shè)計參數(shù)為：裂紋擴展速度200m/s，管道外徑1016mm，壁厚14.7mm，內(nèi)壓10MPa。
外部流場區(qū)域選取范圍為3m×5m。輸送介質(zhì)選用天然氣的主要成分甲烷，初始外流場設(shè)定為空氣。遠(yuǎn)場壓力定義為大氣壓。在底部和側(cè)面應(yīng)用對稱邊界條件。計算段管道全長為35m，裂尖位于距軸向?qū)ΨQ面18.88m處。管壁的擴張速度按照第二章的計算結(jié)果作為流場的運動壁面邊界條件。
計算的結(jié)果揭示了一些重要的現(xiàn)象。
圖5-4給出了XZ對稱面上的馬赫數(shù)分布。管道中的流動在裂尖出口處達(dá)到超聲速，而在對稱面附近發(fā)生了強烈的回流。標(biāo)示為B的位置壓力和密度發(fā)生了強間斷，意味著出現(xiàn)了激波。
為了進(jìn)一步觀察XY截面上的速度，圖5-5繪出了圖5-4中標(biāo)示的A-A截面的速度分布。
矢量箭頭表示的是X-Y面內(nèi)速度分量的方向，等值線表示三維流動區(qū)域中的速率大小。從圖上可以看出外部流場的速度分布沿Y方向從中心到兩側(cè)速度迅速降低。在距管道中心五倍遠(yuǎn)的外流場，氣流速度仍可達(dá)到671m/s。
圖5-5右側(cè)的放大截面上標(biāo)出了A至E的五個點。保持各點的X坐標(biāo)不變，變化Z坐標(biāo)，可在管壁上分別劃出對應(yīng)于截面上不同位置的軸線。
由于管道內(nèi)壁處的氣體靜壓等于該處管道承受的壓力，圖5-6和圖5-7中A-E各點沿Z向的氣體靜壓變化反映了管道承受壓力的軸向衰減規(guī)律。兩圖中的初始內(nèi)壓分別為10MPa和8MPa。
與線性衰減模式（2-47）（某估算結(jié)果在圖5-6中用點劃線標(biāo)出）相比，計算得到的靜壓衰減結(jié)果揭示了幾個重要問題。
其一是回流區(qū)的發(fā)現(xiàn)。在軸向兩側(cè)擴展裂紋的對稱面附近發(fā)現(xiàn)了強烈的回流，這在裂紋擴展初期將導(dǎo)致管壁的裂紋后端的抖動，從而對裂紋驅(qū)動力起到重要影響。隨著裂尖位置的不斷前移，回流區(qū)的壓力將逐漸衰減到大氣壓值。
其二是回流區(qū)的位置。計算表明，回流區(qū)寬度隨管道工作壓力的升高而增長，隨裂尖位置的前移而減小。
其三是衰減區(qū)的壓力變化。氣流通過裂紋尖端時，管壁上作用的氣體壓力沿軸向迅速衰減，經(jīng)驗公式（2-47）中所用的衰減長度一般取為1.5倍管道直徑，最終衰減到零。計算表明衰減長度隨管道工作壓力的升高而增加，并且衰減區(qū)壓力的最低值與管道的工作壓力相關(guān)。由于對裂紋的擴展?fàn)顟B(tài)影響甚微，裂尖前部的壓力不是關(guān)心的對象。
此外，本文計算與公式（2-46）基本相符，這和全尺寸爆裂實驗得到的結(jié)果是一致的。在（2-46）的基礎(chǔ)上，本文提出以下的修正壓力衰減模式以替代（2-47）進(jìn)
通過多種工況驗證，該模式與計算結(jié)果吻合良好。
另外，通過和意大利CSM的實測管理壓力曲線圖5-8相比較，證明圖5-6和圖5-7從趨勢上是成立的。
5.7本章小結(jié)
本章以有限體積法為基礎(chǔ)，對非定常可壓縮的Navier-Stokes方程進(jìn)行了離散和求解。網(wǎng)格的劃分采用非結(jié)構(gòu)性網(wǎng)格，無粘通量采用Roe矢通量分裂，差分格式為二階精度，時間離散采用Guass-Seidel隱式迭代，在人工邊界上采用基于特征線法的無反射邊界條件。
對含動態(tài)裂紋的高壓輸氣管道的流場模擬結(jié)果表明，氣體壓力模式（2-46）和(2-47）式在裂紋充分發(fā)展至不考慮回流區(qū)影響的條件下，對高壓高韌性管道仍然成立
回流區(qū)的發(fā)現(xiàn)對裂紋開始擴展不久的一段時間內(nèi)的裂紋驅(qū)動力起重要影響，對于判定止裂點的位置有著重要的工程意義。（5-29）和（5-30）式可直接應(yīng)用于有限元程序以提高裂紋擴展與止裂的模擬精度。
