浸入式邊界方法的研究及應用進展

王 卓 杜 林 成 龍 孫曉峰

1 北京航空航天大學航空發動機研究院,北京 100191

2 華為技術有限公司,編譯器與編程語言實驗室,北京 100094

3 北京航空航天大學能源與動力工程學院,北京 100191

1 概述

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

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

圖1 用于研究心臟血液流動的模型示意圖

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

盡管IB 方法在處理復雜邊界以及運動邊界問題上具有強大的能力,但是其主要缺陷在于處理高雷諾數流動問題時龐大的計算量.由于湍流邊界層厚度相較于層流而言非常小,通常需要在壁面處布置高分辨率網格來分辨湍流邊界層的線性底層,例如對于Spalart-Allmaras (SA)一方程湍流模型而言(Spalart &Allmaras 1992),通常要求壁面第一層網格的y+小于1.而對于IB 方法來說,由于對笛卡爾的網格加密必須沿空間三個方向同步展開,粗略估計表明,三維情形下網格總量與流動雷諾數的1.5 次方呈正比,而對貼體網格而言僅與雷諾數的0.5 次方呈正比,這也導致對于高雷諾數流動,IB 方法所需的網格量遠遠超出貼體網格.如何在現有的硬件計算能力基礎上拓展IB 方法在高雷諾數流動問題中的應用是進一步擴大該方法的應用范圍以及實現工業規?；瘧玫暮诵奶魬?

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

2 基本思想

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

根據力源建模方式,IB 方法可以進一步細分為連續力源法和離散力源法.

2.1 連續力源法

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

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

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

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

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

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

2.2 離散力源法

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

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

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

3 IB 方法在運動邊界問題中的應用

3.1 生物流動

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

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

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

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

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

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

3.2 流固耦合

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

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

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

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

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

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

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

圖8 長度對剛性/柔性細絲組合體流固耦合響應及尾跡的影響(Ni et al.2023)

從上述研究工作中不難發現,采用IB 方法能夠大幅簡化包含多個運動邊界流固耦合問題的網格復雜程度,特別是對于邊界外形復雜的情形,IB 方法中高質量的笛卡爾網格以及無需網格重構特點相較于傳統貼體網格方法而言,優勢更加明顯.

3.3 運動邊界發聲

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

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

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

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

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

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

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

圖11 開式轉子流致發聲直接數值模擬結果.(a) 瞬時壓力紋影圖,(b) 馬赫數云圖分布

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

4 高雷諾數流動模擬

4.1 實現高雷諾數流動模擬的途徑

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

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

圖12 采用自適應網格加密生成的多層網格系統(Wang et al.2020)

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

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

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

4.2 高雷諾數內流問題的應用

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

圖14 (a) 內流問題示意圖,(b) 采用笛卡爾網格結合自適應加密所生成的網格系統,(c) 流道貼體網格結合自適應加密生成的網格系統(Wang et al.2023)

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

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

圖15 三維亞聲速平面葉柵葉片表面壓力分布數值模擬結果及其與實驗的對比(Wang et al.2023)

5 總結

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

致 謝國家自然科學基金(52022009)資助項目.