舒蔣鵬, 石小濤, 陳小龍, 周柳青, 王思瑩, 胡曉, 張立勝*

(1.三峽大學(xué) 湖北省魚類過壩技術(shù)國際科技合作基地, 湖北 宜昌 443002;2.三峽大學(xué) 水利與環(huán)境學(xué)院, 湖北 宜昌 443002;3.武漢理工大學(xué) 新材料力學(xué)理論與應(yīng)用湖北省重點實驗室, 湖北 武漢 430063)

0 引言

自然界中的魚類經(jīng)過數(shù)十億年的自然選擇,已經(jīng)具備了高機動性、低功耗的游動能力[1],其游動能力是仿魚類水下機器人所無法比擬的。為了探究魚類高效的游動性能,許多學(xué)者對此展開了大量研究。前期研究主要通過理論和實驗的方法對仿生尾鰭開展研究。Lighthill[2]研究了尾鰭的擺動,并提出了大擺幅細長體理論。Freymuth[3]通過煙線法研究擺動水翼在風(fēng)洞中的流場變化。Triantafyllou等[4]研究發(fā)現(xiàn),當(dāng)擺動水翼在斯特勞哈爾數(shù)St為0.25～0.35時,其推進效率最大。Anderson等[5]和Read等[6]通過實驗研究發(fā)現(xiàn),擺動水翼在特定頻率下和特定攻角下的推進效率最高。Lauder等[7]用數(shù)字粒子圖像測速(digital particle image velocimetry, DPIV)技術(shù)測量了擺動水翼的尾部流場,分析了擺動水翼運動參數(shù)對尾渦結(jié)構(gòu)的影響。王肇等[8-9]通過水洞試驗研究了擺動水翼的水動力性能。Schouveiler等[10]通過實驗研究Sr和攻角對擺動水翼推進性能的影響。Iverson等[11]在循環(huán)水道內(nèi)對擺動水翼的力進行直接測量,并結(jié)合粒子圖像測速(particle image velocimetry, PIV)技術(shù)對水翼尾流的渦結(jié)構(gòu)進行評估。近些年,計算流體力學(xué)(computational fluid dynamics,CFD)方法的發(fā)展彌補了理論和實驗的不足,可以對流場結(jié)構(gòu)變化、游動機制以及推進性能等進一步深入研究。張曉慶等[12]采用數(shù)值方法模擬了二維剛性和柔性水翼并與實驗值進行對比,結(jié)果較好。楊亮[13]采用雷諾平均的NS方程(reynolds averaged navier strokes, RANS)求解器和動網(wǎng)格方法研究了運動參數(shù)對擺動水翼水動力性能的影響。文敏華等[14]基于動網(wǎng)格技術(shù)分析了擺動水翼推力中黏性力和壓差力的占比情況。劉煥興[15]采用RANS求解器研究了Sr和最大攻角對擺動水翼推進性能的影響。

上述研究均將擺動水翼放置在均勻流場中,而魚類真實的生存環(huán)境常常存在各種漩渦,其流場較為復(fù)雜,并且擺動水翼作為魚類尾鰭的一種簡化模型,研究流場中渦對擺動水翼的推進性能影響至關(guān)重要。另外有研究指出,魚類可以從漩渦中吸收能量,提高其推進性能。為此本文中基于CFD,利用Fluent的動網(wǎng)格技術(shù)和用戶自定義函數(shù)(user define functions,UDF),通過在擺動水翼前方放置一個D形柱模擬復(fù)雜流場,分析擺動水翼的升沉振幅和擺動頻率對擺動水翼推進性能的影響,并與均勻流場中的擺動水翼對比,所得結(jié)果可為仿魚類水下機器人的研制提供理論參考。

1 數(shù)值計算模型

1.1 數(shù)值計算模型

本文定義2種運動模式:模式1為在均勻流場下的擺動水翼俯仰-升沉耦合運動;模式2為D形柱尾流場(復(fù)雜流場)下的擺動水翼俯仰-升沉耦合運動。擺動水翼前緣與D形柱圓心距離S=4d,d表示D形柱直徑,d=c/2(c為弦長)。模式2示意圖如圖1所示。二維擺動水翼采用NACA0012水翼,樞軸點為c/3處,其滿足的運動規(guī)律[16]可以表示為

圖1 模式2示意圖Fig.1 Schematic of mode 2

y(t)=h0sin(2πft),

(1)

θ(t)=θ0sin(2πft+φ),

(2)

式中:y(t)為升沉速度;θ(t)為俯仰速度;h0為升沉運動振幅;θ0為俯仰運動振幅;φ表示耦合運動相位差;f為水翼運動頻率;t為時間。有研究表明擺動水翼在φ=-90°時推進性能最優(yōu),故本研究φ值恒為-90°。D形柱和擺動水翼的周期性運動圖如圖2所示。

圖2 D形柱和擺動水翼的周期性運動圖Fig.2 Periodic motion diagram of D-cylinder and flapping hydrofoil

1.2 水動力系數(shù)

對于擺動水翼運動,定義水翼斯特勞哈爾數(shù)St=2fh0/U,式中U為來流速度[17]。水翼的瞬時推力系數(shù)Ct、瞬時升力系數(shù)Cy和瞬時力矩系數(shù)Cm分別為

(3)

式中:Fx(t)、Fy(t)、Mz(t)分別是擺動水翼的瞬時推力、瞬時升力、瞬時力矩;ρ為水的密度。

相應(yīng)地,水翼一個運動周期內(nèi)的平均推力系數(shù)CT、平均升力系數(shù)CL、平均力矩系數(shù)CM分別為

(4)

式中T為整個運動周期。

水翼的推進效率η為[18]

(5)

式中CP為一個周期內(nèi)的水翼平均輸入功率系數(shù),其表達式為

(6)

(7)

2 數(shù)值計算方法

2.1 控制方程

本文采用二維不可壓縮黏性納維-斯托克斯(N-S)方程作為流體的控制方程,可以表示為

·u=0,

(8)

(9)

式中:u是流體運動速度;p是流體壓力;ν是流體動力學(xué)黏性系數(shù);湍流模型選擇SSTk-ω兩方程模型。

2.2 計算域網(wǎng)格劃分及參數(shù)設(shè)定

計算域網(wǎng)格劃分及參數(shù)設(shè)定如圖3所示。計算域的尺寸為10c×8c(長度×寬度)。計算域左側(cè)為速度入口,右側(cè)為壓力出口,上下為對稱邊界,擺動水翼邊界設(shè)為無滑移壁面。為了提高計算效率,采用結(jié)構(gòu)化網(wǎng)格跟非結(jié)構(gòu)化網(wǎng)格混合的方法進行網(wǎng)格劃分。計算域分為左、右域,交界處設(shè)為Interface,又將左、右2個流體區(qū)域分為內(nèi)、外域,內(nèi)、外域交界處設(shè)為Interior,內(nèi)域分別為半徑為c和半徑為2c的圓形區(qū)域,采用非結(jié)構(gòu)化四邊形網(wǎng)格劃分,整個計算域除去內(nèi)部計算域的區(qū)域為外域,采用結(jié)構(gòu)化三角形網(wǎng)格劃分,并采用彈簧光順和網(wǎng)格重構(gòu)法對畸變較大的網(wǎng)格進行重新劃分以確保網(wǎng)格的合格性。

圖3 計算域網(wǎng)格劃分及參數(shù)設(shè)定Fig.3 Computing domain grid division and parameter setting

采用SIMPLE算法對連續(xù)方程中的壓力和速度進行耦合,為了提高計算精度,離散方式中選擇時間項采用一階隱式,動量項采用二階迎風(fēng)格式,壓力項選擇二階格式。

2.3 網(wǎng)格無關(guān)性驗證

為了驗證網(wǎng)格的無關(guān)性,以模式1為驗證算例進行網(wǎng)格劃分,劃分3種不同數(shù)量的網(wǎng)格,分別為8.95×104、12.10×104、16.10×104。時間步長選擇Δt=T/200,運動參數(shù)選擇,h0/c=1.0,St=0.25,θ0=25°,f=0.5 Hz,φ=-90°,不同網(wǎng)格數(shù)量時的CT值和相對誤差見表1。從表1中可以看出,網(wǎng)格B的相對誤差最小,因此為了節(jié)省計算時間,本文后續(xù)計算采用網(wǎng)格B。

表1 網(wǎng)格無關(guān)性驗證Tab.1 Mesh agnostic validation

2.4 模型有效性驗證

為了驗證數(shù)值方法的有效性,數(shù)值計算了斯特勞哈爾數(shù)St分別是0.15、0.20、0.25、0.30、0.35、0.40,最大攻角amax=15°,相位差φ=-90°,h0/c=1.0,雷諾數(shù)Re=40 000下的CT值,并將計算值與文獻[6]實驗結(jié)果進行比較,結(jié)果如圖4所示。由圖4可見,平均推力系數(shù)CT的計算值與實驗值吻合良好,具有相同的變化趨勢,說明該數(shù)值方法是有效的。

圖4 平均推力系數(shù)CT計算值與文獻[6]結(jié)果的對比Fig.4 Comparison between the calculated values in this paper and those in Reference[6]

3 計算結(jié)果與分析

本文主要研究擺動頻率f和升沉振幅h0對2種模式下擺動水翼的推進性能影響以及流場變化。

3.1 擺動頻率對擺動水翼推進性能的影響

設(shè)置運動參數(shù):f分別取0.5c、1.0c、1.5c、2.0c、2.5cHz,h0=0.5c,θ0=30°,U=2.0c。

2種模式下的擺動水翼的平均推力系數(shù)CT和增長率δ如圖5所示。由圖5可以看出,2種模式下擺動水翼的CT均隨著f的增加而增大,并且模式2在任何f下始終最大,說明D形柱產(chǎn)生的尾流場可以提升擺動水翼的推進性能。此外,從圖5可以看出,隨著f的增加,δ呈先減小后增大再減小的趨勢,f=0.5cHz時,δ=37.2%,f=1.0cHz時,δ=2.13%,研究發(fā)現(xiàn),與增大f所帶來的對擺動水翼推力改善相比,D形柱引起的尾流場在擺動水翼推力改善中所占的百分比先減小后增大,然后再減小。

圖5 擺動水翼的平均推力系數(shù)CT及增長率δFig.5 Average thrust coefficient CT and growth rate δ of flapping hydrofoil

2種模式下的擺動水翼的平均輸入功率系數(shù)CP和推進效率η如圖6所示。從圖6可以看出,擺動水翼的CP隨著f的增加而增加,而η隨著f的增加而減小。當(dāng)f<1.5cHz時,2種模式下的CP值沒有明顯區(qū)別;當(dāng)f>1.5cHz后,兩者差值逐漸變大,說明D形柱的存在能夠增加擺動水翼的功率消耗。此外,2種模式下擺動水翼的η隨著f的增加而減小,在相同f下,模式2的η始終比模式1的大;在f=0.5cHz時,模式2的η比模式1的η大29.5%;在f=2.5cHz時,模式2的η比模式1的η大4.9%,說明D形柱的存在可以提升擺動水翼的推進效率,但增加效果是逐漸降低。

圖6 擺動水翼的平均輸入功率系數(shù)CP和推進效率ηFig.6 Average input power coefficient CP and propulsive efficiency η of flapping hydrofoil

3.2 升沉振幅對擺動水翼推進性能的影響

設(shè)置運動參數(shù):h0分別取0.5c、1.0c、1.5c、2.0c、2.5c,h0=d,f=1.0 Hz,θ0=30°,U=2.0c。

2種模式下的擺動水翼的平均推力系數(shù)CT和增長率δ如圖7所示。從圖7中可以看出,2種模式下擺動水翼的CT的變化特征,與3.1節(jié)中的變化特征一致,不再進行敘述。此外,從圖7可以看出,隨著h0的增加,δ呈先增大后減小再增大的趨勢,在h0/c=0.5時,有最小增量δ=2.1%,在h0/c=2.5時,有最大增量δ=8.2%,表明與增大h0對擺動水翼推力改善相比,D形柱引起的尾流場在擺動水翼推力改善的增長率是先增大后減小再增大,與3.1節(jié)中的結(jié)果正好相反。

圖7 擺動水翼的平均推力系數(shù)CT及增長率δFig.7 Average thrust coefficient CTof flapping hydrofoil

2種模式下的擺動水翼的平均輸入功率系數(shù)CP和推進效率η如圖8所示。從圖8可以看出,擺動水翼的CP均隨著h0的增加而增加,與3.1節(jié)的變化規(guī)律一致,但值得注意的是,在h0>1.5c時,模式1的CP值均在模式2的上方,說明D形柱的存在能夠減少擺動水翼的功率消耗。從圖8還可以看出,2種模式下擺動水翼的η隨著h0的增加而減小,在相同f下,模式2的η始終比模式1的大;在h0/c=0.5時,模式2的η比模式1的η大0.9%;在h0/c=2.0時,模式2的η比模式1的η大26.1%;在h0/c=2.5時,模式2的η比模式1的η大15.4%,表明D形柱的存在可以提升擺動水翼的推進效率,總體呈先增大后減小的趨勢,在h0/c=2.0時增加最明顯。

圖8 擺動水翼的平均輸入功率系數(shù)CP和推進效率ηFig.8 Average input power coefficient CP and propulsive efficiency η of flapping hydrofoil

3.3 擺動水翼的流場變化

分析對比2種模式下從D形柱脫落于流場中的渦對擺動水翼流場結(jié)構(gòu)的影響。考慮參數(shù)取值分別為f=1.0cHz,h0=1.5c,θ0=30°,U=2.0c。

2種模式下,擺動水翼在一個周期下的流場壓力分布云圖如圖9所示。由圖9可見,當(dāng)t/T=0時,此時擺動水翼的擺角最大,擺動水翼的上、下表面均覆蓋高壓區(qū)和低壓區(qū),值得注意的是,在模式1和模式2中擺動水翼擺角為最大時(t/T=0),擺動水翼下表面靠近前緣的區(qū)域存在一個明顯凸起的強度較大的低壓區(qū),并且模式2的區(qū)域比模式1的更大。隨著擺動水翼運動過程中擺角的減小,這個凸起的低壓區(qū)向擺動水翼的尾緣移動,并在擺動水翼擺角為0時(t/T=1/4)從擺動水翼下表面分離。由于擺動水翼在一個運動周期的對稱性,因此在t/T=1/2時,擺動水翼的上表面靠近前緣的區(qū)域也存在一個明顯凸起的低壓區(qū),此時擺動水翼上表面為低壓區(qū),下表面為高壓區(qū),與擺動水翼在t/T=0時的情況相反。模式1和模式2中擺動水翼上、下兩側(cè)的壓力分布區(qū)域一致,不同之處在于模式2中壓力的幅值更大。擺動水翼兩側(cè)壓力差的變化決定擺動水翼的推進性能,又因為2種模式下擺動水翼的擺角相同,所以模式1和模式2的擺動水翼的平均推力系數(shù)逐漸增大。

(a) 模式1

(b) 模式2圖9 流場壓力分布云圖Fig.9 Flow field pressure distribution cloud

2種模式下,擺動水翼在一個周期下的擺動水翼尾流場的渦結(jié)構(gòu)云圖如圖10所示。設(shè)ωz為流場渦量在z方向的分量,其計算公式[19]為

(a) 模式1

(b) 模式2圖10 擺動水翼尾流場的渦結(jié)構(gòu)云圖Fig.10 Vortex structure cloud image of the wake field of a flapping hydrofoil

(10)

式中:ux、uy分別表示流場在x軸和y軸方向的速度分量。

由圖10可知,擺動水翼在2種運動模式下的尾流場均為反卡門渦街形式。對于模式1,當(dāng)擺動水翼處于最大擺角時(t/T=0),尾緣處存在一個尾緣渦M-1(順時針),并與旋向相同的附著渦LM-1相連接,上端存在一個FM-1渦,與水翼的下表面前緣渦TM-1旋向相同。當(dāng)水翼擺動至平衡位置時(t/T=1/4),此時擺角為0,TM-1和M-1不斷擴大,但M-1并未脫落。此后,擺動水翼反向運動到最大擺角時(t/T=1/2),M-1脫落至尾流場中,形成一個高強度的渦核,TM-1不斷向尾緣移動,并與LM-2相連接,此刻,水翼上、下表面存在旋向相同的TM-2和FM-2渦。當(dāng)水翼反向擺動至平衡位置時(t/T=3/4),TM-1和TM-2不斷擴大,但并未完全從尾緣脫落,最后水翼運動至初始位置,TM-1發(fā)生脫落,與M-1形成旋向相反的渦對。

對于模式2,擺動水翼前方存在D形柱,D形柱不斷脫離2個旋向相反的渦1和渦2,擺動水翼的擺角最大時(t/T=0),D形柱產(chǎn)生的順時針渦1靠近水翼前緣,此時水翼尾緣處存在一個尾緣渦M-1(順時針),上端存在一個FM-1渦,與水翼的下表面前緣渦TM-1旋向相同。當(dāng)水翼擺動至平衡位置時(t/T=1/4),TM-1不斷擴大,此時D形柱產(chǎn)生的順時針渦1與TM-1相互作用,在D形柱產(chǎn)生的順時針渦1的擠壓下,相同時刻下,TM-1渦明顯大于模式1的,M-1在尾緣的帶動下也逐漸拉長擴大,但并未完全從尾緣脫落。此后,擺動水翼反向運動到最大擺角時(t/T=1/2),D形柱產(chǎn)生的逆時針2渦靠近水翼前緣,TM-1不斷擴大,擺動水翼上表面存在逆時針渦TM-2,水翼下表面存在一個FM-2渦,此時M-1完全脫落于尾流場中,形成一個強度高度集中的渦核,與在相同時刻下的模式1中的M-1相比耗散得更快。當(dāng)水翼反向擺動至平衡位置時(t/T=3/4),此時擺角為0,TM-2不斷擴大,TM-1在尾緣的帶動下也逐漸拉長擴大,但并未完全從尾緣脫落,此時D形柱產(chǎn)生的逆時針渦2與TM-2相互作用,在D形柱產(chǎn)生的順時針渦2的擠壓下,相同時刻時,TM-2渦明顯大于模式1的,最后水翼運動至初始位置,TM-1發(fā)生脫落,與M-1形成旋向相反的渦對。D形柱產(chǎn)生的渦1和渦2對擺動水翼的前緣渦TM-1和TM-2產(chǎn)生擠壓作用,使得擺動水翼的尾渦區(qū)域增大,脫落的渦強度更大,渦耗散得更快,渦對之間誘導(dǎo)作用更強,產(chǎn)生更大的推力。

4 結(jié)論

本文基于CFD方法并結(jié)合Fluent的動網(wǎng)格技術(shù),通過改變擺動水翼的擺動頻率和升沉振幅計算擺動水翼在均勻流場下進行俯仰-升沉運動(模式1)和擺動水翼在D形柱尾流場中進行俯仰-升沉運動(模式2)中的推進性能,并結(jié)合流場的壓力分布云圖和渦結(jié)構(gòu)變化分析了D形柱對擺動水翼的水動力性能影響機制,得到了以下結(jié)論:

① 2種模式下,擺動水翼的平均推力系數(shù)和平均輸入功率系數(shù)均隨著擺動頻率和升沉振幅的增加而不斷增加,而推進效率則相反。

② 2種模式下,隨著擺動頻率的增加,D形柱尾流場中的擺動水翼的平均推力系數(shù)始終大于均勻流場中擺動水翼的平均推力系數(shù),并且增長率是先減小后增大再減小。

③ 2種模式下,隨著升沉振幅的增加,D形柱尾流場中的擺動水翼的平均推力系數(shù)始終大于均勻流場中擺動水翼的平均推力系數(shù),并且增長率是先增大后減小再增大。

④ 2種模式下,擺動水翼尾流場均為反卡門渦街形式,但D形柱尾流場中,D形柱脫落的渦對水翼前緣渦發(fā)生擠壓作用,使擺動水翼的尾渦區(qū)域增大,脫落的渦耗散得更快,渦對之間誘導(dǎo)作用更強,產(chǎn)生更大的推力。