基于隱式擴(kuò)散的直接力格式浸沒邊界格子Boltzmann 方法1)

佟 瑩 * 夏 健 * 陳 龍 *, 薛浩天 *

* (南京航空航天大學(xué)航空學(xué)院非定常空氣動(dòng)力學(xué)與流動(dòng)控制工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室,南京 210016)

引言

浸沒邊界(immersed boundary,IB)方法作為一種通用的數(shù)學(xué)框架和方法論在動(dòng)邊界繞流問題中廣泛應(yīng)用[1-4].與常規(guī)貼體網(wǎng)格方法相比,針對(duì)復(fù)雜外形幾何體經(jīng)歷大變形的流固耦合動(dòng)力學(xué)問題,IB 方法在計(jì)算效率和適應(yīng)性等方面具備顯著優(yōu)勢(shì).主要原因有兩個(gè):第一,IB 方法采用兩套彼此獨(dú)立的網(wǎng)格分別離散流體域和浸沒邊界,從而避免了繁瑣復(fù)雜的網(wǎng)格生成和再生過(guò)程;第二,作為一種數(shù)值格式,IB 方法通過(guò)向流體動(dòng)量方程中添加體力項(xiàng)表征界面邊界,無(wú)需對(duì)流動(dòng)求解器代碼進(jìn)行過(guò)多修改,邊界條件施加簡(jiǎn)單.

近年來(lái),基于介觀格子Boltzmann 模型(lattice Boltzmann,LB)發(fā)展流動(dòng)求解器的計(jì)算流體力學(xué)方法(computational fluid dynamics,CFD)被使用處理各類復(fù)雜流動(dòng)問題[5-10].相比于傳統(tǒng)CFD 方法數(shù)值求解宏觀守恒方程,LB 模型通過(guò)介觀流體微團(tuán)的碰撞和遷移描述流體動(dòng)力學(xué)演化過(guò)程,控制方程格式更加簡(jiǎn)單、易于編程、邊界處理方式簡(jiǎn)潔高效、并行特性更好.基于使用笛卡爾網(wǎng)格的共同特征,研究者試圖將LB 與IB 集成為一種靈活處理復(fù)雜邊界問題的CFD 方法.2004 年,Feng 和Michaelides[11]將罰IB 模型與LB 模型耦合,首次提出IB-LB 模型,并將其應(yīng)用于顆粒流.這種方法采用一種反饋機(jī)制描述剛性邊界對(duì)流動(dòng)的影響,算法穩(wěn)定性依賴于恢復(fù)力中引入的自由參數(shù),嚴(yán)重影響了模型計(jì)算效率和數(shù)值精度.隨后,直接力格式的IB 模型被引入LB 模型[12].雖然該模型消除了自由參數(shù)對(duì)時(shí)間步長(zhǎng)的約束,但由于重新引入NS 方程中復(fù)雜的對(duì)流項(xiàng),LB 模型的優(yōu)勢(shì)被破壞.2006 年,Niu 等[13]提出動(dòng)量交換格式的IB-LB 模型,基于分布函數(shù)的動(dòng)量交換準(zhǔn)則計(jì)算邊界力,完整保留了LB 模型的計(jì)算優(yōu)勢(shì).然而,這些傳統(tǒng)IB-LB 模型中,顯式計(jì)算的體積力僅在收斂狀態(tài)下保證速度場(chǎng)滿足無(wú)滑移條件[14],導(dǎo)致模擬結(jié)果呈現(xiàn)流線穿透物面的偽物理現(xiàn)象.

為了避免流線穿透物面,近15 年來(lái)已經(jīng)發(fā)展了各種改進(jìn)的IB-LB 模型.Wu 和Shu[15]提出一種隱式擴(kuò)散界面格式的速度修正IB-LB 模型,通過(guò)向邊界周圍流體節(jié)點(diǎn)施加速度修正強(qiáng)制無(wú)滑移條件.但即使對(duì)于二維流動(dòng),直接求解矩陣方程仍需較大的內(nèi)存占用,且矩陣方程可解性強(qiáng)烈依賴于拉格朗日節(jié)點(diǎn)分布.Luo 等[16]提出多段直接力IB 模型,通過(guò)施加多重直接力強(qiáng)制無(wú)滑移約束,可避免復(fù)雜的矩陣求逆,且穩(wěn)定性不受拉格朗日節(jié)點(diǎn)分布限制.Kang和Hassan[17]將文獻(xiàn)[16]的多段直接力IB 模型與分裂力格式的LB 方程結(jié)合,發(fā)展了隱式擴(kuò)散界面的直接力IB-LB 模型.但是針對(duì)動(dòng)邊界繞流問題,每個(gè)時(shí)間步中經(jīng)固定次數(shù)累積施加的邊界力嚴(yán)重阻礙了大尺度流動(dòng)問題的計(jì)算效率.Dash 等[18]借助附加參數(shù)靈活控制每個(gè)時(shí)間步中邊界力的迭代次數(shù),在保證無(wú)滑移條件的基礎(chǔ)上大幅減少計(jì)算時(shí)間.Yang 等[19]在隱式擴(kuò)散界面背景下改進(jìn)了動(dòng)量交換IB-LB 模型.隨后幾種以分布函數(shù)為操作對(duì)象的IB-LB 模型[20-22]被發(fā)展,相比于以宏觀變量為操作對(duì)象,這些修正方案增加了原有IB 模型的計(jì)算負(fù)擔(dān).

盡管上述IB-LB 模型的數(shù)值預(yù)測(cè)能力在剛性體繞流問題中已經(jīng)取得了較大進(jìn)步,但在可變形動(dòng)邊界繞流問題中仍有待改進(jìn),例如非定常流體力的偽物理震蕩.如圖1 所示,拉格朗日變量描述的邊界力擴(kuò)散向周圍流體節(jié)點(diǎn)施加體積力,歐拉變量描述的體積力重構(gòu)流體速度場(chǎng)使邊界節(jié)點(diǎn)處流體速度和邊界速度之間滿足無(wú)滑移條件.根據(jù)Peskin[1]關(guān)于IB 數(shù)學(xué)框架的工作總結(jié),上述信息交互應(yīng)滿足歐拉/拉格朗日變量同一性,即“擴(kuò)散”應(yīng)保證歐拉變量描述的體積力和拉格朗日變量描述的邊界力等價(jià),“插值”應(yīng)保證歐拉流體速度和拉格朗日邊界速度等價(jià).現(xiàn)有IB-LB 模型中,強(qiáng)制無(wú)滑移保證了速度同一性.然而,體積力和邊界力之間的同一性沒有被強(qiáng)制,這使得預(yù)測(cè)的流體力呈現(xiàn)偽物理震蕩.此外,根據(jù)Huang 和Wu[23]關(guān)于LB 模型的三階Chapman-Enskog分析,恢復(fù)的宏觀方程中出現(xiàn)與體積力二階梯度相關(guān)的附加項(xiàng),其自由參數(shù)與邊界力模型直接相關(guān),即邊界力和體積力之間的關(guān)聯(lián)將直接影響流動(dòng)計(jì)算的數(shù)值精度.

圖1 浸沒邊界方法原理圖.○ 代表拉格朗日邊界節(jié)點(diǎn);●代表歐拉流體節(jié)點(diǎn)Fig.1 Schematic diagram of the immersed boundary method.○ indicate the Lagrangian boundary point;● indicate the Eulerian fluid node

采用隱式擴(kuò)散界面,本文提出一種改進(jìn)的直接力格式IB-LB 模型,其中LB 背景下滿足力同一性的邊界力表達(dá)式基于經(jīng)典IB-NS 模型[1]中的歐拉/拉格朗日變量同一性原則推導(dǎo).關(guān)聯(lián)邊界力與無(wú)滑移約束的線性方程組采用Richardson 迭代方法數(shù)值求解,有效降低矩陣求逆過(guò)程引起的內(nèi)存占用和計(jì)算負(fù)擔(dān).由附加參數(shù)靈活控制子迭代,保證每次時(shí)間推進(jìn)動(dòng)邊界上所有節(jié)點(diǎn)的無(wú)滑移約束.根據(jù)矢量場(chǎng)映射建立描述插值和擴(kuò)散操作的轉(zhuǎn)換矩陣,簡(jiǎn)化信息交互界面上的數(shù)據(jù)交互計(jì)算,有效降低IB 模型的代碼編寫難度.基于改進(jìn)的IB-LB 模型開發(fā)流動(dòng)求解器,數(shù)值計(jì)算能力方面完成兩個(gè)關(guān)鍵改進(jìn):一是數(shù)值解具備完整的二階精度;二是抑制了耦合界面上流體力的非物理震蕩.本文改進(jìn)的IB-LB 模型以期為大變形結(jié)構(gòu)體與黏性流體之間的流固耦合動(dòng)力學(xué)算法開發(fā)提供理論參考.

1 物理模型

如圖1 所示為擴(kuò)散界面背景下IB 方法的原理圖,流體域Ω采用笛卡爾網(wǎng)格離散,浸入流體的物體邊界Γ由獨(dú)立于背景網(wǎng)格的拉格朗日點(diǎn)替代.歐拉流體變量和拉格朗日邊界變量通過(guò)狄拉克函數(shù)關(guān)聯(lián),相互作用方程如下

式中,u,f,F分別為流體速度矢量、體力密度矢量和邊界力密度矢量,δ代表狄拉克函數(shù),x和X為位置矢量(x∈Ω,X∈Γ),s為拉格朗日坐標(biāo),q為邊界有效寬度(dxdy=qds[24]).如圖1 所示,從歐拉流體節(jié)點(diǎn)向拉格朗日邊界節(jié)點(diǎn)插值速度的過(guò)程可理解為:以拉格朗日節(jié)點(diǎn)為中心,影響域內(nèi)所有歐拉節(jié)點(diǎn)速度的加權(quán)平均.同理,從拉格朗日邊界節(jié)點(diǎn)向歐拉流體節(jié)點(diǎn)擴(kuò)散體積力的過(guò)程可以理解為:以歐拉節(jié)點(diǎn)為中心,影響域內(nèi)所有拉格朗日節(jié)點(diǎn)上邊界力的加權(quán)平均.方程式(1)和式(2)描述了浸沒邊界和流體之間的相互作用本質(zhì).

2 計(jì)算模型

2.1 格子Boltzmann (LB)模型

LB 模型使用介觀密度分布函數(shù)的時(shí)空演化描述流體運(yùn)動(dòng),宏觀守恒方程中的基礎(chǔ)流動(dòng)變量(密度ρ,速度u,壓力p)由分布函數(shù)的矩條件簡(jiǎn)單計(jì)算.采用D2Q9 模型[25]離散速度空間 {cα} 定義的速度矢量如下

式中,c=Δx/Δt為格子速度,Δx和Δt分別為格子空間尺度和時(shí)間尺度.對(duì)離散速度空間中考慮體力項(xiàng)映射的Boltzmann 方程沿特征線積分,采用二階梯形法逼近,獲得分布函數(shù)演化方程[26]如下

取Maxwell 分布函數(shù)在離散速度空間 {cα} 中的二階速度截?cái)嗟镁植科胶夥植己瘮?shù)的表達(dá)式

式中,ωα是 {cα} 的加權(quán)系數(shù),ω0=4/9,ω1～4=1/9,ω5～8=1/36,cs是模型的格子聲速,.

Fα表達(dá)式為

基于上述離散速度模型式(3),分布函數(shù)演化方程(4),平衡分布函數(shù)式(5)及體力分布函數(shù)式(6),采用Chapman-Enskog 展開技術(shù)可恢復(fù)IB 形式的NS 方程組,且能保證空間和時(shí)間上具備二階精度[27].其中,流體黏度v與弛豫時(shí)間τ滿足如下關(guān)聯(lián)

分布函數(shù)演化方程(4)右側(cè)碰撞算子描述了流體微團(tuán)碰撞對(duì)分布函數(shù)的影響,Fα則描述了宏觀體積力在離散速度空間中的投影對(duì)分布函數(shù)的貢獻(xiàn),即IB-LB 數(shù)值格式的時(shí)間推進(jìn)可分解為3 個(gè)階段進(jìn)行.

(1) 碰撞

(2) 邊界條件施加

(3)遷移

流體的密度ρ,速度u,壓力p由遷移后更新的fα的矩方程和等溫LB 狀態(tài)方程計(jì)算

本文采用體力密度的半隱式迭代方案構(gòu)造流動(dòng)收斂解.首先給定流場(chǎng)初始條件(ρ0,u0)及體力密度,通過(guò)局部平衡分布函數(shù)式(5)及體力分布函數(shù)式(6)構(gòu)造初始分布函數(shù),并將與相關(guān)的分布函數(shù)、密度和速度表示為,和.由于給定的為預(yù)估,不能準(zhǔn)確滿足無(wú)滑移條件,因此體積力f和速度u需按照滿足無(wú)滑移速度約束的邊界力進(jìn)行修正,這個(gè)過(guò)程由IB 模型來(lái)完成.修正后的體力密度和速度為

使用式(14)左側(cè)的f作為新的重復(fù)上述操作至收斂.

2.2 信息交互界面

歐拉變量和拉格朗日變量之間的信息交互由離散函數(shù)Dlk控制,它被假設(shè)為單值變量函數(shù)φ的張量積[1],表達(dá)式如下

式中,1≤l≤NE,1≤k≤NL.NL為邊界上拉格朗日節(jié)點(diǎn)總數(shù),NE為所有拉格朗日節(jié)點(diǎn)影響域覆蓋的不重合歐拉節(jié)點(diǎn)總數(shù),影響域尺寸由φ確定.

為保證信息交互界面上相互作用方程式(1)和式(2)的離散具備二階精度,本文采用式(17)所示的自身滿足二階連續(xù)、0～2 階矩條件及奇偶條件,一階導(dǎo)數(shù) φ′滿足0～2 階矩條件的分段光順單值函數(shù) φ[28]構(gòu)造Dlk.文獻(xiàn)結(jié)果表明[28],基于式(17)構(gòu)造的Dlk,直接力格式的IB-NS 模型的數(shù)值解具備二階精度.

可變形運(yùn)動(dòng)邊界上每個(gè)拉格朗日節(jié)點(diǎn)具有獨(dú)立的運(yùn)動(dòng)模式,導(dǎo)致大曲率邊界附近影響域重疊瞬時(shí)變化.為保證每個(gè)時(shí)間步中任意邊界節(jié)點(diǎn)處的無(wú)滑移條件能夠被準(zhǔn)確施加,本文發(fā)展了一種考慮影響域重疊的界面構(gòu)造方案,其中心思想是應(yīng)用離散函數(shù)Dlk構(gòu)造描述信息交互的轉(zhuǎn)換矩陣,耦合相同類型節(jié)點(diǎn)間的相互作用.以矢量場(chǎng)映射的方式執(zhí)行插值和擴(kuò)散操作,可保證邊界上所有拉格朗日節(jié)點(diǎn)同時(shí)滿足無(wú)滑移約束,有效抑制影響域重疊引起的震蕩源,從而在物面附近生成逐點(diǎn)光滑的速度場(chǎng)

轉(zhuǎn)換矩陣DE和DL的表達(dá)式如下

對(duì)于剛性靜止壁面,DE與DL只需計(jì)算一次;對(duì)于運(yùn)動(dòng)變形界面,每個(gè)時(shí)間步,由結(jié)構(gòu)運(yùn)動(dòng)更新邊界位置 {X1,X2,···,Xk,···,XNL} 后,需要重新構(gòu)造影響域覆蓋的歐拉節(jié)點(diǎn)集合 {x1,x2,···,xl,···,xNE},并更新NE,然后重新建立DE與DL.

基于上述DE和DL改寫相互作用方程(1)和式(2)如下

式中,uE={u1,u2,···,uNE}和fE={f1,f2,···,fNE} 是存儲(chǔ)在GE中的速度場(chǎng)矢量和體積力場(chǎng)矢量,類似地,UL={U1,U2,···,UNL}和FL={F1,F2,···,FNL}是存儲(chǔ)在GL中的邊界速度場(chǎng)矢量和邊界力密度場(chǎng)矢量.這里,GL表示拉格朗日點(diǎn)構(gòu)成的邊界點(diǎn)集,GE表示所有拉格朗日節(jié)點(diǎn)影響域覆蓋的不重合歐拉流體點(diǎn)集.式(20)和式(21)以場(chǎng)傳遞的方式描述了流體-邊界之間的相互作用.從編程的角度,基于上述表達(dá),可以在代碼中輕松實(shí)現(xiàn)并行計(jì)算并降低內(nèi)存占用.此外,通過(guò)轉(zhuǎn)換矩陣DE和DL,當(dāng)前的擴(kuò)散界面極易拓展至三維流動(dòng).

2.3 邊界力模型

計(jì)算拉格朗日節(jié)點(diǎn)處的邊界力F并將其擴(kuò)散至影響域內(nèi)的歐拉節(jié)點(diǎn)是直接力格式IB 模型的核心.受Peskin[1]工作的啟發(fā),IB 模型的信息交互過(guò)程需滿足變量同一性.因此,本文基于LB 框架推導(dǎo)了滿足變量同一性的邊界力公式,推導(dǎo)過(guò)程如下.

首先,根據(jù)LB 模型的動(dòng)量方程(12)構(gòu)造體積力密度f(wàn)的直接力表達(dá)式

式中,δu=u-u*,δu為歐拉節(jié)點(diǎn)處f引發(fā)的流體速度修正量,為L(zhǎng)B 模型的預(yù)計(jì)算速度,u為重構(gòu)的流體速度.

根據(jù)力同一性的概念,歐拉/拉格朗日變量描述的浸沒物體經(jīng)歷的合流體力等價(jià),即

將體積力密度f(wàn)的直接力表達(dá)式(22)帶入歐拉變量描述的合流體力,可以發(fā)現(xiàn)

綜上得邊界力密度F的表達(dá)式

GL中拉格朗日節(jié)點(diǎn)處流體的質(zhì)量密度矢量ML={M1,M2,···,MNL}來(lái)自歐拉流體節(jié)點(diǎn)的質(zhì)量密度插值

為了保證體積力重構(gòu)的速度場(chǎng)滿足無(wú)滑移約束,建立如下方程組

這里,系數(shù)矩陣A=DE(Δt/2ρE)DL不僅與GL和GE的位置數(shù)據(jù)有關(guān),還與GE中存儲(chǔ)的密度信息有關(guān);FL是線性方程組中待求解的未知量.

采用Richardson 迭代方法[19]數(shù)值求解上述線性方程組

邊界力密度FL由子迭代的累和計(jì)算

基于上述改進(jìn)的IB-LB 模型發(fā)展流動(dòng)求解器,算法流程展示于圖2.IB 模型以運(yùn)動(dòng)界面的位置、速度信息和LB 模型預(yù)計(jì)算的流體密度和速度場(chǎng)信息為輸入變量,任意拉格朗日節(jié)點(diǎn)處的局部流體力H作為邊界力F的作用反力由IB 模型直接輸出,H(Xk,t)=-F(Xk,t)qΔsk,而無(wú)需復(fù)雜的流體應(yīng)力積分計(jì)算,耦合算法穩(wěn)定性更高.針對(duì)變形體繞流問題,非定常局部流體力H在其氣動(dòng)性能分析中十分關(guān)鍵,本文方法在此類問題的流動(dòng)模擬中具有很大的應(yīng)用前景.

圖2 直接力格式IB-LB 模型算法流程圖.淺灰色框指明物面參數(shù)輸出子步驟Fig.2 An overview of one cycle of the IB-LB algorithm with the direct force scheme.The light grey box shows the output sub-step of the surface variables

3 數(shù)值方法驗(yàn)證與分析

本節(jié)以二維Taylor-Green 渦流、靜止圓柱繞流、振蕩圓柱繞流及波動(dòng)翼型繞流4 個(gè)經(jīng)典案例的模擬結(jié)果佐證改進(jìn)的IB-LB 模型的數(shù)值精度、預(yù)測(cè)結(jié)果可靠性、及其在復(fù)雜運(yùn)動(dòng)界面繞流問題中的實(shí)用性.

3.1 Taylor-Green 渦流

作為評(píng)估算法數(shù)值精度的基準(zhǔn)算例,Taylor-Green 渦流已由幾種不同的IB-LB 模型[15,17,19]執(zhí)行模擬.理想的IB 模型應(yīng)具備保持背景LB 求解器固有數(shù)值精度的能力.本節(jié)采用Taylor-Green 渦流模擬結(jié)果來(lái)驗(yàn)證本文改進(jìn)的IB 模型對(duì)背景LB 求解器精度的影響.Taylor-Green 渦流在[-L,L]×[-L,L]的計(jì)算域中解析解表達(dá)式如下

式中,雷諾數(shù)Re根據(jù)計(jì)算域尺寸L和名義速度u0計(jì)算,Re=Lu0/ν .為便于比較,本文使用與文獻(xiàn)[15,17,19]中相同的參數(shù)設(shè)置模擬,其中τ=0.65,Re=10,4 種網(wǎng)格尺度(L/h=10,20,40,80)被使用.取t=0 時(shí)刻速度和密度的解析解為初始條件,外邊界施加周期性條件以確保系統(tǒng)的質(zhì)量和動(dòng)量守恒.以計(jì)算域中是否添加拉格朗日節(jié)點(diǎn)分別執(zhí)行LB 模型模擬和IB-LB 模型模擬,計(jì)算至t=L/u0時(shí)刻流場(chǎng)狀態(tài).拉格朗日節(jié)點(diǎn)處的瞬時(shí)由解析解式(31)構(gòu)造.

式(32) 定義了速度的全局?jǐn)?shù)值誤差L2,這里,N是計(jì)算域中歐拉流體節(jié)點(diǎn)的總數(shù),unum和uexa分別為速度的數(shù)值解和精確解.數(shù)值精度通過(guò)數(shù)值誤差L2與網(wǎng)格尺度h的對(duì)數(shù)斜率來(lái)量化

為了說(shuō)明本文模擬結(jié)果在精度方面的改進(jìn)不針對(duì)擴(kuò)散界面中使用的某一特定單值函數(shù),分別采用4-point cosine[15],4-point piecewise[17]及5-point smoothed piecewise (式(17))單值函數(shù)構(gòu)造滿足力同一性的IB-LB 模型執(zhí)行本算例.圖3 展示了對(duì)數(shù)坐標(biāo)下的數(shù)值誤差L2與網(wǎng)格間距h,3 組IB-LB 模型數(shù)據(jù)與LB 模型數(shù)據(jù)幾乎完全重合.其中,無(wú)拉格朗日節(jié)點(diǎn)嵌入的LB 模型的斜率為2.017,嵌入拉格朗日節(jié)點(diǎn)的IB-LB 模型的斜率分別為2.010,2.013 和2.014.基于相同的LB 背景,速度修正IB-LB 模型[15]采用4-point cosine 單值函數(shù)構(gòu)造Dlk,數(shù)值精度為1.9,多段直接力IB-LB 模型[17]和動(dòng)量交換IB-LB 模型[19]采用4-point piecewise 單值函數(shù)構(gòu)造Dlk,數(shù)值精度分別為1.98 和1.99.當(dāng)前結(jié)果表明,力同一性可改善IB 模型對(duì)LB 結(jié)果精度的弱化,改進(jìn)的IB 模型完整保留了背景LB 模型求解器的數(shù)值精度.同時(shí),當(dāng)前IB-LB 模型流動(dòng)求解器代碼可以為本文計(jì)算模型提供可靠的數(shù)值解,計(jì)算模型具有二階數(shù)值精度.

圖3 對(duì)數(shù)坐標(biāo)下的速度誤差L2 與網(wǎng)格尺度hFig.3 Numerical error L2of velocity versus mesh spacing hfor Taylor-Green flow using the present IB-LB model with different φand the standard LB model

3.2 靜止圓柱繞流

因其流動(dòng)結(jié)構(gòu)由雷諾數(shù)唯一確定,圓柱繞流作為基準(zhǔn)算例在復(fù)雜邊界處理問題中被廣泛應(yīng)用.設(shè)定流動(dòng)參數(shù)U∞=0.1 和D=1.0,通過(guò)改變運(yùn)動(dòng)粘度v調(diào)整雷諾數(shù)Re=DU∞/ν,一系列模擬被執(zhí)行.計(jì)算域尺寸為30D× 15D,圓柱位于計(jì)算域?qū)ΨQ線,圓柱中心到入口邊界的距離為13.5D.入口邊界施加Dirichlet 條件,出口邊界施加Neumann 條件,自由剪切條件應(yīng)用于側(cè)邊界.流場(chǎng)初始條件指定為ρ0=1.0 和u0=(U∞,0).采用網(wǎng)格間距h=1/N均勻劃分背景網(wǎng)格,N為圓柱直徑占據(jù)的格子數(shù),本節(jié)算例N=75,網(wǎng)格收斂指數(shù)(GCI)結(jié)果顯示當(dāng)前網(wǎng)格間距滿足收斂性要求.

首先,通過(guò)與文獻(xiàn)結(jié)果進(jìn)行比較驗(yàn)證了本文改進(jìn)的IB-LB 模型模擬結(jié)果的準(zhǔn)確性.表1 列出Re=20 和Re=40 條件下無(wú)量綱規(guī)則區(qū)長(zhǎng)度2Lw/D和阻力系數(shù)的數(shù)值解,其中FD為圓柱阻力,.對(duì)比結(jié)果表明,本文結(jié)果與現(xiàn)有研究中的數(shù)值模擬結(jié)果[15,17,19-20]具有良好的一致性.我們進(jìn)一步調(diào)整計(jì)算域至40D× 40D(具有相同的網(wǎng)格間距)重復(fù)當(dāng)前模擬,闡明外邊界對(duì)計(jì)算結(jié)果的影響.結(jié)果表明,橫向距離的增加會(huì)導(dǎo)致阻力系數(shù)降低,循環(huán)區(qū)長(zhǎng)度增加,但相對(duì)誤差很小,即30D× 15D的計(jì)算域尺寸下模擬結(jié)果可靠.

表1 Re=20 和40 時(shí)規(guī)則區(qū)長(zhǎng)度和阻力系數(shù)與已有文獻(xiàn)結(jié)果的比較Table 1 Comparison of the obtained recirculating length,the drag coefficient with the previous results at Re=20 and 40

圖4 展示了Re=20 和40 條件下,二維圓柱表面壓強(qiáng)系數(shù)CP與方向角θ的關(guān)系曲線,方向角自前駐點(diǎn)θ=0°經(jīng)圓柱上方指向后基點(diǎn)θ=180°,Cp(θ)=.圖4 中壓強(qiáng)分布曲線光滑無(wú)震蕩,且當(dāng)前模型結(jié)果與前人研究結(jié)果[14-15,29]吻合良好,壓強(qiáng)極大值出現(xiàn)在前滯點(diǎn)(θ=0°),壓強(qiáng)極小值位于分離點(diǎn)附近,圓柱背后分離區(qū)壓強(qiáng)有所回升,以上結(jié)果驗(yàn)證了本文改進(jìn)的IB-LB 模型能夠?yàn)閺?fù)雜結(jié)構(gòu)外形上的拉格朗日節(jié)點(diǎn)提供光滑、可靠的壓強(qiáng)預(yù)測(cè).

圖4 圓柱表面壓強(qiáng)系數(shù)與方向角關(guān)系曲線Fig.4 Pressure coefficient curve of cylinder versus the orientation angle

圖5 為Re=20 和100 時(shí),圓柱附近的流線分布.當(dāng)Re＜ 47,圓柱背后出現(xiàn)穩(wěn)定的對(duì)稱漩渦;當(dāng)Re＞ 47,定常流逐漸發(fā)展為非定常流.數(shù)值結(jié)果表明,子迭代中引入的自由參數(shù)ε不會(huì)影響模型的數(shù)值精度,但其取值直接影響了壁面流線是否泄露.本文所有模擬取 ε=10-8求解的邊界力數(shù)值解能夠保證每個(gè)時(shí)間步中無(wú)滑移邊界條件準(zhǔn)確施加,精確滿足無(wú)滑移約束的速度場(chǎng)抑制了流線穿過(guò)圓柱體表面,如圖5 中所示.圖5(a)展示的定常流中,二維圓柱內(nèi)部包含兩對(duì)對(duì)稱的環(huán)形流線,圓柱背后分布著穩(wěn)定的對(duì)稱規(guī)則區(qū),規(guī)則區(qū)長(zhǎng)度列于表1.圖5(b)展示的非定常流中,可觀察到柱面背后波動(dòng)的漩渦,圓柱體內(nèi)部的環(huán)形流線隨表面脫落的渦旋同步波動(dòng).當(dāng)前數(shù)值結(jié)果驗(yàn)證了本文改進(jìn)的IB-LB 模型可在靈活控制迭代次數(shù)的基礎(chǔ)上保證無(wú)滑移條件準(zhǔn)確施加.

圖5 圓柱周圍流線分布Fig.5 Streamline distribution near the cylinder.

表2 列出了Re=100 時(shí),本文模擬獲得的圓柱表面時(shí)均阻力系數(shù)、升力系數(shù)幅值,及斯特勞哈爾數(shù)S t=Dfvs/U∞數(shù)據(jù)與現(xiàn)有文獻(xiàn)結(jié)果[15,17,30-31]的比較.fvs為圓柱表面渦脫頻率,通過(guò)側(cè)向力系數(shù)CL的快速傅里葉變換獲得.對(duì)比結(jié)果表明,本文模型獲得的數(shù)值解與文獻(xiàn)中實(shí)驗(yàn)、數(shù)值計(jì)算模型的結(jié)果吻合良好,驗(yàn)證了本文改進(jìn)的IB-LB 模型針對(duì)復(fù)雜邊界外形繞流的非定常流動(dòng)演化問題預(yù)測(cè)結(jié)果的有效性.

表2 Re=100 的當(dāng)前模型與其他方法的流動(dòng)特征參數(shù)比較Table 2 Comparison between the present model and other methods at Re=100

3.3 流動(dòng)經(jīng)過(guò)強(qiáng)制振動(dòng)圓柱

在上述Re=40 的定常流動(dòng)中,強(qiáng)制二維圓柱執(zhí)行周期性的側(cè)向振動(dòng)擾動(dòng)周圍流體,圓柱體的側(cè)向振蕩由如下正弦函數(shù)控制

式中,f為振蕩頻率,圓柱振蕩特征數(shù)St=Df/U∞=0.15.在圓柱的受迫振動(dòng)擾動(dòng)下,Re=40 的定常流動(dòng)逐漸發(fā)展為周期性的非定常流動(dòng),渦旋脫落頻率為fvs=0.915f.

圖6 顯示了振蕩圓柱經(jīng)歷的非定常阻力系數(shù)隨圓柱中心位置變化的相圖.采用傳統(tǒng)多段直接力格式的IB-LB 模型[17]執(zhí)行當(dāng)前模擬,獲得的結(jié)果被展示于圖6(a),本文改進(jìn)模型的模擬結(jié)果繪制于圖6(b).如圖6(a)所示,傳統(tǒng)直接力格式的IB-LB 模型的預(yù)測(cè)結(jié)果中,隨圓柱振動(dòng)作用于圓柱表面的流體力合力中呈現(xiàn)偽物理震蕩.相反,本文改進(jìn)模型獲得的非定常流體力隨圓柱振動(dòng)平滑波動(dòng),說(shuō)明滿足力同一性條件的IB-LB 模型能夠有效抑制模擬結(jié)果中非定常流體力的偽物理震蕩.此外,本文改進(jìn)的IB-LB 模型具備預(yù)測(cè)動(dòng)邊界誘導(dǎo)非定常流動(dòng)演化的能力.

圖6 振蕩圓柱阻力系數(shù)與圓柱中心位移的關(guān)系Fig.6 Drag coefficient acting on the oscillation cylinder versus the position of the cylinder center

我們比較了速度修正格式IB-LB 模型[15]和本文改進(jìn)的IB-LB 模型在每個(gè)時(shí)間步推進(jìn)中花費(fèi)的平均CPU 占用時(shí)長(zhǎng)tib-lb,以說(shuō)明模型的計(jì)算效率.這里,tib-lb是流動(dòng)求解器一個(gè)完整的時(shí)間步推進(jìn)所占用的CPU 時(shí)長(zhǎng),其中包含了IB 模型和LB 模型的CPU 占用.由于背景LB 是相同的,tib-lb的比較結(jié)果同樣可以定量評(píng)估兩種IB 模型的效率.兩組模擬在相同的工作站上執(zhí)行,計(jì)算環(huán)境參數(shù)如下,CPU:Intel Xeon(R) Silver 4214 CPU @ 2.2 GHz 2.19 GHz;內(nèi)存:DDR4 2400ZHz 16 GB × 4.針對(duì)Re=40,St=0.15 的振動(dòng)圓柱繞流算例,速度修正IB-LB 模型的tib-lb為0.126 s,本文改進(jìn)模型的tib-lb為0.089 s.當(dāng)前結(jié)果驗(yàn)證了本文改進(jìn)IB-LB 模型的計(jì)算效率優(yōu)勢(shì).

3.4 波動(dòng)翼型繞流

本節(jié)模擬可作為上述振蕩圓柱誘導(dǎo)周圍流體演化的非定常流驗(yàn)證算例的拓展,區(qū)別在于波動(dòng)翼代表了剛-柔運(yùn)動(dòng)疊加的可變形復(fù)雜運(yùn)動(dòng)邊界,此外在翼型前緣和尾緣附近帶有較大的曲率半徑,可驗(yàn)證本文模型對(duì)此類問題的處理能力.均勻流中波動(dòng)的NACA0012 翼型表面的非定常流體力和流場(chǎng)渦量云圖闡明了當(dāng)前模型在大變形結(jié)構(gòu)流固耦合動(dòng)力學(xué)問題中的預(yù)測(cè)結(jié)果可靠性,算法的實(shí)用性及可拓展性.

由平移、旋轉(zhuǎn)和波動(dòng)行波方程的疊加強(qiáng)制翼型振蕩,各行波分量由相同的運(yùn)動(dòng)頻率f控制,且各子運(yùn)動(dòng)間的相位差恒定,翼型波動(dòng)的特征數(shù)S t=L f/U∞=0.54,L為翼型弦長(zhǎng),L/h=100.

側(cè)向平移行波方程如下

式中,側(cè)向平移幅值A(chǔ)1=0.15L.

旋轉(zhuǎn)運(yùn)動(dòng)以距離翼型前緣點(diǎn)0.3L的位置為轉(zhuǎn)動(dòng)中心,旋轉(zhuǎn)運(yùn)動(dòng)行波方程為

式中,旋轉(zhuǎn)幅值θ2=7.75°.

波動(dòng)運(yùn)動(dòng)自距離翼型前緣點(diǎn)0.4L的位置起執(zhí)行,側(cè)向波動(dòng)幅值A(chǔ)3(x/L)為x坐標(biāo)的二次函數(shù)[32],波動(dòng)行波方程為

式中A3(x/L)=0.01L(x/L-0.4)+0.51L(x/L-0.4)2.

經(jīng)過(guò)10 個(gè)翼型運(yùn)動(dòng)周期,流場(chǎng)表現(xiàn)出穩(wěn)定的周期性,其周期特性通過(guò)時(shí)變的流體力和功率參數(shù)描述.圖7 展示了作用于波動(dòng)翼型表面的瞬時(shí)推進(jìn)力系數(shù)隨時(shí)間變化曲線.在波動(dòng)特征數(shù)St=0.54 條件下,周圍流體不能為波動(dòng)翼提供推進(jìn),整個(gè)運(yùn)動(dòng)周期中翼型始終承受流體施加的阻力作用(CT＜ 0).瞬時(shí)變化的復(fù)雜外形承受的流體力合力隨時(shí)間變化趨勢(shì)類似于側(cè)向平移行波曲線的波動(dòng)趨勢(shì),曲線光滑無(wú)震蕩,說(shuō)明當(dāng)前模型能夠準(zhǔn)確預(yù)測(cè)作用在復(fù)雜外形動(dòng)邊界上的非定常流體力,其數(shù)值解中沒有出現(xiàn)偽物理特征.

圖7 波動(dòng)翼型的時(shí)間相關(guān)推力系數(shù)Fig.7 Time-dependent thrust coefficient acting on the undulatory foil

圖8 展示了維持翼型復(fù)合運(yùn)動(dòng)的運(yùn)動(dòng)功率系數(shù)CPS,克服運(yùn)動(dòng)阻力的功率系數(shù)CPD,以及維持運(yùn)動(dòng)的總功率系數(shù)CPT隨時(shí)間變化曲線.其中,運(yùn)動(dòng)功率系數(shù)CPS根據(jù)如下公式計(jì)算

圖8 波動(dòng)翼型的時(shí)間相關(guān)運(yùn)動(dòng)功率系數(shù)、阻力功率系數(shù)和總功率系數(shù)Fig.8 Time-dependent kinematic power coefficient,overcome drag power coefficient,and total power coefficient required by the undulatory foil

式中,H(Xk,t)為流體施加在第k個(gè)拉格朗日節(jié)點(diǎn)的局部流體力,UB(Xk,t)為在第k個(gè)拉格朗日節(jié)點(diǎn)的運(yùn)動(dòng)速度,由復(fù)合運(yùn)動(dòng)行波方程的時(shí)間導(dǎo)數(shù)計(jì)算.克服流體阻力的功率系數(shù),總功率系數(shù)CPT為CPS和CPD的和,CPT=CPS+CPD.圖8所示的數(shù)值結(jié)果表明,維持波動(dòng)翼型運(yùn)動(dòng)所需的CPS,CPD,CPT系數(shù)具備類似于推進(jìn)力系數(shù)CT的周期屬性.與運(yùn)動(dòng)功率系數(shù)CPS相比,抵抗流體阻力所需的功率系數(shù)CPD在總功率系數(shù)CPT中的占比很小.CPT的波動(dòng)趨勢(shì)與CPS的波動(dòng)趨勢(shì)一致,當(dāng)CPT與CPS同時(shí)出現(xiàn)最大值時(shí),CPD為最小值,說(shuō)明波動(dòng)翼通過(guò)運(yùn)動(dòng)功率消耗完成了阻力減.這一結(jié)論與魚游模擬的數(shù)值實(shí)驗(yàn)結(jié)果[32,33]一致.

圖9 展示了波動(dòng)翼型周圍的渦量云圖,根據(jù)一個(gè)運(yùn)動(dòng)周期內(nèi)流場(chǎng)的瞬時(shí)渦量圖可知,波動(dòng)翼型周圍產(chǎn)生了周期性變化的剪切層.剪切層沒有在翼型后緣附近發(fā)展流動(dòng)分離,而是在翼型背后大致1.5 倍弦長(zhǎng)的位置形成核心渦脫落進(jìn)入尾跡.交替的符號(hào)相反的漩渦在波動(dòng)翼尾跡中形成經(jīng)典的卡門渦街.已有研究結(jié)果表明,單一的側(cè)向平移運(yùn)動(dòng)會(huì)在前緣兩側(cè)生成交替的逆向旋轉(zhuǎn)的漩渦,這一結(jié)論與圖9中的渦量云圖結(jié)果相悖.一個(gè)合理的解釋是,柔性波動(dòng)運(yùn)動(dòng)抑制了前緣附近剛性撲翼運(yùn)動(dòng)(平移+旋轉(zhuǎn))誘導(dǎo)的漩渦發(fā)展,在翼型兩側(cè)的剪切層中流體不穩(wěn)定性被削減,不足以支撐交替漩渦的形成.而在波動(dòng)翼的背后,符號(hào)相反的剪切層間相互作用增加了彼此攜帶的流體不穩(wěn)定,從而在波動(dòng)翼型背后形成卡門渦街.

圖9 均勻流動(dòng)中復(fù)合運(yùn)動(dòng)翼型周圍渦量云圖Fig.9 Vortex contour near the swimming foil in a uniform flow

4 結(jié)論

本文針對(duì)可變形動(dòng)邊界繞流的流固耦合動(dòng)力學(xué)問題,在隱式擴(kuò)散界面背景下,提出一種改進(jìn)的直接力格式IB-LB 模型.保留IB 方法在宏觀流場(chǎng)中通過(guò)向邊界附近流體施加體積力干擾速度場(chǎng)分布的思想,浸沒邊界對(duì)流場(chǎng)施加的體積力由使無(wú)滑移速度約束準(zhǔn)確施加的邊界力決定.借助狄拉克函數(shù)滿足的矩條件,從LB 模型的體積力表達(dá)式出發(fā),推導(dǎo)滿足力同一性的邊界力表達(dá)式.LB 背景下,力同一性的數(shù)學(xué)內(nèi)容是,拉格朗日節(jié)點(diǎn)向周圍歐拉節(jié)點(diǎn)的速度修正量擴(kuò)散是歐拉節(jié)點(diǎn)向拉格朗日節(jié)點(diǎn)的質(zhì)量密度插值的伴隨.滿足力同一性的IB 模型不僅完整保留了LB 模型的原有二階精度,而且有效抑制了非定常流體力的偽物理震蕩.另一方面,根據(jù)離散的狄拉克函數(shù)建立描述插值和擴(kuò)散操作的轉(zhuǎn)換矩陣,重疊的影響域中同類型變量被耦合,保證每個(gè)邊界節(jié)點(diǎn)處的無(wú)滑移速度約束條件準(zhǔn)確施加.表征無(wú)滑移條件的線性方程組采用迭代方法求解,無(wú)需矩陣求逆,降低計(jì)算負(fù)擔(dān),算法穩(wěn)定性不受拉格朗日節(jié)點(diǎn)分布約束.