矩陣分裂算法在兩相流方程求解中的應用初步研究

王章立，蘇光輝，秋穗正，田文喜

(1.西安交通大學 核科學與技術學院，陜西 西安 710049；2.上海核工程研究設計院，上海 200233)

?

矩陣分裂算法在兩相流方程求解中的應用初步研究

王章立1,2，蘇光輝1，秋穗正1，田文喜1

(1.西安交通大學 核科學與技術學院，陜西 西安 710049；2.上海核工程研究設計院，上海 200233)

矩陣分裂算法是一種有限差分方法，根據系統特征值的正負，確定對流項差分格式，能描述流場信息傳播的真正方向。本文基于隱式矩陣分裂算法對兩相流均相流模型進行離散，利用矩陣形式追趕法求解離散后的塊三對角矩陣，運用FORTRAN90程序實現數值求解過程。運用噴管內流動和兩相可壓縮性實驗兩個基準題對算法進行校驗，計算結果和解析結果吻合良好，證明求解方法的可靠性。本文成功運用矩陣分裂算法求解兩相流模型，為熱工水力計算程序的開發和改進提供了新的視角。

矩陣分裂算法；塊三對角矩陣；兩相流

熱工水力安全分析程序的發展受到計算流體力學和數值傳熱學的極大影響，這兩個學科中的優秀成果不斷運用到熱工水力安全分析程序中，推動它不斷升級、優化。對流體動力學方程進行求解，主要涉及到方程的空間離散、時間離散、離散代數方程組求解以及求解過程中的穩定性問題。國際上通用的熱工水力程序，對流體動力學方程的空間離散基本上均采用交錯網格，解決常規網格無法檢測出不可壓縮流體求解中不合理的壓力場問題[1-2]，交錯網格技術會給變量存儲、程序編制帶來很大的不便。這些熱工水力程序采用迎風差分格式離散對流[3]，僅根據流場速度的方向判斷上下游。但流場中系統特征值的正負代表流場信息的真正傳播的方向，利用流場速度來進行迎風差分并不能反映流場信息傳播的真正方向。

矩陣分裂算法是一種有限差分方法，根據系統特征值的正負，確定對流項差分格式，能描述流場信息傳播的真正方向。矩陣分裂算法基于原始變量，可很方便地應用于直角坐標系。本文采用隱式矩陣分裂算法對均相流模型進行求解，利用FORTRAN90程序實現整個數值求解過程，運用所編制的程序進行詳細計算和分析。

1 數學物理模型

蒸汽發生器二次側兩相流動采用均相流模型，它假設氣液兩相流速相等，將兩相流體看作均勻混合的單相處理，其控制方程[4-6]如下。

質量守恒：

(1)

式中：ρ為流體密度，kg/m3；u為流速，m/s。

動量守恒：

(2)

式中：p為壓力，Pa；F為單位體積流體所受力，N/m3。

能量方程：

(3)

式中：e為流體內能，J/kg；Q為單位體積流體內熱源，W/m3。

能量方程寫成熵的形式，有：

(4)

式中，s為流體熵，J/(kg·K)。

將式(1)、(2)和(4)展開并寫成矩陣形式，有：

(5)

(6)

(7)

(8)

系數矩陣G具有3個實特征值：

(9)

每個特征值對應的特征向量所構成的矩陣稱為特征矩陣，左特征矩陣和右特征矩陣如下：

(10)

(11)

2 隱式矩陣分裂算法

矩陣分裂算法是一種有限差分方法，根據系統特征值的正負，確定對流項差分格式。它最早運用于氣體動力學，Romstedt[7]首次將該方法拓展到兩相流體動力學中。本文對均相流模型下矩陣分裂算法進行了詳細推導，基于特征值正負方向對分裂后的方程和源項進行離散，并運用矩陣形式追趕法對離散后方程組進行求解。

(12)

首先考慮齊次方程(12)求解過程，以了解矩陣分裂方法基本思路。均相流模型為雙曲線系統，系數矩陣G可對角化，故有：

(13)

(14)

將其帶入式(12)，等式兩邊同乘以R-1，有：

(15)

小時間步長下可將方程線性化，令ΔW=R-1ΔU，有：

(16)

寫成分量形式有：

(17)

式(12)解耦為3個獨立的方程，這3個方程均為對流方程形式。計算流體力學和數值傳熱學中針對對流方程已有非常詳細的討論，一般采用迎風差分格式離散對流項。根據特征值λp的正負，對流項可采用向前差分或向后差分。

當λp≥0時，有：

(18)

當λp<0時，有：

(19)

將式(18)、(19)聯合表示成矩陣形式，有：

(20)

將ΔW=R-1ΔU代入上式，并同乘以矩陣R，有：

(21)

即：

(22)

系數矩陣G分裂為兩部分G+和G-，相應的變量采用向前或向后差分。這時，將矩陣G按特征值分裂為單獨矩陣，將會給G+和G-的求解帶來方便。定義分裂矩陣為：

(23)

有關系式：

(24)

由特征值的正負，將相應的分裂矩陣相加，可得到G+和G-。對于均相流方程組：

(25)

當流體控制方程含有源項時，其源項也應為迎風方向控制體內的源項。

當λp≥0時，有：

(26)

當λp<0時，有：

(27)

系統3個特征值為λ1=u+a，λ2=u-a，λ3=u，在流速u>0的情況下，λ1和λ3為正而λ2為負，將式(26)、(27)聯合寫成矩陣形式，有：

(28)

Si=R[l1Dil2Di+1l3Di]T

(29)

含有源項的均相流方程組，可按式(28)的格式進行離散。如圖1所示網格結構，采用隱式方法離散時間項，可得：

(30)

圖1 網格示意圖Fig.1 Schematic diagram of grid

化簡后有：

(31)

合并后有：

(32)

(33)

離散后的代數方程(32)中，系數矩陣A、未知量U和源項B均為矩陣，所有代數方程構成塊三對角矩陣，可采用矩陣形式追趕法求解[8]，代數方程可表示為：

(34)

系數矩陣可表示為：

(35)

對矩陣K進行LR有：

K=LR

根據矩陣乘法，可得到追趕過程遞推關系式，追的過程為：

L1=B1

(36)

(37)

趕的過程為：

UN=YN

(38)

(39)

3 算法驗證

噴管是熱力系統中重要的工程設備，在工程熱力系統和氣體動力學中均有重要應用，其結構如圖2所示，其結構和工作原理詳見文獻[9-10]。

圖2 漸縮漸擴噴管結構示意圖Fig.2 Schematic diagram of structure for convergent-divergent nozzle

由于噴管內氣體流動參數只與噴管截面積變化率有關，故噴管總截面積可乘以任何比例因子。假設噴管截面積按式(40)規律變化。

(40)

噴管內氣體采用氮氣，駐室內氣體參數為p0=0.5 MPa，s0=6 800 J/(kg·K)。噴管內氣體參數變化范圍內，比熱比γ變化非常小，可取γ=1.4。

不同控制體下噴管喉部參數列于表1。控制體數增加到50后，計算結果基本不再變化，故工況中噴管取50個控制體。

噴管內Ma、壓力和溫度的分布示于圖3。氮氣Ma沿噴管逐漸增加，在噴管喉部處等于1，速度增加到當地聲速。喉部前漸縮噴管內氣體作亞聲速流動，喉部后漸擴噴管內氣體作超聲速流動。氮氣壓力沿噴管逐漸減小，內能逐漸轉化為氣體動能，氣體溫度也逐漸降低。圖3表明，數值計算結果與解析結果符合良好，表明文中所用數值方法合理有效。

表1 控制體敏感性檢查Table 1 Sensitivity study of control volume

圖3 噴管內Ma、壓力和溫度的分布Fig.3 Distributions of Ma, pressure and temperature in nozzle

4 應用舉例

Okawa等[11]研究了搖擺條件下入口質量流速波動對圓管CHF的影響，研究對象如圖4所示，壁面采用均勻熱流密度加熱，入口流體質量流速正弦規律波動。采用一維三流體模型研究圓管內流體流動，非搖擺條件下計算結果如圖5所示，可看出加熱情況下流體質量流速沿流動方向波動振幅不斷衰減。

圖4 單管均勻加熱示意圖Fig.4 Schematic diagram of single tube with uniform heating

采用與Okawa工作中相似的初始條件，模擬加熱段內流體波動情況。圓管被劃分為100個控制體，管道入口質量流速維持2 s恒定后作以2 s為周期的正弦波動。圖6示出了1個周期內不同時刻質量流速變化曲線，流場分布曲線以0.1倍周期為間隔。從圖6可看出，入口流體含有一定含汽率，流經管道逐漸被加熱，流體熵升高，密度變小。兩相流體加熱過程中密度快速變化，其可壓縮性表現為流體質量流速振幅沿管道不斷減小。比較圖5、6，本文計算結果與Okawa計算結果趨勢符合較好，說明本程序能很好地模擬兩相流體流動過程中表現的可壓縮性。

圖5 Okawa計算的不同時刻質量流速變化曲線[11]Fig.5 Transient curve of mass flow rate in different time computed by Okawa[11]

圖6 不同時刻質量流速變化曲線Fig.6 Transient curve of mass flow rate in different time

5 結束語

本文將隱式矩陣分裂算法運用到均相流模型求解中，運用隱式方法離散方程時間項，根據系統特征值構造對流相離散格式，能描述流場信息傳播的真正方向，利用與矩陣形式追趕法求解離散后的塊三對角矩陣。運用噴管內流動和兩相可壓縮性實驗兩個基準題對算法進行校驗，計算結果和解析結果吻合良好，驗證了矩陣分裂算法數值求解的正確性。

[1] 陶文銓. 數值傳熱學[M]. 西安:西安交通大學出版社,2001.

[2] 陶文銓. 計算傳熱學的近代進展[M]. 北京:科學出版社,2005.

[3] STAEDTKE H. Advanced modelling and numerical strategies in nuclear thermal-hydraulics[C]. Ljubljana: Nuclear Society of Slovenia, 2001.

[4] STAEDTKE H. Gasdynamic aspects of two-phase flow[M]. Germany: Federal Republic, 2006.

[5] 徐繼鋆,賈斗南. 沸騰傳熱和汽液兩相流[M]. 北京:原子能出版社,2000.

[6] 閻昌琪. 汽液兩相流[M]. 哈爾濱:哈爾濱工程大學出版社,1995.

[7] ROMSTEDT P. A split-matrix method for the numerical solution of two-phase flow equations[J]. Nuclear Science and Engineering, 1990, 104(1): 1-9.

[8] 鄧建中, 劉之行. 計算方法[M]. 西安:西安交通大學出版社,2001.

[9] JOHN D A. 計算流體力學基礎及其應用[M]. 北京:機械工業出版社,2007.

[10]沈維道,蔣志敏,童鈞耕. 工程熱力學[M]. 北京:高等教育出版社,2001.

[11]OKAWA T, GOTO T, MINAMITANI J, et al. Liquid film dryout in a boiling channel under flow oscillation conditions[J]. International Journal of Heat and Mass Transfer, 2009, 52(15-16): 3 665-3 675.

Preliminary Study on Split Coefficient Matrix Method for Two-phase Flow Equation Solving

WANG Zhang-li1,2, SU Guang-hui1, QIU Sui-zheng1, TIAN Wen-xi1

(1.SchoolofNuclearScienceandTechnology,Xi’anJiaotongUniversity,Xi’an710049,China; 2.ShanghaiNuclearEngineeringResearch&DesignInstitute,Shanghai200233,China)

The split coefficient matrix method is one of finite difference methods, which performs discretization on the convective term according to eigenvalues of the equations. That method can describe the real information propagation direction of flow field. The system equations were integrated based on the split coefficient matrix method in the paper. The generated block triple diagonal algebra matrix was solved with the method of forward elimination and backward substitution. The developed code was assessed by several separated effect tests. These tests were carried out through convergent-divergent nozzles and two-phase flow compression test. Good agreements between the computation and analysis results prove the validity and applicability of the method. The split coefficient matrix method was successfully applied to solve the two-phase flow equations. The results provide a new point of view for the thermal-hydraulic code development.

split coefficient matrix method; block triple diagonal algebra matrix; two-phase flow

2014-08-01;

2014-09-16

王章立(1985—)，男，湖北公安人，工程師，博士，從事核電廠熱工水力安全分析及安全殼方向的設計研究

TL333

A

1000-6931(2015)06-1045-06

10.7538/yzk.2015.49.06.1045