基于IBVCM的二維非定常動邊界數(shù)值模擬研究

盧海天，吳 杰，趙 寧

（南京航空航天大學(xué) 空氣動力學(xué)系，南京 210016）

0 引 言

隨著仿生學(xué)和流體力學(xué)的發(fā)展，魚類的游動機理以及推進力與游動參數(shù)的關(guān)系成為當(dāng)前研究的熱點。根據(jù)魚游動的特性，人們已經(jīng)初步研制出了具有較好性能的仿生魚。但是，就目前而言，仿生魚中仍待解決的問題包括：（1）魚類推進理論模型以及游動算法的研究；（2）在研制上，選擇合適的材料和動力實現(xiàn)仿生魚機構(gòu)的最優(yōu)；（3）仿生魚控制技術(shù)的研究，解決推進與穩(wěn)定的相互制約；（4）數(shù)值模擬上，涉及到單一魚體的縱向前進和對于多個魚的游動及其協(xié)作問題，成果尚不顯著。

相對于實驗研究，數(shù)值計算模擬具有研究周期短、費用少等實驗無法替代的優(yōu)勢。自20世紀60年代開始，許多來自不同領(lǐng)域的學(xué)者對魚游動做出了大量的數(shù)值模擬研究。通過分析數(shù)值結(jié)果，可以獲得相關(guān)參數(shù)對魚游動的影響，從而為仿生魚的設(shè)計和優(yōu)化提供依據(jù)。其中，1960年Lighthill［1］率先對體形瘦長的鰻鱺模式游動的魚提出了細長體模型，得到了定常線性化解。1973年Newman分析了在無粘流體下魚體表面的作用力，計算出魚體表面壓力公式［2］。由于考慮魚體厚度的影響，Newman得出與早期Lighthill，Wu和Newman研究細長型魚體時不同的誘導(dǎo)阻力結(jié)果。80年代中后期童秉綱［3］等建立了模擬魚類游動的三維波動板理論（3DWPT），采用勢流理論的三維非定常渦格法在頻域內(nèi)求解，研究了各類魚的推進性能。1996年Liu和Wassersug應(yīng)用一種新的二維CFD建模方法研究蝌蚪的游動［4］，得出蝌蚪能夠產(chǎn)生高效的推進力，而且其運動特性與自身形狀有關(guān)。

魚游動是典型的非定常動邊界問題，它對傳統(tǒng)的數(shù)值模擬方法提出了巨大的挑戰(zhàn)。目前，已經(jīng)形成的一些數(shù)值方法可以有效地模擬此類問題。其中比較典型的是由Peskin［5］在1977年提出的浸入邊界方法（IBM）。它較好地模擬了動邊界的流動問題，但其主要缺點是在邊界處不能保證無滑移條件。另一種數(shù)值模擬是格子波爾茲曼方法（LBM）。它簡潔且易于實現(xiàn)，沒有涉及到微分方程以及相應(yīng)的代數(shù)方程［6－9］。

本文采用格子波爾茲曼方法（LBM）和浸入邊界速度修正方法（IBVCM）相結(jié)合，僅研究魚體在擺動的模式下影響其推進力的一些基本參數(shù)。涉及到的基本參數(shù)有魚的形狀、頻率、波長和雷諾數(shù)，通過逐個參數(shù)以及參數(shù)組合的形式分析影響魚擺動的基本規(guī)律。此外，從旋渦分布的角度出發(fā)，分析并得到旋渦變化與魚身運動的相互關(guān)系。和傳統(tǒng)的IBM相比，本文的方法更精確地滿足無滑移邊界條件，從而更適合模擬動邊界的流動問題。

1 數(shù)值研究方法與魚擺動的數(shù)學(xué)模型

1.1 浸入邊界速度修正方法（IBVCM）

由于在傳統(tǒng)的IBM中，體積力是預(yù)先設(shè)定的，不能保證由修正的速度場插值得到的邊界速度滿足無滑移條件。為了克服這一難點，本文采用IBVCM對浸入邊界方法進行修正［10］。IBVCM的基本思想是不預(yù)設(shè)體積力，而是先由周圍流場插值得到壁上的速度，再使該速度滿足邊界無滑移條件，從而實現(xiàn)物體上的邊界速度與流體的邊界速度相等。

如圖1所示，粗實線表示物體邊界，可以看出物體邊界與網(wǎng)格的交點Mx與My。根據(jù)動量方程，點Mx僅僅影響點A與點B速度分量u的修正，點My僅僅影響點B與點C速度分量v的修正。點A、點B和點C的速度分量分別表示為（uA，vA）、（uB，vB）、（uC，vC），由流場x方向插值得到點Mx處的速度分量表示為（uix，vix），由流場y方向插值得到點My處的速度分量表示為（uiy，viy）。此外，點Mx處的壁上速度表示為（uMx，vMx），點My處的壁上速度表示為（uMy，vMy）。由線性插值，邊界點Mx的插值速度［11］為：

邊界點My的插值速度為：

則點A與點B之間的速度修正為：

可以看出，當(dāng)點A與點B的速度分量u由式（3）得到修正后，點Mx處的插值速度滿足邊界條件u＝uMx。

類似地，點B與點C之間速度修正為：

某些情形下，比如在一個網(wǎng)格間距內(nèi)網(wǎng)格線與邊界不止有一個交點時，所有的速度分量都應(yīng)進行修正。

圖1 笛卡爾網(wǎng)格線，浸入邊界和它們的交點Fig.1 Cartesian mesh lines，immersed boundary and their intersection points

1.2 格子波爾茲曼方法（LBM）

格子波爾茲曼方法是一種以介觀動理學(xué)模型為基礎(chǔ)的數(shù)值模擬方法。與傳統(tǒng)的數(shù)值方法相比，LBM的主要優(yōu)勢在于計算簡單且容易實現(xiàn)。在本文的數(shù)值模擬中，IBVCM用于處理魚的邊界速度修正，而LBM用于求出笛卡爾網(wǎng)格線上的中間速度場信息。這樣，就把兩種方法進行了耦合。這里，格子波爾茲曼方程為：

其中，τ是單松弛時間，fα是分布函數(shù)，是局部平衡分布函數(shù)，δt是時間步長，eα是粒子速度。在本文的LBM計算中使用D2Q9格子速度模型。這樣，eα可以描述為：

而平衡分布函數(shù)為：

其中，對于α＝0，wα＝4／9；α＝1～4，wα＝1／9；α＝5～8，wα＝1／36。式（7）中的密度ρ與速度V是修正速度場在上一個時間層上的值，可以看出修正速度場對分布函數(shù)fα的影響是通過平衡分布函數(shù)體現(xiàn)的。由質(zhì)量與動量守恒關(guān)系，宏觀密度ρ與流體速度V按照密度分布函數(shù)可計算為：

綜上所述，沿時間方向，當(dāng)魚體每一次變形時，本文LB－IBVCM的基本求解步驟為：

（1）設(shè)體積力為0，應(yīng)用LBM求解流場運動得到中間速度場；

（2）根據(jù)式（3）和式（4）將物體邊界與網(wǎng)格線交點周圍的網(wǎng)格點速度進行一次修正。

1.3 魚擺動的數(shù)學(xué)模型

Videler（1993）在研究中發(fā)現(xiàn)魚體的橫向?qū)ΨQ運動可以由傅里葉級數(shù)中奇數(shù)項的前三項來精確地描述［12］。考慮到第三個之后的頻率項與第一個頻率項相比非常小，在忽略高頻之后，魚體中線的運動可以描述為：

其中ym（x，t）是y方向的中線坐標，a1（x）和b1（x）表示傅里葉因子，f是行波的頻率，t是時間，x是x方向的坐標。將a1（x）和b1（x）表示為：

可以更進一步地推導(dǎo)出公式：

此處，l表示波長，w＝2πf表示圓頻率，am（x）表示波幅。對式（11）關(guān)于時間t求導(dǎo)，得出魚體中線各個點做橫向運動時的速度：

波幅可以選取二階多項式函數(shù)的形式，基于大量的數(shù)值模擬研究，在本文中取為：

2 魚擺動的數(shù)值模擬及參數(shù)分析

本文對影響魚類擺動的基本參數(shù)：形狀，頻率，波長和雷諾數(shù)，采用逐一以及組合的形式進行分析。具體研究的問題是，給定自由來流速度U，模型在流向固定，而在垂直于來流方向魚體中線的各個點做形如式（11）所示的擺動，通過結(jié)果分析從而得出影響魚類擺動的一些基本結(jié)論。魚體表面必須滿足無滑移邊界條件。雷諾數(shù)定義為：Re＝UL／ν，L為魚體長度，ν是運動粘度。升力系數(shù)與阻力系數(shù)定義為：

其中ρ為流場密度。FL和FD分別為升力和阻力，它們分別通過控制體積法求出。對矩形控制體Ω通過積分動量方程得：

其中S是控制體外表面，n＝（n1，n2）是控制體表面邊界的法向量，下標1，2分別表示x方向與y方向。

為了驗證本文的方法，選取Re＝3000，波長l＝0.87，頻率f＝1.3，模擬了NACA0012翼型形狀的魚擺動，所得的阻力系數(shù)為0.062，與文獻［10］中阻力系數(shù)0.06基本一致。這說明本文采用的LB－IBVCM可以很好地模擬魚類的擺動問題。

2.1 網(wǎng)格的劃分

在本文的研究中，計算區(qū)域采用非均勻網(wǎng)格。但在魚體周圍采用均勻網(wǎng)格，尺寸為120×120，其步長為0.01。整個流場區(qū)域的網(wǎng)格為501×351。流場區(qū)域如圖2所示?？紤]到標準的LBM只能應(yīng)用于均勻網(wǎng)格，所以本文采用基于泰勒級數(shù)展開和最小二乘方的LBM（TLLBM）處理非均勻網(wǎng)格［13］。

圖2 流場區(qū)域Fig.2 The flow field

2.2 魚體形狀的影響

本文總共選取五種形狀作為魚體模型進行分析，分別 是 RoboTuna、NACA0006、NACA0012、NACA2412、NACA2424。其中，RoboTuna模型［14］為：

如圖3所示，選取Re＝400，l＝1和w＝4，研究不同形狀魚體的阻力系數(shù)關(guān)于時間的變化關(guān)系。另外固定Re＝400，分別研究不同形狀魚體的阻力系數(shù)關(guān)于頻率和波長的變化關(guān)系，如圖4所示。可以看出，NACA0012和NACA2412翼型的平均阻力相接近，而且均大于NACA0006的阻力，但遠小于Robo－Tuna和NACA 2 4 2 4的阻力。這很好地說明了彎度對阻力的影響不大，而厚度對阻力影響較大。魚的厚度越大，阻力越大。其原因為厚度越大，流經(jīng)翼型后的逆壓梯度越大，流體越容易分離，從而大大加寬了尾流，導(dǎo)致壓差阻力增大。

圖3 Re＝400，l＝1和w＝4時，不同魚體形狀下魚擺動的阻力系數(shù)Fig.3 The drag coefficient of fish swimming for different body shapes at Re＝400，l＝1and w＝4

圖4 不同魚體形狀時的阻力系數(shù)Fig.4 The drag coefficient for different body shapes

2.3 頻率和波長的影響

對于本文中的魚擺動而言，推進力是負的阻力。這主要由于在渦沿著魚體脫落的過程中，上表面的渦從尾尖處脫落至魚體中線的下側(cè)，而下表面的渦從尾尖處脫落至中線的上側(cè)，這就形成了逆卡門渦街，從而產(chǎn)生推進力。由圖4（a）可以看出，隨著w的增加，推進力以及曲線斜率均增大。另外也可看出，Robo－Tuna產(chǎn)生推進力的臨界頻率是w＝14.1，而NACA0012與NACA2412翼型的臨界頻率約在w＝12.8。

如圖4（b）所示，可以看出當(dāng)波長在0.7～1.3的范圍內(nèi)時，阻力系數(shù)隨著波長的增大而減小。對于NACA0012翼型，在Re＝400和w＝4時，圖5給出了不同波長下升力系數(shù)與阻力系數(shù)的關(guān)系。可以看出，不同波長下的相位圖均為關(guān)于CL＝0對稱且類似于蝴蝶形狀的圖形，可知阻力系數(shù)的頻率為升力系數(shù)的兩倍，升力系數(shù)的平均值為0。其原因是魚的形狀是對稱的，來流初始平行于魚體中線，中線上各點亦做橫向?qū)ΨQ運動。而且當(dāng)波長由0.9增大至1.3時，CL的幅度一直在增大，曲線向左移，平均阻力減小。這主要因為在波長的研究中，行波的相速度Cp（Cp＝lf）均大于來流速度U。隨著波長的增大，比值Cp／U增大，魚體周圍的擾動被進一步遏制，流體不易分離，魚體表面的應(yīng)力減小，從而導(dǎo)致阻力減?。?5］。

2.4 頻率、波長、雷諾數(shù)的綜合影響

由于國內(nèi)外很多工作都是關(guān)于大雷諾數(shù)的研究，而低雷諾數(shù)的研究不多見，所以本文僅分析低雷諾數(shù)對推進力的影響。

圖5 Re＝400，w＝4，不同波長下升力系數(shù)與阻力系數(shù)的相位圖Fig.5 Drag and lift coefficient phase diagram varies with wavelength at Re＝400and w＝4

對于NACA0012翼型，為了研究波長以及雷諾數(shù)的影響，如圖6（a）所示，給定w＝14，分析了l＝1和1.7以及不同雷諾數(shù)下的推進力。當(dāng)波長由1增大至1.7后，曲線下移，可知增大雷諾數(shù)與波長可使推進力增大。對于l＝1.7，在Re≥50時即可產(chǎn)生推進力。為了研究雷諾數(shù)和頻率的影響，如圖6（b）所示，給定l＝1，分析了Re＝100和500以及不同頻率下的推進力。當(dāng)雷諾數(shù)由100增大至500后，曲線下移，可知增大雷諾數(shù)與頻率可使推進力增大。在Re＝100時，只有w≥15.6時，才產(chǎn)生推進力。而對于Re＝500，當(dāng)w≥12.5時即可產(chǎn)生推進力，且平均無量綱推進力可達0.63或更多。

圖6 頻率、波長、雷諾數(shù)的綜合影響Fig.6 The combined effects of frequency，wavelength and Reynolds number

2.5 旋渦分析

在圖6（b）中，當(dāng)Re＝500，l＝1和w＝18時，產(chǎn)生的推進力系數(shù)為0.63。圖7顯示了在此運動模式下尾尖沿y方向運動與推進力的變化。對圖中A點與B點進行分析，相應(yīng)時刻的渦分布分別如圖8與圖9所示?？梢园l(fā)現(xiàn)，A點的推進力最大，B點的推進力近似為零。魚的運動與渦的存在，是導(dǎo)致A點與B點推進力差異的主要原因。為了分析尾尖運動對推進力的影響，本文對A點與B點處的旋渦運動進行了分析。

圖7 阻力系數(shù)與魚擺動時相應(yīng)的尾尖的垂直位置Fig.7 The drag coefficient and corresponding vertical location of tail tip of fish swimming

圖8 A點的渦分布Fig.8 The vortex distribution for pointA

圖9 B點的渦分布Fig.9 The vortex distribution for point B

在點A，魚的尾尖向下運動，接近幅度的谷值，而且此時尾尖是在減速的，在改變方向之前是趨于停止的。此刻尾尖處形成順時針方向的渦，魚體下表面有逆時針方向的渦脫落。而且在魚尾后方有相反方向渦的破壞性結(jié)合［16］，形成了新的逆時針方向的渦，所以阻力減小很多。

在點B，魚的尾尖向上運動接近中線的位置，而且尾尖的速度接近于最大值。此刻魚體頭部上表面有順時針方向的渦形成。由于魚的運動，在A點時結(jié)合形成的逆時針方向的渦繼續(xù)向下游運動。同時可以看到，尾尖處不同方向的渦相遇開始產(chǎn)生破壞性的作用，由于還未形成新的渦，所以阻力較A點增加。

3 結(jié) 論

本文通過對魚擺動的一些基本參數(shù)的研究，分析了其對推進力的影響，得到以下結(jié)論：

（1）對于各種模式下魚的擺動，產(chǎn)生推進力時都存在一個臨界頻率、臨界波長和臨界雷諾數(shù)。通過參數(shù)分析，可知在魚體擺動幅度保持不變的情況下，增大頻率、增大波長、增大雷諾數(shù)均可使魚的推進力增大，而相應(yīng)的臨界值也會減小。

（2）魚體中線是以行波的模式振動的。研究表明，魚的形狀不同，產(chǎn)生的推進力也不同。厚度越大，推進力越小。彎度對推進力影響較小。

（3）魚產(chǎn)生的推進力與魚尾尖的位置以及渦有關(guān)。尾尖在不同位置，尾尖的速度不同，因而魚表面渦的分布和強度不同，導(dǎo)致推進力的變化。

［1］LIGHTHILL M J.Note on the swimming of slender fish［J］.JournalofFluidMechanics，1960，9（2）：305－317.doi：10.1017／S0022112060001110

［2］NEWMAN J N.The force on a slender fish－like body［J］.Jour－nalofFluidMechanics，1973，58（4）：689－702.doi：10.1017／S0022112073002429

［3］TONG B G，ZHUANG L X.Hydrodynamic model for fish＇s undulatory motion and its applications［J］.ChineseJournalof Nature，1998，20（1）：1－7.（in Chinese）童秉綱，莊禮賢.描述魚類波狀游動的流體力學(xué)模型及其應(yīng)用［J］.自然雜志，1998，20（1）：1－7.

［4］LIU H，WASSERSUG R，KAWACHI K.A computational fluid dynamics study of tadpole swimming［J］.TheJournalof ExpermentalBiology，1996，199：1245－1260.

［5］PESKIN C S.Numerical analysis of blood flow in the heart［J］.JournalofComputationalPhysics，1977，25（3）：220－252.doi：10.1016／0021－9991（77）90100－0

［6］CHEN H.Volumetric formulation of the lattice Boltzmann method for fluid dynamics：Basic concept［J］.PhysicalReviewE，1998，58（3）：3955－3963.doi：10.1103／PhysRevE.58.3955

［7］INAMURO T，OGATA T，TAJIMA S，et al.A lattice Boltzmann method for incompressible two－phase flows with large density differences［J］.JournalofComputationalPhysics，2004，198（2）：628－644.doi：10.1016／j.jcp.2004.01.019

［8］LALLEMAND P，LUO L S.Theory of the Lattice Boltzmann method：Dispersion，dissipation，isotropy，Galilean invariance，and stability［J］.PhysicalReviewE，2000，61（6）：6546－6562.doi：10.1103／PhysRevE.61.6546

［9］LEE T，LIN C L.A stable discretization of the lattice Boltzmann equation for simulation of incompressible two－phase flows at high density ratio［J］.JournalofComputationalPhysics，2005，206（1）：16－47.doi：10.1016／j.jcp.2004.12.001

［10］SHU C，LIU N，CHEW Y T，et al.Numerical simulation of fish motion by using lattice Boltzmann－immersed boundary velocity correction method［J］.JournalofMechanicalScience andTechnology，2007，21（9）：1352－1358.doi：10.1007／BF03177420

［11］SHU C，LIU N，CHEW Y T.A novel immersed boundary velocity correction－lattice Boltzmann method and its application to simulate flow past a circular cylinder［J］.JournalofComputationalPhysics，2007，226（2）：1607－1622.doi：10.1016／j.jcp.2007.06.002

［12］VIDELER J J.Fish swimming［M］.London：Chapman ＆Hall，1993.

［13］SHU C，NIU X D，CHEW Y T.Taylor series expansion and least squares－based lattice Boltzmann method：three－dimensional formulation and its applications［J］.InternationalJournalof ModernPhysicsC，2003，14（7）：925－944.doi：10.1142／S0129183103005133

［14］ZHU Q，WOLFGANG M J，YUE D K P，et al.Three－dimensional flow structures and vorticity control in fish－like swimming［J］.JournalofFluidMechanics，2002，468：1－28.doi：10.1017／S002211200200143X

［15］BARRETT D S，TRIANTAFYLLOU M S，YUE D K P，et al.Drag reduction in fish－like locomotion［J］.JournalofFluid Mechanics，1999，392：183－212.doi：10.1017／S0022112099005455

［16］WU J，SHU C.Numerical study of flow characteristics behind a stationary circular cylinder with a flapping plate［J］.Physicsof Fluids，2011，23（7）：073601.doi：10.1063／1.3601484