壓力管道裂紋泄漏計算流體力學計算方法

郭超, 魏超, 朱光昱, 唐暉, 劉巧鳳

(生態環境部核與輻射安全中心, 北京 100082)

破前漏(leak-before-break,LBB)技術已經逐漸應用到核電廠壓力管道的設計中,其重要前提是在管道壁形成的穿透裂紋在擴展到發生破壞性斷裂的臨界長度前,泄漏速度足夠大,確保能被檢測到。壓力管道裂紋泄漏率的計算對LBB技術的應用和相關的安全分析具有重要意義。

目前壓力管道裂紋泄漏率計算有兩類方法,一類是基于兩相臨界流相關經驗公式編制泄漏率計算程序進行計算,其從兩相流模型壓降分析出發,判斷是否出現臨界流并計算相應質量流量。國際上編制了相應的貫穿裂紋泄漏率計算程序,例如,美國的SQUIRT、PICEP,法國的ASTEQ、LRCPC,德國的FLORA等;近些年國內的研究院和高校也開展了該類泄漏率計算程序的開發和研究工作[1-6]。這些計算程序主要以早些年代的Moody均勻化平衡態臨界流模型和Henry-Fauske均勻化非平衡態臨界流模型等或在其改進基礎上而開發,利用上述程序可計算獲得泄漏率在數值上的結果。另一類方法是對泄漏裂紋進行建模,利用CFD程序直接計算裂紋泄漏過程。與常見的吸熱蒸發相變不同,壓力管道裂紋泄漏相變屬于壓降導致的蒸發相變,這給CFD計算帶來比較大的挑戰。一些學者利用CFD方法開展了泄漏流動研究,王國濤等[7]利用CFX程序對壓力管道在不同工作壓力下的微小泄漏流場進行了數值模擬;王金生等[8]利用Fluent程序對單相流體在壓力管道中裂紋不同張開方向上的泄漏速率進行了計算研究;岳云飛等[9]采用Fluent程序對非金屬管道微小泄漏的管道內流體與泄漏口的流動狀態進行了分析;高曉輝等[10]對液體泄漏破碎行為進行了數值模擬研究;紀健等[11]對氣液兩相流管道泄漏流場特性進行了數值模擬分析;王瑤等[12]開展了壓力管道孔隙泄漏率計算并采用 CFD 方法獲取了“泄漏率-孔徑”關系;另外,在原油管道泄漏行為研究中CFD方法也被廣泛應用[13-16]。在現階段研究中,基于CFD方法對泄漏過程的研究多側重于泄漏前后流體相態不發生改變的流動,對于泄漏過程發生壓降相變的流動過程研究較為少見。

現建立一種壓力管道裂紋泄漏率CFD計算方法,實現流體壓降相變泄漏過程的瞬態模擬,分析裂紋泄漏過程的流動特性,并對裂紋泄漏率影響因素進行初步分析。

1 壓力管道裂紋泄漏計算方法

核電廠壓力管道出現貫穿裂紋時,管壁內外壓差將驅動管內液體通過裂紋向外泄漏,過冷液體或飽和液體在流動過程中處于快速減壓狀態,在理論上,當靜壓降低到其飽和壓力后,液體開始產生相變(即閃蒸),泄漏液體將從單相狀態過渡到兩相混合的流動狀態。若泄漏壓差足夠大,則會出現兩相臨界流現象。圖1為流體在貫穿裂紋中泄漏流動示意圖[17]。

圖1 通過裂紋的兩相流Fig.1 Two phase flow in crack

本文計算方法主要思想為利用空泡動力學模型來描述泄漏相變的傳質過程。研究對象主要針對大長徑比裂紋,其可忽略泄漏過程中的熱力學非平衡效應影響[18]。計算中考慮氣體可壓縮性,不考慮液體可壓縮性。對于裂紋流道近壁區流動,通常采用近壁模型或壁面函數方法來處理,近壁模型處理方法需修改湍流模型,在邊界層內建立細密的網格,使黏性影響區域(黏性底層和過渡層)能夠通過網格一直到壁面進行解析;而壁面函數處理方法實際上使用半經驗公式“連接”壁面和完全湍流區域之間的黏性影響區,可以使用相對粗糙的網格對壁面附近的高梯度剪切層進行建模,避免了在湍流模型中考慮黏性效應。考慮到裂紋泄漏率計算問題不關注邊界層內部流動細節,也不涉及顯著的流動分離、渦流旋轉等情況,因此不需對邊界層建立十分詳細的網格模型,故使用壁面函數方法處理近壁區的流動。使用壁面函數方法也可以顯著地減少壁面法向的網格數量并節省計算資源。

1.1 相變模型

相變模型采用Rayleigh Plesset Model方程[19-20],該方程可體現流體由壓力變化而導致的相變現象,其表達式為

(1)

式(1)中:RB為汽泡半徑;t為時間;Pv為氣泡的蒸汽壓;P為氣泡周圍液體壓力;ρf為液體密度;σ為液體和蒸汽間的表面張力系數。

忽略二次項(其適用于低頻率震蕩)和表面張力,方程變為Rayleigh Plesset Model簡化模型,即

(2)

氣泡體積變化率為

(3)

氣泡質量變化率為

(4)

式(4)中:mB為氣泡質量;ρg為氣體密度。

若NB為單位體積氣泡數,則氣態相體積分數φg表示為

(5)

單位體積的總質量變化為

(6)

式(6)可用于蒸發相變,將式(6)進行如下改變將同時也適用于冷凝,即

(7)

式(7)中:F為經驗系數,因冷凝速度通常比蒸發速度慢,所以對于冷凝和蒸發F取值不同。通常Fvap=50和Fcond=0.01。

對于蒸發情況式(7)可進一步修正。蒸發一般開始于氣核(nucleation)位置(通常為不可冷凝氣體),隨蒸汽體積分數的增加,氣核點密度必然減小。在式(7)中,讓φnuc(1-φg)替代φg,Rnuc替代RB,則公式變為

(8)

式(8)中:φnuc為氣核的體積分數;Rnuc為氣核的半徑。對于冷凝情況式(8)仍然適用。

1.2 湍流模型

湍流模型可采用具有工程普遍適用性的k-ε模型。該湍流模型可結合Scalable壁面函數來提高近壁計算的魯棒性和準確性。湍流動能和湍流耗散率輸運方程表達式為

ρε+Pkb

(9)

(10)

式中:xj為坐標位置,j=1、2、3,分別表示x、y和z三個方向;Uj為xj方向的速度;k為湍流動能;ε為湍流動能耗散率;μ為分子(動力)黏度;μt為湍流黏度;σk和σε分別為k-ε湍流模型常數;Cε1、Cε2為雷諾應力模型常數;Pkb和Pεb為浮力影響項;Pk為黏性力湍流影響項。

1.3 壁面函數模型

在壁面函數法中,通過經驗公式為平均流和湍流輸運方程提供近壁邊界條件。這些經驗公式把壁面條件(如壁面剪切應力等)與近壁面網格節點的變量聯系起來,假定該節點位于邊界層的完全湍流區域(對數律區)。在該區域,近壁切向速度與壁面剪應力呈對數關系,公式為

(11)

(12)

(13)

式中:u+為無量綱近壁速度;uτ為摩擦速度;Ut表示距離壁面距離為Δy處的切向速度;y+為距離壁面的無量綱距離;τω為壁面切向應力(剪切應力);κ為馮卡門常數;C為依賴于壁面粗糙度的對數層常數。

(14)

式(14)中:y*為對應u*的無量綱距離。

壁面剪切力和無量綱距離表達式變為

τω=ρu*uτ

(15)

(16)

壁面函數方法主要缺點是計算結果取決于最靠近壁的點的位置并且對近壁網格敏感,細化網格并不一定能給出提高精度的解[22]。采用Scalable壁面函數解決該問題,其可以應用在任意精細的網格上,并允許不考慮雷諾數而對模型進行一致性的網格細化。Scalable函數通過y*=max(y*,11.06)來限制對數公式中y*的取值,其中11.06是對數區和線性近壁區的交點處y*的取值。計算出的y*不能低于上述限制,以保證所有網格點都在黏性底層之外(對數區),細化網格導致的結果不一致問題可得到避免。

2 模型建立與計算

2.1 幾何模型與計算模型

本文計算方法研究采用矩形窄縫裂紋簡化真實的裂紋形狀[23]。如圖2所示,L為裂紋深度(即管壁厚度),2a為裂紋長度,δ為裂紋張開位移(crack opening displacement,COD),位置1為入口,位置2為出口。入口與管道內壁的流體相連,出口為管道外部環境。

圖2 矩形流道示意圖Fig.2 Diagram of two phase flow in crack

由于通常裂紋長度2a遠大于COD,因此可建立二維幾何模型進行流道模型簡化。二維模型即沿流動方向的流道矩形截面,其法向為裂紋長度方向。考慮若僅針對流道建立二維模型,則流道出口即為計算模型的邊界,為了保證出口提取的計算結果的準確性,最終建立的二維模型考慮在流道出口位置拓展建立了模擬外部環境的空氣域,避免了可能的計算模型邊界效應對所提取計算結果的影響。對于不支持二維網格算法的程序,可在厚度方向拉伸單層網格并在厚度方向設置對稱邊界條件來進行處理。圖3為裂紋流道簡化模型示例,坐標原點為流道入口中心點,X軸正方向為泄漏流動方向,Y軸為裂紋張開位移方向。

圖3 二維計算模型Fig.3 Two dimensional calculation model

2.2 計算設定

裂紋泄漏CFD計算類型為瞬態計算,根據大量試算經驗,建議計算時間步長設置為1×10-8～1×10-6s,在計算中需要結合具體網格尺寸進行適配。為了更符合泄漏實際過程,計算方法采用了非均相流(相間傳輸)模型,即兩相間可考慮動量、熱量和質量的界面傳遞,其中,質量傳遞模型采用Rayleigh Plesset模型,同時需根據具體工況確定液體飽和蒸氣壓、氣液最大密度比等參數;動量傳遞模型可采用一般的拖曳力模型;對于熱量傳遞,如前所述大長徑比裂紋可忽略熱力學非平衡效應影響。

對于裂紋流道域,泄漏入口可選擇考慮一種流體流入到另一種流體的射流卷吸效應的壓力入口邊界條件;流道模型Y軸方向的上/下邊界設置為無滑移邊界條件,其中粗糙度相關參數根據具體工況流道表面情況確定;流道模型XY平面的法向邊界設定為對稱邊界條件;湍流設置方面采用k-ε湍流模型和Scalable壁面函數的組合;對于出口空氣域,空氣域XY平面法向邊界設定為對稱邊界,其余邊界設定為出口邊界條件。

2.3 網格無關性驗證

分別對網格尺寸為0.05、0.04、0.03、0.025、0.01 mm的裂紋模型(L=36 mm,2a=10 mm,δ=0.5 mm)進行了泄漏率計算,計算過程中殘差處于穩定震蕩狀態且殘差值均在10-4以下,計算過程收斂性良好。待裂紋流道出口流量穩定后,獲得的不同網格尺寸下的出口泄漏流量的結果如圖4所示。當網格尺寸小于0.05 mm后計算結果基本相差不大,兼顧服務器性能和計算效率,網格尺寸最終確定為0.03 mm,局部網格劃分如圖5所示。

圖4 不同網格尺寸下的計算結果對比Fig.4 Comparison of calculation results under different grid sizes

圖5 模型局部網格圖(0.03 mm)Fig.5 Mesh of simulation mode(0.03 mm)

2.4 典型工況計算結果

針對參考文獻[24]中試驗壓力范圍3.85～11.07 MPa下,220 ℃和265 ℃熱態試驗的10個工況開展泄漏率計算。試驗工況的裂紋幾何模型參數如下:試驗段L為80 mm,2a為10 mm,COD為0.6 mm,流道表面粗糙度為6 μm。下文以220 ℃、入口壓力為9.84 MPa工況為代表,從流體組分體積分數、流道壓力、流度以及泄漏率幾方面對裂紋泄漏過程的流動特性進行分析。

2.4.1 裂紋流道內部流體體積分數

圖6為從泄漏開始到泄漏狀態穩定,不同時刻液體體積分數的變化情況。圖7為0.2 ms時刻局部流道內的液體體積分數。

圖6 典型工況泄漏過程液體體積分數云圖Fig.6 Cloud diagram of liquid volume fraction in leakage process

圖7 0.2 ms時刻局部流道內的液體體積分數云圖Fig.7 Liquid volume fraction in the local channel at 0.2 ms

從圖6中可以看到,整個裂紋泄漏過程時間尺度為毫秒量級。初始時刻(0 ms)流道內充滿空氣,液體體積分數為0;泄漏初期泄漏液體逐漸進入裂紋流道,泄漏液體頭部區域由于當地靜壓低于泄漏流體的飽和蒸氣壓,在流動中已伴隨發生著相變現象,0.2 ms時刻流道出口已有少量的氣液兩相流流出;隨著泄漏時間的增加,泄漏流體逐漸充斥在裂紋流道空間內,根據流道出口泄漏率監測值判斷,1 ms時刻泄漏流量已趨于穩定。

為了更直觀地觀察流體內部流動特性情況,通過后處理獲取計算終止時刻泄漏入口到出口中心軸線上的氣、液組分體積分數,如圖8所示。在X=72.5 mm后液體體積分數逐漸開始低于1,即出現了相變。此后,沿流動方向,液體體積分數逐漸震蕩降低,氣體體積分數逐漸震蕩增加,不同位置體積分數的震蕩分布也體現了兩相流動的特點。

圖8 裂紋入口-出口中心點連線上的氣、液體積分數Fig.8 Volume fraction of gas and liquid on the connection line between the center point of crack inlet and outlet

2.4.2 裂紋流道內部壓力及流速

圖9為2 ms時刻裂紋流道靜壓分布情況,流道內壓在相變發生前沿流動方向近似線性降低,在相變后迅速降低,在出口位置降低到外部環境壓力。

圖9 裂紋入口-出口中心點連線上的靜壓(2 ms)Fig.9 Static pressure variation on the connection line between the center point of crack inlet and outlet(2 ms)

圖10是沿流動方向和COD方向的液體速度分布。其中,圖10(a)為流道內部液體流速沿著流動方向的分布,液體流速在相變發生前速度相對穩定,最大流速約86 m/s,在相變發生后迅速增大,最大流速約110 m/s;從該圖中不同Y位置流速分布可看出,流道中心泄漏速度最快,且在Y方向等距離變化下,越接近流道壁面,流速降低越顯著。圖10(b)為相變發生位置后流道內部液體流速沿COD方向(Y軸)的分布,同樣可看出流道中心速度顯著高于靠近壁面速度;從該圖不同X位置流速分布可以看出,在X方向等距離變化下,越接近泄漏出口,流速增大越顯著。

圖10 流動方向和COD方向的液體速度分布Fig.10 Velocity distribution of liquid in the direction of flow direction and COD

2.4.3 裂紋出口泄漏率

在裂紋出口位置設置泄漏率計算監測點(采樣頻率50 kHz)。圖11為泄漏率隨時間的變化。泄漏率呈現出先快速增大然后逐漸降低直至穩定震蕩的狀態。泄漏初期,流道內部空間以空氣為主,流體可以以較高速度流動,隨著泄漏過程的發展,液體逐漸占據流道內部空間,這個過程流道內部壓力分布始終處于變化狀態,從圖6中液體發生相變位置的不斷變化也可以看出這一點,約1 ms后,流道內部壓力分布基本穩定,裂紋出口泄漏率也呈現出的穩定震蕩狀態,這也反映出出口位置泄漏流體時時都處于氣液兩相混合流出的狀態,這與實際物理情況相符。

圖11 裂紋流道出口泄漏率時程Fig.11 Time history of leakage rate at the outlet of crack channel

基于以上對泄漏流體組分體積分數、壓降、流速和泄漏率等參數的分析表明,模擬計算的泄漏流動特性符合相變泄漏過程的物理規律,本文方法可以實現壓力管道裂紋泄漏過程的模擬。

3 泄漏率影響參數初步分析

3.1 計算結果驗證

圖12為2.4節所述220 ℃和265 ℃下10個工況的程序計算結果和試驗結果的對比,橫坐標P0為上游入口壓力,G為單位面積泄漏率。表1為計算結果誤差。

表1 220 ℃和265 ℃工況計算結果誤差Table 1 Calculation result error at 220 ℃ and 265 ℃

圖12 220 ℃和265 ℃下計算結果與試驗結果對比Fig.12 Calculation and test results at 220 ℃ and 265 ℃

從圖12中可以看出,程序計算結果與試驗結果總體上符合性較好,計算值整體上略高于試驗值。從表1中可看出,以試驗值為參考,泄漏率最小誤差3.6%,最大誤差16.3%,平均誤差在15%以內;且隨著裂紋入口滯止壓力增加,計算和試驗結果相對誤差呈現減小的趨勢,高壓下計算結果的符合性優于低壓下計算結果的符合性。對比結果表明本計算方法具有可靠性。

3.2 泄漏率影響參數初步分析

裂紋幾何和流體參數均會影響泄漏率,COD和過冷度分別是代表裂紋幾何和流體特征的典型參數。圖13為矩形裂紋流道下COD對泄漏率的影響,計算工況入口壓力8.92 MPa,流體溫度為223 ℃,L為57.2 mm,2a為63.5 mm,COD分別為0.2、0.25、0.3、0.35、0.4、0.45、0.5、0.55、0.6、0.65、0.7、0.8、0.85 mm(共計13個工況),裂紋表面粗糙度為10.2 μm。從圖13可以看出,當裂紋流道截面尺寸發生變化時,單位面積泄漏率也會變化,即隨著COD的增大而增大。本文計算范圍內單位面積泄漏率隨COD近似線性增長。

圖13 單位面積泄漏率隨COD的變化Fig.13 Variation of leakage rate per unit area with COD

圖14為矩形裂紋流道下過冷度對泄漏率的影響,計算工況入口壓力以及流體溫度變化范圍參考了一般壓水堆的設計壓力及冷卻劑堆芯進出口溫度,入口壓力為15.5 MPa,流體溫度分別為330、320、310、300、290、280 ℃(共計6個工況),即過冷度分別為14.7、24.7、34.7、44.7、54.7、64.7 ℃,其余工況參數L為57.2 mm,2a為63 mm,COD為0.2 mm,裂紋表面粗糙度為10.2 μm。從圖14可知,由于過冷度的增加,流體因溫度降低使得相變產生的蒸汽量減少,引起泄漏率的增加。同時,隨著過冷度的增加,單位面積泄漏率呈現非線性增加(圖14中擬合曲線為EXP指數關系),但其影響程度趨于平緩。

圖14 單位面積泄漏率隨流體過冷度的變化Fig.14 Variation of leakage rate per unit area with subcooling

4 結論

(1)建立了一種壓力管道裂紋泄漏CFD計算方法,可實現帶有相變的裂紋泄漏過程模擬,計算結果與試驗結果符合性良好,平均誤差在15%以內。

(2)對影響泄漏率參數的初步分析表明,單位面積泄漏率隨COD的增加而增加;隨流體過冷度的增加呈現增幅趨緩的增加。

(3)所提出的裂紋泄漏CFD計算方法具有通用性,可用于復雜形狀裂紋的泄漏率計算。