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

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

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

有研究者采用特殊的網(wǎng)格處理技術(shù)模擬流固耦合運(yùn)動(dòng)界面,如動(dòng)網(wǎng)格[5-6]、嵌套網(wǎng)格[7]等.但由于固體運(yùn)動(dòng)的無(wú)規(guī)律性以及運(yùn)動(dòng)過(guò)程中相互間的碰撞,采用動(dòng)網(wǎng)格方法時(shí)將會(huì)產(chǎn)生網(wǎng)格扭曲,甚至網(wǎng)格拓?fù)浣Y(jié)構(gòu)發(fā)生變化的問題.嵌套網(wǎng)格方法雖然適應(yīng)性較好,但是此方法在運(yùn)動(dòng)固體周圍設(shè)置的部件網(wǎng)格限制了其應(yīng)用范圍,當(dāng)計(jì)算大量塊體碰撞的問題時(shí),各個(gè)部件網(wǎng)格將發(fā)生重疊,導(dǎo)致計(jì)算難以順利進(jìn)行.因此,這些網(wǎng)格處理技術(shù)只適用于模擬少量塊體的非碰撞流固耦合問題,而對(duì)于涉及大量離散塊體的流固耦合問題則存在明顯的局限性.

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

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

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

1 CFD-DEM流固耦合方法

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

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

采用基于歐拉框架的CFD描述流體的運(yùn)動(dòng),采用基于拉格朗日框架的DEM描述固體的運(yùn)動(dòng)及碰撞,在離散單元的表面布置浸入邊界點(diǎn),解決固體運(yùn)動(dòng)過(guò)程中與流體間的移動(dòng)且未知的邊界問題.為了實(shí)現(xiàn)流體域與固體域間的強(qiáng)耦合,提出流固耦合系統(tǒng)迭代求解方案.

1.1 流體域數(shù)值模擬

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

(1)

(2)

式中:u為流場(chǎng)速度;t為時(shí)間;p為壓力;ρ為密度;fb為體力項(xiàng);τ為黏性應(yīng)力張量;f為附加體力項(xiàng).

在直接力法中,附加體力項(xiàng)由動(dòng)量方程結(jié)合浸入邊界點(diǎn)上的速度相容條件推導(dǎo)而來(lái).速度相容條件是指在流體固體交界面上需要滿足速度相同條件,即

Un+1=Vn+1

(3)

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

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

(4)

(5)

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

(6)

h取為網(wǎng)格尺寸的倍數(shù);r表示有限元網(wǎng)格節(jié)點(diǎn)與浸入邊界點(diǎn)的距離.結(jié)合式(4)～(6),可實(shí)現(xiàn)物理量的插值與分配.

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

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

(7)

(8)

(9)

在此基礎(chǔ)上,將第n個(gè)時(shí)間步的速度增量分為兩個(gè)部分:

Δu*=Kn+fn+1Δt

(10)

Δu**=Rn+1

(11)

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

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

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

(12)

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

步驟1忽略壓力及附加體力項(xiàng),計(jì)算中間速度場(chǎng)u*,K,u*,K=un+Kn.

步驟2計(jì)算壓力場(chǎng)增量Δpn,k及附加體力項(xiàng)fn+1,k.令pn+1,0=pn,fn+1, 0=0,k=1.

步驟2.1根據(jù)附加體力項(xiàng)的計(jì)算公式計(jì)算fn+1,k.

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

步驟2.2更新中間速度場(chǎng)u*,計(jì)算式為u*=u*, K+fn+1,kΔt.

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

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

1.2 固體域的數(shù)值模擬

根據(jù)牛頓第二定律,固體運(yùn)動(dòng)的控制方程如下:

(13)

(14)

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

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

(15)

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

根據(jù)力-位移法則計(jì)算切向接觸力,具體地,通過(guò)計(jì)算單元所受法向接觸力合力及其作用點(diǎn),并依據(jù)法向接觸力合力的方向確定當(dāng)前時(shí)間步切向力方向,以位移增量法計(jì)算切向接觸力大小[14],詳細(xì)步驟如下:

1.3 流固耦合系統(tǒng)迭代求解方案

為了實(shí)現(xiàn)流體域與固體域間的強(qiáng)耦合,提高計(jì)算的收斂性和穩(wěn)定性,獲得高精度的流場(chǎng)特性,準(zhǔn)確反映固體運(yùn)動(dòng)過(guò)程中的流固耦合界面,提出了流固耦合系統(tǒng)迭代求解方案.為簡(jiǎn)化起見,流體和固體控制方程迭代格式記為

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

(16)

Δδn,k=Rs

(17)

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

假設(shè)第n步計(jì)算已經(jīng)完成,現(xiàn)在進(jìn)入第n+1步的第k次迭代,執(zhí)行過(guò)程如下:

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

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

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

2 數(shù)值算例分析

2.1 圓柱繞流渦激振動(dòng)

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

(18)

(19)

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

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

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

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

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

圖3 圓柱中心振動(dòng)軌跡Fig.3 Centerline displacement phase plot

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

2.2 方塊馳振

為說(shuō)明本文方法計(jì)算固體轉(zhuǎn)動(dòng)耦合問題的準(zhǔn)確性,計(jì)算了經(jīng)典的方塊馳振問題[16-18],參考文獻(xiàn)[17]中采用貼體網(wǎng)格計(jì)算,參考文獻(xiàn)[16,18]中采用基于浸入邊界法的流固耦合數(shù)值模擬方法計(jì)算,提供了強(qiáng)耦合模擬結(jié)果.方塊寬度為L(zhǎng),厚度為Ds,運(yùn)動(dòng)方程如下:

(20)

固定方塊直至下游產(chǎn)生穩(wěn)定脫落的渦,此后釋放方塊.方塊轉(zhuǎn)動(dòng)過(guò)程中的流場(chǎng)渦量示意圖如圖4所示.當(dāng)方塊轉(zhuǎn)動(dòng)角度為正最大值及負(fù)最大值時(shí),渦量場(chǎng)呈反對(duì)稱狀態(tài),如圖4(a)、圖4(c)所示.

圖4 方塊轉(zhuǎn)動(dòng)過(guò)程中的瞬時(shí)渦量場(chǎng)Fig.4 Instantaneous vorticity surrounding rectangular rigid body

表1給出了本文計(jì)算得到的方塊最大轉(zhuǎn)角θmax及轉(zhuǎn)動(dòng)頻率frot,當(dāng)前方法得出的最大轉(zhuǎn)角為15.4°,與文獻(xiàn)[16]中的結(jié)果誤差不超過(guò)4.4%;當(dāng)前方法得出的轉(zhuǎn)動(dòng)頻率為0.019 2,與文獻(xiàn)[18]中的結(jié)果誤差不超過(guò)3.1%.

表1 方塊最大轉(zhuǎn)角θmax及轉(zhuǎn)動(dòng)頻率frot

圖5中繪制了方塊轉(zhuǎn)角變化過(guò)程對(duì)比圖.圖中:U為流體自由流速.由于已有參考文獻(xiàn)中方塊釋放時(shí)間不一致,初始階段的方塊運(yùn)動(dòng)有一定的差異,但 是當(dāng)轉(zhuǎn)動(dòng)逐漸發(fā)展穩(wěn)定之后,各結(jié)果間的吻合度較高,說(shuō)明了本文方法模擬轉(zhuǎn)動(dòng)耦合問題的準(zhǔn)確性.

圖5 方塊轉(zhuǎn)角變化過(guò)程Fig.5 Evolution of rotational angle

2.3 多塊體沉降

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

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

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

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

圖7 塊體運(yùn)動(dòng)速度及流體域速度矢量場(chǎng)變化過(guò)程Fig.7 Instantaneous velocity vector of fluid field and position of rigid bodies

3 結(jié)論

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

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

(2) 計(jì)算了圓柱繞流渦激振動(dòng)、方塊馳振兩個(gè)經(jīng)典算例,并與已有數(shù)值計(jì)算結(jié)果進(jìn)行比較,驗(yàn)證了新方法計(jì)算平動(dòng)耦合問題、轉(zhuǎn)動(dòng)耦合問題的準(zhǔn)確性.

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