三維可變形翼展向變形氣動力特性數值模擬研究

陳 默，張 楠，鄭文濤，李永成

（中國船舶科學研究中心，江蘇 無錫 214082）

0 引 言

隨著現代航空技術的飛速發展，傳統飛行器的設計已經不能滿足人類對飛行性能和飛行品質不斷提升的需求。利用仿鳥類及昆蟲飛行時翅膀展向變形現象而發明的可變翼飛行器一直是各國航空發展領域的研究熱點[1]。可變翼飛行器甚至可利用機翼的扭轉和彎曲輔助替代操縱面，實現飛行器橫滾、偏航以及俯仰等不同姿態變化。

自上世紀80 年代開始，國外就開始了對變形機翼的研究，基于洛克韋爾公司提出的可充分利用飛行器在飛行時產生的扭轉變形概念，美國NASA 蘭利研究中心聯合美國空軍及洛克韋爾公司開展了主動柔性機翼（active flexible wing，簡稱AFW）項目，旨在研究主動柔性機翼技術[2]，且在工程設計中得到初步應用。到1996 年，在“主動柔性機翼”計劃研究基礎上，美國提出了利用機翼扭轉使得結構載荷降低、飛行包線提高的主動氣動彈性機翼（active aeroelastic wing，簡稱AAW）計劃[3]。

2006 年，歐洲主動氣動彈性飛行器（3AS，active aeroelastic aircraft structures）計劃提出了以“主動內部機構”概念來實現機翼結構的主動氣動彈性控制，該結構同樣可以實現對機翼展向變形的控制[4]；2017 年，葡萄牙貝拉國內大學的Santos 等設計了可變展長機翼結構，采用伸縮方式，由齒輪齒條機構帶動外部機翼從做成空殼狀的內部機翼中伸出，機翼的最大展長可達2.5 m，且裝配在無人機上進行了飛行測試[5]。

國內，燕山大學的Guo 等[6]于2016 年提出剛性滑動蒙皮的設計關鍵在于機翼表面蒙皮劃分的概念，并相應研究了兩種有關機翼表面蒙皮劃分的方法；哈爾濱工業大學的王禮佳提出了一種可變翼展的伸縮翼，機翼采取多級伸縮的方式，在內部布置剪叉機構來驅動機翼進行運動，除伸縮外，機構還能夠通過形狀記憶合金絲實現后掠角的變形[7]；南京航空航天大學的宋哲[8]利用形狀記憶合金作為驅動器，設計了一種能夠改變厚度的機翼，其分析結果顯示，厚度增加后的機翼阻力系數降低，升力系數和升阻比均得到了提高。

由于計算流體力學中關于建立復雜變形幾何結構的流場數值模擬技術逐漸趨于成熟，基于數值模擬方法研究變形翼的氣動力特性這一可行途徑已得到廣泛采用。Wang等[9]利用數值模擬方法研究了無人機機翼伸縮對其氣動力特性的影響，結果表明機翼伸縮對升/阻力系數影響不大，但對升阻比影響顯著；程春曉等[10]采用可壓縮Navier-Stokes 方程和Spalart-Allmaras(S-A)湍流模型，對柔性后緣可連續變彎度二維NACA 0012 翼型升/阻力系數、表面壓力分布和流場結構展開了深入的研究，基于該變形方式下的翼型，其升力系數和升阻比均得到較大改善；馬洋等[11]采用RANS 方法，結合SSTk-ω兩方程湍流模型對伸縮/后掠變構型機翼在亞聲速和超聲速階段的氣動力特性開展了研究，相對于固定機翼外形，該構型升阻比優勢明顯。

此外，機動性一直是水中航行體性能研究和設計優化的核心命題之一。舵作為主要操縱裝置，其水動力性能（舵效）的優劣直接影響船舶和水下航行體的活動范圍和機動性能，隨著船舶和水下航行體尺度越來越大，高舵效舵的需求尤其迫切。對于三維可變形機翼的初步研究，可為新型翼舵及新型操縱控制技術的研發和水下航行器的優化設計提供技術支撐。

綜上，國內外學者已針對三維可變形翼的氣動力特性開展了研究。但目前絕大多數變形翼方案均基于改變機翼展長、后掠角和弦長來實現氣動性能的改變，對于機翼展向變形的研究內容較為稀缺。因此，以展向變形作為出發點來開展三維可變形翼氣動力特性的研究具有重要的研究意義。有鑒于此，本文對三維可變形翼定常狀態下給定靜態展向變形和非定常狀態下給定動態展向變形氣動力特性及變化規律進行了數值模擬研究。

1 研究對象

1.1 靜態展向變形模型

靜態變形計算對象為三維NACA 0012 機翼（展弦比b/c=0.5 m/0.2 m=2.5，其中b和c分別表示機翼的展長和弦長），及以其為基礎的一系列給定靜態展向變形的可變形翼（即展向彎曲的剛性機翼），如圖1和圖2所示。靜態展向變形規律按軸線變形表達式（1）來確定，變形起始位置bs為0 mm，變形幅度ha分別為50 mm、100 mm和150 mm。攻角設為α=-18°～18°，間隔2°（抬頭為正），局部角度加密。

1.2 動態展向變形模型

動態變形計算對象基礎原型與靜態展向變形一致，變形規律按軸線變形表達式（2）來確定，變形幅度ha分別為100 mm和150 mm，變形頻率f分別為0.5 Hz、1 Hz、2 Hz、3 Hz和4 Hz。攻角α設為5°。

一個運動周期內動態展向變形翼的運動過程如圖3所示，其中T表示運動周期。

2 數值計算方法

2.1 控制方程

不可壓縮流體連續性方程與RANS方程的張量形式為

2.2 湍流模型

本文采用k-ε湍流模型。湍流動能k和耗散率ε方程為

式中，C1ε= 1.44，C2ε= 1.92，Cμ= 0.09，σk= 1.0，σε= 1.3，Gk、Gb為湍流動能生成項，σk、σε表示湍流Prandtl數，Sk、Sε為源項[12]。

2.3 網格劃分與邊界條件

采用多塊結構化網格劃分方式，計算區域網格的劃分如圖4所示，機翼表面周圍網格如圖5所示。對模型表面附近網格加密，且第一層網格間距根據y+來確定（y+平均值約為1）。動態變形過程利用徑向基函數（radio base function，RBF）進行變形區域網格節點插值置換來實現，即根據邊界位移定義控制點位移插值場，并通過該插值場將網格節點平移到新位置，使網格可以進行非剛性變形。RBF變形基本原理如下[13]：

每個控制節點i的已知位移可展開為

邊界條件定義如下：

速度入口：機翼前緣向前3c，設定來流速度的大小與方向；

壓力出口：機翼尾緣向后6c，設定相對于參考壓力點的流體靜壓值；

壁面：機翼外表面，設定無滑移條件；

外場：距機翼表面3c～5c，速度為未擾動的主流區速度。

2.4 數值求解

控制方程采用有限體積法離散，動量方程中的對流項采用二階迎風差分格式，擴散項采用中心差分格式，湍流模型方程采用二階迎風差分格式，時間項采用二階隱式格式。壓力速度耦合方法選用流場計算中經典的SIMPLE算法。

3 數值計算方法驗證

3.1 網格收斂性驗證

由于網格數量在CFD數值計算中扮演極其重要的角色，本文首先開展網格收斂性驗證工作，以確保數值計算的可靠性。以無變形三維NACA 0012 機翼為計算對象，選取五套網格劃分方案，依照第22 屆ITTC 推薦規程7.5-03-01-02，網格三向加細比均為 2，網格數量分別為7.53 萬、20.55 萬，55.58萬、169.83萬和442.12萬，對應方案代號設為G1、G2、G3、G4和G5。不同數量網格第一層高度相同，y+均為1。五套網格方案不同攻角下，三維NACA 0012機翼升/阻力系數變化曲線如圖6所示。其中，來流速度設為30 m/s，流體介質為空氣，特征雷諾數Re=4.11×105。升力系數CL和阻力系數CD的定義如式（9）所示，其中，D和L分別表示模型受到的阻力和升力，ρ表示空氣密度，本文取為1.225 kg/m3，V為來流速度，A為翼舵平面面積，本文A=0.1 m2。

從圖中可以看出，隨著網格數的增加，計算結果趨于穩定。當網格數增大到第四套網格時，不同攻角下的升/阻力系數近乎不變化。綜合計算精度和計算資源考慮，網格方案G4較為適合本文的計算需求，在下文三維機翼靜態和動態展向變形計算中，網格劃分方案均選用G4。

3.2 計算精度驗證

現存文獻中關于NACA 機翼流體動力的試驗數據一般針對二維（無限展長）情況，本文選用的是三維NACA 0012 機翼，缺乏公開的試驗數據對計算結果進行驗證。因此本文首先考慮雷諾數相同時展弦比變化對機翼數值計算結果的影響，選用的展弦比分別為1（b/c=0.2 m/0.2 m）、2.5（b/c=0.5 m/0.2 m）及無限展長（二維，c=0.2 m），特征雷諾數Re=4.11×105，網格劃分方案一致。

圖7 給出了不同展弦比下NACA 0012 機翼升力系數隨攻角變化的對比曲線。由圖中可以看出，二維翼型計算結果與試驗數據吻合較好，且失速點預報結果一致，從而驗證了本文網格劃分方案和物理建模方法的適用性。此外，展弦比越小，失速角越往后移，對應的最大升力系數也越小，這一趨勢與經典理論吻合。

同時考慮雷諾數變化對三維機翼升力系數數值計算結果的影響，選用的特征雷諾數Re分別為4.11×105、1×106和2×106，計算結果如圖8 所示。從圖中可以看出，在本文后續計算將采用的特征雷諾數Re=4.11×105基礎上，繼續增大雷諾數對升力系數計算結果幾乎沒有影響，且各雷諾數對應的失速角基本一致。

對于三維機翼氣動力計算精度問題，本文選用某舵的舵軸力矩實驗結果進行計算精度驗證，該模型試驗是在中國船舶科學研究中心低速風洞中完成的。其中，舵展弦比b/c=1.595，來流風速為45 m/s，舵角δb=-25°～25°，每間隔5°變一個舵角，特征雷諾數Re=3.42×105。圖9 給出了舵法向力系數隨舵角變化的對比曲線。從圖中可以看出，數值計算結果與試驗數據吻合較好，且失速區間一致，從而驗證了本文網格劃分方案和物理建模方法對三維機翼計算結果的可信度。

3.3 時間步長驗證

考慮到時間步長在非定常計算中起到極其重要的作用，在進行三維機翼動態展向變形研究前開展時間步長選取驗證工作。變形幅度ha為100 mm，變形頻率f=4 Hz，即運動周期T=0.25 s，來流速度設為30 m/s,流體介質為空氣。分別選取時間步長Δt= 0.01 s、0.005 s、0.0025 s 和0.001 25 s 進行計算驗證。計算得到的變形翼升/阻力系數隨時間變化曲線如圖10 所示。從圖中可以看出，隨著時間步長的減小，計算結果趨于穩定。當時間步長Δt≤0.005 s時，除短暫的初始擾動時間段外，升/阻力系數隨時間變化關系幾乎一致。基于計算精度和計算資源的綜合考慮，在下文動態展向變形計算中，時間步長Δt均選用0.005 s。

4 計算結果及分析

4.1 靜態展向變形幅度的影響分析

（1）對于阻力系數，在失速發生之前，負攻角時，不同展向變形幅度下的阻力系數相對于原始無變形三維機翼均略有增大，且相對增幅均在2.8%以下；正攻角時，不同展向變形幅度下的阻力系數均小于原始無變形三維機翼，且下降幅度隨著變形幅度的增大而增大，最大相對降幅達6.0%。因此從阻力產生的能耗考慮，正攻角下展向變形幅度增加對其氣動力性能有利。

（2）對于升力系數，其增減趨勢幾乎與阻力系數一致，且增減幅度更加明顯，負攻角下相對增幅最高可達6.3%，正攻角下最大相對降幅達9.8%。從圖13中橫剖面流線和速度梯度范圍分布對比結果來看，正負攻角下升/阻力系數呈現略微不對稱的現象是由變形翼的三維效果所引起。從提高升力的角度考慮，負攻角下展向變形幅度增加對其氣動力性能有利。此外，不同展向變形幅度下的機翼失速點一致，均在±15o附近；結合圖14 和圖15 可知，展向變形的存在對于機翼迎流面壓力沒有明顯影響，主要是改變背流面負壓值，負攻角下展向變形增加時背流面壓力值會更低。

（3）對于升阻比，不同展向變形幅度下最大升阻比均出現在±6°附近，與原始無變形三維機翼一致，大小略有差異。從升阻比角度考慮，正負攻角下展向變形幅度為50 mm時其氣動力性能均最為優越。

（4）總體而言，靜態展向變形對三維機翼氣動力性能影響不大。

4.2 動態展向變形幅度的影響分析

對三維NACA 0012可變形翼模型以不同變形頻率和不同變形幅度做動態展向變形運動的氣動力特性進行數值模擬，表1 給出了不同變形翼組合的阻力系數和升力系數平均值的比較。圖16 給出了變形幅度一定時不同變形頻率下升/阻力系數時歷變化過程。由以上結果可得給定動態展向變形基本規律：

表1 動態展向可變形翼阻力系數和升力系數平均值比較Tab.1 Comparison of average drag coefficient and lift coefficient of dynamically spanwise deformable airfoil

（1）展向動態變形的存在對阻力系數和升力系數平均值略有影響，阻力系數平均值均隨變形頻率的增大而微幅下降，而對于升力系數，展向變形幅度ha=100 mm 時其平均值隨展向變形頻率的增大而微幅上升，展向變形幅度ha=150 mm 時則無明顯規律可循。因此從兼顧阻力能耗和提高升力角度來考慮，以小幅度高頻做動態展向變形時氣動力性能會略微得到改善。

（2）不同變形頻率下的升力系數和阻力系數時歷變化規律幾乎一致，呈類余弦曲線形式，且阻力系數谷值點均正對應于升力系數峰值點，均出現在變形翼經平衡位置向下做展向變形運動的瞬間。相對應地，變形翼經平衡位置向上做展向變形運動時升力系數達到谷值點，阻力系數峰值點則出現在變形翼向上到達最大變形位置時。此外，升力系數和阻力系數變化幅度均隨變形頻率的增大而增大，且變形幅度越大，變化幅度越明顯。因此在變形翼設計階段，可考慮以不同變形頻率和不同變形幅度來規劃一個展向變形周期內的運動，以達到氣動力性能最優化的結果，例如可使變形翼在由下到上半個周期內低頻變形，而在由上到下半個周期內高頻變形，上下半程變形幅度亦可做相應調整。

圖17 和圖18 分別為變形幅度ha=150 mm 時，截取的縱剖面和橫剖面一個展向變形周期內的速度云圖。從圖中可以看出，變形翼在經平衡位置向下變形到向上變形回至平衡位置半個運動周期內，機翼迎流面附近速度逐漸增大，背流面附近速度逐漸減小，對應的升力系數隨之逐漸減小；同時機翼前緣駐點附近低速區范圍逐漸減小，對應的阻力系數隨之逐漸增大；往后半周期內則正好相反。

5 結 論

本文基于三維NACA 0012 系列展向可變形翼，開展了定常狀態下給定靜態展向變形和非定常狀態下給定動態展向變形氣動力特性及變化規律的數值模擬研究。主要結論如下：

（1）靜態展向變形對三維機翼氣動力性能略有影響，影響程度與攻角有關。正攻角下升力系數略有減小，阻力系數略有增大，負攻角則正好相反。不同靜態展向變形幅度下對應的失速角完全一致。

（2）展向動態變形的存在對三維機翼阻力系數和升力系數平均值略有影響。阻力系數平均值均隨變形頻率的增大而微幅下降，升力系數變化規律與展向變形幅度有關。從兼顧阻力能耗和提高升力角度來考慮，以小幅度高頻做動態展向變形時氣動力性能會略微得到改善。

（3）不同變形頻率下的升力系數和阻力系數時歷變化規律幾乎一致，呈類余弦曲線形式。阻力系數谷（峰）值點對應升力系數峰（谷）值點，出現在變形翼經平衡位置向下（上）做展向變形運動的瞬間；升力系數和阻力系數變化幅度均隨變形頻率的增大而增大，且變形幅度越大，變化幅度越明顯。在變形翼設計階段，可考慮以不同變形頻率和不同變形幅度來規劃一個展向變形周期內的運動，以達到氣動力性能最優化的結果。

上述結論可為展向變形翼的氣動力性能預報以及優化設計提供技術支撐，分析思路適用于變形翼其它變形形式下的展向或弦向變形研究。