，，

(1.華中科技大學(xué) 計算機科學(xué)與技術(shù)學(xué)院，武漢 430074；2.中國艦船研究設(shè)計中心，武漢 430064；3.中國船級社武漢規(guī)范研究所，武漢 430022)

振動翼的推進效率高，能提供較大的瞬時升力[1]，在推進、低速操縱、流體動力機械以及航態(tài)穩(wěn)定控制等領(lǐng)域得到了廣泛關(guān)注[2-3]。振動翼是做垂蕩和縱搖復(fù)合運動的機翼，典型運動見圖1。

圖1 振動翼運動示意

在該坐標(biāo)系下，振動翼以線速度Vx沿x軸做整體平移運動，并且沿y軸以線速度Vy作上下升沉運動，同時振動翼還繞自身的轉(zhuǎn)軸OA以變角速度ω作縱搖運動。振動翼在運動過程中翼面環(huán)量Γ(t)將隨時間變化，時變尾渦面模型及尾渦存儲效應(yīng)對其水動力性能的準(zhǔn)確預(yù)報有很大影響；如何采用合適的計算模型描述這種復(fù)雜現(xiàn)象成為振動翼水動力計算的關(guān)鍵。因此深入研究振動翼非定常水動力性能對于了解水動力產(chǎn)生機理和優(yōu)化振動翼結(jié)構(gòu)與運動參數(shù)[4]意義重大。

1 理論基礎(chǔ)

1.1 理論推導(dǎo)及公式

根據(jù)勢流理論，定義邊界S包圍的流場區(qū)域V，S由機翼表面SB及尾渦面SW組成，見圖2，由于對稱性，只顯示一半邊界。

圖2 計算邊界及坐標(biāo)系示意

當(dāng)場點P(x,y,z)于邊界面S上時，可以得出場點位于邊界面S上的格林定理表述[5]：

(1)

式中：r——場點P(x,y,z)到源點Q(x0,y0,z0)之間的距離；

Δφ——SW兩側(cè)的速度勢差,Δφ=φ+-φ-。

在物體表面SB上，根據(jù)流體不可穿透物面的性質(zhì)，可得固體壁面上的運動學(xué)邊界條件：

(2)

式中：V——物面SB上某點的運動速度。

式(2)右端為該點運動速度的法向投影。

在振動翼尾緣處采用壓力Kutta條件，即：

(p+)T.E=(p-)T.E

(3)

將積分方程(1)離散化可以求解物面上的速度勢，進而可得到振動翼的水動力。在振動翼計算中，振動翼表面SB的網(wǎng)格在弦向可采用余弦分布形式，共劃分為NB個單元；尾渦面SW從尾緣到下游劃分為k個展向條帶，每個條帶沿展向方向的劃分依據(jù)振動翼表面網(wǎng)格展向劃分情況決定，尾渦面共劃分為NW個單元。計算時給每個單元編號為Nj(j=1，2，…，NB+NW)。將場點Pi(xi,yi,zi)置于每個單元的形心處，并把式(1)中的φ(P)不單獨表示，則對每個場點有：

(4)

式中：rij——場點到每個單元的距離。

求解方程(4)得到物面上單元控制點的速度勢φj，再計算得到物面上的速度分布，于是可以得到物面壓力場和受力。在每一時步計算中尾渦面上的Δφj分布可以這樣確定：與尾緣相接的第1個展向條帶中的Δφj由Kutta條件決定，并采用Newton-Raphson迭代方法實現(xiàn)式(3)的滿足；從尾緣第2個展向條帶到下游第k個展向條帶中的Δφj是采用時間步進法，依據(jù)該時步之前的時間歷程t=m·Δt(其中Δt為每次計算的時間間隔，m=1，2，…)內(nèi)計算得到的Δφj并沿尾渦面下泄而確定。對于做周期運動的振動翼，按照上述方法經(jīng)過若干個周期的計算后，可以得到振動翼穩(wěn)定的周期變化的水動力性能。在通常的計算中， 6～8個周期后振動翼水動力性能的變化趨于周期性的穩(wěn)定狀態(tài)。

1.2 粘性摩擦力修正

為了使計算結(jié)果更準(zhǔn)確，需要對在勢流中采用非定常邊界元法得到的結(jié)果進行相應(yīng)的粘性摩擦力修正。將振動翼表面每一單元上的摩擦阻力投影到振動翼推力和升力方向，就可以得到更合乎實際的振動翼水動力性能。摩擦阻力為每一單元上摩擦力的和,近似為

(5)

式中：Cf——摩擦阻力系數(shù)，根據(jù)平板摩擦阻力系數(shù)公式,Cf=0.075/(lgRe-2)2；

ρ——流體密度；

Si——單元面積；

Vti——單元上流體切向速度；

ei——Vti的單位矢量。

1.3 振動翼升力系數(shù)、推力系數(shù)

將振動翼上的合力Ft分別向升力方向eL和推力方向eT投影后，經(jīng)計算就可得到振動翼的升力系數(shù)CL和推力系數(shù)CT：

(6)

2 結(jié)果討論和分析

為了充分檢驗本數(shù)值方法和計算程序的正確與否，首先用非定常邊界元法計算兩種三維振動水翼在純升沉運動下的時均推力系數(shù)曲線，與文獻結(jié)果進行比較；然后計算三維振動翼在不同速度下的升沉縱搖耦合運動的水動力性能，并與相關(guān)試驗結(jié)果進行對比。

2.1 數(shù)值算法的驗證

算例中所選文獻結(jié)果均為勢流結(jié)果，故此處的計算不進行粘性修正以便比較驗證。水翼以零攻角做純升沉運動，運動方程一般為

y(t)=Asin(2πft)

(7)

式中：A——升沉幅度；

f——振動頻率。

翼型剖面為NACA0012。

算例。選用展弦比為lh=5的矩形翼。

圖3 推力系數(shù)曲線的比較

2.2 三維振動翼試驗研究及數(shù)值計算

2.2.1 試驗介紹

振動翼試驗在拖曳水池中進行，見圖4(XOZ平面內(nèi)視圖),其翼面垂直于水面放置，上端面距離水面距離H為1倍弦長。由于水翼振動運動的方向是沿Y方向，所以在此水深下可忽略自由面影響。

圖4 振動翼試驗示意

為了模擬振動翼的運動特性，試驗機構(gòu)見圖5(XOY平面內(nèi)視圖)。電動機通過定比減速器驅(qū)動曲柄OA繞O點以角速度ω勻速旋轉(zhuǎn)，曲柄帶動鉸接在A點的連桿AB運動，其中B點鉸接在沿固定導(dǎo)軌往復(fù)運動的滑塊上，B點的運動軌跡位于圖6中的Y軸上。而水翼CD則固接在連桿AB上，且弦線與連桿垂直，C點為前緣，D點為尾緣。這樣，曲柄OA的轉(zhuǎn)動便帶動水翼做平面運動，即垂蕩和縱搖的復(fù)合運動。當(dāng)整套運動機構(gòu)安裝在以速度V前進的水池拖車上時，便可以進行振動翼的水動力試驗。改變拖車速度的大小，可以得到不同狀態(tài)下的振動翼水動力性能曲線。力傳感器置于拖車與試驗裝置機架之間，采用動態(tài)應(yīng)變儀對信號放大，計算機的數(shù)據(jù)采樣頻率為40 Hz。

圖5 振動翼運動機構(gòu)示意

曲柄OA與該圖中X軸正向之間的夾角為θ，逆時針方向旋轉(zhuǎn)為正。曲柄OA長度為r，連桿AB長度為l。水翼CD在任一時刻的攻角為α，B點的Y軸坐標(biāo)為yB。由圖中幾何關(guān)系可以得到

θ=ω·t

(8)

(9)

(10)

(11)

(12)

由上述各式可以看到，適當(dāng)變化曲柄OA的長度，就可以改變振動翼的運動規(guī)律。在本算例中，確定的各構(gòu)件尺寸為r=0.15 m，l=0.42 m。水翼為矩形，翼型為NACA0021，展長為0.515 m，弦長為0.168 m，曲柄OA轉(zhuǎn)速為62 r/min。

2.2.2 數(shù)值計算結(jié)果

計算中每次迭代的時間間隔ΔT=0.005T。其中：T為曲柄旋轉(zhuǎn)一周的時間。尾渦面長度為振動翼弦長的15倍。用非定常邊界元法分別對拖車速度分別為V=0.5、1.0 m/s時的振動翼的推進性能進行計算，所得的推力系數(shù)瞬時值見圖6。

圖6 推力系數(shù)隨計算時間變化曲線

在一個周期中，振動翼的推力系數(shù)曲線有2個峰。由圖中可以看到，經(jīng)過大約6個周期后，計算數(shù)據(jù)趨于周期性的穩(wěn)定狀態(tài)。一個周期中，車速分別為V=0.5、1.0 m/s的推力系數(shù)試驗值見圖7。

圖7 推力系數(shù)

從數(shù)值計算結(jié)果與試驗數(shù)據(jù)的對比看，計算值在規(guī)律上與試驗值一致，但在數(shù)值上還有一定的誤差。這主要是因為振動翼在粘性流體中做非定常運動，產(chǎn)生的漩渦流動比較復(fù)雜，而目前的理論模型又難以精確定量描述這類情況，所以值得深入研究這些問題和機理。

從圖中可以看到，在一個周期中，前半周期和后半周期的推力系數(shù)不一致，這主要是因為該試驗裝置運動規(guī)律的不完全對稱所造成。從公式(12)中可以看到，dyB/dt在一個周期中變化不完全一樣，所以振動翼的水動力特性也會有少許差異。

3 結(jié)論

用非定常邊界元法對振動翼時變水動力性能進行計算，在大多數(shù)算例中，大約經(jīng)過6個周期后，計算數(shù)據(jù)趨于周期性的穩(wěn)定狀態(tài)。本文數(shù)值結(jié)果與相關(guān)算例、試驗結(jié)果比較吻合，說明所開發(fā)的計算程序合理可行。

[1] READ D A. Oscillating foils for propulsion and maneuvering of ships and underwater vehicles[D]. M. S. Thesis, MIT, 2000.

[2] TRIANTAFYLLOU G S,TRIANTAFYLLOU M S,GROSENBAUGH M A.Optimal thrust development in oscillating foils with application to fish propulsion[J]. Journal of Fluids and Structures, 1993(7)：205-224.

[3] ANTOINE D, JACQUES A A, FRANCOIS D, et al. Computational and experimental investigation of flow over a transient pitching hydrofoil[J]. European Journal of Mechanics B/Fluids, 2009, 28:728-743.

[4] ELLENRIEDER K D, PARKE K R, SORIA J. Fluid mechanics of flapping wings[J]. Experimental Thermal and Fluid Science, 2008,32：1578-1589.

[5] NEWMAN J N. Marine Hydrodynamics[M]. Cambridge: The MIT Press，1977.

[6] ZHU Q, LIU Y, YUE D K. Dynamics of a three-dimensional oscillating foil near the free surface[J]. AIAA J.2006,44(12)：2297-3009.