基于浸潤邊界-格子波爾茲曼通量求解器的柔性結構流固耦合數值模擬

劉 釩,劉 剛,江 雄,舒 昌

(1.中國空氣動力研究與發展中心,綿陽 621000;2.新加坡國立大學 工程學院,新加坡 117576)

0 引 言

具有動態復雜物面邊界物體的繞流問題,是非定常流體力學研究的一個重要課題。特別地,在以仿生微型飛行器為代表的低雷諾數大變形柔性結構流固耦合仿真中具有重要的研究價值[1-2]。在流固耦合數值模擬中,非定常邊界流場與柔性結構變形需要進行耦合計算,無法事先預測的動態變形物面對所使用的非定常流場求解器的動邊界處理能力提出了特殊的要求?；谫N體網格的流場求解器在處理此類問題時需要進行動態網格重構,這往往帶來復雜繁瑣的計算,一般會增加計算時間并降低求解魯棒性。而基于非貼體網格的流場求解方法,如基于笛卡爾網格的浸潤邊界法(Immersed Boundary Method,IBM)有效地規避了網格重構過程。在浸潤邊界法中,物面與流體物質之間的相互作用被轉化為流場方程中的一個體積力項,可對物面兩側區域不加區分地進行流場統一求解。

格子波爾茲曼方法(Lattice Boltzmann Method,LBM)在處理低速不可壓流動問題時,相對于直接求解N-S方程,無需處理壓力泊松方程和應用交錯網格算法,具有相對簡便高效的特性?；谙嗤姆琴N體笛卡爾網格框架,可將LBM與浸潤邊界法相結合,產生了浸潤邊界-格子波爾茲曼方法(IB-LBM)[3,5]。為將其應用有限體積框架內以求解更廣泛的問題,學術界近年來提出了一種新的求解方案,即格子波爾茲曼通量求解器(Lattice Boltzmann Flux Solver,LBFS)[5-9,11]。該方法通過在單元界面處的LBM模型計算流場通量值,且由于無需存儲密度分布函數節約了內存空間。由于該方法基于有限體積法,故可使用非均勻網格求解,從而擺脫了計算時間步長和單元網格尺寸之間的綁定關系,提高了流場數值計算的靈活性和效率。進一步將LBFS與浸潤邊界法相結合以處理復雜動態邊界流動問題,構建了浸潤邊界-格子波爾茲曼通量求解器(IB-LBFS)。同時,引入了隱式速度邊界修正方程,實現了物面無滑移流場邊界條件的精確滿足。

在當前IB-LBFS方法的基礎上,為了將其應用于真實三維大變形柔性結構的流固耦合數值模擬,還需要進行兩方面的工作:首先,需要提升當前IB-LBFS求解器的計算效率,實現可處理大規模網格的并行化計算;同時,需要建立柔性結構動力學求解器,并構建其與IB-LBFS流場求解器之間的流-固交互界面。本文搭建了基于IB-LBFS和絕對節點坐標板殼單元結構力學求解器的流固耦合求解平臺,并通過幾個典型算例確認了該耦合求解器的有效性。

1 浸潤邊界－格子波爾茲曼通量求解器

1.1 控制方程

宏觀流動的Navier-Stokes方程為:

在格子波爾茲曼通量計算中,動量通量ρu和黏性-無黏通量Π的表達式為:

其中:eα為粒子速度,fα為α方向上的密度分布函數,是相應方向上的平衡狀態分布函數為非平衡狀態分布函數。根據LBGK模型的表達式為:

eα和系數wα的選擇取決于不同的格子速度模型。IB-LBFS的控制方程組為:

其中F為浸潤邊界法中物面邊界對應的體積力項。

1.2 離散求解

將式(5)進行時間離散,將求解過程分裂為一個預測步和一個修正步(分別對應式(6)和式(7)):

分別求解式(6)和式(7),可得到中間速度u?和修正速度δu,二者之和為真實流場速度un＋1。

將式(6)離散于有限體積的單元控制體,各個速度方向α上的平衡狀態函數的值由相鄰單元的守恒變量在相應位置插值并轉換得到。經界面重構(圖1),得到界面上的密度ρ、動量通量ρu、無黏-黏性通量Π,使用單元有限體積控制方程解得流場密度場和中間速度場u?。

圖1 相鄰單元界面上的LBM重構(D2Q9模型)Fig.1 LBM reconstruction at cell interface(D2Q9 model)

為了滿足無滑移邊界條件,引入了隱式邊界修正法,即將體積力項F的計算轉化為物面邊界引起的流場修正速度δu的求解。根據無滑移條件,物面速度矢量Un＋1(XlB)的值應與流場在物面相應位置上的流場速度相同,這一關系可由如下的插值關系得到:

其中un＋1＝u?＋δu。作為未知數的物面網格速度修正量假定為,作從物面拉格朗日網格點至流場歐拉網格點的插值,有:

N為構成物面的拉格朗日表面點總數;M為Euler網格點總數;Dij是連續核插值函數,對三維情況有(R為臨界半徑):

2 IB-LBFS：加速優化計算和并行算法

在非定常動邊界流動的IB-LBFS計算中,修正速度方程系數矩陣A的相關計算步驟(矩陣計算、存儲、線性方程組求解)占據了相當的計算時間,為了提高運算效率,有對其進行算法優化的必要。進行了以下三處算法優化:

(1)利用變量物理意義,加速矩陣A生成計算:由于動態邊界的存在,每個非定常時間步中,必須重新計算物面邊界點和背景點的相關系數Dij,再由式(12)得到Aij。Dij的列向量代表了一個物面拉格朗日節點(j＝1,…,N;N為物面點總數)與每個Euler背景網格點(i＝1,…,M;M為背景網格點總數)的相關關系。每個物面拉格朗日點對應的背景歐拉網格點數相對于總網格點數極為有限,這意味著D和A為稀疏矩陣,背景歐拉相關點必須在以該物面點為球心的周圍半徑為R的空間內。因此,在A生成計算中可忽略兩個間距大于2R的物面點之間的相關系數計算。通過添加這一篩選條件,可大大減少A矩陣重構的計算時間。

(2)利用A矩陣的稀疏對稱性使用一維存儲。由于D為M×N的實系數稀疏矩陣,且A可以寫成D·DT的形式,故由矩陣性質可知A亦為正定矩陣。在A生成過程中可掃描其下三角部分,以一維向量形式存儲其下三角非零元素的值及位置信息(圖2)。應用一維存儲可以極大減小稀疏矩陣的內存占用和相關計算時間。

圖2 一維存儲：按行搜索稀疏對稱矩陣A的下三角部分Fig.2 One dimensional storage：search the lower triangle part of the coarse matrix A along each row

(3)利用矩陣A的對稱正定性,使用共軛梯度法對以A為系數矩陣的線性方程組進行迭代求解。對正定對稱矩陣,相對于LU分解求逆等直接解法,共軛梯度法作為一種快速迭代法可以極大加速線性方程組AX＝B的求解。

由于A矩陣的規模取決于所要計算的物面拉格朗日節點數N,故當物面網格規模越大,所使用的加速效率就相對越高。

對典型的動邊界非定常物面繞流問題,在其迭代周期中,單一進程中更為主要的計算時間包括界面單元通量重構計算。該部分的計算量與當前計算進程中處理的背景歐拉網格量成正相關關系。因此,為了實現工程實用的數值模擬能力,有必要將背景流場計算的笛卡爾網格塊劃分為若干子塊進行并行計算。

對單進程IB-LBFS程序,其流場Euler笛卡爾網格由1個塊(Block)構成,構成物面的Lagrange邊界面位于該塊的中心部位(圖3)。為了減小物面邊界和背景網格的插值搜索計算量,該單一塊被分為兩個區域:位于物體周圍的內部區域和位于遠場的外部區域。在內部區域,其Euler網格為均勻網格,在外部區域,Euler網格分布則由于遠場條件逐漸變稀疏。

圖3 歐拉背景網格塊分區：單塊單進程Fig.3 Euler grid domain：single block for serial computation

圖4 歐拉背景網格塊分區：多塊多進程Fig.4 Euler grid domain：multi-block for parallel computation

參考單進程情況下的網格塊結構特性,在并行計算中,對背景網格進行了如圖4所示的分塊處理:原內部均勻區的對應區域被分割為內部子塊(Inner region grid blocks),原外部非均勻區域的對應區域被分割為外部子塊(Outer region grid blocks)。每個子塊對應一個并行進程進行塊內部流場方程的求解。在部分內部子塊中,需要額外進行浸潤邊界的搜索插值和物面邊界修正速度計算。塊與塊之間的內部邊界需要進行必要的數據傳輸,鄰近外邊界的子塊的外邊界面需要給定相應的邊界條件。

在進行修正速度方程組求解時,由于矩陣A具有全局性質,故每一個相關內部子塊進程中均需求解以該矩陣為系數的線性方程組。在計算AX＝B中,A具有全局變量的特性,同時右端項B亦須由各邊界插值進程中的值疊加得到(分裂求解方程組,再進行解X的疊加會導致數值誤差)。因此,限于修正速度的基本求解算法,并行計算的并行效率由于全局變量和增加的數據交換量會受到一定的影響,為此需要在編程時盡量減少涉及的進程數和數據傳遞量。

3 絕對節點坐標板單元

在傳統有限元方法中,對柔性體的描述基于微小變形和轉動量作為其廣義坐標,這一形式難以描述其剛體運動模態,同時對柔性體大變形則需要大量的計算單元求解。絕對節點坐標法(Absolute Nodal Coordinate Formula,ANCF)則基于非增量的在全局坐標系中描述的廣義坐標求解,這一方法可以準確描述柔性體的剛性運動,同時保證了質量陣為常數,可以以較少單元描述柔性體的大變形運動。

對矩形板結構進行了基于ANCF的4節點48自由度薄板單元建模[12-14](見圖5)。在這一模型中假定該板具有如下性質:

(1)板的相對厚度很小,在厚度方向上應力/應變均勻;

(2)初始形狀為矩形。

如圖5所示,該單元節點為4個角點,每個節點包括12個廣義坐標自由度。第i個節點的廣義坐標向量ei包括節點的絕對位置坐標r、節點在兩個物質坐標方向(l1,l2)上的切向梯度?r/?p1、?r/?p2,以及交叉梯度?2r/(?p1?p2)。

板中任一點的坐標由單元插值函數決定:r＝S·e,其中S為二維Hermite廣義坐標插值函數,由兩組一維Hermite插值函數相乘得到。

圖5 48自由度絕對節點坐標薄板單元Fig.5 48 degree of freedom ANCF thin plate element

根據Kirchhoff定律,薄板彈性能可分解為兩部分:一部分為薄板中面產生的拉伸應變能和剪切彈性能,另一部分為板的彎度引起的彎曲彈性能和扭轉彈性能。薄板單元的彈性能U和單元廣義彈性力Qe的表達式為:

其中拉伸和剪切應變ε和側向曲率κ的分量表達式為:

基于絕對節點坐標單元的結構力學求解器合適于與多體動力學求解框架相結合。綜合第2節和第3節理論和算法,搭建了基于IB-LBFS流場求解器和絕對節點坐標單元柔性結構動力學的流固耦合數值計算程序,其計算流程圖和模塊功能圖如圖6所示。

圖6 流固耦合求解器計算流程圖Fig.6 Flowchart of FSI solver

4 數值驗證與分析

4.1 非定常動邊界驗證算例：橫向旋轉球的繞流

旋轉球體的繞流問題是一類典型的動邊界非定常繞流問題,其流動特性由物體外形和運動的邊界條件共同決定。進行橫向旋轉(旋轉角速度與來流方向垂直)的球體繞流作為一個典型問題,可以有效檢驗IB-LBFS求解器的加速優化方法和并行算法。如圖7所示,本算例采用文獻[9]中的參數:來流雷諾數Re＝300,無量綱來流速度U∞＝0.1(該變量與馬赫數的關系有Ma＝U∞)。旋轉無量綱角速度定義為ω－＝ωyD/2U∞的 取值為0.1、0.3、0.5、0.6、0.8、1.0。使用并行IB-LBFS計算程序進行定物面/動態物面計算,并行核數為1152個(其中均勻區內部子塊448個),背景歐拉網格單元數為3.11×107。流場計算域大小為(240D,80D,80D)。

圖8至圖11為ω＝0.1、0.3、0.5、1.0四個角速度下的繞流渦量等值面圖,顯示了相應的脫落渦序列結構。隨著球旋轉角速度的增加,脫落渦的間隔越小,且其脫落軌跡逐漸沿著球面后部切向旋轉速度(＋z方向)方向偏移。圖12和圖13分別給出了球體所受的呈周期性振蕩的阻力系數CD(＋x方向)和側向力系數Cz(＋z)的時均值。計算結果與文獻[9]中的值相符,驗證了流場計算的準確性。

圖14給出了不同角速度下的氣動力St數,可視為渦脫落的無量綱頻率。當ω≤0.3時,St數與角速度ω成正比;0.3≤ω≤0.5為過渡區間;當ω＞0.5后,St數與角速度ω重新具有正比線性關系。

圖7 繞y軸(垂直來流方向)轉動的球體繞流問題Fig.7 Flow around rotating sphere(ωy in y direction)

圖8 無量綱角速度ω＝0.1：繞流渦結構Fig.8 Angular velocityω＝0.1：vortex structure

圖10 無量綱角速度ω＝0.5：繞流渦結構Fig.10 Angular velocityω＝0.5：vortex structure

圖11 無量綱角速度ω＝1.0：繞流渦結構Fig.11 Angular velocityω＝1.0：vortex structure

圖15顯示了本算例中,應用本文第2節中的各種IB-LBFS加速計算方法在其對應計算步驟中所獲得的加速比。其中,使用矩陣生成加速技術可使A矩陣的生成時間減少至直接生成法所需時間的0.05倍;使用共軛梯度法求解線性方程組問題可使計算時間減少為直接解法的0.03倍;在求解方法相同時,應用A矩陣的一維存儲替代其全元素二維存儲可減少1/3的計算時間。

圖13 球體時均側向力系數C Z(對比文獻[8])Fig.13 Lateral force coefficient C Z of rotating sphere

圖14 不同旋轉角速度下的斯特勞哈爾數StFig.14 St Number with different rotating velocities

圖15 加速運算方法加速比Fig.15 Accelerating ratio for the optimization method

表1給出了在本算例進行定常初始場計算中,對相同的整體背景歐拉網格,采用不同的并行進程計算得到的單步迭代時間和相對并行效率。進程數從256增長至576時,其并行效率保持不變;并行進程數進一步增加至1152后,迭代時間進一步減小,但由于單核網格量較少(2.7萬個),數據交換時間占比增加,其并行效率有一定下降。

表1 定物面計算中不同并行進程數下的并行效率Table 1 Parallel efficiency for different parallel process number without body surface updating

表2給出了在1152核并行球體旋轉計算中,增加了物面邊界-背景笛卡爾網格重構計算的非定常迭代步中各部分計算的時間占比。其中物面邊界-背景網格重構計算時間在優化后占比仍在60%以上,這說明了本文相關優化算法的必要性。

表2 非定常并行計算中單步迭代各計算步驟所占時間比例Table 2 Time fraction of computation steps in unsteady surface IB-LBFS parallel computation

4.2 流固耦合驗證算例：矩形旗幟擺動

本節通過研究柔性物面在流場中的變形,考核了基于IB-LBFS和絕對節點坐標法的的流固耦合求解器。矩形柔性旗幟在均勻來流中的擺動運動,是一個典型的三維空間中柔性體的流固耦合問題[15-16]。來流雷諾數為Re＝200,在初始時刻旗幟與均勻來流方向存在夾角α＝0.1π(圖16),厚度為d h＝0.01L,旗幟材料的無量綱彎曲剛度為E－＝1×10-4。柔性旗幟結構使用48自由度絕對結點坐標板單元進行離散。并行分塊進程數為N＝80。

使用中等規模網格對該問題進行并行計算仿真。流場計算域空間范圍為x∈[－12,48],y∈[－15,15],z∈[－15,15],總網格量為240×192×192;均勻區網格區單元尺寸為d h＝0.02,均勻區部分網格量為144×96×96;對矩形旗幟使用41×41個邊界Lagrange點離散,物面網格邊長為d sx＝d sy＝0.025。

圖16 初始狀態：均勻來流中的斜置板Fig.16 Initial configuration：inclined plate in uniform flow

圖17、圖18分別顯示了旗幟擺動中的中間截面z＝0和展向一側橫截面z＝－0.5的流線圖?？梢钥闯?不同于二維索流致擺動算例(可視為展向無限長的旗幟擺動問題),三維有限展向的旗幟擺動的流場顯示出明顯的展向變化和三維效應。

圖17 z＝－0.5截面：密度云圖與流線(t＝50.0)Fig.17 Density contour and streamline at z＝－0.5

圖18 z＝0截面：密度云圖與流線(t＝50.0)Fig.18 Density contour and streamline at z＝0(t＝50.0)

圖19顯示了一個周期(t＝17.36~21.36)之內,旗幟的四個擺動瞬時構型。t＝17.36時,旗幟自由邊擺動至y方向極大值,此時旗幟的整體速度達到極小值,并受到旗幟兩側壓差產生的－y方向側向力,該力使旗幟具有向y負方向擺動的趨勢;t＝18.56時,旗幟處于y中間位置,y方向整體速度絕對值達到極大值,同時由于旗幟構型發生改變,兩側壓差力較之前開始反向,變為＋y方向;在t＝19.36時,旗幟整體－y方向速度在＋y壓差力作用下減至接近0,該方向壓差力在該構型下達到極大;t＝20.56時,旗幟整體達到＋y方向的極速度,繼續運動至圖19(a)構型完成一個完整運動周期。旗幟形成穩定擺動后,其后緣中點B和角點A(見圖16)的y-x方向位移形成了如圖20所示的極限環結構。旗幟擺動的無量綱頻率為St＝Lf/U＝0.252,該值與文獻值(Huang[16],St＝0.26;Tian[15],St＝0.263)相差4%以內。

圖19 1個周期內旗幟擺動的動態構型Fig.19 Flag configuration in a time period

圖20 旗幟后緣點A/B的xy方向振動極限環Fig.20 Limit cycle of point A/B at the afteredge of flag

圖21給出了在相同的背景網格中,不同結構離散單元數得到的后緣點的y方向位移-時間曲線。對5×5 ANCF單元和10×10 ANCF單元描述的旗幟結構,二者結果吻合,說明結構單元離散數達到了網格收斂。

圖22給出了動態物面計算中,串行IB-LBFS計算與并行計算得到的的旗幟阻力系數-時間曲線。兩條曲線重合,表明IB-LBFS并行化程序與串行版本具有良好的一致性。

圖22 旗幟阻力系數-時間曲線：串行與并行Fig.2 Drag coefficient-time curve：comparison between serial and parallel computation

4.3 流固耦合驗證算例：固支板的彎曲變形

一端固支于地面的彈性板在與其垂直方向的均勻來流中,流體載荷將使其產生彎曲變形,并在板的周圍和流動后方形成繞流結構。這一模型與Luhar和Nepf[17]在研究藻類植物在水中的姿態構型中所建立的模型相似。該流固耦合系統示意圖如圖23所示。

來流參數為:均勻來流無量綱速度U0＝0.1,雷諾數Re＝1600。板初始狀態的無量綱幾何參數為:寬度為b＝1,高度L＝5b,厚度h＝0.2b;板密度為ρs,固體與流體的密度比為ρ?＝ρs/ρf＝0.678,無量綱楊氏模量為E?＝＝19 054.9,泊松比νs＝0.4,板受到的浮力方向沿z正方向,其無量綱值為＝0.036 975。計算并行進程數為128,網格均勻區單元尺寸為d h＝0.04;全局計算域尺寸為x∈[ －20,120],y∈[ －40,40],z∈[ －40,40],總網格規模為400×250×200;板表面使用11×51個Lagrange節點模擬,考慮到該板具有一定的薄板特性,在物面邊界中忽略板的厚度。在結構離散中,使用2×10個絕對節點坐標板單元對該彈性板進行描述。

圖23 三維來流中的固支板Fig.23 Clamped plate in three dimensional uniform flow

給定板受到系數為C＝50的結構Rayleigh阻尼力,彎板在均勻來流流場作用下開始逐漸彎曲,并最終使板彈性力與流體載荷力達到平衡,最終收斂于某一穩定構型。圖24給出了該穩定構型的側向投影與Nepf[17]試驗圖像的比較。圖25給出了彎板頂端邊界角點(nd1,nd3)和中點(nd2)在x方向(來流方向)和y方向(浮力方向)的位移-時間曲線,t＞120 s后構型基本達到收斂位置。板端點在垂直方向的無量綱位移為 Δy/b＝0.63(文獻[17]中 Δy＝0.59),水平位移為Δx/b＝2.28(文獻[17]中Δx＝2.14),誤差分別為6.8%和6.5%。造成這一誤差的主要原因在于,本算例中未考慮板的實際厚度,在同樣的彈性模量下,所計算的彎曲剛度小于有限厚度情況彎曲剛度,這造成計算位移略微偏大。

圖26顯示了固支彈性板的彎曲穩定收斂構型下,其三維繞流的渦量等值面結構(|w|＝0.05)。在給定的雷諾數下,彎板后流動顯示出復雜的渦結構條帶系統,其三維效應和相應的分離流動結構值得進一步分析。

圖24 彎曲板穩定構型與Nepf[17]實驗結構對比Fig.24 Steady configuration of the clamped plate compared with the experiment result from Nepf[17]

圖25 柔性板上邊界端點和中點的y-x方向位移Fig.25 Upper boundary points of flexible plate y-x displacement

圖26 彎曲板穩定構型的彎板繞流渦量等值面(|w|＝0.05)Fig.26 Vortex isosurface of the steady configuration of the clamped plate(|w|＝0.05)

5 結 論

本文提出了浸潤邊界－格子波爾茲曼通量求解器加速優化算法,并發展了IB-LBFS的并行算法,構建了基于絕對節點坐標板單元的柔性結構動力學求解器,結合IB-LBFS搭建了流固耦合計算平臺,并使用多個算例驗證了基于IB-LBFS方法的大變形柔性結構流固耦合求解器的并行仿真能力。結果表明:

(1)本文提出將絕對基于絕對結點坐標法(ANCF)的柔性多體動力學框架與格子玻爾茲曼求解器相結合,構成一種新型的流固耦合數值模擬方法。該方法可顯著提高定物面/動態物面迭代計算效率,進一步提高了IB-LBFS解決較大規模網格問題的能力,特別適用于復雜約束柔性體系統的大變形流固耦合數值模擬與仿真。

(2)在當前的流固耦合動態計算中,仍然采用時間松耦合迭代法,為了增大計算時間步長并減小總計算量,可進一步發展柔性體流固耦合隱式迭代算法。并可通過構建更多的絕對節點坐標單元和其他結構單元,實現復雜的工程結構的柔性體流固耦合計算。