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

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

?

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

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

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

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

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

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

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

1 數(shù)學(xué)物理模型

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

質(zhì)量守恒：

(1)

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

動量守恒：

(2)

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

能量方程：

(3)

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

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

(4)

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

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

(5)

(6)

(7)

(8)

系數(shù)矩陣G具有3個實特征值：

(9)

每個特征值對應(yīng)的特征向量所構(gòu)成的矩陣稱為特征矩陣，左特征矩陣和右特征矩陣如下：

(10)

(11)

2 隱式矩陣分裂算法

矩陣分裂算法是一種有限差分方法，根據(jù)系統(tǒng)特征值的正負，確定對流項差分格式。它最早運用于氣體動力學(xué)，Romstedt[7]首次將該方法拓展到兩相流體動力學(xué)中。本文對均相流模型下矩陣分裂算法進行了詳細推導(dǎo)，基于特征值正負方向?qū)Ψ至押蟮姆匠毯驮错椷M行離散，并運用矩陣形式追趕法對離散后方程組進行求解。

(12)

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

(13)

(14)

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

(15)

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

(16)

寫成分量形式有：

(17)

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

當(dāng)λp≥0時，有：

(18)

當(dāng)λp<0時，有：

(19)

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

(20)

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

(21)

即：

(22)

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

(23)

有關(guān)系式：

(24)

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

(25)

當(dāng)流體控制方程含有源項時，其源項也應(yīng)為迎風(fēng)方向控制體內(nèi)的源項。

當(dāng)λp≥0時，有：

(26)

當(dāng)λp<0時，有：

(27)

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

(28)

Si=R[l1Dil2Di+1l3Di]T

(29)

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

(30)

圖1 網(wǎng)格示意圖Fig.1 Schematic diagram of grid

化簡后有：

(31)

合并后有：

(32)

(33)

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

(34)

系數(shù)矩陣可表示為：

(35)

對矩陣K進行LR有：

K=LR

根據(jù)矩陣乘法，可得到追趕過程遞推關(guān)系式，追的過程為：

L1=B1

(36)

(37)

趕的過程為：

UN=YN

(38)

(39)

3 算法驗證

噴管是熱力系統(tǒng)中重要的工程設(shè)備，在工程熱力系統(tǒng)和氣體動力學(xué)中均有重要應(yīng)用，其結(jié)構(gòu)如圖2所示，其結(jié)構(gòu)和工作原理詳見文獻[9-10]。

圖2 漸縮漸擴噴管結(jié)構(gòu)示意圖Fig.2 Schematic diagram of structure for convergent-divergent nozzle

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

(40)

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

不同控制體下噴管喉部參數(shù)列于表1?？刂企w數(shù)增加到50后，計算結(jié)果基本不再變化，故工況中噴管取50個控制體。

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

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

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

4 應(yīng)用舉例

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

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

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

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

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

5 結(jié)束語

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

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

[2] 陶文銓. 計算傳熱學(xué)的近代進展[M]. 北京:科學(xué)出版社,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]. 哈爾濱:哈爾濱工程大學(xué)出版社,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]. 西安:西安交通大學(xué)出版社,2001.

[9] JOHN D A. 計算流體力學(xué)基礎(chǔ)及其應(yīng)用[M]. 北京:機械工業(yè)出版社,2007.

[10]沈維道,蔣志敏,童鈞耕. 工程熱力學(xué)[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—)，男，湖北公安人，工程師，博士，從事核電廠熱工水力安全分析及安全殼方向的設(shè)計研究

TL333

A

1000-6931(2015)06-1045-06

10.7538/yzk.2015.49.06.1045