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

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

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

?

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

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

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

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

矩陣分裂算法；塊三對(duì)角矩陣；兩相流

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

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

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

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

質(zhì)量守恒：

(1)

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

動(dòng)量守恒：

(2)

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

能量方程：

(3)

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

能量方程寫(xiě)成熵的形式，有：

(4)

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

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

(5)

(6)

(7)

(8)

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

(9)

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

(10)

(11)

2 隱式矩陣分裂算法

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

(12)

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

(13)

(14)

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

(15)

小時(shí)間步長(zhǎng)下可將方程線性化，令ΔW=R-1ΔU，有：

(16)

寫(xiě)成分量形式有：

(17)

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

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

(18)

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

(19)

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

(20)

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

(21)

即：

(22)

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

(23)

有關(guān)系式：

(24)

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

(25)

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

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

(26)

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

(27)

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

(28)

Si=R[l1Dil2Di+1l3Di]T

(29)

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

(30)

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

化簡(jiǎn)后有：

(31)

合并后有：

(32)

(33)

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

(34)

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

(35)

對(duì)矩陣K進(jìn)行LR有：

K=LR

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

L1=B1

(36)

(37)

趕的過(guò)程為：

UN=YN

(38)

(39)

3 算法驗(yàn)證

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

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

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

(40)

噴管內(nèi)氣體采用氮?dú)?，駐室內(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后，計(jì)算結(jié)果基本不再變化，故工況中噴管取50個(gè)控制體。

噴管內(nèi)Ma、壓力和溫度的分布示于圖3。氮?dú)釳a沿噴管逐漸增加，在噴管喉部處等于1，速度增加到當(dāng)?shù)芈曀佟：聿壳皾u縮噴管內(nèi)氣體作亞聲速流動(dòng)，喉部后漸擴(kuò)噴管內(nèi)氣體作超聲速流動(dòng)。氮?dú)鈮毫ρ貒姽苤饾u減小，內(nèi)能逐漸轉(zhuǎn)化為氣體動(dòng)能，氣體溫度也逐漸降低。圖3表明，數(shù)值計(jì)算結(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ì)量流速波動(dòng)對(duì)圓管CHF的影響，研究對(duì)象如圖4所示，壁面采用均勻熱流密度加熱，入口流體質(zhì)量流速正弦規(guī)律波動(dòng)。采用一維三流體模型研究圓管內(nèi)流體流動(dòng)，非搖擺條件下計(jì)算結(jié)果如圖5所示，可看出加熱情況下流體質(zhì)量流速沿流動(dòng)方向波動(dòng)振幅不斷衰減。

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

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

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

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

5 結(jié)束語(yǔ)

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

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

[2] 陶文銓. 計(jì)算傳熱學(xué)的近代進(jìn)展[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] 鄧建中, 劉之行. 計(jì)算方法[M]. 西安:西安交通大學(xué)出版社,2001.

[9] JOHN D A. 計(jì)算流體力學(xué)基礎(chǔ)及其應(yīng)用[M]. 北京:機(jī)械工業(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è)計(jì)研究

TL333

A

1000-6931(2015)06-1045-06

10.7538/yzk.2015.49.06.1045