王 卓 杜 林 成 龍 孫曉峰

1 北京航空航天大學(xué)航空發(fā)動(dòng)機(jī)研究院,北京 100191

2 華為技術(shù)有限公司,編譯器與編程語(yǔ)言實(shí)驗(yàn)室,北京 100094

3 北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院,北京 100191

1 概述

計(jì)算流體力學(xué)(computational fluid dynamics,CFD)通過(guò)離散流動(dòng)控制方程并采用計(jì)算機(jī)對(duì)其進(jìn)行迭代求解來(lái)獲得流動(dòng)的時(shí)空特征.流動(dòng)控制方程通常包含多變量且具有非線性,同時(shí)受各種復(fù)雜邊界條件和初始條件的影響,這導(dǎo)致基于解析推導(dǎo)的傳統(tǒng)理論分析手段難以應(yīng)用,此時(shí)CFD 幾乎是對(duì)流動(dòng)進(jìn)行理論分析及預(yù)測(cè)唯一可用的工具(Goldstine 1980).經(jīng)過(guò)數(shù)十年的發(fā)展,CFD 目前廣泛應(yīng)用于工業(yè)流體領(lǐng)域的分析與設(shè)計(jì)等環(huán)節(jié),這也使得實(shí)驗(yàn)環(huán)節(jié)大幅減少,設(shè)計(jì)成本得到了有效地降低.

目前絕大多數(shù)CFD 方法都是基于結(jié)構(gòu)/非結(jié)構(gòu)化的貼體網(wǎng)格展開(kāi),雖然相關(guān)的網(wǎng)格生成技術(shù)在近幾十年來(lái)取得了長(zhǎng)足的進(jìn)步(閻超 等 2011),但這一網(wǎng)格形式能夠處理的復(fù)雜邊界問(wèn)題極其有限,例如對(duì)于圖1 所示的心臟內(nèi)的血液流動(dòng)現(xiàn)象,心臟本身具有復(fù)雜的幾何外形,且外形隨時(shí)間呈現(xiàn)出周期性的改變,而采用貼體網(wǎng)格對(duì)這一現(xiàn)象進(jìn)行研究則會(huì)面臨網(wǎng)格生成及重構(gòu)困難的挑戰(zhàn).實(shí)際上心臟血液流動(dòng)現(xiàn)象是復(fù)雜邊界以及運(yùn)動(dòng)邊界兩類問(wèn)題結(jié)合的典型代表,而如何采用CFD 方法對(duì)這兩類問(wèn)題進(jìn)行準(zhǔn)確高效的求解也是進(jìn)一步拓展CFD 方法應(yīng)用的關(guān)鍵之一.

圖1 用于研究心臟血液流動(dòng)的模型示意圖

浸入式邊界方法(immersed bundary method,IB method)是一種非常適用于復(fù)雜邊界及運(yùn)動(dòng)邊界的CFD 方法,由Peskin (1972,1977)最早提出并且應(yīng)用于圖1 所示的心臟血液流動(dòng)模擬.IB 方法可以通過(guò)在正交的笛卡爾網(wǎng)格上對(duì)邊界進(jìn)行建模從而避免對(duì)于復(fù)雜邊界生成貼體網(wǎng)格,在Peskin 所提出的IB 方法中,固壁邊界對(duì)流動(dòng)的影響在控制方程中引入相應(yīng)的力源項(xiàng)來(lái)進(jìn)行刻畫(huà),如圖2 所示.在邊界位置發(fā)生改變時(shí),只需要對(duì)相應(yīng)的力源項(xiàng)進(jìn)行修改即可,因此能夠避免對(duì)網(wǎng)格進(jìn)行重構(gòu).通過(guò)力源項(xiàng)來(lái)刻畫(huà)邊界的影響早在20 世紀(jì)60 年代就被Sirovich (1967)從理論上證明其數(shù)學(xué)上的理論可行性,因此這一方法在隨后的CFD 發(fā)展中獲得了巨大的進(jìn)步,演化出各種不同的版本形式,并且還與基于雷諾平均的湍流數(shù)值模擬(RANS)、大渦模擬(LES)以及直接數(shù)值模擬(DNS)等各類不同精度的數(shù)值模擬技術(shù)相結(jié)合,被廣泛應(yīng)用于包含復(fù)雜邊界及運(yùn)動(dòng)邊界的各類流動(dòng)問(wèn)題的研究(Mittal &Iaccarino 2005,Verzicco 2023).

盡管IB 方法在處理復(fù)雜邊界以及運(yùn)動(dòng)邊界問(wèn)題上具有強(qiáng)大的能力,但是其主要缺陷在于處理高雷諾數(shù)流動(dòng)問(wèn)題時(shí)龐大的計(jì)算量.由于湍流邊界層厚度相較于層流而言非常小,通常需要在壁面處布置高分辨率網(wǎng)格來(lái)分辨湍流邊界層的線性底層,例如對(duì)于Spalart-Allmaras (SA)一方程湍流模型而言(Spalart &Allmaras 1992),通常要求壁面第一層網(wǎng)格的y+小于1.而對(duì)于IB 方法來(lái)說(shuō),由于對(duì)笛卡爾的網(wǎng)格加密必須沿空間三個(gè)方向同步展開(kāi),粗略估計(jì)表明,三維情形下網(wǎng)格總量與流動(dòng)雷諾數(shù)的1.5 次方呈正比,而對(duì)貼體網(wǎng)格而言僅與雷諾數(shù)的0.5 次方呈正比,這也導(dǎo)致對(duì)于高雷諾數(shù)流動(dòng),IB 方法所需的網(wǎng)格量遠(yuǎn)遠(yuǎn)超出貼體網(wǎng)格.如何在現(xiàn)有的硬件計(jì)算能力基礎(chǔ)上拓展IB 方法在高雷諾數(shù)流動(dòng)問(wèn)題中的應(yīng)用是進(jìn)一步擴(kuò)大該方法的應(yīng)用范圍以及實(shí)現(xiàn)工業(yè)規(guī)?；瘧?yīng)用的核心挑戰(zhàn).

本文將首先對(duì)IB 方法的基本思想進(jìn)行介紹,主要關(guān)注其中刻畫(huà)邊界影響的力源模型;其次將介紹近年來(lái)IB 方法在不同流動(dòng)問(wèn)題中的應(yīng)用,自然界中的生物體繞流通常包含復(fù)雜的固壁邊界,對(duì)其中流動(dòng)機(jī)理的認(rèn)知將有助于工程領(lǐng)域中設(shè)計(jì)更高性能的飛行器、推進(jìn)器等裝置;而流固耦合問(wèn)題則是工程領(lǐng)域典型的一類運(yùn)動(dòng)邊界問(wèn)題,由于涉及固體結(jié)構(gòu)的安全而受到密切關(guān)注.因此本文將以上述問(wèn)題作為復(fù)雜邊界以及運(yùn)動(dòng)邊界問(wèn)題的代表來(lái)具體介紹IB 方法在相關(guān)領(lǐng)域的應(yīng)用進(jìn)展.此外,本文還將介紹IB 方法在運(yùn)動(dòng)邊界發(fā)聲問(wèn)題中的研究進(jìn)展,通常認(rèn)為IB 方法在邊界精度處只有一階精度,這似乎與計(jì)算氣動(dòng)聲學(xué)中的高精度要求相矛盾,但已有的研究進(jìn)展表明邊界處的計(jì)算精度并不會(huì)對(duì)IB 方法在高精度模擬中的應(yīng)用產(chǎn)生根本性的阻礙.最后,本文將對(duì)IB 方法在高雷諾數(shù)流動(dòng)模擬中的研究進(jìn)展進(jìn)行介紹.

2 基本思想

當(dāng)采用與壁面不貼合的笛卡爾網(wǎng)格進(jìn)行流動(dòng)模擬時(shí),一個(gè)關(guān)鍵問(wèn)題是如何刻畫(huà)壁面邊界對(duì)于流動(dòng)的影響,而IB 方法所采取的方式為在控制方程中引入描述邊界影響的力源項(xiàng)從而將壁面內(nèi)部拓展為與外部流場(chǎng)一體化的“偽流場(chǎng)域”.方程(1)給出了適用于不可壓縮流動(dòng)的動(dòng)量方程,其中u為流動(dòng)速度,t為時(shí)間,p和ρ分別為壓力和密度,ν為運(yùn)動(dòng)黏度,f(x,t) 即為刻畫(huà)邊界影響的力源項(xiàng),并且會(huì)隨時(shí)間和空間呈現(xiàn)動(dòng)態(tài)變化.

根據(jù)力源建模方式,IB 方法可以進(jìn)一步細(xì)分為連續(xù)力源法和離散力源法.

2.1 連續(xù)力源法

記xi為壁面拉格朗日坐標(biāo)系下某一離散邊界單元的坐標(biāo),Fi為該單元處壁面對(duì)流體的作用力,由于方程(1)中的流動(dòng)控制方程建立于歐拉坐標(biāo)系,因此f與Fi之間的關(guān)系可以寫(xiě)為方程(2)的形式,其中x為歐拉坐標(biāo)系下的坐標(biāo),δ為Dirac Delta 函數(shù).由于δ函數(shù)具有奇異性,為提高數(shù)值穩(wěn)定性,在實(shí)際處理中通常會(huì)采用具有空間連續(xù)性以及離散守恒性的正則化函數(shù)分布來(lái)對(duì)其進(jìn)行代替從而完成對(duì)力源的坐標(biāo)轉(zhuǎn)換(Peskin 1972,Uhlmann 2005),而這一處理方式也導(dǎo)致連續(xù)力源法通常在邊界處僅有一階精度.

連續(xù)力源法中的另一個(gè)關(guān)鍵點(diǎn)在于如何獲得壁面與流體之間的作用力Fi.早期Peskin 等人所研究的復(fù)雜壁面對(duì)象如心臟等一般為彈性體(Peskin 1977,Saiki &Biringen 1996,Uhlmann 2005),對(duì)于此種類型的彈性壁面邊界可以通過(guò)其物質(zhì)的本構(gòu)方程來(lái)建立壁面內(nèi)部的彈性力分布,而在壁面處則通過(guò)內(nèi)外部作用力的平衡來(lái)間接獲得壁面對(duì)流體的作用力.這一考慮壁面本構(gòu)方程的力源方法并不適用于壁面為剛體的情形,因此為了模擬剛體繞流問(wèn)題,Goldstein 等(1993)通過(guò)雙模態(tài)控制器原理構(gòu)造了反饋力源方法,其表達(dá)式如方程(3)所示,其中ub為壁面處的邊界條件.

反饋力源方法實(shí)際上是利用線性理論控制流動(dòng)方程滿足壁面邊界條件,線性理論能夠很好地達(dá)到控制效果的原因在于力源作用區(qū)域或者控制區(qū)域是一個(gè)很小的區(qū)域,可以近似看成線性區(qū)域.Goldstein 等(1993)將該方法用于不可壓縮槽道湍流模擬以及圓柱繞流問(wèn)題,證明了該方法對(duì)于非定常流動(dòng)模擬的精確性,Zhong 和Sun (2009)、Du 等(2016a,2016b)也利用反饋力源法求解不可壓縮翼型和葉柵流動(dòng),取得了不錯(cuò)的預(yù)測(cè)結(jié)果.反饋力源法的形式極為簡(jiǎn)單,且為顯式力源,可以直接施加到動(dòng)量方程右端求解.然而反饋力源法中的力源形式涉及兩個(gè)經(jīng)驗(yàn)參數(shù)(α、β),這兩個(gè)經(jīng)驗(yàn)參數(shù)用來(lái)調(diào)節(jié)反饋力源的響應(yīng)速度和對(duì)壁面邊界條件的滿足程度,其取值與具體非定常流動(dòng)有關(guān).對(duì)于定常流動(dòng)和不同特征頻率的非定常流動(dòng)而言,α、β一般有所不同,但是當(dāng)其取值過(guò)小時(shí),邊界處的流動(dòng)計(jì)算精度會(huì)有所降低,且收斂較慢,但取值過(guò)大又會(huì)增加系統(tǒng)剛性,引起數(shù)值不穩(wěn)定.此外,反饋力源法允許的 CFL 數(shù)過(guò)小 [O(10-3～10-2)] (Goldstein et al.1993,Lee 2003,Shin et al.2008),為此雖然其實(shí)現(xiàn)簡(jiǎn)單,效率極高,但過(guò)小的 CFL 數(shù)仍然使得其在具體非定常流動(dòng)計(jì)算時(shí)允許的時(shí)間步長(zhǎng)過(guò)小.為了進(jìn)一步放寬CFL 數(shù)的限制,Taira 和Colonius (2007),Li 等(2016)以及Wang 等(2020)提出了一種對(duì)力源進(jìn)行隱式求解的投影式IB方法以實(shí)現(xiàn)對(duì)不可壓縮流動(dòng)的求解.這一方法的核心思想是將力源視作類似于壓力項(xiàng)的拉格朗日乘數(shù)(Lagrange Multiplier),此時(shí)力源項(xiàng)可采用與壓力項(xiàng)求解同樣的數(shù)值思想來(lái)獲得,從而避免了早期連續(xù)力源法中需要人為給定反饋系數(shù)的缺陷,同時(shí)實(shí)現(xiàn)放寬CFL 數(shù)限制的目標(biāo).

從Goldstein 等(1993)構(gòu)造的反饋力源形式來(lái)看,其流體計(jì)算域擴(kuò)展到了固體內(nèi)部區(qū)域,從而形成了流體計(jì)算域和拓展區(qū),拓展區(qū)內(nèi)的流體仍然滿足流動(dòng)控制方程,而力源的強(qiáng)度則同時(shí)取決于內(nèi)外流體的瞬時(shí)流動(dòng),這種思想的可行性實(shí)際上早在 IB 提出之前便已經(jīng)被Sirovich (1967)證明.該研究推導(dǎo)了可壓縮NS 方程的間斷形式,此時(shí)邊界條件完全有間斷形式的源項(xiàng)取代,源項(xiàng)的強(qiáng)度則取決于間斷強(qiáng)度,即邊界內(nèi)外的應(yīng)力或熱傳導(dǎo)等.由此可以看出Sirovich 推導(dǎo)的間斷NS 方程以及利用源項(xiàng)代替固體邊界條件實(shí)際上為Peskin 發(fā)展 IB 方法提供了基本理論依據(jù),同時(shí)也說(shuō)明了 IB 方法通過(guò)延拓流體域處理剛性固體域的理論可靠性.連續(xù)力法通常都只是在浸入式邊界點(diǎn)上施加力源,因此物體內(nèi)部是允許有流動(dòng)的.從物理上來(lái)說(shuō),物體內(nèi)部是固體,是不可能存在流動(dòng)現(xiàn)象的.但是,對(duì)于浸入式邊界方法來(lái)說(shuō),整個(gè)域包括固體和流體是一體離散的,控制方程在兩個(gè)域上沒(méi)有本質(zhì)區(qū)別,將兩個(gè)域隔開(kāi)的界面已經(jīng)不存在了,起到作用的是“力源”.因此,在算法上是允許內(nèi)部的虛擬流動(dòng)的,而且內(nèi)部流動(dòng)可以作為外部真實(shí)流動(dòng)的“光順器”(Goldstein et al.1993).

另外兩種典型的連續(xù)力源法為虛擬彈簧力源法(Lai &Peskin 2000)和罰函數(shù)方法(Kim &Peskin 2007,Huang et al.2011).虛擬彈簧力源法中所構(gòu)造的力源與流體偏離給定位置的程度成正比,類似于一個(gè)具有較大彈性系數(shù)的彈簧作用于流體上,當(dāng)流體稍微偏離預(yù)期位置后便會(huì)產(chǎn)生極大的回復(fù)力將流體“拉”回到固體邊界從而使其滿足無(wú)滑移邊界條件,如圖3(a)所示.在此基礎(chǔ)上Kim 和Peskin 進(jìn)一步提出罰函數(shù)方法,該方法結(jié)合虛擬彈簧力和Peskin 早期利用本構(gòu)方程構(gòu)建邊界力源的方法,將力源分解成無(wú)質(zhì)量邊界的彈性力以及由于無(wú)質(zhì)量邊界隨流體運(yùn)動(dòng)后與有質(zhì)量邊界存在位置差異誘導(dǎo)的彈簧力.該方法的優(yōu)點(diǎn)在于不僅能夠考慮無(wú)質(zhì)量固體邊界的流固耦合問(wèn)題,也能夠考慮固體邊界任意密度分布的流固耦合問(wèn)題.罰函數(shù)方法中力源的給出方式如方程(4)所示,其中Da是達(dá)西系數(shù)(Darcy number),用以描述多孔介質(zhì)滲流特性.Da=K0/L2,K0是特征滲透性,L是特征長(zhǎng)度.從方程(4)的形式不難看出,這也是一種反饋控制系統(tǒng),相當(dāng)于方程(3)中取其中自由參數(shù)是K,如果K趨向于無(wú)窮,力源項(xiàng)就消失了,主控方程就變成了只作用于純粹的流體.如果K趨向于零,那么力源項(xiàng)在主控方程中就處于主導(dǎo)地位,其他的諸如對(duì)流、擴(kuò)散、詳壓力項(xiàng)的值都可以忽略不計(jì).對(duì)于多孔介質(zhì)來(lái)說(shuō),0＜K＜∞,力源的作用是在指定的區(qū)域模擬流體流過(guò)時(shí)的動(dòng)量損失,當(dāng)然,這種情況下主控方程演變?yōu)镹S/Brinkman 方程.通過(guò)調(diào)節(jié)K的取值,可以將固體壁面、流體、多孔介質(zhì)一起求解.同樣的,K的取值與α和β的問(wèn)題類似,造成待求解方程的剛性變大,同時(shí)受到數(shù)值穩(wěn)定性的限制.

圖3 (a) 固壁邊界虛擬彈簧力示意圖,(b) 銳利界面的浸入式邊界方法

2.2 離散力源法

連續(xù)力源法雖然形式簡(jiǎn)單,但是對(duì)于不同壁面邊界條件的刻畫(huà)能力比較受限,大部分方法僅僅只能實(shí)現(xiàn)簡(jiǎn)單的不可壓縮流動(dòng)下無(wú)滑移壁面邊界條件.而離散力源法由于從離散化的 NS 方程出發(fā),因此對(duì)于邊界條件的刻畫(huà)更為靈活,其中最具代表性的是由(Mohd Yusof 1997)所發(fā)展的方法: 將方程(2)寫(xiě)成如方程(5)所示的離散形式,其中Rn+1/2為方程(2)中除力源項(xiàng)之外的所有右端項(xiàng).

Lima E Silva 等(2003)形象地將其發(fā)展的離散力法叫做“physical virtual model (PVM)”,該離散力源法將力源拆分成四項(xiàng),分別叫做加速度力、慣性力、黏性力以及壓力.在離散力法中,力源是在方程離散后加入到方程中,方程求解穩(wěn)定性較好,適合于高雷諾數(shù)計(jì)算,并且可以根據(jù)壁面邊界條件的不同構(gòu)造不同的插值方法,由此適用的壁面邊界條件更加廣泛.Fadlun 等(2000)證明如方程(5)所示的離散力源法能夠在邊界處實(shí)現(xiàn)二階精度,并且可以與大渦模擬方法相結(jié)合,大大地拓寬了IB 方法的適用性.Mohd Yusof (1997)構(gòu)造的離散力源法需要對(duì)壁面位置進(jìn)行判斷和插值操作,判斷笛卡爾網(wǎng)格單元與壁面的內(nèi)外關(guān)系.為此Su 等(2007)利用離散力源法和預(yù)測(cè)-校正算法思想構(gòu)造出了既不顯含經(jīng)驗(yàn)參數(shù),也不需要專門(mén)判斷固體邊界內(nèi)外,同時(shí)理論上可以精確滿足壁面邊界條件的離散力源法,并通過(guò)靜止/運(yùn)動(dòng)圓柱繞流問(wèn)題證明了該方法的時(shí)間精確性,數(shù)值測(cè)試表明該方法相較于連續(xù)力源法能夠?qū)崿F(xiàn)較大的 CFL 數(shù) [O(10-1～1)],但需要求解一個(gè)與力源相關(guān)的線性方程組.基于預(yù)測(cè)-校正的思想,Wu 和Shu (2009)提出了一種與格子玻爾茲曼方法相結(jié)合的IB 方法,其中力源項(xiàng)通過(guò)隱式速度修正的方式來(lái)構(gòu)建,文中采用圓柱及翼型繞流等案例驗(yàn)證該方法的精度.Yang 等(2017)進(jìn)一步將該方法與Yang 等(2015,2016)所發(fā)展的多種通量計(jì)算格式相結(jié)合并應(yīng)用于多種不同類型的三維可壓縮/不可壓縮的流動(dòng)計(jì)算中.Wang 和Zhang (2011)發(fā)展了適用于離散流函數(shù)(discrete stream function)方法的IB 方法,并采用多個(gè)不可壓縮的靜止/運(yùn)動(dòng)邊界案例對(duì)算法進(jìn)行驗(yàn)證,Wang 和Zhang (2011)還結(jié)合了局部網(wǎng)格加密的技術(shù)來(lái)降低網(wǎng)格總量從而提升計(jì)算效率,這一網(wǎng)格處理方式也是目前采用IB方法模擬高雷諾數(shù)流動(dòng)常用網(wǎng)格技術(shù).

在離散力源法思想的指導(dǎo)下,實(shí)際上后來(lái)發(fā)展了無(wú)需力源的浸入式邊界方法,主要代表為混合笛卡爾/浸入式邊界方法(hybrid cartesian/immersed boundary method) (Gilmanov et al.2003,Gilmanov &Sotiropoulos 2005)和鬼點(diǎn)法(ghost-cell method)(Tseng &Ferziger 2003,Mittal et al.2008,Lee &You 2013,Zhang et al.2021)等.這類方法通常又被稱之為銳利界面的IB 方法.在該類方法中,通常需要為每個(gè)邊界外部的流場(chǎng)單元沿壁面外法向構(gòu)造外力點(diǎn)(forcing point),在鬼點(diǎn)法中又被稱之為鏡像點(diǎn)(image point),如圖3(b)所示,其中的P1即為網(wǎng)格點(diǎn)P0所對(duì)應(yīng)的外力點(diǎn),外力點(diǎn)處的流動(dòng)變量一般通過(guò)插值獲得,進(jìn)一步可在壁面附近構(gòu)造出滿足壁面邊界條件的流動(dòng)時(shí)空分布,并代入離散的流動(dòng)控制方程中進(jìn)行求解,以此刻畫(huà)壁面邊界對(duì)流場(chǎng)的影響.值得注意的是,Tseng 和Ferziger (2003)所提出的鬼點(diǎn)法仍然采用的是早期離散力源法的思想,即構(gòu)建鬼點(diǎn)的目的是計(jì)算力源分布,因此與其他類型的鬼點(diǎn)法有著顯著的差異.從上述過(guò)程來(lái)看,銳利界面的IB 方法擁有與貼體網(wǎng)格方法相似的邊界刻畫(huà)思想,但由于網(wǎng)格點(diǎn)并未與壁面完全貼合,這一類方法通常在壁面處并不能滿足守恒性,因此在處理運(yùn)動(dòng)邊界問(wèn)題時(shí),壁面處會(huì)出現(xiàn)不同程度的數(shù)值振蕩(Lee &You 2013,Al-Marouf &Samtaney 2017,Griffith &Leontini 2017).相較于連續(xù)力源法,離散力源法中需要判斷每個(gè)笛卡爾網(wǎng)格單元位于固體內(nèi)部還是外部,并且對(duì)于近壁面的網(wǎng)格單元還需要精確計(jì)算其到固體邊界的距離,該距離在運(yùn)動(dòng)邊界問(wèn)題的模擬中需要根據(jù)固體邊界位置實(shí)時(shí)更新,可能消耗額外的計(jì)算資源,同時(shí)該參數(shù)也是湍流模擬所必須知曉的一個(gè)網(wǎng)格參數(shù),而連續(xù)力源法的施加中并不會(huì)對(duì)這一距離進(jìn)行計(jì)算.因此從這一角度看,離散力源法可能比連續(xù)力源法更加適合用于湍流模擬,而目前已公開(kāi)發(fā)表的關(guān)于IB 方法在湍流模擬中應(yīng)用的文章,全部采用的是離散力源法(Capizzano 2011,2016,Tamaki et al.2017,Wang 等2023),下文章節(jié)將進(jìn)一步的討論IB 方法在湍流模擬中的應(yīng)用.

3 IB 方法在運(yùn)動(dòng)邊界問(wèn)題中的應(yīng)用

3.1 生物流動(dòng)

自然界中的魚(yú)類、鳥(niǎo)類等生物均依靠自身部位 (魚(yú)鰭和鳥(niǎo)的翅膀) 在流體介質(zhì)中的擺動(dòng)來(lái)獲得前進(jìn)的動(dòng)力(Tong et al.2021,Zhang et al.2022a),這一過(guò)程通常表現(xiàn)出較高的推進(jìn)效率以及較低的噪聲水平.因此人類始終希望能夠深入理解自然界生物流動(dòng)中的復(fù)雜機(jī)理,以此來(lái)設(shè)計(jì)具有更高性能的推進(jìn)裝置,CFD 方法是目前廣泛采用的一種研究手段.

以魚(yú)鰭為例,真實(shí)的魚(yú)鰭通常具有復(fù)雜的外形,單個(gè)魚(yú)鰭在擺動(dòng)過(guò)程中的推力狀態(tài)可能與魚(yú)的身體或上下游的其他鰭的擺動(dòng)相互耦合(Park &Sung 2018,Han et al.2020,Mendelson &Techet 2021).與此同時(shí),由于某些魚(yú)類經(jīng)常組成龐大的群體,這種群體行為同樣會(huì)對(duì)單個(gè)個(gè)體的推力狀態(tài)以及流場(chǎng)結(jié)構(gòu)產(chǎn)生重要的影響(Weihs 1973,Gazzola et al.2016,Bao et al.2017).這些基本特征使得傳統(tǒng)的貼體網(wǎng)格方法在處理這一類問(wèn)題上面臨網(wǎng)格生成困難、多體運(yùn)動(dòng)問(wèn)題難以處理的挑戰(zhàn),而這些挑戰(zhàn)正是IB 方法所擅長(zhǎng)之處,因此這一方法也在魚(yú)類游泳等生物流動(dòng)領(lǐng)域獲得了廣泛的應(yīng)用.

為降低預(yù)測(cè)難度同時(shí)把握關(guān)鍵物理因素,早期對(duì)于生物流動(dòng)的數(shù)值模擬研究多采用簡(jiǎn)化固體模型,如用振動(dòng)翼型來(lái)表示魚(yú)鰭或鳥(niǎo)翅(Wang et al.2014,2019,Deng &Caulfield 2015,Deng et al.2015,2016).而近年來(lái)的部分研究工作已經(jīng)開(kāi)始將生物體真實(shí)幾何外形以及擺動(dòng)規(guī)律包含于數(shù)值模擬中.圖4(a)給出了一種外形類似蝙蝠的魚(yú)類示意圖,對(duì)于蝙蝠魚(yú)而言,其推力產(chǎn)生機(jī)制難以通過(guò)簡(jiǎn)化的固體模型進(jìn)行闡述,因此Huang 等(2020)通過(guò)對(duì)真實(shí)蝙蝠魚(yú)身體擺動(dòng)規(guī)律的細(xì)致觀察建立了如圖4(b)所示的固體模型,其中魚(yú)身的擺動(dòng)通過(guò)指定固體邊界在不同時(shí)刻的空間分布來(lái)實(shí)現(xiàn).Zhang 等(2022a)采用IB 方法對(duì)圖4(b)中的蝙蝠魚(yú)模型的推力產(chǎn)生機(jī)制展開(kāi)了系統(tǒng)性的參數(shù)研究,其結(jié)果表明魚(yú)身沿弦向(或流向)的擺動(dòng)波長(zhǎng)對(duì)于推力的產(chǎn)生、推進(jìn)效率以及能量消耗等參數(shù)有著至關(guān)重要的影響,合適的弦向擺動(dòng)波長(zhǎng)有助于獲得更高的推進(jìn)性能,而擺動(dòng)頻率同樣能夠?qū)τ谇熬墱u、翼尖渦等渦系結(jié)構(gòu)產(chǎn)生明顯的影響進(jìn)而改變推進(jìn)狀態(tài).上述研究結(jié)果對(duì)于深入理解包含復(fù)雜三維變形魚(yú)類的推力產(chǎn)生機(jī)理具有重要意義,而這一問(wèn)題對(duì)于采用貼體網(wǎng)格的CFD 方法而言具有極大的挑戰(zhàn)性.

圖4 蝙蝠魚(yú)示意圖(Huang et al.2020,Zhang et al.2022b)

正如前文所提到,魚(yú)類在游泳的過(guò)程中多個(gè)鰭之間的流動(dòng)結(jié)構(gòu)會(huì)產(chǎn)生相互干涉進(jìn)而影響推力狀態(tài),貼體網(wǎng)格在處理多鰭干涉問(wèn)題時(shí)往往面臨不同程度的網(wǎng)格重構(gòu)困難,通常來(lái)講運(yùn)動(dòng)邊界的數(shù)量越多重構(gòu)過(guò)程越復(fù)雜,重構(gòu)過(guò)程消耗的計(jì)算資源也會(huì)越多.而對(duì)于IB 方法而言,由于采用固定的笛卡爾網(wǎng)格,運(yùn)動(dòng)邊界的數(shù)量并不會(huì)對(duì)數(shù)值模擬的過(guò)程產(chǎn)生根本性的影響,因此在處理此類多運(yùn)動(dòng)邊界問(wèn)題上IB 方法同樣具有明顯的優(yōu)勢(shì)(Pan et al.2016).金槍魚(yú)、鱒魚(yú)等是典型的依靠多鰭干涉產(chǎn)生推力的魚(yú)類,其中主要包含腹鰭、臀鰭以及尾鰭等之間的相互干涉,Zhang 等(2022b)采用IB 方法對(duì)金槍魚(yú)中腹鰭與尾鰭之間的相互干涉及推力產(chǎn)生機(jī)理進(jìn)行了數(shù)值模擬研究,結(jié)果表明在腹鰭尾跡的激勵(lì)下,尾鰭所產(chǎn)生的推力能夠出現(xiàn)明顯提高,但腹鰭的推力狀態(tài)基本不受尾鰭的影響,其主要作用機(jī)理在于尾鰭所產(chǎn)生的前緣脫落渦強(qiáng)度會(huì)在于腹鰭尾跡中脫落渦的融合過(guò)程中被強(qiáng)化,這種魚(yú)鰭間的非定常尾跡相互作用機(jī)制在魚(yú)鰭間距較大時(shí)會(huì)表現(xiàn)得更加穩(wěn)定.魚(yú)類群體游泳現(xiàn)象與單個(gè)魚(yú)的多鰭問(wèn)題是相似的,本質(zhì)上也屬于多運(yùn)動(dòng)邊界問(wèn)題,Peng等(2018)、Kelly 和Menzer (2023)均采用IB 方法對(duì)魚(yú)類群體游泳現(xiàn)象中的非定常渦相互作用機(jī)制及其對(duì)推力的影響進(jìn)行了深入研究,其中均使用多個(gè)運(yùn)動(dòng)邊界來(lái)對(duì)魚(yú)群進(jìn)行模擬,不同個(gè)體的尾跡之間相互干涉所形成的流動(dòng)特征如圖5 所示.上述研究也充分體現(xiàn)了IB 方法在處理運(yùn)動(dòng)邊界問(wèn)題上的優(yōu)勢(shì).由于本文主要關(guān)注數(shù)值計(jì)算方法而非物理機(jī)理,因此僅列舉部分IB 方法應(yīng)用的案例,詳細(xì)的物理機(jī)理讀者可參考Zhang 等(2022a)所撰寫(xiě)的文章.

圖5 魚(yú)類群體游泳數(shù)值模擬中不同個(gè)體的尾跡互相干涉(Peng et al.2018)

3.2 流固耦合

流固耦合現(xiàn)象是工程領(lǐng)域備受關(guān)注的問(wèn)題之一,因?yàn)楣腆w結(jié)構(gòu)的振動(dòng)可能會(huì)導(dǎo)致其產(chǎn)生嚴(yán)重的破壞失效進(jìn)而對(duì)生命財(cái)產(chǎn)安全帶來(lái)巨大的損失,一個(gè)典型案例是1940 年美國(guó)Tacoma 大橋在低速自然風(fēng)的激勵(lì)下出現(xiàn)劇烈振動(dòng)進(jìn)而導(dǎo)致橋梁整體坍塌.流固耦合問(wèn)題同樣是一類常見(jiàn)的運(yùn)動(dòng)邊界問(wèn)題,上節(jié)所提到的生物流動(dòng)問(wèn)題中,邊界的運(yùn)動(dòng)軌跡主要是通過(guò)外部輸入的運(yùn)動(dòng)規(guī)律進(jìn)行控制,而對(duì)于流固耦合問(wèn)題,固體邊界的運(yùn)動(dòng)則是由流場(chǎng)激勵(lì)以及固體結(jié)構(gòu)響應(yīng)共同決定的.因此當(dāng)采用貼體網(wǎng)格求解時(shí),由于固體運(yùn)動(dòng)規(guī)律未知,對(duì)于網(wǎng)格重構(gòu)算法的要求也會(huì)相對(duì)更高.對(duì)于具有復(fù)雜邊界或者包含多物體的流固耦合問(wèn)題,應(yīng)用IB 方法同樣具有可以通過(guò)采用正交的笛卡爾網(wǎng)格以及避免由于邊界運(yùn)動(dòng)所引起的網(wǎng)格重構(gòu)來(lái)降低計(jì)算模型的復(fù)雜程度.

圓柱、方柱等鈍體的渦致振動(dòng)(vortex-induced vibration,VIV)是自然界中最為基本的一類流固耦合問(wèn)題,因?yàn)檫@一類問(wèn)題通常具有相對(duì)簡(jiǎn)單的固體幾何,對(duì)應(yīng)的流動(dòng)結(jié)構(gòu)特征更加明顯,有助于準(zhǔn)確地把握流固耦合現(xiàn)象背后的主要物理規(guī)律.Griffith 和Leontini (2017)提出了一種銳利界面的IB 方法并成功將其應(yīng)用于VIV 問(wèn)題,Du 等(2014)、Du 和Sun (2015)將Gilmanov 等(2003)、Gilmanov 和Sotiropoulos (2005)所提出的一種銳利界面IB 方法進(jìn)一步拓展至可壓縮流動(dòng)并利用該方法研究了圓柱VIV 中非定常渦結(jié)構(gòu)的三維特征,同時(shí)證明圓柱旋轉(zhuǎn)對(duì)VIV 的抑制作用,該抑制作用被(Wong et al.2018)的實(shí)驗(yàn)研究所驗(yàn)證.Chen 等(2018)進(jìn)一步采用IB 方法研究了三個(gè)串列圓柱的流致振動(dòng)問(wèn)題,深入分析了圓柱之間距離對(duì)其流固耦合效應(yīng)的影響.Xie等(2019)采用一種基于罰函數(shù)的IB 方法對(duì)圓柱VIV 進(jìn)行了數(shù)值模擬研究,主要關(guān)注附著于圓柱后端的柔性細(xì)絲對(duì)圓柱VIV 響應(yīng)的影響,結(jié)果表明柔性細(xì)絲會(huì)加劇圓柱VIV 的共振幅值,并且會(huì)拓寬出現(xiàn)VIV 鎖定的頻率范圍.總體來(lái)看,目前IB 方法在VIV 等鈍體流固耦合問(wèn)題上的應(yīng)用已經(jīng)比較成熟,相關(guān)的研究極大促進(jìn)了研究人員對(duì)這一現(xiàn)象的物理理解,特別是包含多個(gè)鈍體的情形.

上文中提到橋梁因受自然風(fēng)而產(chǎn)生振動(dòng)變形是建筑領(lǐng)域中一個(gè)重要的流固耦合問(wèn)題,而橋梁本身結(jié)構(gòu)復(fù)雜,包含主梁、欄桿以及懸索等眾多組成部分,每個(gè)部分均可能在風(fēng)中產(chǎn)生非定常脫落渦,進(jìn)而對(duì)橋梁造成激勵(lì)導(dǎo)致其出現(xiàn)流固耦合振動(dòng).Wang 和Cao (2022)采用IB 方法對(duì)主梁上的多條欄桿進(jìn)行了建模,并將其與完全采用貼體非結(jié)構(gòu)網(wǎng)格的方式進(jìn)行了比對(duì),結(jié)果顯示在主梁的基礎(chǔ)上采用IB 方法來(lái)刻畫(huà)側(cè)部欄桿能夠大幅降低網(wǎng)格的復(fù)雜程度.Zhao 等(2020)采用IB 方法對(duì)橋梁甲板與橋墩的組合模型進(jìn)行了數(shù)值模擬,主要關(guān)注在颶風(fēng)、海嘯等極端條件下橋梁的流固耦合響應(yīng),數(shù)值模擬準(zhǔn)確地預(yù)測(cè)實(shí)驗(yàn)結(jié)果,并且表明水面波動(dòng)所產(chǎn)生垂直于橋梁甲板方向的作用力是相鄰橋墩之間甲板上激勵(lì)的主要來(lái)源.上述研究結(jié)果共同表明,對(duì)于橋梁這一類包含多個(gè)不同組件的系統(tǒng)而言,IB 方法能夠有效地降低預(yù)測(cè)其流固耦合響應(yīng)的難度.

航空領(lǐng)域的葉輪機(jī)械葉片顫振是一類典型的包含多運(yùn)動(dòng)邊界的流固耦合問(wèn)題,而目前在對(duì)這一現(xiàn)象進(jìn)行預(yù)測(cè)的時(shí)候,多采用能量法對(duì)問(wèn)題進(jìn)行簡(jiǎn)化,即假定葉片以某一固有特征頻率做小幅振動(dòng),通過(guò)一個(gè)振動(dòng)周期內(nèi)流體對(duì)葉片做的功來(lái)判斷葉片是否穩(wěn)定,即振幅是否會(huì)被放大.能量法的假設(shè)并不能夠準(zhǔn)確描述這一真實(shí)的流固耦合過(guò)程,因此對(duì)于同一問(wèn)題,不同的預(yù)測(cè)代碼之間往往表現(xiàn)出較大的誤差(Holzinger et al.2016),但是對(duì)包含多葉片的流固耦合問(wèn)題進(jìn)行模擬往往又因?yàn)榫W(wǎng)格重構(gòu)等過(guò)程需要消耗巨大的計(jì)算量.為此,Zhong 和Sun (2009)將IB 方法應(yīng)用于葉輪機(jī)械葉片顫振問(wèn)題的模擬預(yù)測(cè)中,大大降低了考慮多葉片流固耦合計(jì)算模型的網(wǎng)格復(fù)雜程度,同時(shí)能夠準(zhǔn)確地捕捉到由于非線性作用所引起的振動(dòng)極限環(huán).Zhong 和Sun (2009)的工作還是在層流條件下開(kāi)展了,而實(shí)際葉輪機(jī)械流動(dòng)通常具有較高的雷諾數(shù),因此Du 等(2016a,2016b)進(jìn)一步將IB 方法拓展至高雷諾數(shù)并對(duì)包含多排葉片的情形展開(kāi)了數(shù)值模擬,深入闡述了由多葉片排干涉所引起的非定常效應(yīng)對(duì)葉片負(fù)荷的影響,其中計(jì)算模型及所得到的流場(chǎng)云圖分布如圖6 所示.Du 等人的研究表明,采用IB 方法可以克服傳統(tǒng)貼體網(wǎng)格方法在處理小軸向間距葉片排問(wèn)題上所面臨的網(wǎng)格生成困難等難題.上述研究工作對(duì)高雷諾數(shù)流動(dòng)的模擬以及IB 方法的工程應(yīng)用進(jìn)行了初步的嘗試,其結(jié)果表明對(duì)于工程領(lǐng)域的高雷諾數(shù)流動(dòng),由于需要布置大量的網(wǎng)格對(duì)流動(dòng)邊界層進(jìn)行分辨,在嚴(yán)格滿足湍流模型要求的網(wǎng)格分辨率下,IB 方法所需的網(wǎng)格量通常會(huì)遠(yuǎn)超傳統(tǒng)貼體網(wǎng)格,因此有必要對(duì)方法進(jìn)行進(jìn)一步改進(jìn)以實(shí)現(xiàn)更加高效的預(yù)測(cè).

圖6 基于IB 方法的葉輪機(jī)械轉(zhuǎn)靜干涉計(jì)算.(a) 轉(zhuǎn)靜葉片排模型,(b) 流場(chǎng)瞬時(shí)渦量分布(Du et al,2016a,2016b)

除了上述工程問(wèn)題的應(yīng)用研究以外,自然界還存在許多其他的流固耦合現(xiàn)象,IB 方法也在其機(jī)理研究中扮演著重要的角色.例如Huang 等(2007)提出了一種適用于細(xì)絲擺動(dòng)[圖7(a)]的IB 方法,通過(guò)構(gòu)建力源的方式來(lái)描述能夠產(chǎn)生柔性變形的細(xì)絲對(duì)流場(chǎng)的影響,進(jìn)而耦合控制細(xì)絲運(yùn)動(dòng)的結(jié)構(gòu)方程來(lái)對(duì)二維流場(chǎng)中細(xì)絲的流固耦合響應(yīng)進(jìn)行求解,Jia 等(2007)、Kim 等(2010)、Uddin 等(2015)都對(duì)兩個(gè)細(xì)絲串列組合的情形進(jìn)行了數(shù)值模擬研究,并對(duì)其中流動(dòng)模態(tài)的耦合特征進(jìn)行了深入分析.細(xì)絲的流固耦合擺動(dòng)問(wèn)題是實(shí)際生活中旗子迎風(fēng)飄揚(yáng)等物理現(xiàn)象的簡(jiǎn)化模型,而這類現(xiàn)象通常具有較強(qiáng)的三維效應(yīng),因此Huang 和Sung (2010)進(jìn)一步采用IB方法對(duì)一個(gè)三維旗子模型的流固耦合響應(yīng)問(wèn)題進(jìn)行了數(shù)值模擬研究,闡明了流動(dòng)雷諾數(shù)對(duì)旗子擺動(dòng)模態(tài)的影響以及相應(yīng)的三維非定常渦結(jié)構(gòu).通過(guò)三維模擬作者發(fā)現(xiàn),某些特征渦結(jié)構(gòu)能夠通過(guò)降低旗子兩側(cè)的壓差來(lái)提升旗子的穩(wěn)定性,降低擺動(dòng)強(qiáng)度.Ni 等(2023)將Huang 等(2007)中的細(xì)絲閉合成一環(huán)形并將局部固定從而組成一剛體與柔性體組合的固壁邊界,如圖7(b)所示,通過(guò)IB 方法求解了柔性體部分的彈性變形,結(jié)果表明通過(guò)改變?nèi)嵝泽w的長(zhǎng)度可以實(shí)現(xiàn)對(duì)其流固耦合運(yùn)動(dòng)模態(tài)的控制，如圖8 所示.Deng 等(2019)采用IB 方法對(duì)附著于圓柱的柔性細(xì)絲進(jìn)行了數(shù)值模擬,結(jié)果表明細(xì)絲的擺動(dòng)能夠有效地降低圓柱的阻力以及升力系數(shù)的波動(dòng).

圖7 (a) 二維柔性細(xì)絲示意圖,(b) 剛性/柔性細(xì)絲組合體

圖8 長(zhǎng)度對(duì)剛性/柔性細(xì)絲組合體流固耦合響應(yīng)及尾跡的影響(Ni et al.2023)

從上述研究工作中不難發(fā)現(xiàn),采用IB 方法能夠大幅簡(jiǎn)化包含多個(gè)運(yùn)動(dòng)邊界流固耦合問(wèn)題的網(wǎng)格復(fù)雜程度,特別是對(duì)于邊界外形復(fù)雜的情形,IB 方法中高質(zhì)量的笛卡爾網(wǎng)格以及無(wú)需網(wǎng)格重構(gòu)特點(diǎn)相較于傳統(tǒng)貼體網(wǎng)格方法而言,優(yōu)勢(shì)更加明顯.

3.3 運(yùn)動(dòng)邊界發(fā)聲

發(fā)展現(xiàn)代高階CFD/CAA 方法是處理運(yùn)動(dòng)邊界復(fù)雜流動(dòng)發(fā)聲問(wèn)題的主要途徑之一(Wang et al.2013,Slotnick et al.2014),但基于貼體網(wǎng)格的高階 CFD/CAA 方法處理涉及流-固-聲多物理場(chǎng)耦合的流致發(fā)聲問(wèn)題時(shí)一般存在網(wǎng)格重構(gòu)質(zhì)量和效率相互矛盾的局限性.為了突破這種局限性,極有希望的解決思路是利用 IB 方法構(gòu)造適合高階 CFD/CAA 方法的可壓縮體積力模型.

事實(shí)上,高階 CFD/CAA 方法結(jié)合 IB 方法的思路早在2011 年就已經(jīng)被Seo 和Mittal(2011)提出并嘗試進(jìn)行靜止物體繞流的直接噪聲模擬 (direct noise computation,DNC),從而實(shí)現(xiàn)流場(chǎng)聲場(chǎng)一體解.Mittal 等(2008)直接改進(jìn)了不可壓縮流動(dòng)下的鬼點(diǎn) IB 方法,未引入體積力模型,通過(guò)構(gòu)造高階插值方法滿足壁面邊界條件,并保持了流固界面尖銳的性質(zhì).然而,不同于傳統(tǒng)低階 CFD 方法,高階 CFD/CAA 方法中的鬼點(diǎn)法必須采用足夠多的網(wǎng)格點(diǎn)構(gòu)造高階插值,使得靠近壁面處笛卡爾網(wǎng)格上的流動(dòng)計(jì)算必須采用高度非對(duì)稱或者偏側(cè)形式的空間離散格式以及濾波或者人工黏性模板,這使得靠近壁面處的網(wǎng)格性質(zhì)同樣高度各向異性.當(dāng)處理不規(guī)則幾何邊界時(shí),這種各向異性極易激發(fā)起數(shù)值偽波,不僅會(huì)使流場(chǎng)及聲場(chǎng)數(shù)值解降階,也極有可能產(chǎn)生數(shù)值不穩(wěn)定性現(xiàn)象.這種由于偏側(cè)格式的使用造成聲場(chǎng)解降階的現(xiàn)象也在Kurbatskii和Tam (1997)關(guān)于笛卡爾網(wǎng)格下涉及曲線邊界時(shí)使用偏側(cè)耗散關(guān)系保持(dispersion relation preserving,DRP)格式時(shí)所發(fā)現(xiàn),并且根據(jù)Trefethen (1982)的理論探討,發(fā)現(xiàn)這種現(xiàn)象在頻散性有限差分方法中是廣泛存在的.Trefethen 通過(guò)分析二維有限差分方法的群速度揭示了對(duì)于波的傳播問(wèn)題,任何各項(xiàng)異性如非對(duì)稱模板數(shù)值格式、不均勻介質(zhì)或者計(jì)算網(wǎng)格等均可能導(dǎo)致雜波或者也稱為寄生波 (parasite waves) 的激發(fā).而無(wú)論是對(duì)于銳利界面IB 方法(Seo &Mittal 2011)還是Kurbatskii 和Tam (1997)針對(duì) CAA 發(fā)展的笛卡爾方法,雖然其采用的笛卡爾網(wǎng)格能夠最大程度上保持各項(xiàng)同性,但由于數(shù)值格式模板的非對(duì)稱性總會(huì)導(dǎo)致不同強(qiáng)度的寄生波,這對(duì)于運(yùn)動(dòng)邊界發(fā)聲問(wèn)題模擬的影響難以評(píng)估.據(jù)此不難發(fā)現(xiàn)對(duì)于涉及聲學(xué)問(wèn)題的模擬,盡量避免各項(xiàng)異性對(duì)數(shù)值準(zhǔn)確性和穩(wěn)定性可能具有顯著好處.

從Sirovich (1967)的流場(chǎng)延拓理論以及Peskin (1972,1977)所提出的體積力方法不難看出,Gilmanov 等(2003)、Gilmanov 和Sotiropoulos (2005)、Mittal 等(2008)、Seo 和Mittal (2011)等研究所采用的銳利界面IB 方法實(shí)際上已經(jīng)放棄了Peskin 所提出的流體域延拓以及原始體積力思想,這直接導(dǎo)致了壁面附近網(wǎng)格上流體計(jì)算必須采用非對(duì)稱數(shù)值格式的結(jié)果,這也是CAA 計(jì)算中所希望避免的.對(duì)比之下,Liu 和Vasilyev (2007)、Sun 等(2012)、Cheng 等(2018,2021b)的工作大致遵循 Peskin 所提出IB 方法中連續(xù)力源的基本思想,將固體域視為流體延拓后的偽流體域統(tǒng)一考慮,并分別針對(duì)高階 CAA 方法發(fā)展了Brinkman 罰函數(shù) IB 方法以及基于Su 等(2007)的半隱式 IB 方法發(fā)展了影響矩陣法 (influence matrix method,IMM).由于體積力模型的采用,在流固界面附近可以采用統(tǒng)一的高階中心差分格式求解.之后,Komatsu 等(2016)分析發(fā)現(xiàn)Liu &Vasilyev (2007)的罰函數(shù) IB 方法不滿足伽利略不變性,從而不適用于運(yùn)動(dòng)邊界問(wèn)題模擬,為此他們修正了能量方程中的源項(xiàng)模型,使之能夠考慮運(yùn)動(dòng)邊界發(fā)聲問(wèn)題,如振蕩圓柱直接噪聲模擬.但根據(jù)罰函數(shù)方法的構(gòu)造原理,可以發(fā)現(xiàn)其仍然存在如下缺點(diǎn): (1) 需要進(jìn)行內(nèi)外流場(chǎng)判斷,對(duì)于三維復(fù)雜幾何邊界問(wèn)題非常耗時(shí);(2) Mask 函數(shù)的采用使得跨越流固界面時(shí)體積力存在跳變;(3) 流固邊界被修改為階梯形,導(dǎo)致壁面邊界條件不能準(zhǔn)確施加在流固邊界;(4) 無(wú)法定義流固界面法線,因此不能處理無(wú)黏問(wèn)題中的無(wú)穿透壁面邊界條件.

相比之下,Sun 等(2012)、Cheng 等(2018,2021b)發(fā)展的 IMM 方法具有一定優(yōu)勢(shì),其體積力模型和原始 IB 方法一樣通過(guò)一個(gè)近似 Dirac 函數(shù)分布到流固界面周圍,因此不需要區(qū)分流固界面內(nèi)外流場(chǎng),并且壁面邊界條件可以精確施加在流固界面處,不需要修改邊界形狀,也可以處理無(wú)黏流動(dòng)問(wèn)題.Sun 等(2012)利用IMM 方法對(duì)復(fù)雜邊界的聲散射問(wèn)題進(jìn)行了數(shù)值模擬,聲場(chǎng)結(jié)果如圖9 所示,通過(guò)將數(shù)值結(jié)果與解析解進(jìn)行對(duì)比來(lái)對(duì)方法的性能進(jìn)行評(píng)估,結(jié)果表明IMM方法能夠?qū)?fù)雜邊界的聲散射問(wèn)題進(jìn)行準(zhǔn)確模擬并獲得與解析解一致的結(jié)果.然而,Sun 等(2012)、Cheng 等(2018)的工作主要針對(duì)均勻背景流動(dòng)下的線性聲散射問(wèn)題,求解線化 Euler 方程,其局限性在于忽略了背景流動(dòng)、非線性以及流體與聲學(xué)相互作用,且其采用的奇異值分解方法 (singular value decomposition,SVD) 求解力源方程組耗時(shí)巨大,難以處理復(fù)雜運(yùn)動(dòng)邊界發(fā)聲問(wèn)題.除此之外,其局限性還體現(xiàn)在方法的收斂性上,和一般貼體網(wǎng)格方法不同,IMM 由于涉及邊界網(wǎng)格和笛卡爾背景網(wǎng)格,其求解的體積力收斂性以及穩(wěn)定性該如何定義仍然模糊不清.這些問(wèn)題直接影響其是否能夠應(yīng)用到實(shí)際復(fù)雜運(yùn)動(dòng)邊界發(fā)聲的多物理場(chǎng)耦合問(wèn)題.

圖9 多個(gè)圓柱散射下的聲場(chǎng)壓力分布(Sun et al.2012)

在IMM 方法的基礎(chǔ)上,Cheng 等(2021b)利用預(yù)測(cè)-校正思想進(jìn)一步推導(dǎo)了非線性主控方程下的可壓縮體積力模型,并且通過(guò)分析體積力模型的物理意義構(gòu)造模型方程實(shí)現(xiàn)該類方法的收斂性和穩(wěn)定性分析,從而據(jù)此開(kāi)發(fā)了魯棒性更強(qiáng)的快速求解算法.由于流場(chǎng)外內(nèi)是采用統(tǒng)一的插值模板進(jìn)行求解,因此邊界處不存在由于采用偏側(cè)差分格式所引起的數(shù)值偽波,數(shù)值穩(wěn)定性明顯提高.利用上述算法,Cheng 等(2021a)研究了二維葉柵非定常流動(dòng)中振動(dòng)誘發(fā)聲共振問(wèn)題的物理機(jī)制,揭示了聲共振工況下葉柵流場(chǎng)分布特征,如圖10 所示.上述基于IB 方法的CAA 算法還被應(yīng)用于更加復(fù)雜的三維運(yùn)動(dòng)邊界流致發(fā)聲問(wèn)題,例如開(kāi)式轉(zhuǎn)子噪聲的產(chǎn)生及傳播,如圖11所示,上述結(jié)果對(duì)于理解氣動(dòng)噪聲背后的流動(dòng)機(jī)理具有重要的工程價(jià)值.

圖10 葉柵聲共振所對(duì)應(yīng)的流場(chǎng)特征.(a) 壓力云圖,(b) 渦量云圖(Cheng et al.2021a)

圖11 開(kāi)式轉(zhuǎn)子流致發(fā)聲直接數(shù)值模擬結(jié)果.(a) 瞬時(shí)壓力紋影圖,(b) 馬赫數(shù)云圖分布

采用IB 方法進(jìn)行CAA 計(jì)算的一個(gè)關(guān)鍵問(wèn)題仍然在于邊界處的處理精度,上文中提到,在連續(xù)力源法中采用連續(xù)分布的函數(shù)來(lái)代替Dirac 函數(shù)導(dǎo)致邊界處通常只能達(dá)到一階精度,這與CAA 計(jì)算中追求的高階精度是相矛盾的,但是從Sun 等(2012)、Cheng 等(2018,2021a)的研究結(jié)果來(lái)看,邊界處的一階精度似乎并不會(huì)對(duì)聲學(xué)模擬的總體準(zhǔn)確性產(chǎn)生特別大的影響,但是進(jìn)一步提高邊界處的處理精度始終是研究IB 方法所追求的目標(biāo)之一.Liang 等(2008) 針對(duì)間斷問(wèn)題的數(shù)值振蕩現(xiàn)象,提出了一種對(duì)間斷函數(shù)的全局描述思想,將間斷函數(shù)表示為光滑函數(shù)和修正項(xiàng)之和,其中修正項(xiàng)由間斷處的躍度條件確定,通過(guò)這樣的轉(zhuǎn)換,在求解微分方程時(shí),一個(gè)存在間斷的問(wèn)題可以轉(zhuǎn)換為光滑問(wèn)題,進(jìn)一步可將譜方法應(yīng)用于間斷問(wèn)題而不會(huì)帶來(lái)Gibbs 振蕩,其數(shù)值結(jié)果表明,采用此種求解方式可有效提高邊界處的求解精度,精度階數(shù)與所構(gòu)造的階躍條件密切相關(guān).這本質(zhì)上屬于一種浸入式界面方法(immersed interface method,IIM) (LeVeque &Li 1994,Xu &Wang 2006a,2006b,Zhong 2007,Gabbard 等 2022),雖然能夠有效地提高間斷處理精度,但是階躍條件的復(fù)雜程度以及相應(yīng)的計(jì)算量都會(huì)隨精度的提升而大幅提高,如何進(jìn)一步提高其數(shù)值穩(wěn)定性并拓展其在CAA 中的應(yīng)用仍有待進(jìn)一步研究.

4 高雷諾數(shù)流動(dòng)模擬

4.1 實(shí)現(xiàn)高雷諾數(shù)流動(dòng)模擬的途徑

雖然IB 方法在運(yùn)動(dòng)邊界及復(fù)雜邊界流動(dòng)模擬等問(wèn)題上展現(xiàn)出了極大的優(yōu)勢(shì),但是其主要劣勢(shì)在于高雷諾數(shù)流動(dòng)模擬時(shí)龐大的計(jì)算量,這一點(diǎn)也直接導(dǎo)致該方法目前在實(shí)際工程領(lǐng)域難以廣泛應(yīng)用.在低雷諾數(shù)層流條件下,笛卡爾網(wǎng)格所需的單元總量與貼體網(wǎng)格的差異并不明顯,但是在高雷諾數(shù)條件下,倘若笛卡爾網(wǎng)格與貼體網(wǎng)格同樣滿足邊界層分辨率的要求,那么笛卡爾網(wǎng)格所需的單元總量會(huì)遠(yuǎn)遠(yuǎn)超出貼體網(wǎng)格,二維及三維情形下,二者網(wǎng)格總量的比值大約分別與雷諾數(shù)及其 1.5 次方呈正比(Mittal &Iaccarino 2005),因此高雷諾數(shù)流動(dòng)模擬對(duì)于浸入式邊界方法來(lái)說(shuō)是一個(gè)巨大的挑戰(zhàn).

采用自適應(yīng)網(wǎng)格加密 (adaptive mesh refinement,AMR) 方式是減少湍流模擬時(shí)網(wǎng)格量的主要方法,這一方法可以直接對(duì)目標(biāo)區(qū)域進(jìn)行局部加密從而避免傳統(tǒng)拉伸網(wǎng)格方式引起的遠(yuǎn)場(chǎng)網(wǎng)格總量上升,如圖12 所示.Berger 和Oliger(1984)、Berger 和Colella(1989)曾采用自適應(yīng)加密的方式來(lái)提升網(wǎng)格在激波間斷處的流場(chǎng)分辨率并說(shuō)明這一方式可以高效地實(shí)現(xiàn)局部網(wǎng)格密度的提升,適合用于準(zhǔn)確描述流場(chǎng)中存在高梯度的區(qū)域.同時(shí),Roma 等(1999)、De Tullio 等(2007)、Angelidis 等(2016)、Al-Marouf 和Samtaney 等(2017)的研究工作都曾在浸入式邊界方法的應(yīng)用中采用自適應(yīng)網(wǎng)格的方法來(lái)提升壁面處的網(wǎng)格分辨率,從而達(dá)到減少網(wǎng)格總量的目的.Wang 等(2020)在有限差分的框架下將自適應(yīng)網(wǎng)格技術(shù)與IB 方法結(jié)合并將其應(yīng)用于運(yùn)動(dòng)邊界問(wèn)題的模擬,結(jié)果表明該數(shù)值工具能夠準(zhǔn)確對(duì)包含運(yùn)動(dòng)邊界的流動(dòng)問(wèn)題進(jìn)行預(yù)測(cè),且計(jì)算量相較于傳統(tǒng)的正交笛卡爾網(wǎng)格大幅降低.而對(duì)于目前已經(jīng)發(fā)表的湍流模擬相關(guān)工作,如Capizzano (2011,2016)、Tamaki 等(2017)、Pu 和Zhou (2018)、Berger 和Aftosmis (2018)、Xu 和Liu (2021)、Cai 等(2021)、Constant 等(2021)、Wang 等(2023)的研究工作,均采用了這種形式的網(wǎng)格來(lái)減少流動(dòng)模擬所需的網(wǎng)格總量.

圖12 采用自適應(yīng)網(wǎng)格加密生成的多層網(wǎng)格系統(tǒng)(Wang et al.2020)

另一方面,即使對(duì)笛卡爾網(wǎng)格進(jìn)行局部加密,想要完全滿足 RANS 模型的y+要求依然需要較大的網(wǎng)格量,因此,上述基于笛卡爾網(wǎng)格的湍流模擬研究工作均采用了不同形式的湍流壁面函數(shù)來(lái)放寬對(duì)于網(wǎng)格y+的要求.壁面函數(shù)將近壁面的流動(dòng)分為兩層,外層與靠近邊界的內(nèi)層分別用RANS 方程和薄邊界層方程來(lái)描述,對(duì)于平板問(wèn)題,如果將薄邊界層方程的壓力項(xiàng)忽略則可以得到速度剖面的解析表達(dá)式,也稱為顯式壁面律,如方程(6)所示(Spalding 1961).Capizzano(2011)在浸入式邊界方法的基礎(chǔ)上引入顯式壁面律來(lái)對(duì)湍流邊界層進(jìn)行建模,在此基礎(chǔ)上對(duì)二維及三維翼型繞流進(jìn)行了測(cè)試,所得到的結(jié)果與實(shí)驗(yàn)結(jié)果較為接近.值得注意的是,顯式壁面律由于忽略了邊界層內(nèi)部的壓力梯度,在實(shí)際應(yīng)用可能導(dǎo)致固壁表面壓力分布等參數(shù)存在一定的數(shù)值的振蕩,同時(shí)對(duì)于摩擦系數(shù)的預(yù)測(cè)也會(huì)在某些壓力梯度劇烈變化的區(qū)域偏離真實(shí)結(jié)果.Wang 等(2023)在銳利界面IB 方法的基礎(chǔ)上,采用方程(6)中的顯式壁面律對(duì)流動(dòng)進(jìn)行模擬,結(jié)果表明即使采用最為簡(jiǎn)單的顯式壁面律也能夠準(zhǔn)確預(yù)測(cè)高雷諾數(shù)流動(dòng)中翼型表面的壓力系數(shù)分布,但是摩擦系數(shù)在局部區(qū)域的收斂性并不完美,并且同樣存在一定的數(shù)值振蕩.Chen 等(2023)將一種基于SA 湍流模型壁面漸進(jìn)解的顯式壁面模型與IB 方法相結(jié)合,計(jì)算所得到的湍流邊界層速度與實(shí)驗(yàn)測(cè)量所得到的壁面律吻合得較好,如圖13 所示.Tamaki 等(2017)、Cai 等(2021)、Chen 等(2021,2023)均嘗試通過(guò)在近壁面采用經(jīng)過(guò)線化后的速度梯度分布并且在壁面處施加滑移邊界條件的方式來(lái)抑制數(shù)值振蕩,數(shù)值結(jié)果表明雖然采用線化的速度梯度分布能夠有效地抑制數(shù)值振蕩,但是局部收斂性仍存在一定的偏離.

圖13 基于SA 顯式壁面模型的湍流邊界層速度分布與壁面律對(duì)比(Chen et al.2023)

從上述研究結(jié)果來(lái)看,如何進(jìn)一步提高壁面模型的準(zhǔn)確性是目前IB 方法所面臨的核心問(wèn)題之一.顯式壁面律的推導(dǎo)過(guò)程中忽略了壓力梯度的影響,這一假設(shè)對(duì)于平板流動(dòng)是適用的,但是對(duì)于某些存在明顯壓力梯度的復(fù)雜流動(dòng),顯式壁面律可能存在較大的誤差.為了提升對(duì)近壁面流動(dòng)的刻畫(huà)精度,Capizzano (2016)采用了一組包含動(dòng)量、內(nèi)能以及湍流輸運(yùn)量的邊界層流動(dòng)控制方程,動(dòng)量方程中包含了壓力梯度的影響,在每個(gè)時(shí)間步對(duì)其進(jìn)行同步離散求解從而給定邊界處網(wǎng)格單元數(shù)值,測(cè)試結(jié)果表明在某些壓力變化劇烈的區(qū)域,包含壓力梯度的處理方式能夠使結(jié)果更接近貼體網(wǎng)格結(jié)果.基于邊界層方程的壁面函數(shù)還被Shi 等(2019)、Ma 等(2019,2021)用于發(fā)展LES 算法以實(shí)現(xiàn)IB 方法在高雷諾數(shù)流動(dòng)問(wèn)題中的應(yīng)用,上述學(xué)者的LES 模擬結(jié)果表明,實(shí)現(xiàn)準(zhǔn)確預(yù)測(cè)需在邊界層附近對(duì)LES 亞格子應(yīng)力模型做適當(dāng)?shù)男拚云胶獗诿婺Ｐ退峁┑哪Σ翍?yīng)力.Berger 和Aftosmis (2018)采用了一組包含壓力梯度的常微分方程對(duì)邊界層流動(dòng)進(jìn)行描述,結(jié)果相較于顯式的壁面律也表現(xiàn)出一定程度的改善.當(dāng)壓力梯度的影響被包含于邊界層方程中時(shí),無(wú)法推導(dǎo)得到解析的速度分布,因此需要對(duì)邊界層控制方程進(jìn)行同步的離散求解,這一求解過(guò)程需要耗費(fèi)一定的計(jì)算時(shí)間,雖然模型的準(zhǔn)確性及適用性有所提升,但是整體計(jì)算效率也會(huì)出現(xiàn)下降.Constant 等(2021)、Xu 和Liu (2021)則通過(guò)增加邊界單元數(shù)量的方式在邊界單元的最外層實(shí)現(xiàn)幾乎一致的y+分布,對(duì)比結(jié)果表明采用此種處理方式,即使對(duì)于方程(6)所給出的顯式壁面律,依然能夠獲得準(zhǔn)確且光滑的壓力及摩擦系數(shù)分布,并且壓力系數(shù)的收斂性同樣出現(xiàn)明顯的改善.如何進(jìn)一步提高IB 方法與壁面模型結(jié)合時(shí)的準(zhǔn)確性以及穩(wěn)定性,同時(shí)兼顧計(jì)算效率,目前仍有待深入研究.

4.2 高雷諾數(shù)內(nèi)流問(wèn)題的應(yīng)用

目前IB 方法主要用于各類外流問(wèn)題,對(duì)于這種情形,通?？梢圆捎眠h(yuǎn)大于固體邊界的方形計(jì)算域并在其中布置正交的笛卡爾網(wǎng)格.而對(duì)于內(nèi)流問(wèn)題,流動(dòng)通常被限制于流道內(nèi)部,與此同時(shí)流道內(nèi)部還可能存在運(yùn)動(dòng)邊界,如對(duì)于上文所提到的航空葉輪機(jī)械領(lǐng)域的轉(zhuǎn)靜干涉、葉片顫振等問(wèn)題,均屬于這一類型.在內(nèi)流條件下,流道邊界所覆蓋的空間范圍可能遠(yuǎn)大于內(nèi)部物體,如圖14(a) 所示,流道邊界通常是靜止的,只有流道內(nèi)部的固體邊界可能是運(yùn)動(dòng)邊界.從圖14(b)可以看出,此時(shí)仍然采用笛卡爾網(wǎng)格對(duì)流道以及內(nèi)部固壁邊界進(jìn)行建模時(shí),用于識(shí)別流道邊界所需的網(wǎng)格量會(huì)遠(yuǎn)超內(nèi)部固壁邊界所需的網(wǎng)格量,即使是在圖14(b)中已經(jīng)采用自適應(yīng)網(wǎng)格的情形下,流道邊界所引入的網(wǎng)格量仍然不可忽視,而在三維高雷諾數(shù)的情形下,這一問(wèn)題會(huì)變得更加突出.對(duì)于圖14 中所示的內(nèi)流問(wèn)題二維邊界模型,Wang 等(2023)給出了定量的網(wǎng)格單元數(shù)量分析,其中用于識(shí)別整個(gè)流道邊界的網(wǎng)格單元數(shù)量占總數(shù)量的70%左右,這一比例與流道邊界覆蓋的范圍密切相關(guān),當(dāng)流道邊界覆蓋的區(qū)域越寬廣,這一比例也會(huì)越高.前文中提到,IB 方法的主要優(yōu)勢(shì)在于處理復(fù)雜邊界和運(yùn)動(dòng)邊界問(wèn)題,而圖14(a)中的流道邊界并不屬于上述任意一種情形,因此想要進(jìn)一步提高IB 方法在處理三維高雷諾數(shù)內(nèi)流問(wèn)題上的計(jì)算效率,則需要對(duì)流道邊界的處理進(jìn)行改進(jìn).

圖14 (a) 內(nèi)流問(wèn)題示意圖,(b) 采用笛卡爾網(wǎng)格結(jié)合自適應(yīng)加密所生成的網(wǎng)格系統(tǒng),(c) 流道貼體網(wǎng)格結(jié)合自適應(yīng)加密生成的網(wǎng)格系統(tǒng)(Wang et al.2023)

針對(duì)內(nèi)流中的運(yùn)動(dòng)邊界問(wèn)題,Ge 和Sotiropoulos (2007)提出了一種適用于曲線網(wǎng)格的浸入式邊界方法,該研究工作首先采用曲線的貼體網(wǎng)格覆蓋整個(gè)流道邊界,然后在內(nèi)部運(yùn)動(dòng)邊界附近構(gòu)造局部正交或近似正交的笛卡爾網(wǎng)格從而使得IB 方法能夠得以應(yīng)用,如圖14(c)所示,該方法避免了采用笛卡爾網(wǎng)格對(duì)流道邊界進(jìn)行處理,從而達(dá)到了降低網(wǎng)格量、提高計(jì)算效率的目的.Ubald 等(2021)采用同樣的思想對(duì)一個(gè)前緣帶有圓柱擾流器的渦輪葉片進(jìn)行了模擬,該研究中流道以及葉片均采用貼體網(wǎng)格進(jìn)行覆蓋,通過(guò)在葉片前緣構(gòu)建局部的笛卡爾網(wǎng)格來(lái)應(yīng)用IB 方法對(duì)圓柱擾流器進(jìn)行建模從而大幅降低網(wǎng)格復(fù)雜程度.類似地,Mochel 等(2014)、Weiss 和Deck(2018)對(duì)火箭和導(dǎo)彈的主體采用貼體網(wǎng)格,在此基礎(chǔ)上對(duì)部分固體型面如鋸齒尾緣、推進(jìn)器上的復(fù)雜結(jié)構(gòu)采用IB 方法進(jìn)行建模,結(jié)果表明此種手段既能降低網(wǎng)格的生成難度,又能夠兼顧計(jì)算效率.上述工作為處理內(nèi)流環(huán)境中的運(yùn)動(dòng)/復(fù)雜邊界問(wèn)題提供了一條可行的思路,但是對(duì)于更加復(fù)雜的運(yùn)動(dòng)邊界問(wèn)題還需要對(duì)上述方法進(jìn)行進(jìn)一步的改進(jìn)以增強(qiáng)其處理工程領(lǐng)域復(fù)雜運(yùn)動(dòng)邊界問(wèn)題的能力.例如對(duì)Du 等(2016a,2016b)、Chen 等(2021)的研究工作中所關(guān)注的葉輪機(jī)械領(lǐng)域的轉(zhuǎn)靜干涉問(wèn)題,此時(shí)運(yùn)動(dòng)邊界靠近流道邊界,難以在流道附近構(gòu)建局部的笛卡爾網(wǎng)格,因此上述的方法在處理此類問(wèn)題時(shí)還需要進(jìn)一步改進(jìn).

Wang 等(2023)的分析表明,當(dāng)流道內(nèi)的運(yùn)動(dòng)邊界靠近流道邊界時(shí),由于識(shí)別流道邊界所采用的貼體網(wǎng)格僅沿壁面法向具有較高的分辨率,因此網(wǎng)格的展弦比較大,此時(shí)由于運(yùn)動(dòng)邊界與網(wǎng)格單元的相對(duì)朝向未知.傳統(tǒng)的銳利界面IB 方法是面向均勻正交的笛卡爾網(wǎng)格發(fā)展而來(lái),通常采用統(tǒng)一的外伸插值距離,而此時(shí)由于網(wǎng)格展弦比較大,統(tǒng)一的外伸插值距離會(huì)出現(xiàn)失效的情形.為解決這一問(wèn)題,Wang 等(2023)針對(duì)銳利界面IB 方法提出了一種自適應(yīng)外伸插值距離的方法,成功將IB 方法拓展至一般的曲線網(wǎng)格,使其能夠處理邊界與不同展弦比切割的情形,結(jié)合上文提到的自適應(yīng)網(wǎng)格加密以及壁面函數(shù),成功對(duì)三維亞聲速平面葉柵以及NASA 轉(zhuǎn)子37 內(nèi)部的跨音流動(dòng)進(jìn)行了模擬.結(jié)果表明,該方法能夠獲得與貼體網(wǎng)格一致的結(jié)果,并且能夠準(zhǔn)確預(yù)測(cè)葉柵葉片表面的壓力分布以及跨音轉(zhuǎn)子的特性曲線,其中葉柵葉片表面的壓力分布數(shù)值模擬結(jié)果及其與實(shí)驗(yàn)的對(duì)比如圖15 所示.在此基礎(chǔ)上,Chen 等(2023)采用上述方法對(duì)一低壓渦輪葉柵以及涵道風(fēng)扇流動(dòng)進(jìn)行了數(shù)值模擬研究,結(jié)果表明該方法同樣能夠準(zhǔn)確預(yù)測(cè)不同工況下葉片表面壓力及下游尾跡的分布狀況,而計(jì)算效率相較于采用傳統(tǒng)的正交笛卡爾網(wǎng)格而言大幅提高.上述研究表明,采用自適應(yīng)外伸插值距離能夠提高IB 方法與不同類型網(wǎng)格的結(jié)合能力從而拓寬方法的適用范圍,但值得注意的是,Tamaki 等(2017)、Berger 和Aftosmis (2018)、Chen 等(2023)、Wang 等(2023)的研究均表明,對(duì)于實(shí)際工程領(lǐng)域的高雷諾數(shù)流動(dòng),流場(chǎng)中可能存在流動(dòng)分離現(xiàn)象,采用湍流壁面模型在這些分離區(qū)域可能出現(xiàn)明顯的預(yù)測(cè)偏差,但是目前IB 方法在高雷諾數(shù)流動(dòng)的應(yīng)用中仍然非常依賴于壁面模型來(lái)提高計(jì)算效率.因此,如何進(jìn)一步提高該方法對(duì)不同類型復(fù)雜流動(dòng)的適應(yīng)能力是將其推向工程應(yīng)用所需要解決的關(guān)鍵問(wèn)題之一,一種可能地解決手段是采用脫體渦方法或大渦模擬等高保真模擬方法來(lái)提高復(fù)雜流動(dòng)的預(yù)測(cè)精度,但其對(duì)計(jì)算資源的需求使其目前還無(wú)法實(shí)現(xiàn)規(guī)?；こ虘?yīng)用.

圖15 三維亞聲速平面葉柵葉片表面壓力分布數(shù)值模擬結(jié)果及其與實(shí)驗(yàn)的對(duì)比(Wang et al.2023)

5 總結(jié)

本文主要對(duì)IB 方法中邊界的建模方式,在復(fù)雜邊界、運(yùn)動(dòng)邊界以及運(yùn)動(dòng)邊界發(fā)聲問(wèn)題等問(wèn)題的應(yīng)用以及高雷諾數(shù)流動(dòng)模擬中的研究進(jìn)展進(jìn)行了綜述.目前IB 方法已廣泛應(yīng)用于諸如生物體流動(dòng)等各類的低雷諾數(shù)流動(dòng)研究中,其核心優(yōu)勢(shì)在于處理復(fù)雜邊界以及運(yùn)動(dòng)邊界的情形.而對(duì)于航空航天、建筑橋梁等領(lǐng)域中包含多運(yùn)動(dòng)邊界的復(fù)雜工程問(wèn)題,IB 方法也有所涉及且同樣能夠有效地降低網(wǎng)格復(fù)雜程度,但由于這些流動(dòng)通常具有較高的雷諾數(shù),而實(shí)現(xiàn)高雷諾數(shù)流動(dòng)的準(zhǔn)確模擬需要投入大量的計(jì)算資源,目前實(shí)現(xiàn)大規(guī)模的應(yīng)用仍需要進(jìn)一步的研究.未來(lái)研究的關(guān)鍵在于進(jìn)一步提高IB 方法針對(duì)高雷諾數(shù)流動(dòng)的計(jì)算效率,同時(shí)提升計(jì)算模型處理三維復(fù)雜高雷諾數(shù)流動(dòng)現(xiàn)象的能力.

致 謝國(guó)家自然科學(xué)基金(52022009)資助項(xiàng)目.