秦如冰，柴 翔，程 旭

(上海交通大學(xué)核科學(xué)與工程學(xué)院，上海 200240)

第四代快堆系統(tǒng)設(shè)計(jì)過(guò)程中關(guān)注的是固有安全性以及被動(dòng)安全性。反應(yīng)堆的固有安全性表示設(shè)計(jì)應(yīng)使其僅依據(jù)自然狀態(tài)即可保持安全狀態(tài)；并確保所有可能的情況下所有性能特征都處于安全范圍內(nèi)。被動(dòng)安全性有更加廣泛的定義，表示無(wú)需任何人為干預(yù)，不用觸動(dòng)信號(hào)，也不需要實(shí)現(xiàn)外部供能，反應(yīng)堆就可以保持安全狀態(tài)。非能動(dòng)停堆組件體現(xiàn)的正是這樣的特性。鈉冷快堆(Sodium-cooled Fast Reactor,SFR)作為第四代反應(yīng)堆系統(tǒng)的重要堆型之一，采用非能動(dòng)停堆組件以保證其安全性已成為國(guó)際上快堆發(fā)展和研究的最主要研究方向之一[1]。

然而，在對(duì)非能動(dòng)停堆組件落棒停堆過(guò)程進(jìn)行模擬分析時(shí)。由于復(fù)雜幾何、存在孔隙結(jié)構(gòu)以及運(yùn)動(dòng)邊界等問(wèn)題的存在，傳統(tǒng)計(jì)算流體力學(xué)程序(Computational fluid dynamics,CFD)所使用的結(jié)構(gòu)化網(wǎng)格或非結(jié)構(gòu)化網(wǎng)格生成方法存在較大的局限性。浸沒(méi)邊界法(Immersed boundary method,IBM)是CFD中的一種方法，該法無(wú)需構(gòu)建貼體網(wǎng)格，而是采用笛卡爾網(wǎng)格，并通過(guò)添加體積力將邊界條件納入控制方程中，簡(jiǎn)單的網(wǎng)格和數(shù)據(jù)結(jié)構(gòu)非常適合處理復(fù)雜幾何和移動(dòng)邊界問(wèn)題[2]。IBM最早于20世紀(jì)70年代由Peskin[3]提出，它最初用于模擬心臟中的血液流動(dòng)模式；經(jīng)過(guò)多年發(fā)展已廣泛應(yīng)用于工程模擬中。

本文針對(duì)鈉冷快堆非能動(dòng)停堆組件落棒過(guò)程模擬困難的問(wèn)題，基于浸沒(méi)邊界法的特性，開(kāi)發(fā)適用于復(fù)雜幾何和移動(dòng)邊界問(wèn)題的數(shù)值模擬程序，并對(duì)程序進(jìn)行驗(yàn)證分析。本文驗(yàn)證了層流狀態(tài)不同雷諾數(shù)的二維(2D)工況，包含固定和振蕩圓柱繞流的模擬結(jié)果。本文工作為后續(xù)三維工況、湍流以及實(shí)際幾何模型的模擬分析驗(yàn)證提供了堅(jiān)實(shí)的基礎(chǔ)。

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

1.1 網(wǎng)格類型

程序中涉及兩種網(wǎng)格。第一種稱為背景網(wǎng)格，為包含整個(gè)計(jì)算域并忽略浸沒(méi)體的笛卡爾網(wǎng)格。第二種網(wǎng)格用來(lái)表示浸沒(méi)體的幾何表面[4]。該表面網(wǎng)格無(wú)孔，因此可將浸沒(méi)體內(nèi)部與流體完全分離。

將表面網(wǎng)格嵌入背景網(wǎng)格中會(huì)把背景網(wǎng)格劃分為三種不同類型。歐拉網(wǎng)格(或流體網(wǎng)格)是指完全在流體區(qū)域中的網(wǎng)格，拉格朗日網(wǎng)格(或固體網(wǎng)格)則完全在浸沒(méi)體內(nèi)，浸沒(méi)邊界網(wǎng)格(IB網(wǎng)格)與浸沒(méi)體表面網(wǎng)格相交。表面網(wǎng)格上的節(jié)點(diǎn)稱為浸沒(méi)邊界節(jié)點(diǎn)(IB節(jié)點(diǎn))如圖1所示。在背景網(wǎng)格的歐拉場(chǎng)的流體網(wǎng)格中求解流動(dòng)的控制方程，并且通過(guò)浸沒(méi)體的拉格朗日網(wǎng)格來(lái)跟蹤表面的位置和形狀。通過(guò)在背景網(wǎng)格中重新定位浸沒(méi)體可以解決邊界的運(yùn)動(dòng)問(wèn)題，并通過(guò)重新生成表面網(wǎng)格可以解決浸沒(méi)體表面可能產(chǎn)生的變形。

圖1 浸沒(méi)邊界法網(wǎng)格劃分示意圖Fig.1 Cell types in IBM

1.2 程序求解過(guò)程

預(yù)估步：

不含重力的不可壓縮流體Navier-Stokes公式為 ：

(1)

連續(xù)性方程如下所示：

·U=0

(2)

使用公式(1)求解不包含浸沒(méi)邊界(IB boundary)區(qū)域的預(yù)估速度u*。在求解矢量時(shí)，公式(1)變?yōu)椋?/p>

Au*=r′u-p

(3)

式中：A——系數(shù)數(shù)組乘以解,u*；

r＇——系數(shù)數(shù)組乘以上一步的速度u。

u*是唯一的未知數(shù)?？蓪的系數(shù)數(shù)組分為對(duì)角矩陣和非對(duì)角矩陣，向量形式化為：

(4)

式中：αp——A的對(duì)角矩陣；

α＇N——A的非對(duì)角矩陣；

r——r＇·u*的矩陣。

公式(4)可以化為：

(5)

式中：αN=r-α＇Nu*。

校正步：

包含校正步速度貢獻(xiàn)(不含壓力)的變量U可通過(guò)下式計(jì)算得到：

(6)

IB網(wǎng)格上的速度通過(guò)公式(7)進(jìn)行插值處理可得：

(7)

式中：x，y——網(wǎng)格坐標(biāo)；

C——插值系數(shù)；

下標(biāo)p——IB網(wǎng)格；

下標(biāo)ib——IB節(jié)點(diǎn)的參數(shù)。

進(jìn)而公式(5)可化為：

(8)

考慮到連續(xù)性方程，公式(8)可變?yōu)椴此煞匠桃郧蠼庑拚龎毫Γ?/p>

(9)

IB網(wǎng)格和流體網(wǎng)格界面的通量可通過(guò)公式(10)計(jì)算得到：

(10)

式中：fib——IB網(wǎng)格和流體網(wǎng)格的界面；

f——靠近IB的流體網(wǎng)格界面；

S——界面面積；

n——界面法向量；

v——界面上的法向速度。

由公式(11)給出計(jì)算：

(11)

通過(guò)以下公式計(jì)算不包括壓力作用的通量：

(12)

根據(jù)公式(9)得到的修正壓力，通量可修正為：

(13)

同樣的，用修正壓力來(lái)修正速度，通過(guò)將公式(6)與壓力梯度相加可得更新后的速度為：

(14)

在程序?qū)嶋H求解過(guò)程中，插值方程系數(shù)C通過(guò)建立未知系數(shù)求解矩陣得到：

C=(MTWM)(-1)MTWA

(15)

A為由ui-uP作為元素構(gòu)成的矢量，矩陣M表示處理狄利克雷邊界條件時(shí)擴(kuò)展模板網(wǎng)格的幾何位置信息，對(duì)于二維工況令Xi=xi-xib，Yi=Yi-Yib，則有：

(16)

對(duì)于固定浸沒(méi)體，包含幾何位置信息的M矩陣只需計(jì)算一次，對(duì)于運(yùn)動(dòng)物體即移動(dòng)邊界問(wèn)題則需要在每個(gè)時(shí)間步上重新進(jìn)行計(jì)算。W則為由權(quán)重wi作為對(duì)角元素的矩陣，權(quán)重與插值網(wǎng)格與浸沒(méi)邊界網(wǎng)格的距離成反比。

對(duì)于諾依曼邊界條件，反應(yīng)插值模板幾何位置信息的矩陣M則應(yīng)寫為：

(17)

與狄利克雷邊界條件類似，諾依曼邊界條件下矢量A包含了插值模板的壓力值以及壓力梯度，即：

(18)

由于邊界條件是通過(guò)插值計(jì)算來(lái)施加的，因此浸沒(méi)邊界網(wǎng)格中速度值的增大或者減小將違反計(jì)算域流體的質(zhì)量守恒條件。因?yàn)榻](méi)邊界網(wǎng)格界面上的速度是通過(guò)浸沒(méi)邊界網(wǎng)格和相鄰網(wǎng)格線性插值得到的，當(dāng)計(jì)算通過(guò)浸沒(méi)邊界網(wǎng)格界面的流量時(shí)，這一影響將變得更加明顯。為了糾正此不足，需要在浸沒(méi)邊界網(wǎng)格界面上應(yīng)用適當(dāng)?shù)目s放因子進(jìn)行流量調(diào)整。流量不平衡的情況通量通過(guò)流動(dòng)比率進(jìn)行衡量，由下式表示流量不平衡情況：

Fimb=Fin-Fout-Fib

(19)

式中：Fin和Fout分別代表流體流入以及流出與浸沒(méi)邊界網(wǎng)格相鄰的流體網(wǎng)格。Fib是由于浸沒(méi)體邊界移動(dòng)(Sib)所產(chǎn)生的流體流動(dòng)速率，并由下式定義：

(20)

由該項(xiàng)衡量彈性所產(chǎn)生的影響，對(duì)于剛體而言，根據(jù)高斯散度定理，該項(xiàng)為零。如果不平衡項(xiàng)Fimb不為零，則通過(guò)調(diào)整并縮小Fin和Fout中較大的值來(lái)進(jìn)行修正。調(diào)整方式為乘以小于1的調(diào)整因子，調(diào)整因子的大小根據(jù)Fin和Fout的比值確定。

2 程序驗(yàn)證分析

為驗(yàn)證本文開(kāi)發(fā)程序正確與否，構(gòu)建以下幾何模型用于驗(yàn)證分析。流體流經(jīng)直徑為D的二維(2D)圓柱，如圖2所示，計(jì)算域是高度為H，長(zhǎng)度為L(zhǎng)i+L0的長(zhǎng)方形。入口處的邊界條件為沿著流動(dòng)方向x施加穩(wěn)定的均勻速度以及零壓力梯度。在出口處則采取質(zhì)量守恒條件和零壓力梯度。計(jì)算域的頂部和底部采用速度的自由流邊界條件。

圖2 計(jì)算域及邊界條件Fig.2 Computational domain and boundary conditions

將圓柱的中心點(diǎn)定義為原點(diǎn)，計(jì)算域的尺寸遵循如下準(zhǔn)則，在x，y方向上滿足[-16D,48D]×[-16D,16D]的幾何尺寸[5]。并將計(jì)算域的網(wǎng)格分不同區(qū)域進(jìn)行細(xì)化，細(xì)化區(qū)域及細(xì)化網(wǎng)格尺寸如圖3所示，網(wǎng)格在圓柱的周圍是均勻的，對(duì)應(yīng)區(qū)域范圍為-D≤x≤D，-D≤y≤D。根據(jù)以上原則，在程序中生成用于模擬計(jì)算的網(wǎng)格如圖4所示，為節(jié)省計(jì)算資源，在模擬固定圓柱工況時(shí)圓柱內(nèi)部的背景網(wǎng)格因不參與計(jì)算不需要進(jìn)行加密。通過(guò)圖4圓柱所在位置我們可以看出程序計(jì)算時(shí)無(wú)需生成貼體網(wǎng)格。

圖3 計(jì)算域劃分及網(wǎng)格尺寸Fig.3 Computational domain decomposition and grid spacings

圖4 程序計(jì)算所生成的網(wǎng)格Fig.4 Mesh in the code

程序驗(yàn)證是基于流體流經(jīng)直徑為D的固定或振蕩圓柱的二維(2D)模擬，并在各種雷諾數(shù)(Re=U∞D(zhuǎn)/μ)下進(jìn)行的，并將計(jì)算結(jié)果與文獻(xiàn)中的數(shù)據(jù)進(jìn)行對(duì)比。斯特勞哈爾數(shù)(St)，阻力(CD)和升力(CL)系數(shù)由式(21)定義,其中fv表示脫落頻率，F(xiàn)d和Fl分別表示單位長(zhǎng)度的阻力和升力。各種算例驗(yàn)證過(guò)程中所涉及的幾何模型、雷諾數(shù)以及所需比較的參數(shù)如表1所示。

(21)

式中：L——再循環(huán)長(zhǎng)度;

a——圓柱與再循環(huán)中心之間的距離;

b——兩個(gè)再循環(huán)中心的垂直距離;

θ——后駐點(diǎn)的分離角。

穩(wěn)態(tài)流動(dòng)下的渦旋特征參數(shù)如圖5所示。

表1 程序驗(yàn)證分析工況匯總

圖5 渦旋特征參數(shù)Fig.5 Characteristic geometrical parameters in the steady regime

Re=30時(shí)的流動(dòng)為穩(wěn)態(tài)，特征在于位于圓柱后穩(wěn)定的再循環(huán)區(qū)域。圓柱體周圍的網(wǎng)格尺寸為Δx=Δy=0.01D，經(jīng)過(guò)模擬后得到的數(shù)據(jù)，圖5中定義的所有特征幾何參數(shù)與表2中的文獻(xiàn)數(shù)據(jù)相當(dāng)，經(jīng)細(xì)化網(wǎng)格后模擬的結(jié)果誤差小于10%。

表2 Re=30時(shí)模擬結(jié)果與文獻(xiàn)對(duì)比

隨著雷諾數(shù)的增大，根據(jù)Williamson[9]和Norberg[10]的研究，Rec=40為流動(dòng)轉(zhuǎn)為不穩(wěn)定流動(dòng)的臨界值，因此，模擬了在Re=100和185時(shí)的工況，即具有渦旋脫落的2D非穩(wěn)態(tài)區(qū)域模擬，圓柱體周圍的網(wǎng)格尺寸為Δx=Δy=0.01D，圖6展示的渦旋輪廓為著名的卡門渦街，其特征在于渦旋的周期性脫落，對(duì)流并從圓柱體上向后散開(kāi)，將Re=185時(shí)的渦脫落結(jié)構(gòu)與Guilmineau[12]和Pinelli[5]的工作進(jìn)行對(duì)比，結(jié)果符合較好。

圖6 Re=185時(shí)的渦脫落Fig.6 Vorticity countours at Re=185

CD和CL隨時(shí)間t的變化，其中時(shí)間做無(wú)量綱處理作為自變量，如圖7所示，表明升力和阻力波動(dòng)的幅度隨著雷諾數(shù)的增加而增加，符合Guilmineau和Queutey[12]的研究結(jié)論。

對(duì)于不穩(wěn)定流動(dòng)兩個(gè)不同雷諾數(shù)下得到的模擬數(shù)據(jù)：斯特勞哈爾數(shù)，平均阻力(在10個(gè)時(shí)間段內(nèi)計(jì)算)和均方根升力系數(shù)與表3和表4中總結(jié)的文獻(xiàn)數(shù)據(jù)相當(dāng)。

表3 Re=100時(shí)模擬結(jié)果與文獻(xiàn)對(duì)比Table 3 2D flow past a fixed cylinder at Re=100. Numerical and experimental data from the literature are provided for comparison

表4 Re=185時(shí)模擬結(jié)果與文獻(xiàn)對(duì)比 Table 4 2D flow past a fixed cylinder at Re=185. Numerical and experimental data from the literature are provided for comparison

圖7 CD和CL隨時(shí)間變化圖Fig.7 Temporal evolutions of CD and CL for the 2D flow at Re=100 (a);Re=185 (b)

針對(duì)Re=185時(shí)做網(wǎng)格敏感性分析，得到阻力系數(shù)在不同網(wǎng)格尺寸下的變化情況。在上述模擬分析驗(yàn)證過(guò)程中，共進(jìn)行了5次加密，笛卡爾網(wǎng)格的最小尺寸為Δx=Δy=0.01D。為了進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證，本文分別計(jì)算了加密次數(shù)為1至4次時(shí)的阻力系數(shù)與升力系數(shù)變化情況，對(duì)應(yīng)的笛卡爾網(wǎng)格最小尺寸分別為：Δx=Δy=0.16D、0.08D、0.04D、0.02D。不同網(wǎng)格尺寸下得到的阻力系數(shù)隨時(shí)間變化如圖8所示，可以看出，當(dāng)最小網(wǎng)格尺寸為0.02D時(shí)與0.01D下的計(jì)算結(jié)果相近。

在Δx=Δy=0.01D以及0.02D下分別對(duì)阻力系數(shù)曲線做快速傅里葉變換，阻力系數(shù)的頻譜分析圖9如所示，當(dāng)網(wǎng)格最小尺寸為0.01D時(shí)，頻率等于0時(shí)振幅最大為1.28，當(dāng)網(wǎng)格最小尺寸為0.02D時(shí)，頻率等于0時(shí)振幅最大為1.29。由此可知當(dāng)網(wǎng)格尺寸小于0.02D時(shí)即可滿足網(wǎng)格無(wú)關(guān)性要求。

圖8 不同網(wǎng)格尺寸下，阻力系數(shù)CD隨時(shí)間的變化Fig.8 Drag coefficient CD as a function of the time for the 2D flow past an oscillating cylinder at different cell size

圖9 阻力系數(shù)頻譜分析圖Fig.9 Frequency-amplitude plot of drag coefficients

為了驗(yàn)證程序處理移動(dòng)物體的能力，在Re=500下，對(duì)流體流經(jīng)正弦運(yùn)動(dòng)的圓柱體工況進(jìn)行了二維模擬。圓柱體在垂直方向上以固定的振幅比A(=ymax/D)=0.25以及頻率比F=f0/fv=0.975(f0為圓柱振蕩頻率)。計(jì)算域與計(jì)算固定圓柱體時(shí)模擬的域相同，圓柱體周圍的網(wǎng)格尺寸為Δx=Δy=0.02D。脫落頻率fv為Re=500時(shí)，基于二維固定圓柱繞流的模擬結(jié)果。一旦流動(dòng)穩(wěn)定為2D渦脫落狀態(tài)，我們便開(kāi)始移動(dòng)圓柱體。

在圖10中詳細(xì)描述了圓柱做正弦振蕩的流程。 左邊的五張圖顯示了在整個(gè)渦旋脫落周期一半的過(guò)程中五個(gè)瞬間繪制的渦旋輪廓，在此期間，升力沿向上的方向起作用。與右邊所示的Blackburn[6]的模擬結(jié)果進(jìn)行比較，對(duì)比結(jié)果顯示計(jì)算較為準(zhǔn)確，進(jìn)而可以很好地預(yù)測(cè)渦脫落過(guò)程的空間動(dòng)態(tài)以及附著點(diǎn)和分離點(diǎn)隨著時(shí)間的演變。圖11顯示了在最后10個(gè)振蕩周期內(nèi)取升力系數(shù)的平均值隨時(shí)間的變化，其中時(shí)間做無(wú)量綱化處理2πt/T，T為圓柱振蕩周期。同樣，結(jié)果可以認(rèn)為與Blackburn[6]的參考數(shù)據(jù)非常匹配。

圖10 在Re=500時(shí)，二維振蕩圓柱繞流的瞬時(shí)渦輪廓Fig.10 Instantaneous contours of vorticity for the 2D flow past an oscillating cylinder at Re=500

圖11 在Re=500時(shí)，升力系數(shù)CL隨時(shí)間的變化Fig.11 Lift coefficient CLas a function of the time for the 2D flow past an oscillating cylinder at Re=500

3 結(jié)論

本研究開(kāi)發(fā)了基于浸沒(méi)邊界法的CFD程序求解器。對(duì)移動(dòng)和固定的二維圓柱繞流工況在不同雷諾數(shù)下進(jìn)行模擬驗(yàn)證分析，雷諾數(shù)變化范圍為30到500，包括穩(wěn)態(tài)流動(dòng)和不穩(wěn)定流動(dòng)共4種不同工況。將計(jì)算所得的特征幾何參數(shù)、阻力系數(shù)、升力系數(shù)以及渦脫落情況等參數(shù)與文獻(xiàn)進(jìn)行對(duì)比，模擬結(jié)果顯示，使用本程序計(jì)算具有良好的效率和準(zhǔn)確性。并針對(duì)Re=185的工況進(jìn)行網(wǎng)格敏感性分析，通過(guò)阻力系數(shù)完成了網(wǎng)格無(wú)關(guān)性驗(yàn)證。但目前的模擬范圍還只限于2D工況以及雷諾數(shù)較小的層流，本文為未來(lái)工作聚焦于驗(yàn)證3D以及湍流下模擬結(jié)果的準(zhǔn)確性，并針對(duì)鈉冷快堆非能動(dòng)停堆組件停堆落棒實(shí)際情況中面臨的核反應(yīng)堆熱工水力復(fù)雜幾何問(wèn)題進(jìn)行模擬打下了基礎(chǔ)。