大平板上液膜沖擊附板行為的數值模擬研究

胡 真，胡 珀，楊小杰，杜卡帥

(1.上海交通大學 核科學與工程學院，上海 200240；2.上海核工程研究設計院有限公司，上海 200233)

在第3代自主先進壓水堆CAP1400的設計中，非能動余熱排出(PRHR)系統在事故狀態下的長期冷卻過程中將發揮至關重要的作用。該系統中的安全殼內置換料水箱(IRWST)在此過程中將蒸發產生的水蒸氣通過冷凝回流系統形成冷凝液并最終返回IRWST，可保證IRWST的水量可利用率，從而支持一回路的長期冷卻。在此過程中，安全殼垂直壁面上的附板和返回槽會造成一定量的冷凝液損失，而該損失將直接影響冷凝回流的份額，因此預測冷凝回流的損失量對于PRHR長期冷卻分析具有重要意義。

降液膜在安全殼外部壁面上的研究目前已有顯著成效[1-7]，但關于降液膜在安全殼內壁面上沖擊附板時的流動行為卻鮮有研究，而對于流體沖擊板面行為的研究也多為研究單個液滴的濺射和鋪展行為[8-9]，因此當前國內外尚無可行的計算方法來有效估算該損失。商用計算軟件FLUENT中的EWF模型(歐拉壁面模型)具有計算效率高、能很好模擬液膜在壁面的流動等優勢[4,10]，在模擬全尺寸的模型時計算成本相對較小，故本文將以EWF模型為工具，結合相關實驗數據，通過構建合理的數值模型，研究大平板上液膜沖擊附板時的濺射行為，旨在通過研究液膜沖擊附板的過程為冷凝回流損失率做出相關預測。

1 數學模型

1.1 模型建立及邊界條件

圖1 WAFT實驗臺架[2] Fig.1 WAFT experiment bench[2]

液膜實驗裝置(WAFT)如圖1所示。實驗環境為常溫常壓，實驗所用平板涂層為Carbozinc 11 HS有機富鋅涂料。實驗過程中液膜分配箱將水均勻分配到平板上端，水以漫溢的方式沿平板平鋪而下，沖擊附板，發生濺射后的部分水流經返回槽收集后進行稱重以評估實際的損失。

實驗所用實驗板及其焊接方式如圖2所示，板長5 m、寬1.2 m，板端末設有高0.1 m的返回槽，距返回槽2 m處置有附板。

圖2 實驗板設置及附板焊接示意圖Fig.2 Plate setting and welding diagram of attached plate

根據實驗要求選用不同厚度和寬度的附板，附板通過角焊的方式固定在板面上，焊腳高度為l，附板高度為H，焊腳高比β定義為：

β=l/H

(1)

在實驗過程中，通過控制變量法改變液膜分配箱流出水的流量、溫度、附板尺寸以及焊腳高比來分析各因素對液膜濺射的影響。

在數值研究過程中，選擇實驗板及其上方0.1 m高處氣-液流動的三維空間區域作為計算域，模擬計算中各尺寸采用與實驗1∶1的比例進行建模計算，模型切面示于圖3。

圖3 三維模型切面示意圖Fig.3 Section diagram of 3-D model

具體模型尺寸列于表1。由于實驗板上方未加玻璃罩遮擋，故將上方區域邊界設為壓力出口。在計算區域內，液膜由上方入口處流入后沿壁面鋪展開，待其發展充分后沖擊位于板面下部的附板，飛濺損失的液體由上方壓力出口射出，其余液體通過高0.1 m的返回槽流出。具體邊界條件設置列于表2。

表1 三維模型的幾何尺寸Table 1 Geometry dimension of 3-D model

表2 邊界條件Table 2 Boundary condition

為方便計算，本文在建模過程中做了如下假設：1) 計算過程中不考慮液膜的波動特性；2) 不考慮液膜與空氣之間的質量能量交換，即不考慮液膜在沖擊附板和流動過程中的蒸發損失；3) 不考慮實驗中返回槽上的濾網損失以及鋼板失配影響。

1.2 控制方程

由于本文只研究大平板上降液膜的流動和濺射行為，因此傳熱情況和能量方程不在考慮之中。假設流體不可壓縮，其穩態下連續性方程為：

?(ρui)/?xi=Smax

(2)

式中：ρ為流體密度；ui為流體在xi方向的速度；Smax為源項。

動量方程可描述為：

(3)

式中：i,j為卡笛爾坐標系中的不同坐標方向；u為速度，m/s；p為壓力，Pa；μ為動力黏度，Pa·s；δ為Kronecker函數；Smon-i為源項；ρu′iu′j為基于Boussinesq假設的雷諾應力。

本文采用Realizablek-ε湍流模型描述液膜的湍流流動，該湍流模型的模擬輸運方程k和ε方程分別如式(4)、(5)所示。

Gk+Gb-ρε+Sk

(4)

式中：k為湍流動能；Gk為由主流速度梯度產生的湍流動能；Gb為由浮升力產生的湍流動能；σk為k方程的普朗特數，σk=1.0；ε為湍流耗散率；Sk為用戶自定義的源項。

(5)

EWF模型可用來預測壁面上薄液膜的產生和流動情況，此外，該模型還可用來模擬液滴飛濺、剝離、液膜分流等現象。該模型適用于液膜厚度較薄、膜內液體流動與壁面平行的情況[12]。

液膜的質量守恒方程為：

(6)

液膜的連續性方程為：

(7)

式(7)右邊第1項表示動壓氣體壓力在壁面方向上的重力分量以及液體表面張力對液膜的影響，第2～4項分別表示在液膜流動方向上的重力作用、氣膜交界面處的黏性剪切力作用和膜內的黏性作用。

本文所用的液膜分離模型[13-14]為O’Rouke模型。FLUENT中對液膜分離的判定準則為：

(8)

在底板和附板的邊界設置中均需開啟EWF模型，其中附板和焊縫的DPM模型選項設置為reflect，而底板的該選項設置為trap，確保液膜在沖擊附板前貼板流動。

本文模型選用基于壓力的穩態求解器，由于底面設有附板，為更好地實現計算收斂，選用標準近壁面模型進行計算。

1.3 網格無關性分析

考慮到計算的準確性，采用結構化網格劃分的方式對計算區域進行網格劃分。由于應用了標準近壁面模型，貼近實驗板的第1層網格的y+值保持在20左右，其他壁面附近的第1層網格的y+值保持在20以上。本文選取24萬到96萬不同數目的網格模型，對每10 min的損失水量進行檢測，網格無關性驗證結果如圖4所示。為綜合計算效率和計算準確性，本文最終選用48萬網格的模型。

圖4 網格無關性驗證Fig.4 Grid independence validation

2 結果與分析

2.1 冷凝回流損失比較

圖5 冷凝回流損失結果比較Fig.5 Comparison of condensation reflux loss results

為驗證數值模型在預測冷凝回流損失結果上的正確性，本文選取焊腳高比為1.0時不同工況下的模擬計算值與實驗值進行比較。圖5為不同流量的降液膜沖擊不同厚度和長度附板的計算值與實驗值對比。從圖5可明顯看出，小流量下液膜沖擊附板后損失率很小，且計算值與實驗值極為接近；中流量(500 m3/h)和大流量(1 000 m3/h)情況下，液膜損失率顯著增大，且計算值和實驗值間的誤差也有明顯提高，個別工況的計算值與實驗值間的誤差超過5%?？傮w而言，通過本文所建數值模型得到的冷凝回流損失量與實驗誤差能控制在8%以內，即本文所建數值模型能較好地預測冷凝回流損失量。

2.2 降液膜沖擊附板時的濺射行為分析

降液膜沖擊附板的整個行為過程如圖6所示，降液膜在平板上經過充分發展，沖擊附板后的行為過程可分為4個階段，即液膜沖擊附板階段、第1次液膜再聚合階段、液膜分離階段、第2次液膜再聚合階段。

液膜厚度和液膜流動的速度場如圖7所示，圖中前3階段液膜行為的速度場速度范圍均為0～1 m/s，第2次液膜再聚合階段的速度場速度范圍為0～2 m/s，且速度方向均為沿板面方向。

在液膜沖擊附板階段，液膜沖擊附板時，其流動方向發生改變，并沿焊縫方向繼續流動。在液膜沖擊附板后第1次發生分離時，液膜主要分成3個方向流動：根據無滑移條件假設，黏性底層因與壁面保持相對靜止狀態，因此仍緊貼壁面流動，湍流核心區的液體受表面張力和液體內部作用力影響，大部分液體在沿焊縫方向流動一定距離后繼續回到附板表面貼壁面流動，少部分液體克服表面張力影響脫離液膜飛濺射出。

圖6 液膜沖擊附板的行為過程分析Fig.6 Analysis of film impinging on attached plate

圖7 液膜厚度和速度場分析Fig.7 Analysis of falling film thickness and velocity field

在液膜第1次聚合階段，當附板寬度足夠短、焊縫高度足夠高時，液膜在沖擊附板后可能無法在附板上完成第1次再聚合而直接以射流形式流入返回槽中，這從附板上液膜厚度小于附板前的液膜厚度分布可看出。

在降液膜沖擊附板的過程中，液膜沖擊附板后部分克服表面張力作用飛濺射出的液體是造成液膜損失，即冷凝回流損失的主要原因。

(9)

在液膜的第2次聚合階段，由于液膜在脫離附板的過程中速度存在差異，因此脫離壁面射出的液體回到主板進行聚合所需時間不同，導致附板后的平板上液膜厚度分布呈現不均勻的現象。而從附板后平板上液膜厚度大于附板上液膜厚度可發現，參與第2次聚合的液滴中有部分液滴在沖擊附板后直接越過附板飛回到平板上。

2.3 液膜流量對液膜損失的影響

圖8 液膜流量對液膜損失的影響Fig.8 Influence of film flow rate on film loss

圖8為液膜流量對液膜損失的影響?？煽闯?，在模擬范圍內，小流量下液膜損失率較小，隨著流量增大，2、3號附板的損失率先增大后略有減小，而1號附板的損失率在大流量下出現“突增”現象。附板尺寸列于表3。

表3 附板尺寸Table 3 Size of different attached plates

通過增加不同流量的數值模擬點，發現隨著流量的逐漸增大，1號附板出現損失率變化的突增點。分析其原因為：由于1號附板的寬度最小、厚度最大，從而液膜在到達附板時具備的速度最大，且沿附板流動時獲得的垂直于附板方向的速度最大。

根據O’Rourke等[14]和Wu等[15]的研究，液膜的分離過程和液滴從液膜中剝離主要受We影響，而上述過程是影響液膜在沖擊附板后損失的主要原因。由此，定義附板與返回槽間的高度差和返回槽高度的比值為液膜濺射空間數He：

He=(L-H)/L

(10)

式中：L為返回槽高度(100 mm)；H為附板高度。

從而通過比較液膜損失率ηloss和We以及He之間的關系，得到關系式(11)：

ηloss∝ΔWex·Hey0

(11)

式中：ΔWe=|We-Wecr|，Wecr為液膜沖擊附板時導致液膜損失率出現“突增”現象時的We，在本文模擬范圍內取1.5；x和y均為常數。由于小流量下液膜損失率較小，且在模擬過程中未表現出明顯規律，預測價值較小，故本文主要針對較大流量(Q=500 kg/h和1 000 kg/h)分析液膜損失率和以上二者的關系。

在本文模擬范圍內，分別取y=-1.2，x=0.82。ηloss與Wex·Hey的關系如圖9所示，誤差線為10%的相對誤差，可看出關系式與數據點之間擬合關系良好。

2.4 附板焊縫高度對液膜損失的影響

圖10為附板焊腳高比對液膜損失的影響?？煽闯?，在實驗和模擬范圍內，隨著焊腳高比的不斷減小，模擬計算獲得的液膜損失率與實驗值間的誤差越來越大，本文所用計算模型的預測準確性越來越低。

圖9 不同流量下損失率和Wex·Hey之間的關系Fig.9 Relationship between film loss rate and Wex·Hey under different flow rates

圖10 附板焊腳高比對液膜損失的影響Fig.10 Influence of attached plate’s welding height on film loss

圖11a、b為液膜沖擊焊腳高比β為1.0和0.3的附板時垂直于板面方向的速度分布場，圖11c、d為液膜的速度矢量場。從圖11a可看到，EWF模型模擬β=1.0的附板時對速度場的模擬與實際情況較接近，能較好地模擬液膜沖擊時飛濺射出的現象。圖11b顯示，液膜沖擊附板時液體在拐角處獲得較高焊腳附板工況更大的z方向速度，而由于EWF模型基于液膜內速度與壁面平行的基本假設，因此在經過拐角處時，可明顯看到液膜z方向速度以極高的速度梯度減小。

從圖11c、d中圓圈標記處能更加清晰地看到EWF模型對不同于壁面方向速度的修正。對于焊腳高比較大的工況，由于液膜垂直于運動方向的速度不大，EWF模型的這種修正對液膜的流動影響較??；而對于焊腳高比較小的工況，EWF模型對液膜垂直于運動方向的速度進行修正，使液膜運動方向強行貼合壁面，導致液膜的流動不符合慣性和實際作用，最后導致模型的計算結果與實驗結果偏差較大。

圖11 液膜沖擊不同焊腳高度的附板時速度流場Fig.11 Velocity field of falling film impinging on attached plate with different welding heights

3 總結與討論

本文對大平板上降液膜沖擊附板的流動行為進行分析，建立了相關幾何模型，利用商用CFD程序FLUENT建立了降液膜沖擊附板行為的相關數值模擬方案并進行了數值模擬研究。經過對計算結果的分析和與相關實驗結果的比較，發現大平板上降液膜沖擊附板的流動行為過程主要可分為4個階段：液膜沖擊附板階段、液膜在附板上第1次再聚合階段、液膜脫離附板分離階段、液膜在平板上第2次再聚合階段。通過對4個階段研究，主要得到以下結論。

1) 液膜在小流量下沖擊附板時普遍損失率較小，而在中大流量下沖擊附板有較大損失，且損失率受附板尺寸影響；損失率和韋伯數、受附板高度影響的液膜空間數之間的關系為ηloss∝Wex·Hey(0

2) FLUENT內EWF模型對于液膜沖擊高焊腳的附板能實現較好的預測；而對于低焊腳的附板，EWF模型對液膜沖擊附板后的速度有較大修正，導致EWF模型對液膜沖擊低焊腳的附板后的損失無法很好預測。