小河道撲動水翼裝置推水流動特性研究

華爾天 陳萬前 湯守偉 謝榮盛 郭曉梅 徐高歡

(1.浙江工業(yè)大學機械工程學院， 杭州 310023； 2.先進水利裝備浙江省工程研究中心， 杭州 310018；3.浙江水利水電學院機械與汽車工程學院， 杭州 310018)

0 引言

平原河網(wǎng)地區(qū)地勢平緩、水體流速較小，水動力不足引發(fā)的水體污染對居民生活質(zhì)量與身體健康造成了巨大的影響[1]。調(diào)水引流作為常用的治理措施之一，能夠提升河網(wǎng)水動力，有效改善水環(huán)境質(zhì)量[2-4],但對分流能力小的支杈的水動力狀況幫助較小[5]。而且降低轉(zhuǎn)速、增大葉輪直徑雖然可以改善泵站在超低揚程下運行效率低、穩(wěn)定性差[6-8]等問題，卻依舊無法適應平原小河道揚程幾乎為零的情況。文獻[9-10]模擬魚尾擺動形式研制的仿尾鰭式微型壓電泵，具有結(jié)構(gòu)簡單、流量大、揚程低等優(yōu)點，驗證了魚尾擺動形式對水體的推動作用。受此啟發(fā)，本文提出利用撲動水翼進行水體推動的方法，以提升平原小河道的水動力條件。

撲動水翼是由魚類游動方式簡化而來的，具有阻力小、噪聲低、效率高、機動性好等特點。自文獻[11]提出撲動水翼可以作為一種新式螺旋槳以來，撲動水翼的推進性能成為了研究熱點[12-15]，撲動水翼的應用范圍也被不斷擴展，比如利用串聯(lián)水翼從非恒定場中獲得推進力的波浪滑翔機[16-17]，以及利用水翼從波浪中獲取能量的海洋能量采集裝置等[18-19]。文獻[20]對比了撲動、擺動以及純俯仰3種運動形式，發(fā)現(xiàn)3種方式在合適的參數(shù)下均能在尾流中形成雙列反卡門渦街；文獻[21]則對不同運動參數(shù)下的撲動水翼進行了研究，分析了升沉以及俯仰之間的相位差、攻角的偏置以及運動方程中高次項對撲動水翼推進性能的影響；文獻[22]研究了斯特勞哈爾數(shù)和最大攻角的變化對推力和流體力學效率的影響，確定了在一定的參數(shù)范圍內(nèi)能達成高效率與高推力的共存，并且試驗中測得的效率峰值超過了70%。文獻[23]對常規(guī)正弦運動的4種組合方式進行了模擬分析，通過分析各形式下的尾流結(jié)構(gòu)以及升阻系數(shù)，提出俯仰的相位比升沉超前90°時推進效果最佳。文獻[24]利用改變撲動水翼的運動形式的方法增加了撲動水翼的平均迎角，有效提升了撲動水翼在較大斯特勞哈爾數(shù)St時的推力以及效率。上述研究均表明在合適的運動參數(shù)下，撲動水翼能夠從尾流的雙列反卡門渦街獲得強大的推進力，有著優(yōu)異的水動力性能。

以上對撲動水翼的研究主要集中于撲動水翼的推進性能，關于撲動水翼對水體的推動效果的研究較少。為了研究撲動水翼的推水性能，本文利用CFD方法建立撲動水翼的二維計算模型，對不同頻率以及不同來流速度情況下?lián)鋭铀淼耐扑阅苓M行模擬研究，并搭建撲動水翼試驗裝置對模擬結(jié)果進行驗證。

1 運動模型

1.1 運動描述

水翼的研究分為柔性水翼以及剛性水翼兩種，其中柔性水翼自由度較多，運動復雜，且其優(yōu)勢在較高撲動頻率區(qū)域；而剛性水翼的運動更為簡單、便于控制、成本低廉，有著更為實際的工程意義[25]，因此本文選擇剛性水翼為研究對象。撲動水翼的運動由升沉與俯仰兩種運動結(jié)合而成，根據(jù)文獻[23]選擇常規(guī)正弦運動，如圖1所示。圖中，Amax表示撲動水翼的升沉幅值；θmax表示撲動水翼的俯仰幅值；T表示運動周期。

圖1 撲動水翼的撲動形式示意圖Fig.1 Schematic of flapping form of hydrofoil

撲動水翼的基本運動方程為

(1)

式中y(t)——撲動水翼升沉位移，m

θ(t)——撲動水翼俯仰位移，rad

ω——撲動角頻率，rad/s

φ——升沉與俯仰運動的相位差，rad

對式(1)求導，得撲動水翼的速度方程為

(2)

1.2 運動參數(shù)與力學模型

在撲動水翼的研究中，瞬時推力系數(shù)Ct以及瞬時升力系數(shù)Cy是衡量撲動水翼水動力性能的關鍵參數(shù)，計算公式分別為

(3)

式中Ft(t)——水平方向瞬時推力，N

Fy(t)——豎直方向瞬時升力，N

ρ——流體密度，kg/m3

s——翼型的展長，m

平均推力系數(shù)及平均升力系數(shù)定義為

(4)

式中n——計算總周期數(shù)

k——取值周期數(shù)，取2

除此以外，為表征撲動水翼的推水性能，還需要計算撲動水翼裝置的流量、揚程以及效率。

流量公式定義為

(5)

b——流道寬度，取0.8 m

平均揚程可由出入口的壓力差換算獲得，計算公式為

(6)

g——重力加速度，m/s2

撲動水翼運動的平均輸入功率計算公式為

(7)

M(t)——翼型繞轉(zhuǎn)軸的瞬時扭矩，N·m

(8)

水力推進效率ηpr則可以表示為

(9)

式中CP——瞬時功率系數(shù)

(10)

由此，可得裝置泵效率ηpu計算公式為

(11)

2 數(shù)值方法

2.1 控制方程與湍流模型

本文利用商業(yè)CFD軟件Fluent進行數(shù)值計算，考慮黏性二維不可壓縮流動的雷諾平均納維-斯托克斯方程(RANS)和時均連續(xù)性方程，可知其運動控制方程為

(12)

(13)

xi、xj——i、j方向空間坐標，m

p——流體壓強，Pa

υ——層流運動粘性系數(shù)，Pa·s

υt——湍流粘性系數(shù)，Pa·s

本文采用Realizablek-ε湍流模型求解N-S方程的模型，相應方程參見文獻[26]。

2.2 計算域與網(wǎng)格生成

計算域尺寸以及遠場邊界條件對計算結(jié)果的影響不可忽略，為使翼型后續(xù)流場能夠充分發(fā)展，常取翼型后流域長為20c，故本文中設置x方向總長為8 m，并設置x正方向長度為6 m，翼型旋轉(zhuǎn)中心置于原點位置處。為了模擬河流推水時河道兩岸對流場的影響，設置流道寬度為0.8 m。計算域及初始時刻無俯仰角的網(wǎng)格劃分如圖2所示。

圖2 計算域及網(wǎng)格劃分Fig.2 Computational domain and grid partition

由于撲動水翼需要進行起伏以及俯仰結(jié)合的運動，還需要結(jié)合動網(wǎng)格技術實現(xiàn)不同時刻的網(wǎng)格重構(gòu)。為了提高計算效率，利用非結(jié)構(gòu)化三角形網(wǎng)格劃分以翼型起始位置為中心的長1 m、寬0.8 m的長方形變形域，利用結(jié)構(gòu)化網(wǎng)格劃分其余計算域。為了更好地捕獲翼型以及壁面的邊界流場，在其周圍劃分邊界層，第1層網(wǎng)格厚度為0.000 1 m(y+<2.71)。此外，為保持翼型周圍網(wǎng)格的質(zhì)量，將翼型邊界層及其向外數(shù)層網(wǎng)格單獨分為一個固定域，該區(qū)域隨翼型運動，不參與網(wǎng)格重構(gòu)。

為驗證網(wǎng)格數(shù)量的無關性，使用不同網(wǎng)格數(shù)量進行數(shù)值模擬計算。結(jié)果表明，在網(wǎng)格數(shù)量達到8.3萬后，其對計算結(jié)果的影響可以忽略，因此為減少計算的時間，選擇網(wǎng)格數(shù)為82 994。

2.3 數(shù)值計算

本文基于ANSYS Fluent軟件，結(jié)合Realizablek-ε湍流模型以及動網(wǎng)格技術進行撲動水翼推水性能的數(shù)值研究。設置近翼型側(cè)為壓力入口，遠翼型側(cè)為壓力出口，兩側(cè)壁面與翼型表面為無滑移壁面，3個計算域的交界面設為Interface，近壁計算采用增強壁面處理。

動網(wǎng)格模型采用彈性光順法以及局部重構(gòu)法。對于彈性光順法，為不影響較遠區(qū)域的網(wǎng)格，設定彈性系數(shù)為0.8，為不影響邊界處網(wǎng)格節(jié)點分布，設定邊界點松弛因子為0.000 6。對于局部重構(gòu)法，設置最大網(wǎng)格扭曲率為0.7以及合適的最大最小網(wǎng)格長度標準。最后利用用戶自定義指令(UDF)進行運動域運動形式的設定，得到不同時刻的動網(wǎng)格如圖3所示。求解器設置中，選用Coupled算法耦合壓力場與速度場，采用二階迎風格式來離散時間。而步長設置應小于最小網(wǎng)格尺度與流場速度的比值，需要根據(jù)不同工況的計算情況進行調(diào)整。

圖3 1個周期內(nèi)4個主要時刻的動網(wǎng)格更新Fig.3 Dynamic mesh updating at four main points in a cycle

2.4 方法驗證

為驗證數(shù)值計算方法的有效性，選取NACA0012翼型，在最大攻角相同的情況下，對St不同的多種工況進行計算，并與文獻[21]的試驗結(jié)果進行對比。基本參數(shù)設置與文獻[21]相同，U=0.4 m/s，φ=π/2，c=Amax=0.1 m，θmax=π/12，Re=4×104，其中U表示計算域入口的流速。

將撲動水翼的平均推力系數(shù)隨St變化的數(shù)值計算結(jié)果與文獻[21]試驗數(shù)據(jù)進行對比，結(jié)果如圖4所示。

圖4 數(shù)值模擬結(jié)果與文獻[21]試驗數(shù)據(jù)的對比Fig.4 Comparison of numerical simulation results with experimental data in reference[21]

由圖4可知，平均推力系數(shù)的計算值與試驗值相符合，具有較好的一致性，本文使用的數(shù)值計算方法是正確有效的。

3 數(shù)值結(jié)果與分析

計算模型的各參數(shù)設置如下：弦長c為0.3 m，撲動水翼的旋轉(zhuǎn)中心與前緣距離l為0.2c，起伏幅值Amax為0.5c，俯仰幅值θmax為π/6，相位角φ為-π/2，計算時展長s為0.3 m。

3.1 撲動頻率對撲動水翼裝置水動力性能的影響

為了研究撲動頻率對撲動水翼裝置水動力性能的影響，本文設置撲動水翼擺動頻率從0.1 Hz逐漸增加到1 Hz，每0.1 Hz變化一次；此后每0.5 Hz增加一次，最高頻率為5 Hz。從其中選擇較為典型的4組進行對比分析。

圖5為不同撲動頻率下?lián)鋭铀淼乃矔r推力系數(shù)隨時間變化的曲線。為了便于對比，將各自周期的無量綱時間的相對值作為橫坐標。由圖5可以看出，當撲動水翼從平衡位置開始向上撲動，瞬時推力系數(shù)先增大再減小，撲動至最大撲動幅值時瞬時推力系數(shù)到達谷點，隨后撲動水翼開始反向撲動，瞬時推力系數(shù)再次開始增加，過平衡位置后迅速達到另一波峰。而且在一個運動周期內(nèi)，瞬時推力系數(shù)總是大于零，即在撲動的全過程撲動水翼總是對水流起到推動作用。

圖5 不同頻率下?lián)鋭铀淼乃矔r推力系數(shù)變化曲線Fig.5 Variation of instantaneous thrust coefficient of flapping hydrofoil at different frequencies

圖6為不同撲動頻率下?lián)鋭铀淼乃矔r升力系數(shù)隨時間變化的曲線。由圖6可以看出，瞬時升力系數(shù)曲線基本沿軸線對稱分布，所有工況下的平均升力系數(shù)均接近零。對比圖5以及圖6可以發(fā)現(xiàn)，在瞬時推力系數(shù)的極大值處(撲動平衡位置)瞬時升力系數(shù)取極值，而在瞬時推力系數(shù)的極小值處(即最大撲動位置處)瞬時升力系數(shù)為零。

圖6 不同頻率下?lián)鋭铀淼乃矔r升力系數(shù)變化曲線Fig.6 Variation of instantaneous lift coefficient of flapping hydrofoil at different frequencies

為分析撲動水翼尾流的結(jié)構(gòu)形態(tài)對其性能的影響，選擇渦量在z方向分量Ωz為特征量，并由此繪制各工況下的渦量等值線圖，z向渦量分量的計算公式為

(14)

式中ux、uy——x、y方向的速度，m/s

圖7為不同頻率下各時刻的渦量等值線圖。可以看出，撲動水翼在有限通道內(nèi)進行撲動推水時，尾流中的反卡門渦街可能會發(fā)生偏移，被壓縮側(cè)受到的來自壁面的作用更為強烈，消散更為迅速，導致雙列反卡門渦街兩側(cè)的旋渦強度差異較大。綜合觀察圖5～7可以發(fā)現(xiàn)，瞬時推力系數(shù)與瞬時升力系數(shù)的大小與反卡門渦街的旋渦強度密切相關。以f=1 Hz為例，該工況下的尾部射流偏向斜上方，導致反卡門渦街下方的旋渦強度明顯強于上方旋渦。當翼型下?lián)鋾r，下方強旋渦對翼型的作用遠大于翼型上撲時上方弱旋渦對翼型的作用，導致翼型下?lián)鋾r對水流的推動力明顯強于上撲時的推動力，故圖5中f=1 Hz時瞬時推力系數(shù)曲線中每周期的第2個波峰遠高于第1個波峰。此外，尾渦作用也是導致瞬時推力系數(shù)極大值偏離平衡位置的原因之一。

圖7 頻率為0.1、0.5、1 Hz不同時刻的渦量等值線圖Fig.7 Contour maps of vorticity at different frequencies of 0.1 Hz, 0.5 Hz and 1 Hz

圖8為撲動水翼裝置的平均流量以及平均揚程隨著撲動頻率變化的曲線。從圖8可以看出，撲動水翼裝置的流量與其撲動頻率基本成正比，而揚程與撲動頻率的平方成正比。結(jié)合圖7以及圖8可以發(fā)現(xiàn)，在模擬過程中尾流的偏移并不會對流量以及揚程產(chǎn)生明顯的影響，即尾渦的偏移雖然增大了推力的瞬時變化，但是對平均推力影響較小。

圖8 流量以及揚程隨撲動頻率的變化曲線Fig.8 Variation of flow rate and lift with flutter frequency

圖9為撲動水翼裝置的推進效率以及泵水效率隨撲動頻率變化的曲線。可以看出，撲動水翼的推進效率與泵水效率隨著撲動頻率的增加不斷增加，但是增加的速率隨頻率增加不斷降低。對比相同頻率下的推水效率與泵水效率可以發(fā)現(xiàn)，推水效率總是高于泵水效率，而且兩者之間的差值隨著頻率的增加不斷減小，這也許是由計算段的水力損失所占的比重隨著撲動頻率的增加有所降低導致的。

圖9 推進效率以及泵水效率隨撲動頻率的變化曲線Fig.9 Variation of propulsion efficiency and pump efficiency with flutter frequency

3.2 撲動水翼裝置的推水特性

為了探究不同流量下?lián)鋭铀硌b置的揚程以及泵水效率的變化規(guī)律，選取出入口流速為0～1 m/s，每0.1 m/s變化一次進行模擬計算研究。計算過程中固定撲動頻率為1 Hz。圖10(圖中Q0.7表示來流速度為0.7 m/s時的流量)給出的是以最大效率點流量為流量標度時，撲動水翼裝置的揚程以及效率隨著流量變化的曲線。

由圖10可知，撲動水翼裝置的揚程隨著流量的增加不斷減小，但在流量較小時有著明顯的波動；而撲動水翼裝置的泵水效率隨著流量的增加先增大再減小。隨著流量的增大，撲動水翼與水流之間的相對速度不斷減小，翼型對水流的推動力越來越小，故裝置的揚程隨著流量的增加不斷減小。除此之外，在最佳效率點時翼型在撲動過程中首次對水流起到了阻礙作用。此外，該裝置具有超低的揚程，當流量為0即封閉計算域出入口時獲得的平均揚程也僅有0.024 m，在河流推水這種低揚程、大流量的應用場合具有優(yōu)勢。

圖10 頻率為1 Hz時，揚程與效率隨流量的變化曲線Fig.10 Variation of lift and efficiency with flow rate when frequency was 1 Hz

圖11為頻率1 Hz時，不同流速下3種流場的渦量等值線圖。可以發(fā)現(xiàn)：在低流速下，尾渦上下并列并且快速消散，這是因為新尾渦的生成過程會對已經(jīng)形成的舊渦造成影響，甚至是直接破壞來不及流向下游的舊渦，尾渦的能量直接損失，能量使用效率低下。而隨著流速的不斷增加，每個周期中生成的兩個旋渦逐漸分離，排列形成雙列反卡門渦街，在渦街中心形成射流加速水流的流動，尾渦能量逐漸得到利用，裝置效率不斷增加。當流速超過最佳效率點后，尾渦中每個周期中生成的兩個反方向旋渦遠離的過程中，尾渦由雙列反卡門渦街轉(zhuǎn)變?yōu)閱瘟蟹纯ㄩT渦街，相互作用減小，無法形成有效射流推動水流的流動，裝置效率開始下降。

圖11 頻率為1 Hz、不同流速時不同時刻的渦量等值線圖Fig.11 Vorticity contours at different time with frequency of 1 Hz and different flow velocities

圖12為撲動水翼在最佳效率點平衡位置的渦量等值線圖，此時撲動水翼裝置的尾跡呈現(xiàn)為如圖所示的雙列反卡門渦街形態(tài)。圖13為撲動水翼在最佳效率點平衡位置的壓力云圖，綜合觀察圖12和圖13可以發(fā)現(xiàn)，尾渦在向后移動的過程中逐漸耗散，尾渦能量中的一部分被轉(zhuǎn)化為壓力勢能，壓力逐漸增加。此外觀察其他工況下的渦量等值線圖可知，揚程在0.010～0.020 m范圍內(nèi)時，裝置尾跡均為相似的雙列反卡門渦街，且具有較好的推水性能，能適應揚程波動較大的應用場合。

圖12 撲動水翼在最佳效率點平衡位置的渦量等值線圖Fig.12 Vorticity contour map of flapping hydrofoil at equilibrium position of optimal efficiency point

圖13 撲動水翼在最佳效率點平衡位置的壓力云圖Fig.13 Pressure contour of flapping hydrofoil at equilibrium position of optimal efficiency point

4 性能試驗

為了驗證撲動水翼裝置的推水性能，設計了撲動水翼試驗裝置，如圖14所示。伺服電機與減速箱提供裝置運行所需的動力；曲柄滑塊機構(gòu)將旋轉(zhuǎn)運動轉(zhuǎn)化為直線運動，帶動葉片進行升沉運動的同時，配合兩組曲柄滑塊機構(gòu)之間的相位差實現(xiàn)葉片的俯仰運動，從而實現(xiàn)葉片的正弦撲動。

圖14 撲動水翼裝置示意圖Fig.14 Schematic of flapping hydrofoil device1.伺服電機與減速箱 2.曲柄滑塊機構(gòu) 3.導軌 4.背板 5.撲動葉片 6、7.支架

圖15為撲動水翼裝置實物圖。試驗裝置架設在4 m×0.5 m×0.7 m的開放流道上方。通過聲學多普勒流速儀測量流道出口處各點的流速，并利用高速攝像機進行流場的捕捉。其中聲學多普勒流速儀的量程為0～3 m/s，精度為0.005 m/s，采樣頻率為25 Hz；高速攝像機選用Phantom VEO 340L，幀速為100 f/s。

圖15 撲動水翼裝置試驗臺Fig.15 Flapping hydrofoil device test bench1.撲動水翼裝置 2.水槽 3.方形流道 4.葉片 5.伺服驅(qū)動器

試驗葉片展長為0.3 m，試驗頻率為0.1～1 Hz，每0.1 Hz變化一次。利用聲學多普勒流速儀測量方形流道出口處均勻設置的6個測點上流速，單次測量不小于12個周期，多次測量取其平均值為測點流速。以6個測點的平均流速作為平面流速。將試驗流速與仿真結(jié)果進行對比，如圖16所示。頻率為1 Hz時，試驗流速為0.264 m/s，總不確定度u為0.020 1 m/s。

圖16 試驗與仿真的流速對比Fig.16 Comparison of flow velocity between test and simulation

由圖16可知，試驗與模擬的結(jié)果整體趨勢相同，裝置的流量隨著運動頻率的增加而增大。但模擬結(jié)果明顯高于試驗結(jié)果，造成該差異的主要原因是試驗為保證葉片運動順利進行，在葉片兩側(cè)留有較大的空隙，導致試驗過程中部分水體通過兩側(cè)縫隙從壓力側(cè)轉(zhuǎn)移至吸力側(cè)，大大降低了裝置對水體的推動作用。此外，葉片固定塊以及泄漏也會對結(jié)果造成一定影響。

利用高速攝像機拍攝試驗流場結(jié)果如圖17所示。由于油墨的擴散較快，每個時刻僅有一個尾渦存留，可以發(fā)現(xiàn)試驗流場與仿真流場基本相同，從側(cè)面印證了仿真結(jié)果的正確性。

圖17 試驗與仿真流場圖Fig.17 Flow field diagrams of test and simulation

5 結(jié)論

(1)撲動水翼裝置的尾渦結(jié)構(gòu)形態(tài)是性能影響的關鍵因素，當尾渦處于雙列反卡門渦街時，推水性能較高。

(2)在不同頻率下，撲動水翼裝置的流量與撲動頻率成正比，揚程與撲動頻率的平方成正比。而裝置的推進效率以及泵水效率均隨著頻率的增加而增大，但增加的速率不斷降低。

(3)在固定的頻率(1 Hz)下，撲動水翼裝置的揚程隨著頻率的增加而減小，而效率隨著流量的增加先增大后減小。

(4)撲動水翼裝置具有超低的揚程，當頻率為1 Hz時，該裝置的最大平均揚程僅有0.024 m；而且揚程在0.010～0.020 m范圍內(nèi)，該裝置尾流均呈現(xiàn)雙列反卡門渦街，具有較好的推水性能，在河流推水這種低揚程、大流量的應用場合具有優(yōu)勢。