基于浸入邊界法的高解析度CFD-DEM流固耦合方法

毛 佳, 肖景文, 趙蘭浩, 底瑛棠

(河海大學 水利水電學院,南京 210098)

流固耦合問題常見于不同工程領域,根據固體的物理性質,可劃分為流體與連續介質、流體與非連續介質間的流固耦合問題.目前,處理流體與連續介質間的流固耦合問題的數值模擬技術已較為成熟,廣泛應用于水利工程[1]、船舶工程[2]等領域;而流體與非連續介質間的流固耦合,例如滑坡涌浪[3]、冰區航行[4]等問題的數值模擬仍具有較大的挑戰性,主要難點在于固體運動過程中流固耦合界面的準確描述.

有研究者采用特殊的網格處理技術模擬流固耦合運動界面,如動網格[5-6]、嵌套網格[7]等.但由于固體運動的無規律性以及運動過程中相互間的碰撞,采用動網格方法時將會產生網格扭曲,甚至網格拓撲結構發生變化的問題.嵌套網格方法雖然適應性較好,但是此方法在運動固體周圍設置的部件網格限制了其應用范圍,當計算大量塊體碰撞的問題時,各個部件網格將發生重疊,導致計算難以順利進行.因此,這些網格處理技術只適用于模擬少量塊體的非碰撞流固耦合問題,而對于涉及大量離散塊體的流固耦合問題則存在明顯的局限性.

考慮到流固耦合運動界面的模擬難度較高,有研究者提出簡化的歐拉-歐拉多相流模型[8],采用流體本構模型模擬離散塊體,使得固體和流體間的界面成為內部界面,避免了采用特殊的措施處理流體域的移動且未知的邊界.然而,此類方法不能體現固體的非連續性,較適用于解決固體呈現出顯著流體化特征的問題.實際上,流體與固體具有截然不同的物理力學性質,需要采用不同的方程進行描述[9],因此,有研究者提出了屬于歐拉-拉格朗日方法的離散相模型(DPM)[10].此方法能夠考慮固體的非連續性,但由于無法模擬固體間的相互碰撞,只能用于模擬稀疏兩相流問題.

近年來,國內外研究者普遍采用耦合數值方法模擬流固耦合問題,即將基于歐拉框架的計算流體力學(CFD)方法與基于拉格朗日框架的離散元法(DEM)耦合,建立CFDEM流固耦合數值模擬方法.目前的耦合方法大多基于經驗公式計算流固耦合作用力,實際上可以將這些耦合數值模擬方法看作DPM的一種改進,通常將此類方法稱為非解析的(Unresolved)CFDEM耦合方法[11-12].非解析的CFDEM耦合方法是將計算流體力學方法與非連續介質方法相結合進行流固耦合問題研究的一次有益嘗試,以簡化的方式處理固體與流體的相互作用力,回避了兩種介質之間邊界移動且未知的問題,能夠有效描述非連續介質在流場中的宏觀運動規律,具有計算量小的優點.但與實際物理過程相比,該方法仍存在一些明顯的不足:① 目前的經驗公式僅限于圓形顆粒;② 流固耦合邊界條件無法準確滿足;③ 多顆粒流固耦合問題中,尾流對顆粒運動的影響難以充分反映.

為了解決非解析流固耦合方法的上述缺點,本文提出高解析度(Resolved)CFD-DEM流固耦合方法.采用Navier-Stokes方程組作為控制方程模擬流體的運動,并在固定網格上采用投影法求解流場運動變量;采用離散單元法模擬非連續固體的任意運動;引入浸入邊界法(IBM)表示固體和流體間的移動邊界,并計算兩者間的相互作用力.通過強耦合迭代求解算法,提高計算的收斂性和穩定性,獲得高精度的流場特性,準確反映固體運動過程中的流固耦合界面.

1 CFD-DEM流固耦合方法

相比于基于經驗公式的非解析流固耦合方法,本文提出的高解析度方法能夠準確反映流固耦合作用力,同時,在流固耦合界面上嚴格滿足速度邊界條件,能夠在塊體周圍獲得高解析度的流場.本文提出的基于浸入邊界法的CFD-DEM流固耦合方法示意圖如圖1所示.

圖1 基于浸入邊界法的CFD-DEM流固耦合方法示意圖Fig.1 Sketch of resolved CFD-DEM algorithm based on immersed boundary method

采用基于歐拉框架的CFD描述流體的運動,采用基于拉格朗日框架的DEM描述固體的運動及碰撞,在離散單元的表面布置浸入邊界點,解決固體運動過程中與流體間的移動且未知的邊界問題.為了實現流體域與固體域間的強耦合,提出流固耦合系統迭代求解方案.

1.1 流體域數值模擬

在黏性不可壓縮Navier-Stokes方程組右端項中引入附加體力項,得到如下控制方程組:

(1)

(2)

式中:u為流場速度;t為時間;p為壓力;ρ為密度;fb為體力項;τ為黏性應力張量;f為附加體力項.

在直接力法中,附加體力項由動量方程結合浸入邊界點上的速度相容條件推導而來.速度相容條件是指在流體固體交界面上需要滿足速度相同條件,即

Un+1=Vn+1

(3)

式中:Un+1表示將笛卡爾網格點上的流體速度插值到浸入邊界點上得到的速度;Vn+1表示同一浸入邊界點上的固體速度;n表示計算步.

引入插值函數I(φ,Xi),將定義在笛卡爾網格上的物理量插值到浸入邊界點上,同時,引入分配函數D(Φ,x),將定義在浸入邊界點上的物理量分配到笛卡爾網格上,插值函數與分配函數分別滿足如下條件:

(4)

(5)

式中:φ(x)表示定義在笛卡爾網格節點x上的物理量,包括u,p,f等;Φ(Xi)表示定義在浸入邊界點Xi上的物理量;gb表示笛卡爾網格點的集合;NIBP為插值計算某笛卡爾網格節點上物理量時所涉及到的浸入邊界點數目;ΔSi表示浸入邊界點的控制域;δ表示狄拉克函數,在二維情況下滿足δ(x)=d(x)d(y),

(6)

h取為網格尺寸的倍數;r表示有限元網格節點與浸入邊界點的距離.結合式(4)～(6),可實現物理量的插值與分配.

基于特征線法,式(2)的時間離散格式如下:

Δu=Kn+Rn+1+fn+1Δt

(7)

(8)

(9)

在此基礎上,將第n個時間步的速度增量分為兩個部分:

Δu*=Kn+fn+1Δt

(10)

Δu**=Rn+1

(11)

將中間速度場增量式(10)與速度場增量的修正量式(11)代入速度邊界條件Vn+1=Un+1中,整理可得:

fn+1Δt=D[I(fn+1Δt)]=

D[Vn+1-I(un+Kn+Rn+1)]

(12)

由式(12)可見,附加體力項與速度場、壓力場耦合,因此需要通過迭代計算求解上述三類未知量,具體步驟如下.

步驟1忽略壓力及附加體力項,計算中間速度場u*,K,u*,K=un+Kn.

步驟2計算壓力場增量Δpn,k及附加體力項fn+1,k.令pn+1,0=pn,fn+1, 0=0,k=1.

步驟2.1根據附加體力項的計算公式計算fn+1,k.

值得注意的是,當浸入邊界點與網格點重合,式(12)中的fn+1Δt=D[I(fn+1Δt)]成立,附加體力項的計算無需迭代.然而,通常情況下,浸入邊界點與網格點不重合,通過單次計算附加體力項并不能滿足速度邊界條件,因此需要迭代求解體力項.令fn+1,k,0=fn+1, k-1,i=1,根據fn+1,k,iΔt=fn+1,k,i-1×Δt+D[Vn+1-I(un+Kn+Rn+1+fn+1,k,i-1Δt)]計算fn+1, k, i.在此基礎上,計算‖Vn+1-I(un+Kn+Rn+1+fn+1, k, i-1Δt)‖,進行收斂性判斷,其中‖·‖表示任一范數,如果范數小于容差ε,則結束步驟2.1,否則令i=i+1,迭代計算直到滿足收斂條件.最終,fn+1, k=fn+1, k, i.

步驟2.2更新中間速度場u*,計算式為u*=u*, K+fn+1,kΔt.

步驟2.4收斂性檢查.計算‖I(un+Kn+Rn+1,k)-I(un+Kn+Rn+1,k-1)‖,如果小于容差ε,則滿足收斂性檢查,結束步驟2的迭代計算,否則,令k=k+1,重復步驟2.1～2.4,直到滿足收斂條件.

步驟3更新壓力場pn+1=pn+Δpn, k,校正速度場un+1=u*+Rn+1,k.

1.2 固體域的數值模擬

根據牛頓第二定律,固體運動的控制方程如下:

(13)

(14)

Fc可分為法向接觸力Fcn和切向接觸力Fct,Fcn的計算式如下:

Fcn=pf∮Sβc∩βtnSβc∩βt(φc-φt)dS

(15)

式中:pf為罰函數;βc為接觸單元;βt為目標單元;n為重疊域邊界的單位外法向向量;φc、φt分別為定義在接觸單元、目標單元上任一點的勢函數[13];S為重疊區域.

根據力-位移法則計算切向接觸力,具體地,通過計算單元所受法向接觸力合力及其作用點,并依據法向接觸力合力的方向確定當前時間步切向力方向,以位移增量法計算切向接觸力大小[14],詳細步驟如下:

1.3 流固耦合系統迭代求解方案

為了實現流體域與固體域間的強耦合,提高計算的收斂性和穩定性,獲得高精度的流場特性,準確反映固體運動過程中的流固耦合界面,提出了流固耦合系統迭代求解方案.為簡化起見,流體和固體控制方程迭代格式記為

(Δun, k, Δpn, k)=Rf

(16)

Δδn,k=Rs

(17)

式中:Rf表示流體方程右端項,與un、pn、δn+1, k-1相關;Rs表示固體方程右端項,與δn、un+1, k、pn+1, k相關.

假設第n步計算已經完成,現在進入第n+1步的第k次迭代,執行過程如下:

步驟1求解流體方程,更新速度場un+1, k與壓力場pn+1, k.

步驟2求解固體方程,更新固體的位置δn+1, k.

步驟3收斂性檢查.計算‖Δζn, k-Δζn, k-1‖,如果小于容差εζ,則滿足收斂性檢查,結束當前時間步的流固耦合系統迭代計算,否則,令k=k+1,重復步驟1～3,直到滿足收斂條件.其中,ζ可取為變量Δu、Δp以及Δδ,εζ為對應于各變量ζ的給定的小值.

2 數值算例分析

2.1 圓柱繞流渦激振動

圓柱繞流渦激振動[15-16]被廣泛用于驗證流固耦合算法解決固體平動耦合問題的準確性.圓柱具有水平向x及豎直向y兩個方向的自由度,運動方程如下:

(18)

(19)

式中:ξ為阻尼比;U*為折合流速;D為圓柱直徑;m*=ms/mf為質量比,ms、mf分別為圓柱質量與同等體積下的流體質量;CL、CD分別為升力系數、阻力系數;X、Y為圓柱無量綱位移,分別滿足X=x/D、Y=y/D.

計算域大小為[-15D, 45D]×[-20D, 20D],初始時刻圓柱圓心坐標為(0, 0).雷諾數Re=200,ξ=0.01,U*=5.0,m*=4/π.計算域左側邊界為均勻入流邊界,右側邊界為自由出流邊界,上、下邊界為滑移固壁.為分析網格敏感性,采用了3種網格尺寸,分別記為L1、L2及L3,對應地,計算域離散為204×112、408×224及816×448個非結構化四邊形單元,圓柱周圍的網格尺寸分別為D/20、D/40及D/80,在圓柱表面布置753個浸入邊界點.

固定圓柱直至下游產生穩定的卡門渦街,此后釋放圓柱,可在圓柱下游觀察到如圖2所示的渦(圖中Ω為渦量).圓柱振動頻率為0.187,與漩渦脫落頻率接近.

圖2 圓柱下游渦量場(L3)Fig.2 Vortexes field downstream of cylinder (L3)

由于流場拖曳力的作用,圓柱的振動中心偏離(0, 0),本文方法計算得到的圓柱振動中心為(0.632D, 0),與參考文獻[15]中的(0.62D, 0)接近.為方便對比,將圓柱振動中心平移至(0, 0),在此基礎上繪制了圓柱“8字形”振動軌跡對比圖,如圖3所示.由于L1網格較粗,無法精確捕捉下游脫落的漩渦,所以使得x方向與y方向的位移較小,當采用較密網格L2、L3時,本文計算結果與前人的貼體網格計算結果[15]、基于浸入邊界法的強耦合數值模擬結果[16]極為接近,一方面說明了加密網格可顯著提高數值計算的精度,另一方面說明了本文方法模擬固體平動耦合問題的準確性.

圖3 圓柱中心振動軌跡Fig.3 Centerline displacement phase plot

為說明流固耦合迭代對計算效率的影響,比較了計算時間t=3 s時,迭代次數取1和3的耗時.迭代1次的工況耗時為 7 899.38 s,迭代3次的工況耗時為 23 150.77 s,約為迭代1次工況耗時的2.9倍.可見,增加流固耦合迭代次數對計算效率有顯著影響.然而,對于流場的變化或流固耦合作用劇烈的強耦合問題,若采用弱耦合求解格式,即用上一個時間步的力來求解當前步的運動,即使流體、固體的計算收斂,仍不能代表流固耦合系統一定收斂.因此,對于此類問題,即使會損失計算效率,也必須采用強耦合求解格式來保證計算精度.

2.2 方塊馳振

為說明本文方法計算固體轉動耦合問題的準確性,計算了經典的方塊馳振問題[16-18],參考文獻[17]中采用貼體網格計算,參考文獻[16,18]中采用基于浸入邊界法的流固耦合數值模擬方法計算,提供了強耦合模擬結果.方塊寬度為L,厚度為Ds,運動方程如下:

(20)

固定方塊直至下游產生穩定脫落的渦,此后釋放方塊.方塊轉動過程中的流場渦量示意圖如圖4所示.當方塊轉動角度為正最大值及負最大值時,渦量場呈反對稱狀態,如圖4(a)、圖4(c)所示.

圖4 方塊轉動過程中的瞬時渦量場Fig.4 Instantaneous vorticity surrounding rectangular rigid body

表1給出了本文計算得到的方塊最大轉角θmax及轉動頻率frot,當前方法得出的最大轉角為15.4°,與文獻[16]中的結果誤差不超過4.4%;當前方法得出的轉動頻率為0.019 2,與文獻[18]中的結果誤差不超過3.1%.

表1 方塊最大轉角θmax及轉動頻率frot

圖5中繪制了方塊轉角變化過程對比圖.圖中:U為流體自由流速.由于已有參考文獻中方塊釋放時間不一致,初始階段的方塊運動有一定的差異,但 是當轉動逐漸發展穩定之后,各結果間的吻合度較高,說明了本文方法模擬轉動耦合問題的準確性.

圖5 方塊轉角變化過程Fig.5 Evolution of rotational angle

2.3 多塊體沉降

為說明本文方法模擬流體與多固體間的流固耦合作用及固體間碰撞作用的能力,模擬了多塊體沉降問題.

計算域尺寸為[0, 10 m]×[0, 10 m],采用兩類正方形塊體,共計28個,邊長分別為0.698 m、0.548 m,在重力作用下自由下落.初始時刻位置如圖6所示.A點坐標為(0.802 m, 7.945 m),B點坐標為(1.427 m, 6.750 m),同層塊體間的x向距離均為1.100 m,相鄰兩層塊體間y向距離為1.120 m.水體密度為1 000 kg/m3,動力黏性為0.01 kg/(m·s),塊體密度為 2 600 kg/m3.計算域離散為480×480個結構化四邊形單元,在邊長為 0.698 m、0.548 m 的塊體表面分別布置404、320個浸入邊界點.本算例中所有邊界均設為無滑固壁.

圖6 塊體初始時刻位置Fig.6 Initial positions of solids

圖7中給出了塊體運動速度(v)及流體域速度矢量場變化過程.初始階段(見圖7(a)),水體受塊體運動擾動的同時,向塊體施加了作用力.隨著塊體的進一步沉降,水體和塊體間的相互作用逐漸增強(見圖7(b)～圖7(d)),此階段在塊體周圍可以觀察到大量的渦,受渦的影響,塊體在沉降過程中伴隨著轉動.部分塊體在2.8 s左右開始接觸計算域底部,此后流場的擾動顯著降低(見圖7(e)～圖7(f)).此算例說明了本文方法不限于模擬流體和圓形顆粒的相互作用,能夠描述流場的復雜變化及固體間的碰撞,可實現復雜流固耦合問題中固體運動界面的準確描述.

圖7 塊體運動速度及流體域速度矢量場變化過程Fig.7 Instantaneous velocity vector of fluid field and position of rigid bodies

3 結論

提出了基于浸入邊界法的高解析度CFD-DEM流固耦合方法,主要結論如下:

(1) 區別于非解析CFDEM耦合方法依賴于經驗公式計算流固耦合作用力,新方法采用浸入邊界法表示固體和流體間的移動邊界,并計算兩者間的相互作用力.

(2) 計算了圓柱繞流渦激振動、方塊馳振兩個經典算例,并與已有數值計算結果進行比較,驗證了新方法計算平動耦合問題、轉動耦合問題的準確性.

(3) 模擬了多塊體沉降問題,說明了新方法不局限于圓形顆粒,可以有效描述復雜流固耦合問題中的運動界面,獲得高解析度的流場.