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

劉 釩,劉 剛,江 雄,舒 昌

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

0 引 言

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

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

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

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

1.1 控制方程

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

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

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

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

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

1.2 離散求解

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

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

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

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

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

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

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

2 IB-LBFS：加速優(yōu)化計算和并行算法

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

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

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

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

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

由于A矩陣的規(guī)模取決于所要計算的物面拉格朗日節(jié)點數(shù)N,故當物面網(wǎng)格規(guī)模越大,所使用的加速效率就相對越高。

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

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

圖3 歐拉背景網(wǎng)格塊分區(qū)：單塊單進程Fig.3 Euler grid domain：single block for serial computation

圖4 歐拉背景網(wǎng)格塊分區(qū)：多塊多進程Fig.4 Euler grid domain：multi-block for parallel computation

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

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

3 絕對節(jié)點坐標板單元

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

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

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

(2)初始形狀為矩形。

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

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

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

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

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

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

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

4 數(shù)值驗證與分析

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

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

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

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

圖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節(jié)中的各種IB-LBFS加速計算方法在其對應計算步驟中所獲得的加速比。其中,使用矩陣生成加速技術可使A矩陣的生成時間減少至直接生成法所需時間的0.05倍;使用共軛梯度法求解線性方程組問題可使計算時間減少為直接解法的0.03倍;在求解方法相同時,應用A矩陣的一維存儲替代其全元素二維存儲可減少1/3的計算時間。

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

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

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

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

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

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

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

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

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

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

圖16 初始狀態(tài)：均勻來流中的斜置板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)之內(nèi),旗幟的四個擺動瞬時構型。t＝17.36時,旗幟自由邊擺動至y方向極大值,此時旗幟的整體速度達到極小值,并受到旗幟兩側壓差產(chǎn)生的－y方向側向力,該力使旗幟具有向y負方向擺動的趨勢;t＝18.56時,旗幟處于y中間位置,y方向整體速度絕對值達到極大值,同時由于旗幟構型發(fā)生改變,兩側壓差力較之前開始反向,變?yōu)椋珁方向;在t＝19.36時,旗幟整體－y方向速度在＋y壓差力作用下減至接近0,該方向壓差力在該構型下達到極大;t＝20.56時,旗幟整體達到＋y方向的極速度,繼續(xù)運動至圖19(a)構型完成一個完整運動周期。旗幟形成穩(wěn)定擺動后,其后緣中點B和角點A(見圖16)的y-x方向位移形成了如圖20所示的極限環(huán)結構。旗幟擺動的無量綱頻率為St＝Lf/U＝0.252,該值與文獻值(Huang[16],St＝0.26;Tian[15],St＝0.263)相差4%以內(nèi)。

圖19 1個周期內(nèi)旗幟擺動的動態(tài)構型Fig.19 Flag configuration in a time period

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

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

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

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

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

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

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

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

給定板受到系數(shù)為C＝50的結構Rayleigh阻尼力,彎板在均勻來流流場作用下開始逐漸彎曲,并最終使板彈性力與流體載荷力達到平衡,最終收斂于某一穩(wěn)定構型。圖24給出了該穩(wěn)定構型的側向投影與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ěn)定收斂構型下,其三維繞流的渦量等值面結構(|w|＝0.05)。在給定的雷諾數(shù)下,彎板后流動顯示出復雜的渦結構條帶系統(tǒng),其三維效應和相應的分離流動結構值得進一步分析。

圖24 彎曲板穩(wěn)定構型與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ěn)定構型的彎板繞流渦量等值面(|w|＝0.05)Fig.26 Vortex isosurface of the steady configuration of the clamped plate(|w|＝0.05)

5 結 論

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

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

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