基于浸沒邊界法的流固耦合模擬分析

秦如冰，柴 翔，程 旭

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

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

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

本文針對鈉冷快堆非能動停堆組件落棒過程模擬困難的問題，基于浸沒邊界法的特性，開發適用于復雜幾何和移動邊界問題的數值模擬程序，并對程序進行驗證分析。本文驗證了層流狀態不同雷諾數的二維(2D)工況，包含固定和振蕩圓柱繞流的模擬結果。本文工作為后續三維工況、湍流以及實際幾何模型的模擬分析驗證提供了堅實的基礎。

1 數學模型

1.1 網格類型

程序中涉及兩種網格。第一種稱為背景網格，為包含整個計算域并忽略浸沒體的笛卡爾網格。第二種網格用來表示浸沒體的幾何表面[4]。該表面網格無孔，因此可將浸沒體內部與流體完全分離。

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

圖1 浸沒邊界法網格劃分示意圖Fig.1 Cell types in IBM

1.2 程序求解過程

預估步：

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

(1)

連續性方程如下所示：

·U=0

(2)

使用公式(1)求解不包含浸沒邊界(IB boundary)區域的預估速度u*。在求解矢量時，公式(1)變為：

Au*=r′u-p

(3)

式中：A——系數數組乘以解,u*；

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

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

(4)

式中：αp——A的對角矩陣；

α＇N——A的非對角矩陣；

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

公式(4)可以化為：

(5)

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

校正步：

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

(6)

IB網格上的速度通過公式(7)進行插值處理可得：

(7)

式中：x，y——網格坐標；

C——插值系數；

下標p——IB網格；

下標ib——IB節點的參數。

進而公式(5)可化為：

(8)

考慮到連續性方程，公式(8)可變為泊松方程以求解修正壓力：

(9)

IB網格和流體網格界面的通量可通過公式(10)計算得到：

(10)

式中：fib——IB網格和流體網格的界面；

f——靠近IB的流體網格界面；

S——界面面積；

n——界面法向量；

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

由公式(11)給出計算：

(11)

通過以下公式計算不包括壓力作用的通量：

(12)

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

(13)

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

(14)

在程序實際求解過程中，插值方程系數C通過建立未知系數求解矩陣得到：

C=(MTWM)(-1)MTWA

(15)

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

(16)

對于固定浸沒體，包含幾何位置信息的M矩陣只需計算一次，對于運動物體即移動邊界問題則需要在每個時間步上重新進行計算。W則為由權重wi作為對角元素的矩陣，權重與插值網格與浸沒邊界網格的距離成反比。

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

(17)

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

(18)

由于邊界條件是通過插值計算來施加的，因此浸沒邊界網格中速度值的增大或者減小將違反計算域流體的質量守恒條件。因為浸沒邊界網格界面上的速度是通過浸沒邊界網格和相鄰網格線性插值得到的，當計算通過浸沒邊界網格界面的流量時，這一影響將變得更加明顯。為了糾正此不足，需要在浸沒邊界網格界面上應用適當的縮放因子進行流量調整。流量不平衡的情況通量通過流動比率進行衡量，由下式表示流量不平衡情況：

Fimb=Fin-Fout-Fib

(19)

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

(20)

由該項衡量彈性所產生的影響，對于剛體而言，根據高斯散度定理，該項為零。如果不平衡項Fimb不為零，則通過調整并縮小Fin和Fout中較大的值來進行修正。調整方式為乘以小于1的調整因子，調整因子的大小根據Fin和Fout的比值確定。

2 程序驗證分析

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

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

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

圖3 計算域劃分及網格尺寸Fig.3 Computational domain decomposition and grid spacings

圖4 程序計算所生成的網格Fig.4 Mesh in the code

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

(21)

式中：L——再循環長度;

a——圓柱與再循環中心之間的距離;

b——兩個再循環中心的垂直距離;

θ——后駐點的分離角。

穩態流動下的渦旋特征參數如圖5所示。

表1 程序驗證分析工況匯總

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

Re=30時的流動為穩態，特征在于位于圓柱后穩定的再循環區域。圓柱體周圍的網格尺寸為Δx=Δy=0.01D，經過模擬后得到的數據，圖5中定義的所有特征幾何參數與表2中的文獻數據相當，經細化網格后模擬的結果誤差小于10%。

表2 Re=30時模擬結果與文獻對比

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

圖6 Re=185時的渦脫落Fig.6 Vorticity countours at Re=185

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

對于不穩定流動兩個不同雷諾數下得到的模擬數據：斯特勞哈爾數，平均阻力(在10個時間段內計算)和均方根升力系數與表3和表4中總結的文獻數據相當。

表3 Re=100時模擬結果與文獻對比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時模擬結果與文獻對比 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隨時間變化圖Fig.7 Temporal evolutions of CD and CL for the 2D flow at Re=100 (a);Re=185 (b)

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

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

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

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

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

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

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

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

3 結論

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