一種求解低速跨流域流固耦合問題的離散統(tǒng)一氣體動理學(xué)格式

王 勇，劉 沙，2，卓叢山，2，鐘誠文，2，*

(1. 西北工業(yè)大學(xué) 航空學(xué)院，西安 710072；2. 西北工業(yè)大學(xué) 翼型、葉柵空氣動力學(xué)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

0 引 言

運(yùn)動邊界和流固耦合（Fluid-structure interaction，F(xiàn)SI）問題是航空航天工程、船舶和海洋工程、土木和交通工程等領(lǐng)域經(jīng)常面臨的多場耦合問題[1-2]。數(shù)值模擬FSI 問題目前已發(fā)展大量方法。近些年，隨著技術(shù)進(jìn)步，以微納米機(jī)電系統(tǒng)（micro-electro-mechanicalsystems，MEMS）和超低軌道航天器為代表的應(yīng)用領(lǐng)域面臨跨流域運(yùn)動邊界和流固耦合問題，如MEMS中微尺度梁的振動問題[3]、航天器返回過程[4]和超低軌道航天器與外層稀薄大氣相互作用問題[5]。在這類問題中，由于特征尺度的顯著減小或氣體密度顯著降低，流體的稀薄氣體效應(yīng)是這類流固耦合問題需要考慮的重要影響因素之一。根據(jù)克努森數(shù)（Kn）定義[6]，流動可劃分為連續(xù)流、滑移流、過渡流和自由分子流四個流域。由于上述新問題中流體流域通常處于滑移流域或過渡流域，而傳統(tǒng)的流固耦合求解方法中流體部分通常采用基于Navier-Stokes 方程的適用于連續(xù)流域的宏觀方法，因此有必要發(fā)展可考慮稀薄氣體效應(yīng)影響的流固耦合求解新方法。

近年來，在跨流域流動數(shù)值模擬方法方面，離散統(tǒng)一氣體動理學(xué)格式（discrete unified gas kinetic scheme，DUGKS）[7]這一新興數(shù)值方法得到了廣泛關(guān)注。該方法在格式通量計算部分，借鑒了格子Boltzmann 方法（lattice Boltzmann method，LBM）[8]的特征線構(gòu)造思想，較統(tǒng)一氣體動理學(xué)格式（unified gas kinetic scheme，UGKS）[9]顯著簡化，計算效率進(jìn)一步提升。在網(wǎng)格實(shí)施方面，由于標(biāo)準(zhǔn)LBM 采用遷移-碰撞模式，Boltzmann 方程的對流項(xiàng)隱藏在分布函數(shù)沿特征線的遷移過程中[10]，因此計算網(wǎng)格需采用均勻笛卡爾網(wǎng)格；而DUGKS 在實(shí)施過程中由于采用對流項(xiàng)顯式求解思想，盡管計算復(fù)雜度有所提升，數(shù)值耗散也有所增加，但可采用的網(wǎng)格也更加靈活，結(jié)構(gòu)、非結(jié)構(gòu)及混合網(wǎng)格均可，對不規(guī)則計算域的適用性顯著提高。在數(shù)值耗散方面，與傳統(tǒng)的有限體積LBM[11-12]相比，由于對流項(xiàng)計算過程中引入遷移-碰撞思想，因此數(shù)值耗散小，對網(wǎng)格尺度和時間步長約束小。目前該類方法已具有全流域流動計算能力，并在快速發(fā)展中[13-14]?？紤]到運(yùn)動邊界問題的重要性，目前已有開展少量工作[15]，仍需進(jìn)一步完善。

在運(yùn)動邊界數(shù)值模擬方面，網(wǎng)格技術(shù)通常分為基于動網(wǎng)格技術(shù)的任意拉格朗日-歐拉法（arbitrary Lagrangian-Eulerian，ALE）[16]和基于笛卡爾網(wǎng)格的浸入邊界法（immersed boundary method，IBM）[17]。ALE可采用貼體網(wǎng)格，在模擬高雷諾數(shù)流動及結(jié)構(gòu)運(yùn)動模式相對簡單的FSI 問題時有較大優(yōu)勢，計算效率較高，如航空工程中的顫振數(shù)值模擬。IBM 由于壁面附近采用笛卡爾網(wǎng)格，適用于雷諾數(shù)較低的流動問題，模擬大變形動邊界問題有顯著優(yōu)勢。目前，基于DUGKS的ALE[15]和IBM[18]均有發(fā)展，其中ALE-DUGKS 在模擬考慮稀薄氣體效應(yīng)影響的運(yùn)動邊界問題時有良好表現(xiàn)。而在低速連續(xù)流域范圍，以LBM、DUGKS為代表的介觀數(shù)值方法，由于無需求解宏觀方法中的壓力泊松方程[19]，方法實(shí)施也相對簡單。ALEDUGKS 在低速連續(xù)流和稀薄流的獨(dú)特優(yōu)勢是為進(jìn)一步發(fā)展低速跨流域運(yùn)動邊界和流固耦合數(shù)值方法奠定基礎(chǔ)。在流固耦合數(shù)值模擬方面，發(fā)展較為成熟的是松耦合計算方法[20]，即在每個數(shù)值模擬時間步長內(nèi)采用計算流體力學(xué)（computational fluid dynamics，CFD）和計算結(jié)構(gòu)動力學(xué)（computational structural dynamics，CSD）各自成熟的求解方法，并在流-固耦合界面進(jìn)行數(shù)據(jù)交換實(shí)現(xiàn)流體域和固體域的耦合推進(jìn)。

綜上，本文的主要工作是針對MEMS 和航天工程面臨的跨流域流動運(yùn)動邊界和流固耦合問題，采用新興跨流域數(shù)值方法DUGKS，將已發(fā)展的運(yùn)動邊界ALE-DUGKS 擴(kuò)展至低速，同時兼顧高速流動的流固耦合計算框架，利用網(wǎng)格的貼體性采用較少的網(wǎng)格量更高效地模擬邊界層流動。針對連續(xù)流域圓柱繞流強(qiáng)迫振動和自由振動兩個典型算例進(jìn)行模擬，驗(yàn)證方法的計算能力，并進(jìn)一步采用跨流域求解方法驗(yàn)證方法的跨流域計算能力和拓展性。

1 基于DUGKS 的流固耦合求解方法

1.1 ALE-DUGKS 格式

ALE-DUGKS 是王勇等[15]在原始DUGKS[7]基礎(chǔ)上構(gòu)造的跨流域流動運(yùn)動邊界求解格式。原始DUGKS 格式構(gòu)造基于Boltzmann-BGK 方程：

式中，下標(biāo)k為控制體界面編號，下標(biāo)b為界面中點(diǎn)處相應(yīng)變量值。式（8）和式（9）中的V和S分別為控

1.2 流固耦合計算方法

ALE-DUGKS 具備運(yùn)動邊界問題求解能力，耦合計算結(jié)構(gòu)動力學(xué)并采用松耦合求解方法后，可將DUGKS發(fā)展成可考慮稀薄氣體效應(yīng)的跨流域流固耦合計算框架。在實(shí)施過程中，流體域網(wǎng)格變形采用Laplace光順法[23]，結(jié)構(gòu)域控制方程為：

2 算例驗(yàn)證

2.1 靜止圓柱繞流算例

為了驗(yàn)證DUGKS 方法并為后續(xù)算例提供對比數(shù)據(jù)，首先進(jìn)行靜止圓柱繞流模擬。離散速度模型采用D2Q9 模型。參考文獻(xiàn)[25，26]的算例設(shè)置，計算域見圖1，此處D為圓柱直徑，根據(jù)阻塞率5%布置上、下邊界遠(yuǎn)場距離。計算域左側(cè)邊界和上、下邊界設(shè)置為入口邊界條件，邊界外法線為nb；對于進(jìn)入計算域的分布函數(shù)方向eα·nb＜0，分布函數(shù)給定為平衡態(tài)；離開計算域的分布函數(shù)方向由內(nèi)場分布函數(shù)插值獲得。計算域右側(cè)邊界設(shè)置為出口邊界條件，分布函數(shù)由內(nèi)場插值獲得。圖2 為計算網(wǎng)格，總網(wǎng)格單元為37 506，其中壁面附近采用四邊形貼體網(wǎng)格，最小網(wǎng)格尺度為D/256，其余為三角形網(wǎng)格。初始化來流速度為u=U∞=0.05，v=0； 初始密度 ρ∞=1.0。雷諾數(shù)范圍為Re= 60～150。升力系數(shù)CL定義為：

圖1 圓柱繞流計算域和邊界條件Fig. 1 The configuration and boundary conditions for the flow around a circular cylinder

圖2 圓柱繞流計算網(wǎng)格Fig. 2 A mesh for the flow around a circular cylinder

式中f為自然渦脫落頻率。本文與Prasanth[25]和He[26]的數(shù)值模擬結(jié)果以及Williamson[27]的平行渦脫落模型吻合良好。在雷諾數(shù)Re= 60～80 范圍內(nèi)數(shù)值結(jié)果略高于平行渦脫落模型，主要原因是隨著雷諾數(shù)降低，黏性影響范圍增大，遠(yuǎn)場邊界條件對結(jié)果影響顯著。減小阻塞率后，可明顯改善數(shù)值結(jié)果和模型之間的差異。

圖3 三組雷諾數(shù)下圓柱繞流升力系數(shù)時間發(fā)展歷程Fig. 3 Time evolutions of lift coefficients for flows with three Reynolds numbers around a circular cylinder

圖4 Re = 60～150 靜止圓柱繞流St 對比Fig. 4 A comparison of St for flows around a stationary cylinder at Re = 60～150

2.2 Re = 100 圓柱橫向強(qiáng)迫振動算例

在靜止圓柱繞流算例基礎(chǔ)上，進(jìn)一步模擬圓柱橫向強(qiáng)迫振動算例以驗(yàn)證算法的運(yùn)動邊界模擬能力。計算參數(shù)與算例1 相同，雷諾數(shù)Re= 100。圓柱y方向運(yùn)動方程為：

式中，A為圓柱振幅，fy為圓柱結(jié)構(gòu)頻率，fRe=100為靜止圓柱繞流自然渦脫落頻率，A*為圓柱無量綱振幅，f*為與自然渦脫落頻率的比值。在文獻(xiàn)[28]中，Kumar 等將圓柱渦脫落分為非鎖頻區(qū)、過渡區(qū)和鎖頻區(qū)。圖5 為各區(qū)分界線圖，與Kumar 等結(jié)果整體十分吻合。圖6 為五組A*和f*組合下的升阻力系數(shù)時間發(fā)展歷程和升力系數(shù)功率譜圖，其中阻力系數(shù)CD為將式（20）中FL替 換成FD。當(dāng)無量綱頻率f*＜1 時，非鎖頻區(qū)（圖6（a））和鎖頻區(qū)（圖6（b））有清晰的分界線。在非鎖頻區(qū)，由于流體自然渦脫落頻率和結(jié)構(gòu)頻率接近，升阻力系數(shù)出現(xiàn)類似“拍”現(xiàn)象，功率譜圖中流動自然渦脫落頻率占優(yōu)且包含其他無規(guī)律頻率信息；在鎖頻區(qū)，渦脫落頻率鎖定在結(jié)構(gòu)頻率，同時其他頻率與主頻呈現(xiàn)清晰的倍頻現(xiàn)象。當(dāng)無量綱頻率f*＞ 1 時，非鎖頻區(qū)和鎖頻區(qū)之間存在過渡區(qū)。在非鎖頻區(qū)和過渡區(qū)，升阻力系數(shù)同樣出現(xiàn)“拍”現(xiàn)象（圖6（c）和圖6（d））。當(dāng)流動自然渦脫落頻率占優(yōu)時，流動呈現(xiàn)非鎖頻現(xiàn)象；當(dāng)結(jié)構(gòu)自然頻率占優(yōu)時，流動則處在過渡區(qū)范圍。對于功率譜圖，非鎖頻區(qū)和過渡區(qū)除主頻外，還包含其他無規(guī)律頻率信息，鎖頻區(qū)其他頻率信息則呈現(xiàn)清晰的倍頻現(xiàn)象（圖6（e））。

圖5 Re = 100 圓柱橫向強(qiáng)迫振動非鎖頻區(qū)、過渡區(qū)和鎖頻區(qū)Fig. 5 No lock-in, transitional, and lock-in regimes for flows around a forced-oscillating circular cylinder at Re = 100

圖6 五組不同振幅A＊和頻率f ＊下升阻力系數(shù)時間發(fā)展歷程和升力系數(shù)功率譜對比Fig. 6 Time evolutions of lift and drag coefficients, and the power spectra of lift coefficient at five different amplitudes and frequencies of oscillation

2.3 Re = 100 圓柱自由振動算例

為驗(yàn)證ALE-DUGKS 模擬流固耦合問題的能力，首先對Re= 100 圓柱自由振動進(jìn)行模擬。結(jié)構(gòu)自由振動控制方程為x和y方向兩自由度無阻尼振動方程，對應(yīng)于式（18）中的質(zhì)量矩陣M、剛度矩陣K和外力矩陣Fext分別為：時間內(nèi)，圓柱處在鎖定狀態(tài)，流場發(fā)展成周期性流動后釋放圓柱。圖7 為U*= 6.25 時兩個方向位移時間發(fā)展歷程和收斂后的耦合位移，其中x和y為圓柱在兩個方向的位移值。可以發(fā)現(xiàn)周期性非定常升力使得圓柱在y方向產(chǎn)生較大位移，非定常阻力則使圓柱在x方向運(yùn)動到平衡位置后以較小振幅振動，且兩方向耦合位移呈現(xiàn)類似“8”的結(jié)構(gòu)；圓柱在其他折減速度時有類似現(xiàn)象。圖8 為圓柱y方向最大位移和渦致振動頻率對比，其中ymax/D為無量綱最大位移，St= 0.166 為靜止圓柱繞流的無量綱渦脫落頻率。本文與Singh[29]、Zhang[30]等數(shù)值結(jié)果吻合良好，主要差異在鎖頻和非鎖頻的過渡區(qū)間。圖9 為U*= 4.0 時y方向位移時間歷程，與圖7 結(jié)果不同，此時流動處在非鎖頻狀態(tài)，圓柱位移很小。圖10 則為流動處在鎖頻和非鎖頻的過渡狀態(tài)，需要較長的數(shù)值模擬時間獲得周期性振動結(jié)果。此外，數(shù)值模擬發(fā)現(xiàn)，阻塞率較小時，圓柱最大位移偏高，即遠(yuǎn)場距離對圓柱位移有壓縮效應(yīng)。這也是本文與文獻(xiàn)[29]相同，取5%阻塞率的原因。

圖7 圓柱自由振動折減速度U＊ = 6.25 結(jié)構(gòu)位移Fig. 7 Displacements of the free-oscillating circular cylinder at U＊ = 6.25

圖8 Re = 100 圓柱自由振動y 方向最大位移和渦脫落頻率對比Fig. 8 Comparisons of the maximum y-direction displacement and vortex shedding frequency for the flow around a freeoscillating circular cylinder at Re = 100

圖9 圓柱自由振動U＊ = 4.0 時y 方向位移時間發(fā)展歷程Fig. 9 Time evolution of the y-direction displacement of a freeoscillating circular cylinder at U＊ = 4.0

圖10 圓柱自由振動U＊ = 7.8 時y 方向位移時間發(fā)展歷程Fig. 10 Time evolution of the y-direction displacement of a freeoscillating circular cylinder at U＊ = 7.8

2.4 Re = 60～150 圓柱自由振動算例

與Re= 100 圓柱自由振動算例不同，該算例主要驗(yàn)證結(jié)構(gòu)自然頻率固定時來流風(fēng)速變化對結(jié)構(gòu)位移和頻率的影響。結(jié)構(gòu)參數(shù)仍為式（23）～式（25），此時無量綱頻率給定為Fs= 16.6/Re，其中16.6 使得Re=100 時結(jié)構(gòu)無量綱頻率約等于靜止圓柱繞流無量綱渦脫落頻率。圖11 為y方向最大位移和渦脫落頻率對比，與Prasanth 等[25]的數(shù)值結(jié)果十分吻合，僅在鎖頻和非鎖頻的過渡區(qū)間有差異。主要原因是Prasanth 等在計算過程中通過動態(tài)調(diào)節(jié)雷諾數(shù)（增加或減?。┮苑治鲲L(fēng)速的動態(tài)變化對過渡區(qū)流動機(jī)理的影響；本文主要為方法驗(yàn)證，忽略動態(tài)影響機(jī)理分析，計算中每給定一個雷諾數(shù)獲得收斂結(jié)果后完成計算。對于y方向位移和升阻力系數(shù)，在非鎖頻和鎖頻區(qū)與算例2.3 的Re= 100 自由振動算例類似，但在過渡區(qū)較為復(fù)雜。圖12 為Re= 130 的時間歷程，經(jīng)過頻譜分析，當(dāng)無量綱時間小于1 800 時，y方向位移不斷增加，此時渦脫落頻率中圓柱的自然渦脫落頻率占優(yōu)；當(dāng)時間大于1 800 后，渦脫落頻率鎖定在結(jié)構(gòu)自然頻率，圓柱位移也收斂到穩(wěn)定的周期性振動狀態(tài)。

圖11 Re = 60～150 圓柱自由振動y 方向最大位移和渦脫落頻率對比Fig. 11 Comparisons of the maximum y-direction displacement and vortex shedding frequency for flows around a free-oscillating circular cylinder at Re = 60～150

圖12 Re = 130 圓柱自由振動y 方向位移和升阻力系數(shù)時間發(fā)展歷程Fig. 12 Time evolutions of the y-direction displacement, and lift and drag coefficients for the flow around a free-oscillating circular cylinder at Re = 130

2.5 Re = 60 圓柱橫向強(qiáng)迫振動連續(xù)流和稀薄流對比

為驗(yàn)證算法的擴(kuò)展能力，采用跨流域計算框架對Re= 60 圓柱橫向強(qiáng)迫振動進(jìn)行數(shù)值模擬?？缌饔蛴嬎阒?，克努森數(shù)Kn、馬赫數(shù)Ma和雷諾數(shù)Re關(guān)系為[31]：

圖13 Ma = 0.086 6、Re = 60 圓柱強(qiáng)迫振動升力系數(shù)時間發(fā)展歷程Fig. 13 Time evolutions of lift coefficients for the flow around a forced-oscillating circular cylinder at Re = 60 and Ma = 0.086 6

針對弱可壓縮流動，式（7）和式（27）求解見文獻(xiàn)[15]。離散速度采用28×28 高斯點(diǎn)。Ma= 0.4 時圓柱自然渦脫落頻率St= 0.132，與Canuto 等[32]結(jié)果吻合，略低于Ma= 0.086 6 時的結(jié)果。圖14 為兩組馬赫數(shù)下升力系數(shù)時間發(fā)展歷程對比，隨著馬赫數(shù)增加，稀薄程度增加，升力系數(shù)呈減小趨勢。Canuto 等[32]采用宏觀方法和無滑移邊界分析了壓縮性對靜止圓柱繞流的流動特性影響，因方法的局限性無法考慮稀薄氣體效應(yīng)影響，數(shù)值結(jié)果在低雷諾數(shù)時存在一定偏差。但與本文靜止圓柱和振動圓柱趨勢一致，即Re=60 時馬赫數(shù)對St影響較小，但最大升力系數(shù)降低。

圖14 兩組馬赫數(shù)下圓柱強(qiáng)迫振動升力系數(shù)時間發(fā)展歷程Fig. 14 Time evolutions of lift coefficients for flows around a forced-oscillating circular cylinder at two different Ma

2.6 Ma = 5.0 圓柱橫向強(qiáng)迫振動稀薄流算例

2.2～2.5 節(jié)的算例驗(yàn)證了算法的低速跨流域運(yùn)動

圖15 Ma = 5.0 圓柱強(qiáng)迫振動物理網(wǎng)格和速度網(wǎng)格Fig. 15 The meshes in the (a) physical- and (b) velocityspace[21]for the flow around a forced-oscillating circular cylinder at Ma = 5.0

其中，A*為圓柱無量綱振幅，U*= 100 為折減速度。圖17 為A*= 0.05 時的升力系數(shù)時間發(fā)展歷程。由于振幅較小，振動頻率較低，圓柱所受升力系數(shù)較小。此外，與靜止圓柱繞流相同，馬赫數(shù)固定時，來流越稀薄，升力越大。圖18 為Kn= 0.1、A*分 別取0.05 和0.1 時的升力系數(shù)時間發(fā)展歷程對比。由于來流為稀薄流，圓柱表面無分離渦和渦脫落現(xiàn)象，流動呈現(xiàn)為附著流的未失穩(wěn)狀態(tài)，此時升力系數(shù)與振幅表現(xiàn)為線性增長關(guān)系。考慮到不同馬赫數(shù)、克努森數(shù)、振幅A*和折減速度U*組合下的圓柱強(qiáng)迫振動流動分析超出了本文的討論范圍，課題組將會在接下來的論文中進(jìn)行系統(tǒng)討論。

圖17 兩組努森數(shù)下圓柱強(qiáng)迫振動升力系數(shù)時間發(fā)展歷程Fig. 17 Time evolutions of lift coefficients for flows around a forced-oscillating circular cylinder at two different Kn numbers

圖 18 Kn = 0.1 時兩組無量綱振幅A＊的升力系數(shù)時間發(fā)展歷程Fig. 18 Time evolutions of lift coefficients for flows around a forced-oscillating circular cylinder at two different A＊ values

3 結(jié) 論

本文在原始DUGKS 基礎(chǔ)上，耦合任意拉格朗日-歐拉方法和計算結(jié)構(gòu)動力學(xué)，實(shí)現(xiàn)了可考慮稀薄氣體效應(yīng)的低速兼顧高速流動的流固耦合問題數(shù)值模擬框架。針對低速連續(xù)流域的圓柱強(qiáng)迫振動和自由振動兩個典型算例，采用LBM 中高效的D2Q9 離散速度模型，該框架可獲得與宏觀方法一致的數(shù)值結(jié)果。鑒于DUGKS 的全流域、全速域計算能力，更改離散速度模型后，即可具備模擬跨流域運(yùn)動邊界模擬能力。在算法實(shí)施方面，本文初步驗(yàn)證了顯式ALEDUGKS 模擬流固耦合問題能力，未來可發(fā)展隱式格式，進(jìn)一步提高其計算效率。在應(yīng)用方面，考慮到稀薄氣體效應(yīng)對新型跨流域飛行器低雷諾數(shù)可壓縮繞流場的顯著影響[35]，以及國內(nèi)外較少開展低雷諾數(shù)可壓縮圓柱繞流分析，本文后續(xù)將深入研究圓柱靜、動態(tài)氣動特性數(shù)值模擬。