一種求解低速跨流域流固耦合問題的離散統一氣體動理學格式

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

(1. 西北工業大學 航空學院，西安 710072；2. 西北工業大學 翼型、葉柵空氣動力學重點實驗室，西安 710072)

0 引 言

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

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

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

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

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

1.1 ALE-DUGKS 格式

ALE-DUGKS 是王勇等[15]在原始DUGKS[7]基礎上構造的跨流域流動運動邊界求解格式。原始DUGKS 格式構造基于Boltzmann-BGK 方程：

式中，下標k為控制體界面編號，下標b為界面中點處相應變量值。式（8）和式（9）中的V和S分別為控

1.2 流固耦合計算方法

ALE-DUGKS 具備運動邊界問題求解能力，耦合計算結構動力學并采用松耦合求解方法后，可將DUGKS發展成可考慮稀薄氣體效應的跨流域流固耦合計算框架。在實施過程中，流體域網格變形采用Laplace光順法[23]，結構域控制方程為：

2 算例驗證

2.1 靜止圓柱繞流算例

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

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

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

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

圖3 三組雷諾數下圓柱繞流升力系數時間發展歷程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 圓柱橫向強迫振動算例

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

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

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

圖6 五組不同振幅A＊和頻率f ＊下升阻力系數時間發展歷程和升力系數功率譜對比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 圓柱自由振動算例

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

圖7 圓柱自由振動折減速度U＊ = 6.25 結構位移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 方向位移時間發展歷程Fig. 9 Time evolution of the y-direction displacement of a freeoscillating circular cylinder at U＊ = 4.0

圖10 圓柱自由振動U＊ = 7.8 時y 方向位移時間發展歷程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 圓柱自由振動算例不同，該算例主要驗證結構自然頻率固定時來流風速變化對結構位移和頻率的影響。結構參數仍為式（23）～式（25），此時無量綱頻率給定為Fs= 16.6/Re，其中16.6 使得Re=100 時結構無量綱頻率約等于靜止圓柱繞流無量綱渦脫落頻率。圖11 為y方向最大位移和渦脫落頻率對比，與Prasanth 等[25]的數值結果十分吻合，僅在鎖頻和非鎖頻的過渡區間有差異。主要原因是Prasanth 等在計算過程中通過動態調節雷諾數（增加或減小）以分析風速的動態變化對過渡區流動機理的影響；本文主要為方法驗證，忽略動態影響機理分析，計算中每給定一個雷諾數獲得收斂結果后完成計算。對于y方向位移和升阻力系數，在非鎖頻和鎖頻區與算例2.3 的Re= 100 自由振動算例類似，但在過渡區較為復雜。圖12 為Re= 130 的時間歷程，經過頻譜分析，當無量綱時間小于1 800 時，y方向位移不斷增加，此時渦脫落頻率中圓柱的自然渦脫落頻率占優；當時間大于1 800 后，渦脫落頻率鎖定在結構自然頻率，圓柱位移也收斂到穩定的周期性振動狀態。

圖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 方向位移和升阻力系數時間發展歷程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 圓柱橫向強迫振動連續流和稀薄流對比

為驗證算法的擴展能力，采用跨流域計算框架對Re= 60 圓柱橫向強迫振動進行數值模擬?？缌饔蛴嬎阒?，克努森數Kn、馬赫數Ma和雷諾數Re關系為[31]：

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

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

2.6 Ma = 5.0 圓柱橫向強迫振動稀薄流算例

2.2～2.5 節的算例驗證了算法的低速跨流域運動

圖15 Ma = 5.0 圓柱強迫振動物理網格和速度網格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 時的升力系數時間發展歷程。由于振幅較小，振動頻率較低，圓柱所受升力系數較小。此外，與靜止圓柱繞流相同，馬赫數固定時，來流越稀薄，升力越大。圖18 為Kn= 0.1、A*分 別取0.05 和0.1 時的升力系數時間發展歷程對比。由于來流為稀薄流，圓柱表面無分離渦和渦脫落現象，流動呈現為附著流的未失穩狀態，此時升力系數與振幅表現為線性增長關系?？紤]到不同馬赫數、克努森數、振幅A*和折減速度U*組合下的圓柱強迫振動流動分析超出了本文的討論范圍，課題組將會在接下來的論文中進行系統討論。

圖17 兩組努森數下圓柱強迫振動升力系數時間發展歷程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＊的升力系數時間發展歷程Fig. 18 Time evolutions of lift coefficients for flows around a forced-oscillating circular cylinder at two different A＊ values

3 結 論

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