王國才，李春華，陳平劍，趙小全

（1.海軍駐景德鎮(zhèn)地區(qū)軍事代表室；2.中國直升機設(shè)計研究所，江西 景德鎮(zhèn) 333001；3.總參謀部陸航部駐景德鎮(zhèn)地區(qū)軍事代表室，江西 景德鎮(zhèn) 333002）

直升機依靠旋翼來產(chǎn)生升力和操縱力，因此旋翼對操縱輸入的動態(tài)響應(yīng)特性在很大程度上決定了直升機的飛行性能[1]。許多機動飛行科目，如垂直躍上/躍下、快速側(cè)移等，其操縱輸入是階躍式的快速輸入，會導(dǎo)致旋翼拉力、誘導(dǎo)速度、揮舞角等氣動參數(shù)和旋翼性能出現(xiàn)動態(tài)的非定常變化，尤其是關(guān)鍵的總距輸入改變會對旋翼的操縱性能有較大的影響。因此，開展旋翼槳葉在總距突增時的氣動響應(yīng)特性分析具有重要的理論和實際意義[2]。

槳葉的氣動響應(yīng)與旋翼入流分布密切相關(guān)，目前旋翼入流的求解方法主要有動量理論、廣義動態(tài)入流理論、渦流理論方法和CFD方法等。動量理論只給出了最簡單的均勻誘導(dǎo)速度分布，不適合于槳葉氣動載荷及動態(tài)響應(yīng)特性計算。廣義動態(tài)入流理論是20世紀80年代由Peters、He等人[3]發(fā)展起來的旋翼入流計算方法，該方法有較好的精度，但難以準確模擬旋翼尾跡畸變對入流分布的瞬時影響。而CFD方法分析操縱響應(yīng)在目前還不成熟。渦流理論方法[4-5]，尤其是基于自由尾跡的渦方法計入了尾跡的收縮和渦線的畸變，允許渦線自由移動，考慮了尾跡自誘導(dǎo)作用，可以更好地模擬實際尾跡形狀及入流變化，目前廣泛用于模擬槳葉尾跡和分析旋翼入流特性，具有快捷高效、計算精度高的優(yōu)點。而且，近年來Leishman[6]等人開始提出了一種時間精確的自由尾跡分析模型，該模型利用隨時間變化的尾跡迭代算法，精確預(yù)測尾跡的瞬時變化過程，能夠較好地計算各時間點旋翼入流的非定常影響。國內(nèi)對時間精確自由尾跡計算方法的研究很少，而關(guān)于自由尾跡方法對旋翼槳葉動態(tài)響應(yīng)特性的研究目前還未見相關(guān)的文獻發(fā)表。

本文在課題組已有的旋翼準定常松弛尾跡研究的基礎(chǔ)[7]上，通過建立隨時間變化的旋翼自由尾跡模型、槳葉氣動模型及動力學模型和非定常翼型模型，并結(jié)合旋翼配平模型，嘗試建立一種新的時間精確自由尾跡計算方法。并以該方法分析槳葉在總距突增時，旋翼的力和力矩系數(shù)非定常變化的時間歷程，以及槳葉的動態(tài)揮舞響應(yīng)特性。

1 計算模型

1.1 尾跡模型

隨著旋翼旋轉(zhuǎn)，槳葉后緣脫出螺旋形尾跡，且無論在懸停還是前飛，各片槳葉的尾跡渦均交織在一起[5]。因此，要利用B-S定理來預(yù)測旋翼附近的誘導(dǎo)速度，需先計算尾跡渦線的幾何形狀。

旋翼渦線以當?shù)厮俣龋ê杂闪魉俣群托碚T導(dǎo)速度 vind）自由運動，尾跡渦線的控制方程[7]可描述為：

對于尾跡的離散采用PC2B方法[6]，即用一個二階向后差分方法將網(wǎng)格[ψ,ψ+Δψ]×[?,?+Δ?]中心點(ψ+Δψ/2,?+Δ?/2)的差分表達式近似代替控制方程中的偏微分項，計算點與周圍節(jié)點的依賴關(guān)系如圖1所示。

圖1 尾跡離散點的計算示意圖

利用上述有限差分格式，尾跡控制方程離散為：

1.2 槳葉模型

本文槳葉模型采用Weissinger-L 簡單升力面模型[6-7]，該模型對槳葉以定環(huán)量的渦段來模擬，即沿槳葉展向分成若干小段，每一小段中附著渦位于1/4弦線，環(huán)量Γ為常數(shù)，相同數(shù)目的控制點布置在3/4弦線中點處。由于載荷沿槳葉展向變化，不同強度的渦從槳葉后緣拖出，構(gòu)成旋翼的尾跡。

旋翼尾跡分為近尾跡和遠尾跡，近尾跡由沿每片槳葉后緣脫出的自由尾隨渦和脫體渦構(gòu)成槳葉渦面，其各段渦線的強度為槳葉相鄰兩展向微段的附著環(huán)量之差，而每段渦線環(huán)量的弦向變化形成脫體渦。沿槳葉后緣脫出的尾跡在轉(zhuǎn)過一定角度后聚合并卷起為槳尖渦線，構(gòu)成遠尾跡。

Weissinger-L 氣動模型的難點在于計算各渦線段的環(huán)量值。本文采用邊界條件即控制點的法向誘導(dǎo)速度為零來進行求解，如式(4)：

為了求解方程(4)，需先求出影響系數(shù)，對于任意控制點I，影響系數(shù)I(u,v,w)ij是指所有渦線，環(huán)量為1時對控制點的誘導(dǎo)速度，在獲得影響系數(shù)后，邊界條件(4)可表示為：

由于渦線位置難以確定，因此氣動模型中對環(huán)量的求解需結(jié)合1.1的自由尾跡模型。

槳葉的揮舞動力方程通過離心力矩、慣性力矩、氣動力矩、彈簧力矩的平衡來建立，方程如下：

1.3 非定常翼型模型

要準確計算槳葉的拉力，需獲得詳細的翼型升、阻特性參數(shù)。而對于大多數(shù)翼型而言，該參數(shù)難以獲得。為此，本文采用了Beddoes 翼型模型[8]計算翼型的升、阻特性。其迎角表示的翼型升力系數(shù)計算公式為：

式中：M為馬赫數(shù)。

對式中的f，Beddoes 提出的一個經(jīng)驗表達式[8]：

式中：0α為翼型的零升迎角；1α為氣流分離點為0.7時失速角；而1S和S2定義了翼型的靜態(tài)失速特性。1α、1S和S2的數(shù)值可表示成馬赫數(shù)M的指數(shù)關(guān)系式，具體取值參見文獻[8]。

翼型的阻力系數(shù)表示為：

式中：

Cd0為零升阻力系數(shù)，αDD為阻力發(fā)散角，也可分別表示成馬赫數(shù)M的指數(shù)關(guān)系式[8]。

1.4 配平模型

利用式(6)獲得了槳葉揮舞角后，需要進一步求出旋翼的總距和周期變距[5]。計算思路為：固定旋翼的軸傾角和前飛速度，調(diào)整旋翼總距角 θ0.7和周期變距角 θ1c、θ1s以滿足給定的旋翼拉力系數(shù)CT和槳盤后倒角 a1、側(cè)傾角 b1。

平衡求解通過將非線性輸出矢量對操縱矢量進行泰勒級數(shù)展開至一階并化簡得：

于是輸入和輸出之間存在如下關(guān)系：

式中，J為雅可比矩陣，定義為：

在求解過程中，初始值不一定在真實解附近，為保證穩(wěn)定性，引入一阻尼因子λ（0＜λ≤1.0）

當滿足 Δ x＜RMS（RMS為指定小量）時，認為獲得收斂的配平操縱量。

2 計算步驟

利用上節(jié)給出的各分析模型，本文建立了新的時間精確自由尾跡方法，其分析旋翼氣動響應(yīng)特性的計算步驟為：

步驟一，給出飛行條件、旋翼參數(shù)，計算旋翼平均入流和尾跡初始位置；

步驟二，沿方位角步進，迭代計算出定常飛行條件的尾跡幾何形狀，并完成配平分析；

步驟三，對旋翼的總距操縱輸入做突增改變；

步驟四，方位角步進，計算尾跡節(jié)點的誘導(dǎo)速度，使用PC2B算法更新尾跡的幾何形狀；

步驟五，計算槳葉的環(huán)量和氣動力變化，分析槳葉的揮舞響應(yīng)特性，返回步驟四，計算下一個時間步的尾跡形狀和動態(tài)響應(yīng)特性，直到各氣動響應(yīng)值達到相對穩(wěn)定。

3 驗證算例

旋翼入流計算的準確性是載荷響應(yīng)和揮舞響應(yīng)計算的關(guān)鍵。因此，利用本文的方法首先計算了旋翼入流，來驗證方法的正確性。以馬里蘭大學的2MRTS 旋翼[9]為算例，圖2計算了模型旋翼在0.15和0.3 前進比的旋翼縱向入流分布，圖中橫坐標表示槳盤縱向位置（來流方向）與旋翼半徑的比值。可以看出，入流分布為典型前大后小，入流梯度與前進比有關(guān)，計算結(jié)果與實驗值[9]吻合較好。

圖2 前飛時旋翼縱向入流的計算值與實驗值對比

圖3以全尺寸直升機旋翼[10]為算例，計算了懸停時槳葉總距以200 (°)/s的速度增加到12°，之后保持總距不變時，旋翼拉力系數(shù)變化的時間歷程，圖中橫坐標表示方位角的變化。可以看出懸停時的槳葉的拉力沒有周期性變化，拉力系數(shù)迅速增大且有過沖，然后降低并在一周以后逐步達到穩(wěn)定狀態(tài)，計算結(jié)果與實驗結(jié)果[10]吻合較好。

圖3 懸停時計算的旋翼拉力系數(shù)響應(yīng)與實驗值對比

通過對入流和總距增加時拉力系數(shù)響應(yīng)曲線的實驗值與計算值對比，可以知道本文的方法有較好的準確性，可用來進行旋翼的響應(yīng)特性分析。

4 總距突增響應(yīng)分析

本節(jié)以表1的模型旋翼為算例，計算總距階躍突增時旋翼的力、力矩系數(shù)變化及揮舞響應(yīng)特性。輸入條件為：旋翼穩(wěn)定飛行后，旋轉(zhuǎn)一周，然后總距作階躍式突增。

表1 模型旋翼參數(shù)

4.1 懸停

拉力系數(shù)為0.01時，計算獲得的懸停旋翼總距值為0.233 rad，周期變距為0。此時旋翼揮舞形成錐度角為2.78°，側(cè)倒角和后倒角為0°。由于懸停時槳葉的氣動環(huán)境相同，因而旋翼的拉力系數(shù)和錐度角均無周期性變化。圖4為將槳葉總距增加到0.245 rad和0.25 rad時，其旋翼拉力系數(shù)變化的時間歷程，受階躍總距輸入影響，拉力系數(shù)也近乎階躍式突增，隨著拉力系數(shù)的增大，旋翼入流也增大，于是槳葉的有效迎角又降低了。因此，拉力系數(shù)又迅速回落，之后振蕩上升，經(jīng)過約兩個周期后，達到穩(wěn)定狀態(tài)。穩(wěn)定后的拉力系數(shù)比總距增加前分別增加了20%和26%。

圖4 總距突增時旋翼拉力系數(shù)變化的時間歷程（懸停）

圖5為總距增加到0.245 rad和0.25 rad時，旋翼錐度角的響應(yīng)特性，初始錐度角2.78°。相比拉力系數(shù)變化趨勢，總距突增后錐度角的增大較為平緩，之后呈振蕩趨勢，也在約兩個周期后在更大的錐度角位置達到穩(wěn)定，穩(wěn)定后的錐度角分別為3.12°和3.27°。

圖5 總距突增時的錐度角變化的時間歷程（懸停）

4.2 前飛

圖6分別為前進比為0.1和0.2時的總距突增時的旋翼拉力系數(shù)響應(yīng)歷程。受來流影響，前飛時槳葉的氣動環(huán)境有周期性變化，因此，旋翼拉力系數(shù)也呈現(xiàn)周期性變化。前進比為0.1和0.2，拉力系數(shù)為0.01時計算的槳葉總距分別是0.219 rad和0.204 rad。與懸停類似，總距階躍突增到0.25 rad后，拉力系數(shù)迅速增大，然后回落，并振蕩達到穩(wěn)定。槳葉總距增加到0.25 rad后，前飛前進比0.2的槳葉拉力系數(shù)比0.1 大。

圖6 總距突增時旋翼拉力系數(shù)變化的時間歷程（前飛）

圖7為前飛時旋翼錐度角的響應(yīng)歷程?？偩嘣黾忧?，前進比0.2時的錐度角比0.1的略大?？偩嗤辉龊螅F度角較快增大且過沖回落，然后達到穩(wěn)定，前飛時錐度角也呈現(xiàn)周期性變化，但幅度并不大。對比圖5可知，相同拉力系數(shù)時，懸停時錐度角最小。前進比越大，錐度角也越大。

圖7 總距突增時的錐度角變化的時間歷程（前飛）

圖8為前進比為0.1，總距突增時，俯仰和滾轉(zhuǎn)力矩變化的時間歷程。與拉力系數(shù)類似，前飛時受周期氣動環(huán)境的影響，俯仰和滾轉(zhuǎn)力矩均有周期性變化趨勢。在前進比0.1時，模型旋翼的俯仰力矩和滾轉(zhuǎn)力矩的峰值相位差約為32°?？偩嘣黾雍?，力矩系數(shù)的幅值都增加了，但俯仰力矩的均值相比總距變化前變化不明顯，而滾轉(zhuǎn)力矩略有增加。

圖8 總距突增時旋翼力矩系數(shù)變化的時間歷程（μ=0.1）

圖9為前進比為0.2，俯仰和滾轉(zhuǎn)力矩變化的時間歷程。

圖9 總距突增時旋翼力矩系數(shù)變化的時間歷程（μ=0.2）

隨著前進比的增加，相同拉力系數(shù)時的俯仰和滾轉(zhuǎn)力矩均有所增大。在前進比為0.2時，模型旋翼的俯仰力矩和滾轉(zhuǎn)力矩的峰值相位差約為41°，比前進比為0.1時大。俯仰力矩仍然是均值基本不變，幅值增大。而滾轉(zhuǎn)力矩系數(shù)是先減小后增大。在前進比為0.2時，算例計算的俯仰力矩略大于滾轉(zhuǎn)力矩系數(shù)。

5 結(jié)論

1）無論是懸停還是前飛，槳葉總距的階躍突增，會使旋翼拉力系數(shù)產(chǎn)生階躍式的增大，然后迅速回落且振蕩變化，約在2個周期后達到穩(wěn)定。而錐度角的響應(yīng)相對平滑，先增加過沖后再回落，并振蕩穩(wěn)定。拉力系數(shù)和錐度角的穩(wěn)態(tài)值主要由旋翼總距和飛行速度決定。

2）前飛時，受周期氣動環(huán)境的影響，在總距改變之前和之后，旋翼的拉力系數(shù)、俯仰力矩系數(shù)、滾轉(zhuǎn)力矩系數(shù)和錐度角都有周期性變化，但錐度角的周期性變化較小。

3）錐度角與拉力系數(shù)和前進比有關(guān)。在相同拉力系數(shù)時，懸停時錐度角最小，隨著前進比增大，錐度角也越大。

4）前飛時，總距突增后俯仰力矩幅值增大，均值幾乎不變。而滾轉(zhuǎn)力矩系數(shù)不但幅值增大，且均值也略有增加。前進比的大小會改變俯仰和滾轉(zhuǎn)力矩系數(shù)峰值之間的相位差。

5）本文建立的新的時間精確自由尾跡計算方法，能夠有效計算總距突增時的旋翼拉力、力矩系數(shù)和揮舞響應(yīng)等旋翼的非定常氣動特性。

[1]徐國華,王適存.具有后掠槳尖的旋翼氣動特性計算方法[J].空氣動力學學報,1999,17(3)∶111-112.

[2]PETERS D A,NINH H Q.Dynamic inflow for practical applications[J].Journal of the American Helicopter Society,1988,33(4)∶64-68.

[3]HE C J.Development and application of a generalized dynamic wake theory for lifting rotors[D].Georgia Institute of Technology,1989.

[4]MILLER W O,BLISS D B.Direct periodic solutions of rotor free wake calculations[J].Journal of the American Helicopter Society,1993,38(2)∶53-60.

[5]BAGAI A,LEISHMAN J G.Rotor free wake modeling using a pseudo implicit technique-including comparisons with experimental data[J].Journal of the American Helicopter Society,1995,40(3)∶1276-1285.

[6]BHAGWAT M J,LEISHMAN J G.Time-accurate free-vortex wake model for dynamic rotor response[C]//American Helicopter Society Specialist Meeting.2000.

[7]李春華,徐國華.懸停和前飛狀態(tài)傾轉(zhuǎn)旋翼機的旋翼自由尾跡計算方法[J].空氣動力學學報,2005,23(2)∶152-156.

[8]BEDDOES T S.A wake model for high resolution airloads[C]//Proceedings of the 2ndInternational Conference on Basic Rotorcraft Research.1985.

[9]ELLIOT J W,ALTHOFF S L,SAILEY R H.Inflow measurements made with a laser velocimeter on a helicopter model in forward flight[R].NASA TM 100545,1988.

[10]CARPENTER P J,FRIDOVICH B.Effect of a rapid blade-pitch increase on the thrust and induced-velocity response of a full-scale helicopter rotor[R].NACA TN 3044,1953.