張 璜，薄涵亮

(清華大學(xué) 核能與新能源技術(shù)研究院，北京 100084)

液滴運(yùn)動(dòng)現(xiàn)象廣泛存在于各種核能動(dòng)力裝置中。例如，蒸汽發(fā)生器內(nèi)的汽水分離裝置中就存在大量液滴運(yùn)動(dòng)。汽水分離裝置的作用是將蒸汽發(fā)生器產(chǎn)生的濕蒸汽流進(jìn)行干燥。歐美各國(guó)規(guī)定經(jīng)汽水分離器干燥后的蒸汽濕度要小于0.25%，否則高速蒸汽流夾帶著這些液滴會(huì)造成二回路管道、關(guān)閉件腐蝕，以及汽輪機(jī)的汽蝕，影響汽輪機(jī)壽命。因此，汽水分離裝置的分離效率關(guān)系到核電系統(tǒng)運(yùn)行的經(jīng)濟(jì)性和安全性。

常見(jiàn)的汽水分離裝置有旋葉式分離器、重力式分離器和波形板式分離器，它們內(nèi)部的液滴在蒸汽場(chǎng)中運(yùn)動(dòng)，流出或被分離器捕獲。通過(guò)建立汽水分離器內(nèi)液滴運(yùn)動(dòng)模型并對(duì)其進(jìn)行求解，可得到液滴在分離器內(nèi)部的運(yùn)動(dòng)軌跡和數(shù)目分布，進(jìn)而能較為細(xì)致地研究分離器對(duì)液滴的分離機(jī)理，從而為設(shè)計(jì)高性能的汽水分離器提供理論依據(jù)。

不同學(xué)者采用的液滴運(yùn)動(dòng)模型雖不同，但該模型均可表述為1組一階常微分方程。由于該組常微分方程通常無(wú)解析解，所以一般采用數(shù)值方法對(duì)其進(jìn)行求解。龐鳳閣等[1]采用高階方程的Runge-Kutta(RK)算法對(duì)波形板內(nèi)的液滴運(yùn)動(dòng)模型進(jìn)行求解，得到了直徑為7～10 μm的液滴運(yùn)動(dòng)軌跡。陳韶華等[2-3]采用定(單)步長(zhǎng)的RK算法求解液滴運(yùn)動(dòng)模型，對(duì)重力分離空間內(nèi)直徑為0～350 μm的液滴和波形板分離器內(nèi)直徑為10.5～12.5 μm的液滴的運(yùn)動(dòng)行為均作了研究。王曉墨等[4]采用向前Euler算法求解液滴運(yùn)動(dòng)模型，得到了直徑為10、20和900 μm的液滴在波形板內(nèi)的運(yùn)動(dòng)軌跡。張謹(jǐn)奕等[5-6]考慮了液滴在三維空間中受到的Magnus升力和Saffman升力，采用4級(jí)4階顯式單步長(zhǎng)RK算法對(duì)該模型進(jìn)行數(shù)值求解，但只對(duì)波形板汽水分離器內(nèi)直徑大于50 μm的液滴的運(yùn)動(dòng)軌跡和分離效率進(jìn)行了研究。

向前Euler算法雖易編程實(shí)施，但它是顯式算法且僅有一階精度，穩(wěn)定性和精度均不能保證[7]。單步長(zhǎng)RK算法雖能保證數(shù)值求解結(jié)果具有較高精度，但它依然是顯式算法，數(shù)值穩(wěn)定性不夠好[7]。

針對(duì)以上不足，本文擬提出三維空間內(nèi)液滴運(yùn)動(dòng)模型數(shù)值求解的一種算法，即算子分裂算法(OS算法)。本文首先簡(jiǎn)述適用于三維空間的液滴運(yùn)動(dòng)模型，隨后詳細(xì)介紹OS算法的數(shù)值離散格式，最后分別采用RK算法和OS算法對(duì)波形板汽水分離器內(nèi)的液滴運(yùn)動(dòng)模型進(jìn)行求解，從相容性、穩(wěn)定性、計(jì)算精度和計(jì)算效率等方面對(duì)這兩種算法進(jìn)行比較。

1 三維空間液滴運(yùn)動(dòng)模型的數(shù)學(xué)表述

在三維空間中的單個(gè)液滴，不僅有相對(duì)于固定坐標(biāo)系的平動(dòng)，還有繞自身對(duì)稱(chēng)軸的轉(zhuǎn)動(dòng)和液滴表面的變形運(yùn)動(dòng)。在上述運(yùn)動(dòng)的同時(shí)，除了受到流場(chǎng)的作用力，還可能會(huì)對(duì)流場(chǎng)的變化產(chǎn)生影響。

因此，考慮到三維空間液滴運(yùn)動(dòng)的復(fù)雜性，本文提出如下假設(shè)：

1) 液滴和流場(chǎng)間無(wú)質(zhì)量和熱量交換；

2) 液滴為剛性圓球，在運(yùn)動(dòng)過(guò)程中不發(fā)生形變，因此液滴可看作剛體，其運(yùn)動(dòng)可分解為平動(dòng)和轉(zhuǎn)動(dòng)；

3) 流場(chǎng)對(duì)液滴有作用力，而液滴在運(yùn)動(dòng)過(guò)程中對(duì)流場(chǎng)無(wú)影響。

基于以上假設(shè)，液滴的轉(zhuǎn)動(dòng)方程為：

(1)

液滴運(yùn)動(dòng)方程為：

(2)

其中：md=πρdd3/6，ρd為液滴密度；vd為液滴速度；液滴所受曵力FD=πCDρf|vf-vd|(vf-vd)d2/8，CD為曵力系數(shù)；液滴所受慣性質(zhì)量力FA=πρfd3[d(vf-vd)/dt]/12；液滴所受體積力FV=πd3(ρd-ρf)g/6，g為重力加速度；液滴所受Magnus升力FM=πCMaρfd3(vf-vd)×(ωd-ωf/2)/8，CMa為Magnus升力系數(shù)；液滴所受Saffman升力FS=1.615CSa(Rμ)2(ρfμf)0.5d2×|ωf|-0.5[(vf-vd)×ωf]，CSa為Saffman升力系數(shù)，(Rμ)2表示液滴內(nèi)部環(huán)流量對(duì)升力的影響，μf為流體動(dòng)力黏性系數(shù)。

將以上各表達(dá)式代入式(1)和(2)，并將方程左邊系數(shù)歸一化得：

(3)

其中：λ1=-15ρf/16π、λ2=3ρf/(4ρdd+2ρfd)、λ3=3ρf/(4ρd+2ρf)、λ4=[1.615(μd+2μf/3)2/(μd+μf)2(μfρf)0.5]/(ρdπd/6+ρfπd/12)、λ5=2(ρd-ρf)/(2ρd+ρf)為歸一化系數(shù)；CM、CD、CMa、CSa等系數(shù)的值詳見(jiàn)文獻(xiàn)[5]。

式(3)即為本文需求解的三維空間液滴運(yùn)動(dòng)模型的數(shù)學(xué)表述式。

2 OS算法的離散格式

算子LA為：

(4)

算子LB為：

(5)

LA求解液滴的角速度變化和由曵力加速度引起的液滴速度變化；LB求解由Magnus升力、Saffman升力加速度和體積力加速度引起的液滴速度變化。這樣就將求解式(3)的過(guò)程轉(zhuǎn)換成分別求解式(4)和(5)，這一步稱(chēng)為算子分裂。

假設(shè)時(shí)間步長(zhǎng)τ=ti+1-ti，將(ωd，vd)記為Yd，從時(shí)刻ti到時(shí)刻ti+1，OS算法實(shí)施步驟為：

Yi+1=LAi(τ/2)LBi(τ)LAi(τ/2)·Yi

(6)

其中，LAi和LBi為離散的LA和LB算子。

OS算法需將式(6)表示的過(guò)程對(duì)時(shí)間離散。本文采用隱式算法求解算子LA，選取具有二階精度的隱式梯形算法[7]離散式(4)：

(7)

同時(shí)采用RK算法離散式(5)，得：

(8)

引入定義：

(9)

那么，K1、K2、K3、K4為：

(10)

OS算法針對(duì)式(3)的離散格式為式(7)和(8)，對(duì)每一個(gè)時(shí)間步長(zhǎng)的計(jì)算步驟為式(6)。

3 OS算法的驗(yàn)證

本文將RK算法[5]也用于求解式(3)，以模擬直徑為1～50 μm的液滴在波形板中的運(yùn)動(dòng)軌跡，從兩種算法的相容性、穩(wěn)定性、計(jì)算精度和計(jì)算效率等方面來(lái)驗(yàn)證OS算法的性能。

本算例中，采用無(wú)鉤型波形板[1]，其幾何形狀示于圖1a。波形板的波數(shù)為2，入口及出口端長(zhǎng)度均為L(zhǎng)=10 mm，波形板間距H=12.5 mm，曲折角α=45°，折邊B=25 mm，并定義w=B/H。假設(shè)波形板在冷態(tài)工況下工作，用不可壓縮空氣代替蒸汽[1]，其密度ρf=1.225 kg/m3，動(dòng)力黏性系數(shù)μf=1.789 4×10-5Pa·s；液滴密度ρd=998 kg/m3，動(dòng)力黏性系數(shù)μd=9.98×10-4Pa·s。波形板內(nèi)液滴和氣體的運(yùn)動(dòng)可簡(jiǎn)化為二維，所以可忽略液滴受到的體積力，可令式(3)中λ5=0。

波形板內(nèi)氣體的二維流動(dòng)滿(mǎn)足Navier-Stokes方程，其連續(xù)性方程和動(dòng)量守恒方程為：

(11)

(12)

其中，υf=μf/ρf為流體運(yùn)動(dòng)黏性系數(shù)。利用FLUENT 6.3.26對(duì)波形板內(nèi)氣相流場(chǎng)進(jìn)行模擬，流動(dòng)設(shè)置為穩(wěn)態(tài)，入口速度分別設(shè)為3、5和7 m/s，黏性模式選用標(biāo)準(zhǔn)κ-ε湍流模型，固壁設(shè)置為無(wú)滑移邊界，出口設(shè)置為充分發(fā)展流動(dòng)，選用SIMPLE算法求解，各項(xiàng)殘差設(shè)置為10-4 [8]。當(dāng)入口速度為3 m/s時(shí)，計(jì)算結(jié)果采用Tecplot處理得到的波形板內(nèi)流場(chǎng)的流線(xiàn)示于圖1b。

為了模擬液滴在波形板內(nèi)流場(chǎng)中的運(yùn)動(dòng)軌跡，需將FLUENT 6.3.26計(jì)算所得的流場(chǎng)速度和速度梯度保存為輸出文件，作為求解液滴運(yùn)動(dòng)模型的已知流場(chǎng)速度和旋度。

3.1 OS算法相容性分析

隱式梯形法的局部截?cái)嗾`差為o(τ2)[7]，RK算法的局部截?cái)嗾`差為o(τ5)[7]。

采用隱式梯形法和RK算法構(gòu)造的OS算法的局部截?cái)嗾`差階數(shù)推導(dǎo)過(guò)程為：采用OS算法求解式(3)，當(dāng)計(jì)算到第i步時(shí)，液滴的角速度和速度可寫(xiě)為yi=(ωdi，vdi)，采用OS算法計(jì)算得到第i+1步的數(shù)值解為：

yi+1=[(yieLAτ/2+o(τ2))eLBτ+

o(τ5)]eLAτ/2+o(τ2)

(13)

利用式(3)得到第i+1步的精確解為：

(14)

式(14)減去式(13)所得結(jié)果即為OS算法的局部截?cái)嗾`差，最終有：

ο(τ2)

(15)

由式(15)可知，OS算法構(gòu)造的離散格式的局部截?cái)嗾`差為o(τ2)，與式(3)相容。

3.2 OS算法穩(wěn)定性分析和RK算法時(shí)間步長(zhǎng)比較

提出一個(gè)算法，首先必須保證它構(gòu)造的離散格式和原方程相容。除此之外，更重要的是保證算法的收斂性。證明算法的收斂性往往非常困難[9]，但只要對(duì)實(shí)際問(wèn)題建立合適的物理與數(shù)學(xué)模型，則假定由相容與穩(wěn)定的離散格式就能獲得收斂的數(shù)值解[10]。

由于所求解的式(3)為非線(xiàn)性常微分方程組，因此本文僅討論RK算法和OS算法的局部初值穩(wěn)定性。在本文計(jì)算的波形板流場(chǎng)中，直徑小于50 μm的液滴的平動(dòng)雷諾數(shù)Red(Red=ρf|vf-vd|d/μf)和轉(zhuǎn)動(dòng)雷諾數(shù)Reω(Reω=ρf|ωd-ωf/2|d2/4μf)總體上分別小于6.2和1，且隨著直徑的減小，Red和Reω更加趨近于滿(mǎn)足上述要求。此時(shí)可取CM=16π/Reω、CD=24/Red、CMa=1、CSa=1，在二維情況下，式(3)變?yōu)椋?/p>

(16)

其中：α=4μf/d2；β=36μf/(2ρd+ρf)d2；γ=2ρf/(4ρd+2ρf)；η=19.38(μfρf)0.5/πd(2ρd+ρf)。本文提出的OS算法主要針對(duì)直徑小于50 μm的液滴，所以比較RK算法和OS算法穩(wěn)定性時(shí)，可直接針對(duì)式(16)進(jìn)行分析。

式(16)的系數(shù)矩陣Gcreep為：

圖1 波形板的幾何形狀(a)和內(nèi)部流場(chǎng)流線(xiàn)(b)

(17)

其中：ωd=(0，0，ωd)；ωf=(0，0，ωf)；vd=(vdx，vdy，0)；vf=(vfx，vfy，0)。系數(shù)矩陣Gcreep包含要求解的未知數(shù)ωd、vdx和vdy，采用量綱分析估計(jì)Gcreep的大小。本文中，液滴速度和流場(chǎng)速度的量級(jí)為10 m/s，液滴角速度的量級(jí)為102rad/s，流場(chǎng)旋度的量級(jí)為103rad/s，液滴直徑的量級(jí)為10-5m，結(jié)合已給出的液滴和流場(chǎng)的物性參數(shù)，系數(shù)矩陣Gcreep約為：

(18)

Gcreep特征值的實(shí)部為ζ1=-1.08×107、ζ2=ζ3=-3.24×103，可見(jiàn)式(16)是剛性方程。

用RK算法求解式(3)時(shí)，根據(jù)絕對(duì)穩(wěn)定性條件[11]理論預(yù)測(cè)出計(jì)算不同直徑液滴所需的最大時(shí)間步長(zhǎng)τmaxTRK=-2.78/ζmin，ζmin為系數(shù)矩陣Gcreep的最小特征值(若特征值為復(fù)數(shù)，取特征值的實(shí)部)。通過(guò)數(shù)值實(shí)驗(yàn)，得到不同直徑液滴情況下保持RK算法穩(wěn)定的最大時(shí)間步長(zhǎng)τmaxNRK，將兩者同時(shí)列于表1。由表1可見(jiàn)，理論值和數(shù)值實(shí)驗(yàn)值很吻合。

OS算法的實(shí)施過(guò)程見(jiàn)式(6)，ti+1時(shí)刻，yi+1=LAi(τ/2)LBi(τ)LAi(τ/2)·yi，令GOS=LAi(τ/2)LBi(τ)LAi(τ/2)，則GOS的具體表達(dá)式為：

(19)

其中，χ=γ(ωf/2-ωd)-ηωf/|ωf|0.5。GOS的特征值ζ1=(1+ατ/2)2/(1-ατ/2)2，ζ2=ζ3=(1-βτ/2)2(1-τ2χ2/2+τ2χ4/24)/(1+βτ/2)2，要使OS算法穩(wěn)定，只需保證GOS特征值的實(shí)部的絕對(duì)值小于1。已知α<0，β>0，那么必有|ζ1|<1，|ζ2|=|ζ3|=((1-βτ/2)2/(1+βτ/2)2)|1-τ2χ2/2+τ2χ4/24|<|1-τ2χ2/2+τ2χ4/24|，因此只需滿(mǎn)足|1-τ2χ2/2+τ2χ4/24|<1就能保證OS算法穩(wěn)定，其解為：

(20)

令τTOS=(12/χ2)0.5，不同液滴直徑下的τTOS列于表1。由表1可知，τTOS的量級(jí)為10-1～10-3，而τmaxTRK和τmaxNRK的量級(jí)為10-6～10-9。值得注意的是，τTOS僅是保證OS算法穩(wěn)定的充分條件，所以使OS算法穩(wěn)定的最大時(shí)間步長(zhǎng)τmaxOS≥τTOS。通過(guò)數(shù)值實(shí)驗(yàn)也發(fā)現(xiàn)，對(duì)于不同直徑的液滴，τmaxOS均可放寬到1 s，因此也證實(shí)了上述觀(guān)點(diǎn)。通過(guò)以上分析可知，OS算法可取的時(shí)間步長(zhǎng)大于RK算法時(shí)間步長(zhǎng)的取值，因此OS算法的穩(wěn)定性?xún)?yōu)于RK算法。

3.3 OS算法計(jì)算精度和計(jì)算效率分析

由3.1節(jié)可知，RK算法的局部截?cái)嗾`差為o(τ5)，而OS算法的局部截?cái)嗾`差為o(τ2)，為保證OS算法依然具有較好的精度，對(duì)其時(shí)間步長(zhǎng)選取如下：

(21)

RK算法的時(shí)間步長(zhǎng)τRK取表1中給出的τmaxNRK，注意此時(shí)τOS仍比τRK大。

比較了兩種算法計(jì)算得出的不同直徑液滴在波形板內(nèi)的軌跡形狀。兩種算法中，液滴進(jìn)入波形板的速度和位置均相同。液滴進(jìn)入波形板的速度與波形板內(nèi)流場(chǎng)入口速度相等，液滴入射方向垂直于波形板入口，液滴初始位置在波形板入口處均勻分布。同一直徑的液滴計(jì)算10個(gè)。圖2為4種不同直徑液滴以不同初始速度在波形板內(nèi)的運(yùn)動(dòng)軌跡示意圖?？梢?jiàn)，兩種算法得到的軌跡形狀相似。

表1 τmaxTRK、τmaxNRK和τTOS的比較

記RK算法計(jì)算得到液滴恰好流出或被波形板捕獲的位置(終端位置)為(xRK，yRK)，OS算法得到液滴的終端位置為(xOS，yOS)。為了定量比較兩種算法所得液滴終端位置的精度，分別定義x軸和y軸終端位置的相對(duì)誤差為εx和εy：

(22)

兩種算法計(jì)算所得不同直徑液滴終端位置平均值及其相對(duì)誤差列于表2～4。由表2～4可看出，εx在0.6%以?xún)?nèi)，εy在2.7%以?xún)?nèi)，可見(jiàn)OS算法和RK算法精度相當(dāng)。

在計(jì)算表2的同時(shí)，記錄每個(gè)液滴在波形板內(nèi)經(jīng)過(guò)完整軌跡所需的計(jì)算時(shí)間(取計(jì)算機(jī)的CPU耗時(shí))并取平均值，定義RK算法計(jì)算單個(gè)液滴的平均CPU耗時(shí)除以O(shè)S算法的平均CPU耗時(shí)為CPU耗時(shí)倍數(shù)，最后所得結(jié)果列于表5。

由表5可見(jiàn)，3種工況的平均CPU耗時(shí)倍數(shù)均大于35。對(duì)于直徑相對(duì)較小的液滴，如1 μm液滴，CPU耗時(shí)倍數(shù)可達(dá)300以上。因此，表2～5結(jié)果可證實(shí)，在保證OS算法和RK算法精度相當(dāng)?shù)那疤嵯?，OS算法的計(jì)算效率優(yōu)于RK算法。

液滴初始速度：a——3 m/s；b——5 m/s；c——7 m/s

表2 液滴初始速度為3 m/s時(shí)兩種算法計(jì)算所得不同直徑液滴終端位置平均值及其相對(duì)誤差

表3 液滴初始速度為5 m/s時(shí)兩種算法計(jì)算所得不同直徑液滴終端位置平均值及其相對(duì)誤差

表4 液滴初始速度為7 m/s時(shí)兩種算法計(jì)算所得不同直徑液滴終端位置平均值及其相對(duì)誤差

表5 兩種算法對(duì)不同直徑單個(gè)液滴的平均CPU耗時(shí)及CPU耗時(shí)倍數(shù)

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

本文針對(duì)三維空間液滴運(yùn)動(dòng)模型提出了一種新的數(shù)值求解方法——OS算法，并采用該算法求解波形板內(nèi)液滴運(yùn)動(dòng)軌跡，同時(shí)與RK算法所得結(jié)果進(jìn)行了對(duì)比。得到主要結(jié)論如下：

1) 本文構(gòu)造的OS算法具有局部二階精度；

2) 保持OS算法穩(wěn)定的最大時(shí)間步長(zhǎng)的量級(jí)可取10-1～10-3s，而保持RK算法穩(wěn)定的最大時(shí)間步長(zhǎng)的量級(jí)為10-6～10-9s，因此OS算法的穩(wěn)定性?xún)?yōu)于RK算法；

3) 在計(jì)算精度相同的情況下，OS算法的計(jì)算效率優(yōu)于RK算法。

因此OS算法的提出為數(shù)值求解三維空間液滴運(yùn)動(dòng)模型，進(jìn)而計(jì)算波形板的分離效率提供了一種相對(duì)穩(wěn)定、高效的算法。

參考文獻(xiàn)：

[1] 龐鳳閣，于瑞俠，張志儉. 波形板汽水分離器的機(jī)理研究[J]. 核動(dòng)力工程，1992，13(3)：9-13.

PANG Fengge, YU Ruixia, ZHANG Zhijian. The study of mechanism of corrugated baffles separator[J]. Nuclear Power Engineering, 1992, 13(3): 9-13(in Chinese).

[2] 陳韶華，黃素逸，趙緒新. PWR蒸汽發(fā)生器水滴重力分離研究[J]. 華中理工大學(xué)學(xué)報(bào)：自然科學(xué)版，1997，25(1)：69-71.

CHEN Shaohua, HUANG Suyi, ZHAO Xuxin. A study on the gravity separation of water droplets in steam generators for PWR[J]. Journal of Huazhong University of Science and Technology: Nature Science Edition, 1997, 25(1): 69-71(in Chinese).

[3] 陳韶華，黃素逸，趙緒新. 波形板汽水分離器汽水兩相分離機(jī)理研究[J]. 華中理工大學(xué)學(xué)報(bào)：自然科學(xué)版，1998，26(1)：5-7.

CHEN Shaohua, HUANG Suyi, ZHAO Xuxin. On separation mechanism of steam water two phase flows in a steam water separator with corrugated baffle[J]. Journal of Huazhong University of Science and Technology: Nature Science Edition, 1998, 26(1): 5-7(in Chinese).

[4] 王曉墨，黃素逸. 波形板汽水分離器中液滴軌跡的數(shù)值模擬[J]. 核動(dòng)力工程，2003，24(6)：582-585.

WANG Xiaomo, HUANG Suyi. Numerical simulation of droplets track in separator with wave-shape vanes[J]. Nuclear Power Engineering, 2003, 24(6): 582-585(in Chinese).

[5] 張謹(jǐn)奕，薄涵亮，孫玉良，等. 三維空間液滴運(yùn)動(dòng)模型[J]. 清華大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，53(1)：96-101.

ZHANG Jinyi, BO Hanliang, SUN Yuliang, et al. Three-dimensional droplet motion model[J]. Journal of Tsinghua University: Science and Technology, 2013, 53(1): 96-101(in Chinese).

[6] 張謹(jǐn)奕，薄涵亮. 液滴模型在波形板汽水分離器中的應(yīng)用[J]. 原子能科學(xué)技術(shù)，2012，46(增刊)：776-781.

ZHANG Jinyi, BO Hanliang. Application of droplet model in wave-type vane steam separator[J]. Atomic Energy Science and Technology, 2012, 46(Suppl.): 776-781(in Chinese).

[7] 李大侃. 常微分方程數(shù)值解[M]. 杭州：浙江大學(xué)出版社，1994：21-42.

[8] 李嘉. 波形板汽水分離器的理論和實(shí)驗(yàn)研究[D]. 武漢：華中科技大學(xué)，2007.

[9] 陸金甫，關(guān)治. 偏微分方程數(shù)值解法[M]. 2版. 北京：清華大學(xué)出版社，2004：19-28.

[10] 陶文銓. 數(shù)值傳熱學(xué)[M]. 2版. 西安：西安交通大學(xué)出版社，2001：48-52.

[11] 李慶揚(yáng)，關(guān)治，白峰杉. 數(shù)值計(jì)算原理[M]. 北京：清華大學(xué)出版社，2000：372-376.