一種模擬動(dòng)邊界繞流的銳利界面浸入邊界法*

郭 濤，張晉銘，張紋惠,3，王文全

（1. 昆明理工大學(xué)建筑工程學(xué)院，云南 昆明 650500；2. 四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065；3. 中國(guó)能源建設(shè)集團(tuán)云南省電力設(shè)計(jì)院有限公司，云南 昆明 650051）

繞流是自然界以及實(shí)際工程運(yùn)用中一個(gè)十分普遍的物理現(xiàn)象，也是一個(gè)經(jīng)典的流體力學(xué)問題。例如：機(jī)翼在流動(dòng)空氣中的振動(dòng)問題；導(dǎo)彈無側(cè)滑大攻角飛行時(shí)，背風(fēng)面分離渦的不對(duì)稱運(yùn)動(dòng)導(dǎo)致的附加偏航力問題；水輪機(jī)導(dǎo)葉、葉片在水流作用下的擺動(dòng)和振動(dòng)問題；魚在水中通過魚鰭的擺尾提供推力；昆蟲利用翅膀的拍動(dòng)使其在空中飛行等等。這些運(yùn)動(dòng)及工程都有涉及到翼型、橢圓在流場(chǎng)中的繞流問題。當(dāng)流體繞過不具有流線型的鈍體結(jié)構(gòu)時(shí)，在鈍體尾部會(huì)發(fā)生周期性的漩渦脫落，產(chǎn)生周期性的上下拖曳力，也就是升、阻力。隨著來流情況的變化，當(dāng)渦脫頻率接近結(jié)構(gòu)的自振頻率時(shí)，漩渦脫落和結(jié)構(gòu)振動(dòng)相互鎖定，形成共振，也就是渦激振動(dòng)，其在實(shí)際生活中產(chǎn)生的危害十分常見并且影響也十分廣泛。例如：1940 年剛通車數(shù)月的塔科馬大橋發(fā)生顫振風(fēng)毀事件，震驚世界；以及近期的虎門大橋渦激振動(dòng)事件等。因此，流固耦合問題一直是計(jì)算流體力學(xué)領(lǐng)域較為關(guān)注的一類重要問題。

一般地，研究流固耦合問題時(shí)，通常采用貼體網(wǎng)格，這類網(wǎng)格質(zhì)量良好且生成速度快，但是在實(shí)際計(jì)算中，為了適應(yīng)不斷變化的邊界條件，在每一個(gè)時(shí)間步中都要重新劃分網(wǎng)格，這樣就大大增加了計(jì)算成本。鑒于這些原因，Peskin提出了浸入邊界法，并應(yīng)用到模擬人類心臟血液流動(dòng)中的瓣膜運(yùn)動(dòng)中。浸入邊界法（immersed boundary method，IBM）在整個(gè)物理計(jì)算區(qū)域中采用一套笛卡爾網(wǎng)格，避免了在每一個(gè)時(shí)間步更新網(wǎng)格的繁瑣，使用力場(chǎng)來代替固體邊界，使流固耦合問題中的動(dòng)邊界處理起來更加簡(jiǎn)單和精確。因此，力源項(xiàng)的獲得是浸入邊界法求解流體問題的關(guān)鍵。Peskin 方法，早期是將固體邊界離散為一組由彈簧連接的單元，當(dāng)邊界發(fā)生變形時(shí)，就會(huì)產(chǎn)生一個(gè)意圖使其回到原位置的回復(fù)力，通過邊界上的回復(fù)力結(jié)合delta 函數(shù)求得力源項(xiàng)。但該方法涉及到固體的彈性變形，主要適用于柔性體，而且該變形不易通過流場(chǎng)物理量求得，不具有一般性。Goldstein 等在此基礎(chǔ)上發(fā)展了虛擬邊界法來求解剛性邊界問題，主要使用反饋力法，配合譜方法來漸進(jìn)地施加所需的邊界條件，該法也稱反饋力法。在浸入邊界法的發(fā)展過程中，Jamaludin 等提出了直接力法，即通過直接構(gòu)造受力區(qū)的速度場(chǎng)來推導(dǎo)力源項(xiàng)，從而確保滿足邊界條件，基本實(shí)現(xiàn)了對(duì)邊界的準(zhǔn)確描述，成為學(xué)者們的主要研究方向之一。如Fadlum 等將其應(yīng)用于流固耦合研究，獲得了與貼體網(wǎng)格下計(jì)算結(jié)果一致性較好的結(jié)果；Le 等提出了一種隱式直接力浸入邊界法，Wu 等、王文全等提出了一種對(duì)速度進(jìn)行二次修正的隱式直接力浸入邊界法；Uhlmann等提出了顯式直接力法。直接力法相比于Peskin 的經(jīng)典IBM 法，無須通過繁瑣的delta 函數(shù)完成拉格朗日點(diǎn)與歐拉點(diǎn)上信息交換，通過直接施加浸入邊界處的速度邊界，由此推導(dǎo)出力源項(xiàng)，不依賴于固體的本構(gòu)關(guān)系，也不進(jìn)行反饋調(diào)整，故名直接力法。

按照是否使用delta 函數(shù)或能否對(duì)界面進(jìn)行清晰描述，可將IBM 法分為擴(kuò)散界面（diffused-interface）浸入邊界法和銳利界面（sharp-interface）浸入邊界法。理論上來講，經(jīng)典浸入邊界法、反饋力法和直接力法均屬于擴(kuò)散界面浸入邊界法。銳利界面法不直接計(jì)算力源項(xiàng)來施加邊界條件，主要是在界面處通過局部插值，重構(gòu)邊界網(wǎng)格點(diǎn)附近的流動(dòng)參數(shù)（速度），以此作為邊界條件將其直接施加到計(jì)算域進(jìn)行求解。理論上可以具有高階精度，在可壓縮流體以及激波沖擊問題中均得到了利用。本文將亦采用一種改進(jìn)的銳利界面浸入邊界法，將整個(gè)物理區(qū)域劃分成純流體區(qū)域以及包含固體的次流體區(qū)域。在流體次區(qū)域邊界的法向，構(gòu)造“虛擬點(diǎn)—受力點(diǎn)—垂足點(diǎn)”的計(jì)算結(jié)構(gòu)。受力點(diǎn)速度由流體域內(nèi)某一虛擬點(diǎn)和邊界上垂足點(diǎn)之間的速度線性插值確定,實(shí)現(xiàn)復(fù)雜動(dòng)邊界的流動(dòng)數(shù)值模擬。本文利用C++編程，通過圓柱繞流經(jīng)典算例與其他方法的文獻(xiàn)結(jié)果及試驗(yàn)結(jié)果進(jìn)行對(duì)比，驗(yàn)證該方法的有效性和準(zhǔn)確性?；谠摲椒ǎM(jìn)一步探究不同軸長(zhǎng)比、不同攻角 θ 下的橢圓柱繞流的渦結(jié)構(gòu)分布特征、流線、壓力等流場(chǎng)現(xiàn)象和升阻力系數(shù)、渦脫頻率等水力不穩(wěn)定現(xiàn)象。

1 銳利界面浸入邊界法求解過程

1.1 控制方程

將整個(gè)物理區(qū)域（包括流體、固體）看作不可壓縮牛頓流體的粘性流動(dòng)。在整個(gè)計(jì)算域內(nèi)，流體采用Euler 描述，固體采用Lagrange 描述。其連續(xù)性方程和動(dòng)量方程可表示為

1.2 時(shí)間推進(jìn)

采用二階時(shí)間分布投影法對(duì)時(shí)間進(jìn)行離散，來求解控制方程（1）、（2），具體步驟如下。

（1）忽略力源項(xiàng)，計(jì)算考慮初始?jí)毫ψ饔孟碌闹虚g速度，即：

1.3 空間離散

2 驗(yàn)證算例

如圖2 所示，計(jì)算域?yàn)橐痪匦螀^(qū)域，長(zhǎng)×寬=15×10，流體次區(qū)域?yàn)檫呴L(zhǎng)1.5的正方形，次區(qū)域左側(cè)距流場(chǎng)入口的距離為4.25，次區(qū)域中央浸沒一個(gè)剛性圓柱，其約束方式為全固定邊界，圓柱直徑為，圓心坐標(biāo)為(0,0)。計(jì)算域左側(cè)為均勻來流入口邊界，采用Dirichlet 邊界條件，即=,=0；上、下兩側(cè)均為無穿透邊界；右側(cè)為自由出流邊界，采用Neumenn 邊界，即?/?=0，?/?=0 。整個(gè)流場(chǎng)采用一套間距為Δ=Δ=0.025的均勻四邊形網(wǎng)格，固體邊界采取與流體網(wǎng)格尺度相等的等間距離散，時(shí)間步長(zhǎng)為Δ=0.002。

圖1 固體邊界處理Fig. 1 Treatment of the solid boundary

圖2 整體計(jì)算域及邊界條件Fig. 2 Computational domain and boundary conditions

圖3 為=160 s、雷諾數(shù)=300 時(shí)圓柱繞流的流場(chǎng)速度分布（范圍： - 1≤/≤1 ）。從圖中可看出，速度等值線表現(xiàn)光滑，說明計(jì)算過程中對(duì)速度的更新、修正并沒有造成流場(chǎng)的震蕩。由此可知，采用設(shè)置中間速度并不斷更新速度的方法是適用的。同時(shí)，固體邊界周圍的速度分布也符合真實(shí)的圓柱繞流流場(chǎng)分布規(guī)律。說明利用這種銳利界面法描述固體對(duì)流體的耦合作用和運(yùn)用雙線性插值方法計(jì)算虛擬點(diǎn)的速度是可行的。

圖3 流向速度分布Fig. 3 Isolines of velocity

圖4 為一個(gè)渦脫落周期內(nèi)圓柱尾跡渦隨時(shí)間的演化過程（范圍： - 8≤wD/≤8 ，為渦量；實(shí)線為正值，虛線為負(fù)值）。從中可看出柱體表面的邊界層分布、分離剪切層的卷起和圓柱后方的分離區(qū)域清晰可見。在一個(gè)周期內(nèi)圓柱尾部出現(xiàn)明顯的正、負(fù)漩渦交替脫落，向下游發(fā)展，呈現(xiàn)出典型的卡門渦街現(xiàn)象，與試驗(yàn)得到的卡門渦街現(xiàn)象一致。說明本文對(duì)邊界附近的流動(dòng)參數(shù)插值的計(jì)算方法及程序是可行的。

圖4 一個(gè)周期內(nèi)尾跡渦的演化Fig. 4 Isolines of vorticity against time

圖5 為=300 時(shí)尾渦脫落穩(wěn)定后的升力系數(shù)及阻力系數(shù)的時(shí)程曲線，可以看出，隨著時(shí)間的推進(jìn)，流場(chǎng)逐漸趨于穩(wěn)定，在流場(chǎng)穩(wěn)定后升阻力系數(shù)都隨時(shí)間的發(fā)展而呈現(xiàn)出穩(wěn)定的周期擺動(dòng)。其阻力系數(shù)、Strouhal 數(shù)()的全時(shí)域統(tǒng)計(jì)結(jié)果見表1 所示。其中升力系數(shù)=2F/ρ，阻力系數(shù)=2F/ρ，Strouhal 數(shù)=/,為渦的脫落頻率。對(duì)比表中結(jié)果可知，相比于文獻(xiàn)[16,29]的數(shù)值解，本文的計(jì)算結(jié)果與文獻(xiàn)[30]的實(shí)驗(yàn)值更為接近。驗(yàn)證了本方法及程序的準(zhǔn)確性。

圖5 升力、阻力系數(shù)時(shí)程曲線Fig. 5 Variations of the lift and drag coefficients with time

表1 本文結(jié)果與其他文獻(xiàn)結(jié)果的對(duì)比（ R e=300 ）Table 1 The results comparison of average drag coefficients and Strouhal number at Re=300

3 橢圓柱主動(dòng)運(yùn)動(dòng)算例

3.1 計(jì)算模型

目前，以圓柱為對(duì)象的鈍體繞流的研究已較為充分，對(duì)橢圓柱的關(guān)注度相對(duì)要少一些，但橢圓柱作為介于圓柱和平板之間更具有代表性的鈍體，其形狀更具流線型，流動(dòng)阻力低，具有較好的對(duì)抗流體的作用，而且剪切邊界層的卷起、尾跡渦的脫落、轉(zhuǎn)捩和耗散等也很復(fù)雜，其研究對(duì)工程實(shí)際具有重要的現(xiàn)實(shí)意義。取計(jì)算域尺寸、邊界條件、網(wǎng)格尺寸、時(shí)間步長(zhǎng)等均與驗(yàn)證算例一致，雷諾數(shù)=800。不同之處在于：流體次區(qū)域中央浸沒的固體為一系列不同軸長(zhǎng)比=/的橢圓柱（見圖6），其繞樞軸點(diǎn)（原點(diǎn)）旋轉(zhuǎn)擺動(dòng)時(shí)攻角 θ （擺角）的變化規(guī)律由下式控制：

圖6 橢圓柱繞流計(jì)算模型Fig. 6 Computational model of flow around an elliptical cylinder

式中： θ為最大攻角值，取80°； ω =2π為振蕩角頻率，為運(yùn)動(dòng)頻率， β 為曲線形狀參數(shù)。該橢圓柱的擺動(dòng)形式由水翼在正弦振蕩的基礎(chǔ)上演變而來（當(dāng) β=1 時(shí)為正弦波的半個(gè)波長(zhǎng)）。 β=2.5 代表非正弦振蕩，前后0.2 s 為勻速開、關(guān)時(shí)間段，中間持續(xù)0.6 s，其1 個(gè)開、關(guān)周期內(nèi)的方波曲線如圖7 所示。

圖7 一個(gè)周期內(nèi)橢圓柱的瞬時(shí)攻角變化曲線Fig. 7 Changes of the instantaneous angle of attack of the elliptical cylinder within a period

3.2 不同軸長(zhǎng)比對(duì)橢圓柱繞流的影響

由圖8 可以看出，隨著的不斷增大，平均阻力系數(shù)以及最大升力系數(shù) ||均逐漸增大。由于橢圓柱垂直于水流方向的投影面積不斷增大，當(dāng)流體流過橢圓柱時(shí)，對(duì)橢圓柱的沖擊力增大，因此橢圓柱對(duì)流體的反作用力也會(huì)增大，橢圓柱所受到來自流體的阻力必然會(huì)增大。但是，相比阻力而言，升力的增大幅度比較小。由圖9 可以看出，渦脫頻率隨著軸長(zhǎng)比的變化趨勢(shì)可以分成3 個(gè)階段。初期：隨著軸長(zhǎng)比的增大而降低，當(dāng)=0.2 時(shí)，橢圓柱的渦脫頻率最大，約為0.6 Hz；第二階段：當(dāng)進(jìn)入某一區(qū)域（=0.5～0.7）后，開始穩(wěn)定下來，約為0.38 Hz；后期：當(dāng)超出這一區(qū)間后，渦脫頻率將再次隨著的增大而遞減，“穩(wěn)定現(xiàn)象”消失，在=0.9 時(shí)，達(dá)到最小值。

圖8 平均阻力系數(shù)、最大升力系數(shù)隨軸長(zhǎng)比的變化Fig. 8 Variation of lift and drag coefficients with axis ratio

圖9 渦脫頻率隨軸長(zhǎng)比的變化Fig. 9 Variation of vortex shedding frequency with axis ratio

由圖10 可知，當(dāng)=0.2 時(shí)，橢圓柱更具流線型，流動(dòng)阻力低，尾部并沒有明顯的分離渦，流場(chǎng)比較平穩(wěn)；當(dāng)=0.5 時(shí)，尾部流線波動(dòng)開始逐漸明顯；隨著軸長(zhǎng)比的增大，流態(tài)波動(dòng)更加劇烈；當(dāng)=0.7 時(shí)，在后方開始產(chǎn)生一個(gè)尺度較小的分離渦，其位置向上偏離橢圓柱的中心。

圖10 不同軸長(zhǎng)比下的流線圖Fig. 10 The instantaneous streamlines with different axis ratio

圖11 為=800 時(shí)不同軸長(zhǎng)比工況下，同一時(shí)刻下的瞬時(shí)渦量圖（ - 8≤wD/≤8 ，w為渦量；實(shí)線為正值，虛線為負(fù)值）?？梢钥闯觯寒?dāng)=0.2 時(shí)，尾跡渦以S 型模態(tài)脫落，呈現(xiàn)出上側(cè)負(fù)漩渦與下側(cè)正漩渦交替脫落向下游演化發(fā)展的卡門渦街現(xiàn)象。脫落的正負(fù)漩渦之間的縱向間距、周期較小，形成明顯的2 列，尾渦強(qiáng)度（尺寸）也較小，漩渦的分離點(diǎn)基本位于軸上；隨著軸長(zhǎng)比的增大，卡門渦街規(guī)律依然明顯，但漩渦強(qiáng)度（尺寸）和周期逐漸增大，漩渦分離點(diǎn)偏離軸，且脫落位置增大，正、負(fù)漩渦中心間距減小逐漸趨于一列。

圖11 不同軸長(zhǎng)比工況下的瞬時(shí)渦量Fig. 11 Variation of instantaneous vorticity with axis ratio

3.3 不同攻角對(duì)橢圓柱繞流的影響

圖12 給出了=0.2，=800 時(shí)，不同攻角 θ 下，升、阻力系數(shù)的變化規(guī)律。

圖12 不同攻角對(duì)升、阻力系數(shù)的影響Fig. 12 Variations of lift and drag coefficients with angle of attack

從圖13 中也可看出，當(dāng) θ =0時(shí)橢圓柱上下表面的壓力對(duì)稱，壓差接近零，升力系數(shù)很?。浑S著攻角增大，橢圓柱的上表面的壓力，不管從數(shù)值和承壓面積上來看都明顯比下表面大，即上下表面的壓差增大，故升力系數(shù)隨攻角的增大而增大。當(dāng)攻角發(fā)展到臨界攻角 θ =45時(shí)，壓差達(dá)到最大，之后升力系數(shù)的值隨著攻角的增大而逐漸減小。

圖13 不同攻角下的瞬時(shí)壓力場(chǎng)Fig. 13 Instantaneous pressure fields at different angles of attack

由圖14 漩渦脫頻率曲線可以看到，渦脫頻率隨著攻角 θ 的不斷增大，整體呈現(xiàn)先增加再逐漸遞減的趨勢(shì)。其變化趨勢(shì)可以分成3 個(gè)階段：初期， θ =2°～5°階段，隨著攻角 θ 的增大而增加， θ =2時(shí)，渦脫頻率達(dá)到最大，約為=0.61 Hz；第二階段，當(dāng)攻角進(jìn)入 θ =2°～10°區(qū)域后，有所穩(wěn)定，遞減速率較小，在0.6 Hz 左右浮動(dòng)；后期，當(dāng)攻角 θ ＞10后，渦脫頻率將隨著攻角 θ 增大而逐漸減小，后期曲線逐漸平緩，減小的速率逐漸減慢。

圖14 不同攻角對(duì)渦脫頻率的影響Fig. 14 Variation of vortex shedding frequency with angle of attack

圖15 為=800 時(shí)不同攻角工況下橢圓柱尾跡渦的對(duì)比，由圖可見，固體表面分布的邊界層、分離剪切層的卷起和橢圓柱后側(cè)的分離區(qū)，以及渦的脫落等現(xiàn)象都清晰可見。與圓柱繞流一樣，橢圓柱繞流尾跡渦也是呈現(xiàn)出明顯的一對(duì)正負(fù)漩渦從上下兩側(cè)交替脫落，向下游發(fā)展的卡門渦街現(xiàn)象。但是，橢圓柱擺動(dòng)的攻角幅值對(duì)流態(tài)起決定性作用。隨著攻角的增加，渦的脫落周期、尺寸、形狀和演化形式等均有較大差異，渦脫落的反對(duì)稱性越發(fā)明顯。類似于水輪機(jī)活動(dòng)導(dǎo)葉大攻角運(yùn)行或?qū)棢o側(cè)滑大攻角飛行時(shí)，背流（風(fēng)）面分離渦的不對(duì)稱運(yùn)動(dòng)。攻角越大越容易產(chǎn)生不希望的附加偏航力，或加大水輪機(jī)活動(dòng)導(dǎo)葉轉(zhuǎn)動(dòng)軸折斷的風(fēng)險(xiǎn)。

圖15 不同攻角下的瞬時(shí)渦量圖Fig. 15 Variation of the instantaneous vorticity with the angle of attack

隨著橢圓柱擺動(dòng)攻角的增大，漩渦脫落頻率降低，周期加長(zhǎng)，漩渦尺寸加大。當(dāng)攻角 θ ≤30時(shí)，在每個(gè)周期內(nèi)，有兩個(gè)符號(hào)相反的渦從橢圓柱上下側(cè)脫落，在下游形成交替出現(xiàn)的渦對(duì)，是典型的反對(duì)稱S 型脫落模態(tài)，如圖16(a)所示，類似于靜止圓柱的渦脫落方式。當(dāng) θ =45°～60°時(shí)，上下側(cè)脫落的漩渦尾部均拖著一個(gè)方向相同的子渦（2 個(gè)渦核）；在向下游演化發(fā)展的過程中，上側(cè)脫落的負(fù)漩渦的子渦與母渦融合，而下側(cè)脫落的正漩渦的子渦則產(chǎn)生分裂，分裂后的次正漩渦位于最上側(cè)，靠近負(fù)漩渦，最終尾跡渦表現(xiàn)為3 排，在下游形成了一對(duì)渦和一個(gè)單渦交替脫落的“P+S”Ⅰ型模態(tài)，如圖16(b) 所示。當(dāng)θ≥75時(shí)，則情況相反，上側(cè)脫落的負(fù)漩渦的子渦發(fā)生分裂,分裂后的次負(fù)漩渦位于最下側(cè)，靠近正漩渦，而橢圓柱下側(cè)脫落的正漩渦的子渦則與母渦則融合，在下游形成了新的“P+S”Ⅱ型模態(tài)，如圖16(c)所示。

圖16 不同攻角下的渦脫落模態(tài)示意圖Fig. 16 The diagram of vortex modes at different angles of attack

4 結(jié) 論

采用一種銳利界面（sharp-interface）浸入邊界法對(duì)二維橢圓柱繞流問題進(jìn)行了數(shù)值計(jì)算，并討論了不同軸長(zhǎng)比、攻角 θ 對(duì)渦結(jié)構(gòu)分布特征及升阻力系數(shù)、渦脫頻率等水力不穩(wěn)定現(xiàn)象的影響，通過對(duì)比分析，得到以下結(jié)論：

(1)在模擬鈍體繞流時(shí)，通過與文獻(xiàn)和試驗(yàn)結(jié)果的對(duì)比，本文的銳利界面浸入邊界法結(jié)果與試驗(yàn)結(jié)果吻合較好，驗(yàn)證了本方法的準(zhǔn)確性和有效性；

(2)對(duì)于橢圓柱繞流問題，隨著不同軸長(zhǎng)比的增加，升阻力系數(shù)呈遞增趨勢(shì)，而流場(chǎng)流態(tài)波動(dòng)也越來越劇烈，尾跡渦脫落位置、強(qiáng)度（尺寸）、間距、周期也越來越大。周期加長(zhǎng)，渦脫頻率則降低；

(3)當(dāng)攻角發(fā)生變化時(shí)，平均阻力系數(shù)隨攻角的增大而呈遞增趨勢(shì)，而最大升力系數(shù)，以及最大升力系數(shù)和平均阻力比（升阻比），均表現(xiàn)為前期隨著攻角的增大而增大，超過臨界攻角后，再降低的趨勢(shì)；升阻比的臨界攻角為25°，而最大升力系數(shù)的臨界角為45°；

(4)當(dāng)攻角發(fā)生變化時(shí)，渦脫頻率總體呈隨著攻角增大而逐漸減小的趨勢(shì)。漩渦脫落頻率降低，周期加長(zhǎng)，漩渦尺寸加大；當(dāng)攻角 θ ≤30時(shí)，尾渦呈反對(duì)稱S 型脫落模態(tài)，當(dāng) 4 5≤θ≤60時(shí)，尾渦呈反對(duì)稱“P+S”Ⅰ型脫落模態(tài)，當(dāng) θ ≥75時(shí)，尾渦呈反對(duì)稱“P+S”Ⅱ型脫落模態(tài)。