柔性翼仿生推進(jìn)特性數(shù)值模擬

李永成，鄭文濤

(中國(guó)船舶科學(xué)研究中心，江蘇 無(wú)錫 214082)

0 引 言

自然界中的魚類具備超凡的運(yùn)動(dòng)特性，如快速性、高機(jī)動(dòng)性等。柔性波動(dòng)方式更是魚類諸多游動(dòng)方式中推進(jìn)特性最優(yōu)的運(yùn)動(dòng)方式之一[1]。魚類采取柔性波動(dòng)方式運(yùn)動(dòng)時(shí)，整個(gè)身體采取橫向振動(dòng)變形，其典型代表如幔魚等。關(guān)于魚形體柔性波動(dòng)推進(jìn)的研究最早可以追溯到1970年英國(guó)Lighthill[2-3]提出的著名細(xì)長(zhǎng)體理論，該理論使得行波推進(jìn)運(yùn)動(dòng)方程在勢(shì)流領(lǐng)域首次擁有了數(shù)值解。童秉綱院士、陸夕云院士[4-6]等在此基礎(chǔ)上開創(chuàng)了三維行波板理論。隨著計(jì)算機(jī)的發(fā)展，邵雪明[7]、王志東[8]等開展了大量的魚形體粘性行波推進(jìn)研究工作。

上述研究工作大多針對(duì)于魚形體在低雷諾數(shù)下的推進(jìn)特性，而關(guān)于高雷諾數(shù)魚形體的推進(jìn)特性則少有涉及。不同雷諾數(shù)下魚形體的柔性推進(jìn)特性勢(shì)必存在差異，鑒于此，本文在前期研究[9]的基礎(chǔ)上開展非定常狀態(tài)下可變形翼高雷諾數(shù)推進(jìn)特性的數(shù)值模擬研究。對(duì)不同來(lái)流速度、變形位置及不同變形幅值組合形成的翼型進(jìn)行數(shù)值模擬，給出了推力特性以及周圍流場(chǎng)情況，為下一步進(jìn)行連續(xù)變形翼流體動(dòng)力特性試驗(yàn)研究提供初步思路。

1 數(shù)值計(jì)算方法

考慮不可壓縮流動(dòng)湍流問(wèn)題，其運(yùn)動(dòng)控制方程為：

本文的柔性翼行波推進(jìn)運(yùn)動(dòng)屬于典型的動(dòng)邊界問(wèn)題，變形翼擺動(dòng)推進(jìn)過(guò)程中，翼型本身時(shí)刻都在做扭曲變形運(yùn)動(dòng)，壁面附近計(jì)算網(wǎng)格的位置和形狀也必須時(shí)刻進(jìn)行相應(yīng)的調(diào)整，因此需要用到Fluent中的動(dòng)網(wǎng)格技術(shù)來(lái)捕捉變形翼的運(yùn)動(dòng)邊界。根據(jù)前期工作經(jīng)驗(yàn)，本文選取動(dòng)網(wǎng)格技術(shù)中的網(wǎng)格重構(gòu)技術(shù)對(duì)動(dòng)邊界進(jìn)行捕捉，具體動(dòng)網(wǎng)格技術(shù)見(jiàn)參考文獻(xiàn)[10-11]。

2 計(jì)算模型、網(wǎng)格及數(shù)值方法驗(yàn)證

2.1 計(jì)算模型

計(jì)算對(duì)象為二維的NACA0012翼型，具體如圖1所示。

圖1 計(jì)算模型示意圖Fig. 1 Diagram of computation model

可變形翼舵中弧線長(zhǎng)c= 0.2 m, 后緣變形起始位置xs，以翼型尾部任意一點(diǎn)a（初始橫縱坐標(biāo)分別為x0和y0，其中x0＞xs）為例，翼型變形運(yùn)動(dòng)方程如下式：

其中f為運(yùn)動(dòng)頻率，A為運(yùn)動(dòng)幅值。

一個(gè)運(yùn)動(dòng)周期內(nèi)變形翼的運(yùn)動(dòng)過(guò)程如圖2所示。其中T為運(yùn)動(dòng)周期，其與頻率f的關(guān)系為T= 1/f。

2.2 網(wǎng)格劃分與邊界條件

本文選取計(jì)算域大小如圖3所示。翼型前緣向前3c，設(shè)定來(lái)流速度的大小與方向；翼型向后6c，設(shè)定相對(duì)于參考?jí)毫c(diǎn)的流體靜壓值；翼型外表面，設(shè)定無(wú)滑移條件，同時(shí)加載UDF使其按照給定運(yùn)動(dòng)方程變形運(yùn)動(dòng)；上下邊界距離計(jì)算模型表面約4c，邊界條件亦設(shè)為速度入口條件，以模擬無(wú)界流場(chǎng)。

本文計(jì)算中使用的網(wǎng)格為結(jié)構(gòu)化網(wǎng)格。網(wǎng)格劃分基本原則為在模型表面附近網(wǎng)格加密，其中第一層網(wǎng)格間距根據(jù)y+確定（y+ 平均值約為1左右）。全局網(wǎng)格大約45萬(wàn)，網(wǎng)格無(wú)關(guān)性驗(yàn)證工作見(jiàn)參考本文獻(xiàn)[10]。

2.3 數(shù)值方法驗(yàn)證

開展數(shù)值方法精度驗(yàn)證工作，計(jì)算模型與參考文獻(xiàn)[12]一致，即二維NACA0012翼型。該翼型做拍動(dòng)運(yùn)動(dòng)，相關(guān)參數(shù)設(shè)置為：來(lái)流速度0.4 m/s，升沉幅值0.75C（C為弦長(zhǎng)）。計(jì)算翼型在頻率為 0.3下的推力系數(shù)和推進(jìn)效率，結(jié)果如圖4所示。

圖4 數(shù)值方法驗(yàn)證結(jié)果Fig. 4 Validation results of numerical method

由圖4可知，本文計(jì)算得到的拍動(dòng)翼推力系數(shù)和推進(jìn)效率與參考文獻(xiàn)相應(yīng)的結(jié)果近乎一致，最大誤差不超過(guò)4%，符合數(shù)值計(jì)算精度的要求。因此本文采用的定網(wǎng)格技術(shù)可用于數(shù)值求解柔性翼行波推進(jìn)運(yùn)動(dòng)。

2.4 時(shí)間步長(zhǎng)選取

考慮到時(shí)間步長(zhǎng)在非定常計(jì)算中起到極其重要的作用，本部分開展時(shí)間步長(zhǎng)選取驗(yàn)證工作。以變形起始位置s= 0.4c, 來(lái)流速度V= 1 m/s, 運(yùn)動(dòng)周期T= 0.5 s，運(yùn)動(dòng)幅值A(chǔ)= 0.021（對(duì)應(yīng)的θ值為10°）為例分別選取時(shí)間步長(zhǎng)Δt= 0.000 1 s, 0.001 s和 0.005 s進(jìn)行計(jì)算驗(yàn)證。計(jì)算得到的變形翼阻力系數(shù)隨時(shí)間變化曲線如圖5所示。

圖5 時(shí)間步長(zhǎng)驗(yàn)證結(jié)果Fig. 5 Validation results of time step

其中，CD表示阻力系數(shù)，其表達(dá)式如下式：

其中：D表示變形翼所受到的阻力值，ρ為水的密度，值為998.2 kg/m3。

從圖5可以看出，當(dāng)時(shí)間步長(zhǎng)減小至0.001 s后，計(jì)算得到的變形翼阻力系數(shù)與時(shí)間步長(zhǎng)為0.000 1 s下計(jì)算得到的阻力系數(shù)變化曲線近乎吻合。因此在后續(xù)計(jì)算中時(shí)間步長(zhǎng)Δt設(shè)為0.001 s。

3 計(jì)算結(jié)果及分析

3.1 柔性變形翼推進(jìn)特性分析

對(duì)不同來(lái)流速度、運(yùn)動(dòng)頻率和運(yùn)動(dòng)幅值下的變形翼推進(jìn)特性進(jìn)行數(shù)值模擬，共有45個(gè)不同組合的工況。具體工況為：來(lái)流速度設(shè)為V= 0.1 m/s, 0.5 m/s 和1.0 m/s運(yùn)動(dòng)幅值設(shè)為A= 0.021 m，0.043 m和0.069 m（對(duì)應(yīng)的θ值分別為10°，20°和30°），運(yùn)動(dòng)頻率f的變化范圍隨來(lái)流速度不同而相應(yīng)變化。具體表現(xiàn)為：V= 0.1 m/s時(shí)，f取值分別為0.4, 0.5, 0.6, 1.0和2.0 Hz；V= 0.5 m/s時(shí)，f取值分別為1.0, 1.5, 2.0, 2.5和3.0 Hz；V= 1.0 m/s時(shí)，f取值分別為2.5, 3.0, 3.5, 4.0和4.5 Hz。其原因在于，通過(guò)查閱文獻(xiàn)[1,4]可知，來(lái)流速度越大，克服流體阻力做功所需的運(yùn)動(dòng)頻率亦越大。因此為了達(dá)到研究運(yùn)動(dòng)參數(shù)對(duì)變形翼推進(jìn)性能的影響同時(shí)又兼顧計(jì)算效率，本文在低來(lái)流速度時(shí)，選取較低的運(yùn)動(dòng)頻率與之匹配，反之亦然。

在進(jìn)行計(jì)算結(jié)果分析前，首先對(duì)變形翼推力系數(shù)和推進(jìn)效率進(jìn)行定義。變形翼推進(jìn)效率實(shí)質(zhì)上為輸出功率與輸入功率的比值，輸出功率可以表征為平均推力與來(lái)流速度的乘積，輸入功率則可以表示為翼型表面速度壓力乘積沿著翼型邊界的積分，具體如下式：

其中：CT和?分別為推力系數(shù)和推進(jìn)效率；T為變形翼產(chǎn)生的推力；p和v分別為變形翼表面速度和壓力分布；Γ為變形翼邊界；為某一運(yùn)動(dòng)時(shí)刻。

上述計(jì)算工況下變形翼推力系數(shù)值如圖6所示。

圖6 不同計(jì)算工況下變形翼推力系數(shù)值Fig. 6 Thrust coefficient values of deformed wing under different calculation conditions

由圖6可以看出，變形翼的推力系數(shù)值隨著運(yùn)動(dòng)頻率、變形幅值的增大相應(yīng)增加，此外來(lái)流速度越大，系統(tǒng)所需克服流體阻力做的功也越大。具體來(lái)說(shuō)，以運(yùn)動(dòng)幅值A(chǔ)= 0.021 m為例，當(dāng)來(lái)流速度從0.1 m/s逐漸增大到0.5 m/s和1.0 m/s時(shí)，變形翼可產(chǎn)生正推力對(duì)應(yīng)的運(yùn)動(dòng)頻率分別為0.5 Hz, 2.5 Hz和4.0 Hz。

為了對(duì)不同運(yùn)動(dòng)參數(shù)組合下變形翼的推進(jìn)效率進(jìn)行準(zhǔn)確客觀的評(píng)價(jià)，引入一個(gè)無(wú)量綱參數(shù)斯特哈爾數(shù)，。具體表達(dá)式如下：

斯特哈爾數(shù)將變形翼的主要運(yùn)動(dòng)參數(shù)結(jié)合在一起，可用于準(zhǔn)確評(píng)價(jià)變形翼的推進(jìn)特性。不同斯特哈爾數(shù)下變形翼的推進(jìn)效率值如圖7所示。

圖7 不同St數(shù)下變形翼推進(jìn)效率Fig. 7 Propulsion efficiency of deformable wing at different St numbers

由圖7可以看出，不同來(lái)流速度下變形翼的推進(jìn)效率值隨著St數(shù)的增加呈現(xiàn)先增大后減小的趨勢(shì)，存在最佳St數(shù)對(duì)應(yīng)最大推進(jìn)效率。不同來(lái)流速度下最大推進(jìn)效率對(duì)應(yīng)的最佳St值位于0.25～0.40之間。這一結(jié)論與參考文獻(xiàn)中[12]給出的最佳St=0.30較為吻合。

3.2 柔性變形翼流場(chǎng)結(jié)構(gòu)分析

上文對(duì)不同計(jì)算工況下的變形翼運(yùn)動(dòng)參數(shù)推進(jìn)特性進(jìn)行了研究, 本部分將對(duì)變形翼的流場(chǎng)結(jié)構(gòu)進(jìn)行分析，對(duì)上述計(jì)算結(jié)果進(jìn)行解釋。圖8為不同St數(shù)下變形翼周圍的流場(chǎng)結(jié)構(gòu)以及周圍壓力分布圖。其中Z為Z方向的渦量，其表達(dá)式如下式：

圖8 不同St數(shù)下變形翼尾流場(chǎng)渦量圖Fig. 8 Vorticity diagram of the wake flow field of the deformed wing at different St numbers

其中，u和v分別為水平方向和垂直方向的速度分量。

由圖8可以看出，不同St數(shù)變形翼尾流場(chǎng)中呈現(xiàn)交錯(cuò)分布的渦環(huán)結(jié)構(gòu)。其產(chǎn)生的原因在于，由于變形翼上下邊界的運(yùn)動(dòng)使得翼型表面的流動(dòng)分離，并且在來(lái)流的作用下向后移動(dòng)，從而呈現(xiàn)出交錯(cuò)分布的渦環(huán)結(jié)構(gòu)。注意到不同St數(shù)下變形翼流場(chǎng)渦環(huán)結(jié)構(gòu)不同之處在于在低St數(shù)下（St=0.10）,變形翼尾流場(chǎng)中的渦環(huán)結(jié)構(gòu)呈現(xiàn)卡門渦街形式（BvK）, 即翼型上緣脫落的渦環(huán)在尾流場(chǎng)中呈現(xiàn)順時(shí)針旋轉(zhuǎn)，而下緣脫落的渦環(huán)在尾流場(chǎng)中則呈現(xiàn)逆時(shí)針旋轉(zhuǎn)；隨著St數(shù)的增加（St=0.35, 0.80），變形翼尾流場(chǎng)中則呈現(xiàn)反卡門渦街（rB-vK），即即翼型上緣脫落的渦環(huán)在尾流場(chǎng)中呈現(xiàn)逆時(shí)針旋轉(zhuǎn)，而下緣脫落的渦環(huán)在尾流場(chǎng)中則呈現(xiàn)順時(shí)針旋轉(zhuǎn)。如此二者之間相互作用形成一股向后的射流，根據(jù)動(dòng)量定理，向后的射流亦給一個(gè)作用于向前的推力，如此變形翼便可以向前運(yùn)動(dòng)[17]。

4 結(jié) 語(yǔ)

本文通過(guò)開展不同來(lái)流速度、不同運(yùn)動(dòng)幅值以及不同運(yùn)動(dòng)頻率的數(shù)值計(jì)算，獲得變形翼推力系數(shù)和推進(jìn)效率隨運(yùn)動(dòng)參數(shù)的變化規(guī)律。此外，通過(guò)比較分析不同St數(shù)下變形翼周圍渦流場(chǎng)結(jié)構(gòu)，定性地從流動(dòng)機(jī)理的角度對(duì)上述計(jì)算結(jié)果進(jìn)行闡釋。本文工作可為后續(xù)試驗(yàn)開展提供技術(shù)依據(jù)和理論指導(dǎo)。