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

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

* (南京航空航天大學航空學院非定常空氣動力學與流動控制工業和信息化部重點實驗室,南京 210016)

引言

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

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

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

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

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

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

1 物理模型

如圖1 所示為擴散界面背景下IB 方法的原理圖,流體域Ω采用笛卡爾網格離散,浸入流體的物體邊界Γ由獨立于背景網格的拉格朗日點替代.歐拉流體變量和拉格朗日邊界變量通過狄拉克函數關聯,相互作用方程如下

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

2 計算模型

2.1 格子Boltzmann (LB)模型

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

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

取Maxwell 分布函數在離散速度空間 {cα} 中的二階速度截斷得局部平衡分布函數的表達式

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

Fα表達式為

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

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

(1) 碰撞

(2) 邊界條件施加

(3)遷移

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

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

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

2.2 信息交互界面

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

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

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

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

轉換矩陣DE和DL的表達式如下

對于剛性靜止壁面,DE與DL只需計算一次;對于運動變形界面,每個時間步,由結構運動更新邊界位置 {X1,X2,···,Xk,···,XNL} 后,需要重新構造影響域覆蓋的歐拉節點集合 {x1,x2,···,xl,···,xNE},并更新NE,然后重新建立DE與DL.

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

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

2.3 邊界力模型

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

首先,根據LB 模型的動量方程(12)構造體積力密度f的直接力表達式

式中,δu=u-u*,δu為歐拉節點處f引發的流體速度修正量,為LB 模型的預計算速度,u為重構的流體速度.

根據力同一性的概念,歐拉/拉格朗日變量描述的浸沒物體經歷的合流體力等價,即

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

綜上得邊界力密度F的表達式

GL中拉格朗日節點處流體的質量密度矢量ML={M1,M2,···,MNL}來自歐拉流體節點的質量密度插值

為了保證體積力重構的速度場滿足無滑移約束,建立如下方程組

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

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

邊界力密度FL由子迭代的累和計算

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

圖2 直接力格式IB-LB 模型算法流程圖.淺灰色框指明物面參數輸出子步驟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 數值方法驗證與分析

本節以二維Taylor-Green 渦流、靜止圓柱繞流、振蕩圓柱繞流及波動翼型繞流4 個經典案例的模擬結果佐證改進的IB-LB 模型的數值精度、預測結果可靠性、及其在復雜運動界面繞流問題中的實用性.

3.1 Taylor-Green 渦流

作為評估算法數值精度的基準算例,Taylor-Green 渦流已由幾種不同的IB-LB 模型[15,17,19]執行模擬.理想的IB 模型應具備保持背景LB 求解器固有數值精度的能力.本節采用Taylor-Green 渦流模擬結果來驗證本文改進的IB 模型對背景LB 求解器精度的影響.Taylor-Green 渦流在[-L,L]×[-L,L]的計算域中解析解表達式如下

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

式(32) 定義了速度的全局數值誤差L2,這里,N是計算域中歐拉流體節點的總數,unum和uexa分別為速度的數值解和精確解.數值精度通過數值誤差L2與網格尺度h的對數斜率來量化

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

圖3 對數坐標下的速度誤差L2 與網格尺度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 靜止圓柱繞流

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

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

表1 Re=20 和40 時規則區長度和阻力系數與已有文獻結果的比較Table 1 Comparison of the obtained recirculating length,the drag coefficient with the previous results at Re=20 and 40

圖4 展示了Re=20 和40 條件下,二維圓柱表面壓強系數CP與方向角θ的關系曲線,方向角自前駐點θ=0°經圓柱上方指向后基點θ=180°,Cp(θ)=.圖4 中壓強分布曲線光滑無震蕩,且當前模型結果與前人研究結果[14-15,29]吻合良好,壓強極大值出現在前滯點(θ=0°),壓強極小值位于分離點附近,圓柱背后分離區壓強有所回升,以上結果驗證了本文改進的IB-LB 模型能夠為復雜結構外形上的拉格朗日節點提供光滑、可靠的壓強預測.

圖4 圓柱表面壓強系數與方向角關系曲線Fig.4 Pressure coefficient curve of cylinder versus the orientation angle

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

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

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

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

3.3 流動經過強制振動圓柱

在上述Re=40 的定常流動中,強制二維圓柱執行周期性的側向振動擾動周圍流體,圓柱體的側向振蕩由如下正弦函數控制

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

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

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

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

3.4 波動翼型繞流

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

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

側向平移行波方程如下

式中,側向平移幅值A1=0.15L.

旋轉運動以距離翼型前緣點0.3L的位置為轉動中心,旋轉運動行波方程為

式中,旋轉幅值θ2=7.75°.

波動運動自距離翼型前緣點0.4L的位置起執行,側向波動幅值A3(x/L)為x坐標的二次函數[32],波動行波方程為

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

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

圖7 波動翼型的時間相關推力系數Fig.7 Time-dependent thrust coefficient acting on the undulatory foil

圖8 展示了維持翼型復合運動的運動功率系數CPS,克服運動阻力的功率系數CPD,以及維持運動的總功率系數CPT隨時間變化曲線.其中,運動功率系數CPS根據如下公式計算

圖8 波動翼型的時間相關運動功率系數、阻力功率系數和總功率系數Fig.8 Time-dependent kinematic power coefficient,overcome drag power coefficient,and total power coefficient required by the undulatory foil

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

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

圖9 均勻流動中復合運動翼型周圍渦量云圖Fig.9 Vortex contour near the swimming foil in a uniform flow

4 結論

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