振動翼推進性能計算及試驗研究

，，

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

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

圖1 振動翼運動示意

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

1 理論基礎

1.1 理論推導及公式

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

圖2 計算邊界及坐標系示意

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

(1)

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

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

在物體表面SB上，根據流體不可穿透物面的性質，可得固體壁面上的運動學邊界條件：

(2)

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

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

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

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

(3)

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

(4)

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

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

1.2 粘性摩擦力修正

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

(5)

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

ρ——流體密度；

Si——單元面積；

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

ei——Vti的單位矢量。

1.3 振動翼升力系數、推力系數

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

(6)

2 結果討論和分析

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

2.1 數值算法的驗證

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

y(t)=Asin(2πft)

(7)

式中：A——升沉幅度；

f——振動頻率。

翼型剖面為NACA0012。

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

圖3 推力系數曲線的比較

2.2 三維振動翼試驗研究及數值計算

2.2.1 試驗介紹

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

圖4 振動翼試驗示意

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

圖5 振動翼運動機構示意

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

θ=ω·t

(8)

(9)

(10)

(11)

(12)

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

2.2.2 數值計算結果

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

圖6 推力系數隨計算時間變化曲線

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

圖7 推力系數

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

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

3 結論

用非定常邊界元法對振動翼時變水動力性能進行計算，在大多數算例中，大約經過6個周期后，計算數據趨于周期性的穩定狀態。本文數值結果與相關算例、試驗結果比較吻合，說明所開發的計算程序合理可行。

[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.