劉 釩,劉 剛,江 雄,舒 昌

(1.中國(guó)空氣動(dòng)力研究與發(fā)展中心,綿陽(yáng) 621000;2.新加坡國(guó)立大學(xué) 工程學(xué)院,新加坡 117576)

0 引 言

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

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

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

1 浸潤(rùn)邊界－格子波爾茲曼通量求解器

1.1 控制方程

宏觀流動(dòng)的Navier-Stokes方程為:

在格子波爾茲曼通量計(jì)算中,動(dòng)量通量ρu和黏性-無黏通量Π的表達(dá)式為:

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

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

其中F為浸潤(rùn)邊界法中物面邊界對(duì)應(yīng)的體積力項(xiàng)。

1.2 離散求解

將式(5)進(jìn)行時(shí)間離散,將求解過程分裂為一個(gè)預(yù)測(cè)步和一個(gè)修正步(分別對(duì)應(yīng)式(6)和式(7)):

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 絕對(duì)節(jié)點(diǎn)坐標(biāo)板單元

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

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

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

(2)初始形狀為矩形。

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

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

圖5 48自由度絕對(duì)節(jié)點(diǎn)坐標(biāo)薄板單元Fig.5 48 degree of freedom ANCF thin plate element

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

其中拉伸和剪切應(yīng)變?chǔ)藕蛡?cè)向曲率κ的分量表達(dá)式為:

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

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

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

4.1 非定常動(dòng)邊界驗(yàn)證算例：橫向旋轉(zhuǎn)球的繞流

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

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

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

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

圖8 無量綱角速度ω＝0.1：繞流渦結(jié)構(gòu)Fig.8 Angular velocityω＝0.1：vortex structure

圖10 無量綱角速度ω＝0.5：繞流渦結(jié)構(gòu)Fig.10 Angular velocityω＝0.5：vortex structure

圖11 無量綱角速度ω＝1.0：繞流渦結(jié)構(gòu)Fig.11 Angular velocityω＝1.0：vortex structure

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

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

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

圖15 加速運(yùn)算方法加速比Fig.15 Accelerating ratio for the optimization method

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

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

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

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

4.2 流固耦合驗(yàn)證算例：矩形旗幟擺動(dòng)

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

使用中等規(guī)模網(wǎng)格對(duì)該問題進(jìn)行并行計(jì)算仿真。流場(chǎng)計(jì)算域空間范圍為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;對(duì)矩形旗幟使用41×41個(gè)邊界Lagrange點(diǎn)離散,物面網(wǎng)格邊長(zhǎng)為d sx＝d sy＝0.025。

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

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

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

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

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

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

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

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

4.3 流固耦合驗(yàn)證算例：固支板的彎曲變形

一端固支于地面的彈性板在與其垂直方向的均勻來流中,流體載荷將使其產(chǎn)生彎曲變形,并在板的周圍和流動(dòng)后方形成繞流結(jié)構(gòu)。這一模型與Luhar和Nepf[17]在研究藻類植物在水中的姿態(tài)構(gòu)型中所建立的模型相似。該流固耦合系統(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。計(jì)算并行進(jìn)程數(shù)為128,網(wǎng)格均勻區(qū)單元尺寸為d h＝0.04;全局計(jì)算域尺寸為x∈[ －20,120],y∈[ －40,40],z∈[ －40,40],總網(wǎng)格規(guī)模為400×250×200;板表面使用11×51個(gè)Lagrange節(jié)點(diǎn)模擬,考慮到該板具有一定的薄板特性,在物面邊界中忽略板的厚度。在結(jié)構(gòu)離散中,使用2×10個(gè)絕對(duì)節(jié)點(diǎn)坐標(biāo)板單元對(duì)該彈性板進(jìn)行描述。

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

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

圖26顯示了固支彈性板的彎曲穩(wěn)定收斂構(gòu)型下,其三維繞流的渦量等值面結(jié)構(gòu)(|w|＝0.05)。在給定的雷諾數(shù)下,彎板后流動(dòng)顯示出復(fù)雜的渦結(jié)構(gòu)條帶系統(tǒng),其三維效應(yīng)和相應(yīng)的分離流動(dòng)結(jié)構(gòu)值得進(jìn)一步分析。

圖24 彎曲板穩(wěn)定構(gòu)型與Nepf[17]實(shí)驗(yàn)結(jié)構(gòu)對(duì)比Fig.24 Steady configuration of the clamped plate compared with the experiment result from Nepf[17]

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

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

5 結(jié) 論

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

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

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