基于直角網(wǎng)格法的單個和陣列布置下柔性水翼繞流數(shù)值模擬*

辛建建 陳振雷? 石凡 石伏龍

1) (寧波大學(xué)海運(yùn)學(xué)院, 寧波 315211)

2) (長沙理工大學(xué)水利工程學(xué)院, 長沙 410114)

研究柔性水翼在不可壓縮流體中的水動力特性, 對于船舵和減搖鰭等海洋結(jié)構(gòu)物的設(shè)計和性能優(yōu)化具有重要意義.本文將自主開發(fā)的徑向基函數(shù)虛擬網(wǎng)格法求解器擴(kuò)展到模擬繞單個或多個柔性水翼的不可壓縮流動問題.數(shù)值模型基于虛擬網(wǎng)格有限差分法考慮浸入邊界對流場的影響, 引入緊支徑向基函數(shù)(compact supported radial basis function, CSRBF)以物面 Lagrangian 質(zhì)點(diǎn)追蹤復(fù)雜的柔性動邊界.基于該方法, 首先模擬了均勻流中主動拍動的柔性水翼, 升阻力系數(shù)良好的網(wǎng)格收斂性結(jié)果驗(yàn)證了本文方法的精度和可靠性.并研究了柔性水翼在不同振蕩頻率下的水動力特性, 闡述了柔性水翼的推力生成機(jī)制.然后模擬了繞陣列布置柔性水翼的流動現(xiàn)象, 研究了不同間距和不同振蕩頻率下水翼表面的升阻力系數(shù)變化規(guī)律和尾渦特性, 觀察到緊密布置的柔性水翼在高頻振蕩下推力系數(shù)存在顯著的放大效應(yīng), 同時推力為零時的臨界頻率提前.

1 引 言

魚類游動是學(xué)者和工程師開發(fā)先進(jìn)高效推進(jìn)設(shè)備的重要靈感來源.例如, 仿生機(jī)器魚以類似游魚擺動尾鰭的方式獲得推力, 相比傳統(tǒng)螺旋槳推進(jìn)器, 具有機(jī)動性強(qiáng)、噪聲低及隱蔽性好等優(yōu)點(diǎn), 在海洋勘測和水下偵查等民用和軍用等領(lǐng)域具有重要的應(yīng)用價值.另外, 許多魚類趨向于成群游動,通過從鄰近同伴產(chǎn)生的渦流中獲得能量, 以提高游動效率和節(jié)省體力.因此, 研究游魚或魚群在水中的游動特性, 闡明魚尾鰭橫向擺動的推力產(chǎn)生機(jī)制和水動力相互作用機(jī)制, 有助于船舵、減搖鰭和仿生推進(jìn)器等海洋結(jié)構(gòu)物的設(shè)計和性能優(yōu)化, 并為小型自主潛航器或無人機(jī)群的最優(yōu)排布提供技術(shù)指導(dǎo).然而, 游魚或魚群游動是一個復(fù)雜的水動力現(xiàn)象, 涉及周期性的旋渦脫落, 渦流與魚體之間的能量交換以及游魚之間的水動力相互干涉作用.

隨著計算機(jī)硬件和數(shù)值算法的迅速發(fā)展, 計算流體力學(xué) (computational fluid dynamics, CFD)技術(shù)為魚類等柔性翼型結(jié)構(gòu)物的水動力特性研究提供了強(qiáng)有力的手段.然而, 數(shù)值模擬此流固耦合現(xiàn)象對于CFD研究是極具挑戰(zhàn)性的, 由于涉及復(fù)雜幾何形狀、界面大變形等難點(diǎn).典型的求解技術(shù)是貼體網(wǎng)格方法[1], 界面附近網(wǎng)格隨物體一起運(yùn)動以精確地描述運(yùn)動界面, 適合于模擬高雷諾數(shù)流動.然而, 貼體網(wǎng)格生成對經(jīng)驗(yàn)要求高, 擬合柔性多體運(yùn)動邊界涉及復(fù)雜的網(wǎng)格生成和重構(gòu)過程, 也會極大增加計算量和引入額外數(shù)值誤差.另一方面, 浸入邊界法[2]在近些年受到極大的關(guān)注.浸入邊界法在固定背景(一般為直角)網(wǎng)格上求解Navier-Stokes方程, 通過在流體動量方程中附加力源項(xiàng)[3]或重構(gòu)浸入邊界附近插值模板[4]以考慮浸入邊界對流場的影響.其網(wǎng)格生成簡單, 無需網(wǎng)格重構(gòu), 適于模擬復(fù)雜幾何形狀、多體動邊界流動問題.例如, Huang 等[5]采用懲罰浸入邊界法,Tullio 和Pascazio[6]采用直接力浸入邊界法研究柔性絲帶的流致振蕩特性, 吳曉笛等[7]采用基于浸入邊界的格子玻爾茲曼求解器模擬了雙圓柱兩自由度渦激振動現(xiàn)象.

針對魚類游動的水動力機(jī)理研究, 許多學(xué)者基于直角網(wǎng)格方法對繞單個或多個魚形水翼的流動問題進(jìn)行了數(shù)值模擬, 并取得卓有成效的進(jìn)展.Yeo等[8]基于混合直角網(wǎng)格方法研究了單個魚體的自主推進(jìn)性能和魚體大變形時的操縱性, 他們展示了仿生鯉魚在體長范圍內(nèi)執(zhí)行70°急轉(zhuǎn)彎的能力.Tian等[9]基于貼體動網(wǎng)格技術(shù)的空間時間有限元法模擬了一對母子魚在均勻流中的振蕩特性,研究了魚體間相對距離、雷諾數(shù)、相位差和相速度對力系數(shù)、能效和渦結(jié)構(gòu)的影響.Bergmann等[10]基于虛擬區(qū)域法模擬了單個和并列布置兩魚體的游動特性.Khalid 等[11,12]采用 Mittal等[13]的多維虛擬網(wǎng)格法研究了串列布置下非同步振蕩的兩游魚的水動力性能, 他們發(fā)現(xiàn)不同振蕩頻率的兩游魚產(chǎn)生非線性相互作用, 生成了新的水動力分量, 而且, 在下游魚體的導(dǎo)邊處尾渦會發(fā)生分裂, 導(dǎo)致上游魚體的阻力降低.國內(nèi)方面, 王亮和吳錘結(jié)[14]研究了三魚群游動以揭示“槽道效應(yīng)”在魚群游動中的節(jié)能機(jī)制, 王文全[15]研究了魚體自主游動的尾渦模式和魚體運(yùn)動特性.

目前, 學(xué)者對單魚體、串列或并列布置的雙魚體的流體特性進(jìn)行了大量研究.然而, 尚未知串列或并列布置時形成的力模式和尾渦結(jié)構(gòu)是否適用于魚群游動, 當(dāng)并列和串列布置同時存在時是否有新的力模式和尾渦結(jié)構(gòu)出現(xiàn).因此, 研究流體中以陣列布置的魚群游動能更好地再現(xiàn)魚群游動的水動力特性.鑒于此, 本文擬數(shù)值研究單個或陣列布置下魚形水翼在流場中的水動力特性, 以揭示魚體游動的推力產(chǎn)生機(jī)制和魚群游動的水動力優(yōu)勢.數(shù)值模型基于自主開發(fā)的徑向基函數(shù)虛擬網(wǎng)格法[16,17],在固定直角交錯網(wǎng)格上求解不可壓縮Navier-Stokes方程, 引入緊支徑向基函數(shù)到虛擬網(wǎng)格法以處理任意復(fù)雜的柔性邊界, 一個面積分?jǐn)?shù)表示方法通過提高浸入界面的局部質(zhì)量守恒性以抑制動邊界的壓力振蕩問題.采用該方法模擬了不同激勵頻率下以魚形游動的水翼在均勻流中的水動力特性,并進(jìn)行了網(wǎng)格收斂性的驗(yàn)證分析.然后, 研究了9個以陣列布置的魚形水翼在不同間距和不同激勵頻率下的力系數(shù)變化規(guī)律和尾渦形態(tài).

2 數(shù)學(xué)模型和數(shù)值方法

2.1 控制方程和數(shù)值離散

非定常黏性不可壓縮流體的控制方程如下:

其中 u 是速度向量, 具有直角網(wǎng)格分量 (u,v) , ρ 是密度, t 是時間, p 是壓力, ν 是流體的運(yùn)動黏性系數(shù), V和a是固體的速度和加速度; n是固體表面的外法向向量; Γ 表示固體表面.

流體控制方程(1)和方程(2)采用有限差分法在交錯直角網(wǎng)格上離散, 采用虛擬網(wǎng)格法施加浸入動邊界的無滑移邊界條件(3)式.在本文求解器中,以分步法解耦速度和壓力, 以二階Total Variation Diminishing Runge-Kutta(TVD-RK2)格式進(jìn)行時間積分, 采用時間顯式格式處理對流項(xiàng), 黏性項(xiàng)以 Crank–Nicolson 格式半隱式處理.對于空間離散格式, 采用TVD-MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)[18]格式離散對流項(xiàng), 黏性項(xiàng)以中心差分格式離散.采用不完全Cholesksy預(yù)處理的雙共軛梯度法 (Incomplete Cholesksy Biconjugate gradient stabilized method, BiCGSTAB)[19]求解壓力泊松方程離散形成的大型稀疏線性方程組.

2.2 虛擬網(wǎng)格法

為了考慮靜止或運(yùn)動物體對流場的影響, 本文采用虛擬網(wǎng)格法.該方法的基本理念是, 通過在浸入物體內(nèi)部布置有限數(shù)量的虛擬網(wǎng)格以表示靜止或運(yùn)動邊界的存在, 以適當(dāng)?shù)牟逯的K及技術(shù)(如雙線性插值方法), 施加浸入動邊界條件.

圖1 虛擬網(wǎng)格法的二維插值格式 (a)無虛擬網(wǎng)格; (b) 一個虛擬網(wǎng)格; (c) 兩個虛擬網(wǎng)格Fig.1.Two-dimensional interpolation stencil of the ghost cell method: (a) No ghost cells; (b) one ghost cell; (c) two ghost cells.

虛擬網(wǎng)格的插值重構(gòu)過程可細(xì)分為三步: 標(biāo)記鏡像點(diǎn)、插值重構(gòu)和滿足鏡像條件.首先, 對于每個虛擬網(wǎng)格(GC), 沿著物體表面的外法向從虛擬網(wǎng)格單元中心延伸到流體域以確定唯一的一個鏡像點(diǎn)(MP), 如圖1所示.然后是采用雙線性插值格式計算每個鏡像點(diǎn)的流體變量, 為了避免內(nèi)插模板點(diǎn)是虛擬網(wǎng)格自身的情況, 鏡像點(diǎn)的計算需要分如圖1中的三種情況進(jìn)行討論.最后是根據(jù)鏡像點(diǎn)變量通過滿足鏡像條件得到虛擬網(wǎng)格變量值, 更多細(xì)節(jié)參見文獻(xiàn)[17].需要注意的是由于交錯網(wǎng)格布置, 兩個方向速度分量或壓力的離散單元位置并不重合, 因此, 需要對每個方向的速度分量或壓力分別應(yīng)用以上的重構(gòu)過程.在整個計算域(液體和固體)求解Navier-Stokes方程, 邊界條件便得到隱式滿足.

2.3 CSRBF界面追蹤方法

浸入邊界法的網(wǎng)格生成過程簡單, 但由于浸入邊界本身并未包含在計算網(wǎng)格中, 精確表示背景網(wǎng)格中浸入邊界的存在是個難點(diǎn).本文采用基于CSRBF的隱式等值面方法重構(gòu)任意界面, 以光滑的基函數(shù)擬合離散的物面控制點(diǎn)構(gòu)造任意復(fù)雜物體的等值曲面, 根據(jù)等值面的特性識別場點(diǎn)屬性.本文CSRBF界面描述方法簡單、以有限的物面控制點(diǎn)獲得較高的精度, 適于描述不規(guī)則、大變形甚至凹形界面.

針對任意浸入界面, 基于有限數(shù)量的界面位置信息的離散控制點(diǎn), 引入徑向基函數(shù)構(gòu)造如下場函數(shù):

其中M為物面控制點(diǎn)的數(shù)量, βj為權(quán)重系數(shù),P(xi)為線性多項(xiàng)式.在文獻(xiàn) [16]中, 我們采用MQ (Multi-quadrics)核函數(shù)擬合復(fù)雜剛性界面.然而, 當(dāng)涉及復(fù)雜薄壁邊界時, 由于矩陣病態(tài), 擬合過程可能失敗.為了解決該問題, 本研究采用緊支 Wendland’s C2 函數(shù)[20], 得到的矩陣為正定稀疏矩陣, 矩陣的病態(tài)性得到緩解, 收斂性得到提高.Wendland’s C2 函數(shù)表達(dá)如下:

3 數(shù)值結(jié)果和討論

本文徑向基函數(shù)虛擬網(wǎng)格法的計算程序采用Fortran 90語言編寫, 數(shù)值算法的流程圖見圖2.以圓柱繞流、翼型繞流和圓柱振蕩等算例驗(yàn)證了本文方法的精度、可靠性和模擬剛性動邊界的能力[16,17].

3.1 魚形水翼在均勻流中的振蕩運(yùn)動

這里模擬了均勻流中柔性水翼的橫向振蕩運(yùn)動以驗(yàn)證本文方法模擬柔性動邊界流動的能力.水翼剖面為NACA0012機(jī)翼, 采用徑向基函數(shù)擬合141個樣品點(diǎn)進(jìn)行追蹤.為了模擬水中游魚的運(yùn)動,柔性水翼的橫向振蕩采用一個波浪函數(shù)描述:

圖2 靜止或運(yùn)動邊界流動求解算法的流程圖Fig.2.Flowchart of numerical algorithm for stationary or moving boundary flows.

y(x,t)是水翼表面質(zhì)點(diǎn)在y方向的位移, ω 是振蕩角頻率, 波數(shù) k =2π .振蕩幅值 A (x) 表示為A(x)=A1+A2x+A3x2, 其中系數(shù)分別為 A1=0.04 , A2=?0.16和 A3=0.32 .因此, 水翼運(yùn)動在導(dǎo)邊和隨邊的最大幅值分別為 0.04 和 0.2 m.雷諾數(shù) (Re)和Strouhal數(shù)(St)是兩個重要的流動參數(shù), Strouhal數(shù)定義為 St=2Amaxf/U∞, 其中 Amax是隨邊的最大幅值, f 是振蕩幅值 (f = ω/2π).根據(jù) Sui等[21]的工作, 雷諾數(shù)設(shè)置為 Re = 400, 選擇了五個振蕩角頻率 (ω = 0.15, 0.30, 0.45, 0.60 和 0.75)以研究振蕩頻率對力系數(shù)的影響, 分別對應(yīng)的Strouhal數(shù) 為 St= 0.191, 0.382, 0.573, 0.764 和0.955.計算區(qū)域設(shè)置為 [–3b, 9b] × [–3b, 3b], 翼弦長度 b = 1 m.入口邊界為均勻流, 均勻流速度 U∞=0.05 m/s, 在出口采用出流邊界條件, 上下邊界采用輻射邊界條件.對于壓力, 在出口指定壓力為零,其余邊界采用Neumann邊界條件, 選擇固定的時間步長 ? t = 0.002 s.

首先對激勵頻率 ω = 0.45時魚體振蕩運(yùn)動進(jìn)行了網(wǎng)格收斂性研究, 采用非均勻網(wǎng)格對局部區(qū)域 [–0.2b, 1.3b] × [–0.3b, 0.3b]進(jìn)行網(wǎng)格細(xì)化, 細(xì)化區(qū)域分別采用 150 × 80, 300 × 150, 600 × 300,對應(yīng)在x-方向的最小網(wǎng)格間距 ? x 分別為0.01b,0.005b和0.0025b, 對應(yīng)在y-方向的最小網(wǎng)格間距?y分別為 0.0075b, 0.004b 和 0.002b.表1 比較了本文三套網(wǎng)格下計算得到的力系數(shù)與Sui等[21]的參考結(jié)果.從表1知, 隨著網(wǎng)格細(xì)化, 平均阻力系數(shù)和均方根升力系數(shù)逐漸趨近Sui等[21]的Lattice Boltzmann方法結(jié)果, 而且, 兩套細(xì)網(wǎng)格上的結(jié)果基本重合, 說明本文方法獲得收斂解, 其中中間網(wǎng)格(300 × 150)上計算的平均阻力系數(shù)相比參考結(jié)果略大2.22%, 均方根升力系數(shù)略小4.17%.

鑒于300 × 150的中間網(wǎng)格上已計算得到相對合理的收斂解, 將300 × 150的中間網(wǎng)格用于其他工況的計算.圖3將五個頻率下的力系數(shù)與Sui等[21]的Lattice Boltzmann方法的計算結(jié)果進(jìn)行了比較.平均阻力系數(shù)和均方根升力系數(shù)與參考結(jié)果合理地符合.對于角頻率 ω ≤0.45 , 本文力系數(shù)與參考結(jié)果較好地符合, 當(dāng) ω >0.45 , 可觀察到力系數(shù)與參考結(jié)果的微小偏差.以角頻率ω=0.75為例, 本文平均阻力系數(shù)略大于參考結(jié)果7.1%, 均方根升力系數(shù)略小于參考結(jié)果3.3%.從圖3可觀察到, 隨著振蕩頻率的增加, 平均阻力系數(shù)迅速降低.臨界頻率為 ω =0.6 , 此時平均阻力系數(shù)從正值變?yōu)樨?fù)值, 即水翼在流體中運(yùn)動所遭受的阻力變?yōu)橥屏? 這是由于高頻振蕩產(chǎn)生的推力即為水中魚類向前運(yùn)動的動力來源.當(dāng)角頻率從0.6增加到0.75時, 推力系數(shù)從大約為0顯著增加到0.21.另外, 當(dāng)角頻率從0.15增加到0.75時, 均方根升力系數(shù)呈現(xiàn)近乎線性的增加趨勢(從0.21 到 3.1).

表1 本文力系數(shù)結(jié)果與 Sui等[21]的參考結(jié)果比較Table 1.Comparison of the force coefficients between the present method and Sui’s method[21].

圖3 五 個 頻 率 下 的 力 系 數(shù) 比 較 (a) 平 均 阻 力 系 數(shù) ;(b) 均方根升力系數(shù)Fig.3.Comparison of force coefficients under five oscillation frequencies: (a) Average drag coefficient; (b) root mean square lift coefficient.

圖4展示了三個角頻率下 (ω =0.15 , 0.45,0.75)升阻力系數(shù)隨無因次時間的變化趨勢.初始階段經(jīng)歷一些不穩(wěn)定的波動后, 三個頻率下的升力和阻力系數(shù)都呈現(xiàn)出周期性的變化趨勢.升力系數(shù)的振蕩周期是阻力系數(shù)的兩倍, 該特性與圓柱繞流的卡門渦街現(xiàn)象類似.隨著振蕩頻率從 ω =0.15 增加到 ω =0.75 , 平均阻力系數(shù)的平衡位置從0.36降低到–0.3左右, 但是升力系數(shù)的振蕩幅值從0.4顯著增加到4.3左右, 與圖3的平均阻力系數(shù)和均方根升力系數(shù)的變化趨勢一致.對于 ω =0.15 ,阻力和升力系數(shù)呈現(xiàn)小幅值和規(guī)則的正弦曲線變化.當(dāng)頻率增加到 ω =0.45 , 阻力系數(shù)曲線的波峰變得略微平坦, 波谷變得較為尖銳, 同時, 升力系數(shù)曲線的正弦曲線特性發(fā)生輕微偏斜.當(dāng)頻率進(jìn)一步增加到 ω =0.75 , 阻力系數(shù)觀察到更加平坦的波峰和更加尖銳的波谷, 類似于脈沖形狀, 而且, 升力系數(shù)曲線觀察到更明顯的波峰波谷偏斜現(xiàn)象.

圖4 三個頻率下升阻力系數(shù)隨無因次時間的變化趨勢(a) 阻力系數(shù); (b) 升力系數(shù)Fig.4.Variation of drag and lift coefficients with the dimensionless time under three oscillation frequencies: (a)Drag coefficient; (b) lift coefficient.

進(jìn)一步分析了柔性水翼振蕩的水動力特性.圖5—圖7分別展示了在三個振蕩頻率下相位角θ= π/2和 3 π/2 時的局部瞬時渦等值線, 根據(jù)(5)式,相位角 θ = π/2 和 3 π/2 分別對應(yīng)隨邊位移處于正的和負(fù)的最大位移時.由于柔性水翼周期性的擺動,沿著水翼表面的剪切流變得不穩(wěn)定, 一對渦從水翼表面交替脫落并沿著流動方向流向尾流遠(yuǎn)方.當(dāng)水翼尾部在前半個周期擺動到正的最大位移時, 一個順時針旋轉(zhuǎn)的渦從水翼上表面生成.當(dāng)水翼尾部在后半個周期擺動到負(fù)的最大位移時, 一個逆時針旋轉(zhuǎn)的渦從水翼下表面脫落.這對從水翼尾部交替脫落的反向旋轉(zhuǎn)的渦列就是卡門渦街, 因而產(chǎn)生周期性的升力系數(shù)變化.在這個渦列的尾流形成一個能量缺損區(qū), 由此產(chǎn)生一個正向阻力.隨著振蕩頻率的 增 加, 觀 察 到 ω =0.45 時尾渦 變 換為 逆 von-Karman渦街形態(tài), 此時, 在水翼尾流區(qū)產(chǎn)生一個能量增強(qiáng)區(qū), 阻力降低, 當(dāng)振蕩頻率增加到ω=0.6時, 阻力甚至降低到負(fù)值, 即為推力, 如圖3所示.在更高的振蕩頻率如 ω =0.75 時, 水翼快速的拍動運(yùn)動導(dǎo)致逆卡門渦街過渡為偏斜尾渦, 如圖5所示.另外, 隨著振蕩頻率增加, 附著在水翼尾部的渦列隨著水翼的加速運(yùn)動以更快的速度脫落, 如圖7 所示.在 ω =0.15 , 0.45, 0.75 時分別可觀察到兩對、四對和七對渦列, 當(dāng)水翼以高頻振蕩時(ω=0.75), 迅速脫落的渦列在尾流區(qū)緊致排列, 而且,渦之間發(fā)生相互排擠, 導(dǎo)致渦形狀的變形, 在尾流遠(yuǎn)方, 渦列甚至發(fā)生偏斜.

圖5 ω = 0.15 時在兩個時刻的局部瞬時渦等值線 (a) θ=/2; (b) θ =3π/2Fig.5.Local instantaneous vortex contours for ω = 0.15 at wo moments: (a) θ = π/2 ; (b) θ =3π/2 .

圖6 ω = 0.45 時在兩個時刻的局部瞬時渦等值線 (a) θ=π/2; (b) θ =3π/2Fig.6.Local instantaneous vortex contours for ω = 0.45 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

圖7 ω = 0.75 時在兩個時刻的局部瞬時渦等值線 (a) θ=π/2; (b) θ =3π/2Fig.7.Local instantaneous vortex contours for ω = 0.75 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

3.2 陣列布置魚形水翼在均勻流中的振蕩運(yùn)動

本節(jié)模擬了9個以陣列布置水翼在均勻流中的振蕩運(yùn)動以研究魚群游動中振蕩水翼表面的力模式和尾渦形態(tài).為了簡化計算模型, 研究聚焦于中心位置的水翼, 該水翼以(6)式指定的模式進(jìn)行振蕩運(yùn)動, 其余水翼保持靜止.幾何模型如圖8所示, 計算區(qū)域設(shè)置為 [–4b, 12b] × [–4b, 4b].水翼之間的橫向間距和縱向間距相同, 研究了三個間距下魚群的力系數(shù)特性, 間距 G分別為 G = 0.25,0.5 和 1 m.根據(jù) 3.1 節(jié)的網(wǎng)格收斂性結(jié)果, 采用均勻網(wǎng)格對水翼附近 [–2.5b, 2.5b] × [–1.5b, 1.5b]的區(qū)域進(jìn)行細(xì)化, 在x和y方向的最小網(wǎng)格間距分別為 0.005 和 0.004 m, 對應(yīng)于 3.1 節(jié)的收斂網(wǎng)格.邊界條件和計算參數(shù)的設(shè)置類似于3.1節(jié), 本節(jié)研究了五個振蕩角頻率 (ω = 0.15, 0.30, 0.45, 0.60和0.75)下中間位置水翼的力模式和尾渦形態(tài).

圖8 陣列布置水翼游動的幾何模型Fig.8.Geometrical model of the undulating hydrofoil in array arrangement.

首先對中間位置水翼表面的升阻力系數(shù)進(jìn)行了定量分析.圖9比較了三個位置下群游中間位置水翼和單水翼游動中振蕩頻率和水翼表面力系數(shù)之間的關(guān)系.在水翼間距 G = 1 m 時, 中間位置水翼的平均阻力系數(shù)和升力均方根系數(shù)與3.1節(jié)單水翼繞流算例近似, 由于距離遠(yuǎn), 鄰近水翼對中間位置水翼的影響較為微弱, 中間水翼保持相對獨(dú)立的游動模式.當(dāng)間距縮小為 G = 0.5 m 時, 水翼之間存在顯著的水動力相互作用, 力系數(shù)對振蕩頻率的敏感性增加, 當(dāng)振蕩頻率從ω = 0.15增加到ω =0.75 時, 平均阻力系數(shù)從 0.32 降低到–0.38, 升力均方根系數(shù)從0.41增加到4.11.當(dāng)縮小到最窄間距G = 0.25 m時, 此時水動力相互作用最為顯著,在振蕩頻率從 ω =0.15 增加到 ω =0.75 時, 平均阻力系數(shù)從0.32降低到–0.67, 升力均方根系數(shù)從0.75增加到10.10.值得注意的是, 相比單水翼繞流, 水翼表面阻力變?yōu)橥屏Φ呐R界頻率從大約0.60減小到0.45左右, 說明水翼群游時, 水翼可以以更少的能量消耗獲得同樣推力.在低振蕩頻率ω =0.15時, 三個間距下的水動力相互作用力系數(shù)變化不大, 因?yàn)榇祟l率下脫落的尾渦結(jié)構(gòu)相對穩(wěn)定獨(dú)立, 受到周圍水翼的影響小.在振蕩頻率 ω =0.75 時, G = 0.25 m 時群游中間水翼相比單水翼繞流可獲得2.5倍的平均推力增加, 升力均方根系數(shù)可獲得3倍增加, 說明增大激勵頻率可放大水翼之間的水動力相互作用, 也說明群游水翼在相同的運(yùn)動形式下能獲得更高的游動效率, 這反映出魚類群游的水動力學(xué)優(yōu)勢.

圖9 三個間距下群游中間位置水翼和單水翼在五個頻率下的力系數(shù)比較 (a)平均阻力系數(shù); (b)均方根升力系數(shù)Fig.9.Comparison of force coefficients between the central hydrofoil in array arrangement for three gap distances and the single hydrofoil under five frequencies: (a) Average drag coefficient; (b) root mean square lift coefficient.

圖10展示了G = 0.5 m時群游中間位置水翼游動中三個頻率下升力和阻力系數(shù)隨無因次時間的變化趨勢.群游水翼的升阻力系數(shù)變化趨勢與單水翼類似, 然而相比單水翼游動, 隨著激勵頻率的增加, 升力系數(shù)以更大的幅值和更大的偏斜呈現(xiàn)周期性的振蕩特性, 阻力系數(shù)呈現(xiàn)更加明顯的波谷平坦和波峰尖銳的現(xiàn)象.

圖10 群游中間位置水翼在三個頻率下升阻力系數(shù)隨無因次時間的變化趨勢 (a)阻力系數(shù); (b)升力系數(shù)Fig.10.Variation of drag and lift coefficients on the central hydrofoil in array arrangement with the dimensionless time under three oscillation frequencies: (a) Drag coefficient; (b) lift coefficient.

圖11—圖13展示了群游水翼在間距G =0.5 m時三個頻率下兩個時刻的局部渦場.中間位置水翼的游動特性和力系數(shù)受兩個因素影響.首先, 中間水翼的頭部位于上游水翼尾部的低壓區(qū),其尾部位于下游水翼頭部的高壓區(qū), 此壓差增加了中間水翼的推力; 其次, 中間水翼與其并列的上下兩水翼形成一個半封閉區(qū)域, 兩水翼之間的流體受到限制而沿著水翼游動方向流動, 減少了能量向兩側(cè)的耗散, 提高了推進(jìn)效率.當(dāng)中間水翼尾部向上擺動到最大幅值時(如圖11(a)), 中間水翼上表面的流體由于伯努利效應(yīng)而得到加速, 導(dǎo)致壓力降低, 下表面的流體減速, 導(dǎo)致壓力升高, 因此, 上下表面之間壓差形成向上的升力, 相比單水翼游動時增大了升力幅值, 如圖4所示; 當(dāng)中間水翼尾部向下擺動到最大幅值時(如圖11(b)), 中間魚體形成向下的升力.隨著Sr數(shù)增加, 魚體之間的劇烈水動力相互作用, 扭曲了周期性的對稱交替渦泄, 導(dǎo)致升力系數(shù)的正弦振蕩特性發(fā)生顯著的偏離, 如圖10中的升力系數(shù)曲線.而且, 由于激勵頻率增加, 相速度增加, 渦列以更快的速度脫落, 在下游魚體的阻礙作用下, 渦列沿著下游魚體上下表面發(fā)生分叉, 形成兩列緊密布置的渦列向遠(yuǎn)處, 如圖12和圖13所示.

圖11 ω = 0.15 時在兩時刻的水翼群游局部渦場圖(a) θ = π/2 ; (b) θ =3π/2Fig.11.Local instantaneous vortex contours of undulating hydrofoils in array arrangement for ω = 0.15 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

圖12 ω = 0.45 時在兩時刻的水翼群游局部渦場圖 (a) θ=π/2; (b) θ =3π/2Fig.12.Local instantaneous vortex contours of undulating hydrofoils in array arrangement for ω = 0.45 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

圖13 ω = 0.75 時在兩時刻的水翼群游局部渦場圖 (a) θ=π/2; (b) θ =3π/2Fig.13.Local instantaneous vortex contours of undulating hydrofoils in array arrangement for ω = 0.75 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

4 結(jié) 論

本文采用自主開發(fā)的徑向基函數(shù)虛擬網(wǎng)格法研究了繞單個或陣列布置柔性水翼的不可壓縮流動問題.采用虛擬網(wǎng)格法在浸入邊界內(nèi)部布置有限數(shù)量的虛擬網(wǎng)格以施加物面邊界條件, 引入緊支徑向基函數(shù)通過物面Lagrangian質(zhì)點(diǎn)擬合物面等值面函數(shù)以追蹤任意復(fù)雜的物面邊界.在繞單柔性水翼流動計算中, 網(wǎng)格收斂性研究表明當(dāng)網(wǎng)格細(xì)化時升阻力系數(shù)迅速趨近于參考結(jié)果, 同時, 不同激勵頻率下的升阻力系數(shù)與參考結(jié)果符合較好, 驗(yàn)證了本文方法的精度和可靠性.而且, 平均阻力系數(shù)隨振蕩頻率的增加迅速降低, 臨界頻率為 ω =0.6 , 此時水翼在流體中運(yùn)動所遭受的阻力變?yōu)橥屏? 即水翼高頻橫向振蕩產(chǎn)生了自身向前運(yùn)動的推力.同時, 在流場中觀察到尾渦交替脫落的卡門渦街現(xiàn)象, 隨著振蕩頻率增加, 渦列以更快的速度脫落,渦之間甚至發(fā)生相互排擠和變形.然后, 分析了不同間距和不同激勵頻率下陣列布置水翼的水動力特性, 隨著間距縮小, 力系數(shù)呈現(xiàn)顯著的放大效果,在間距 G = 0.25 m 時, 推力系數(shù)的臨界頻率減小為 ω =0.45 , 說明緊密布置的水翼在相同的運(yùn)動形式下能獲得更高的游動效率.