孫青林，梁 煒，陳增強，賀應平

(1.南開大學 計算機與控制工程學院，天津 300071；2.航宇救生裝備有限公司實驗部，湖北 襄陽，441022)

襟翼偏轉翼傘氣動性能數(shù)值模擬分析

孫青林1，梁 煒1，陳增強1，賀應平2

(1.南開大學 計算機與控制工程學院，天津 300071；2.航宇救生裝備有限公司實驗部，湖北 襄陽，441022)

為研究襟翼偏轉對翼傘氣動性能的影響，對不同襟翼偏轉情況分別建立CFD模型，通過有限體積法進行空間離散并求解RANS方程，模擬翼傘在轉向與雀降階段的氣動性能，進而結合最小二乘法進行參數(shù)辨識，實現(xiàn)翼傘氣動模型的修正.模擬結果表明:襟翼偏轉會引起翼傘壓強分布改變，失速迎角減小，升阻力系數(shù)突增，對翼傘氣動性能造成復雜的影響;修正的翼傘氣動模型可以較好描述翼傘氣動性能與襟翼偏轉的變化規(guī)律，相比傳統(tǒng)氣動模型有效地提高了計算精度，為翼傘在轉向與雀降階段的精確建模提供參考.

翼傘；襟翼偏轉；氣動特性；數(shù)值模擬；氣動模型

翼傘是一種前緣有切口，利用沖壓空氣保持一定形狀的柔性飛行器.翼傘由柔性不透氣的材料包裹，通過肋片氣孔平衡內(nèi)部壓力，增強飛行性能的穩(wěn)定性，已廣泛應用于飛行器回收、物資運輸、港口消霧等領域[1-2].目前翼傘動力學模型中較少提及襟翼偏轉的計算方法，為實現(xiàn)翼傘精確建模與自主歸航，研究翼傘的氣動性能有重要的意義.

本文基于上述研究成果，重點研究襟翼偏轉對翼傘氣動性能的影響，首先引入基于CFD的襟翼偏轉氣動計算方法，將襟翼偏轉下的翼傘氣動模型的修正問題轉為參數(shù)辨識問題；然后針對轉向與雀降對應的不同襟翼偏轉方式進行CFD數(shù)值模擬，研究翼傘氣動性能的變化，并通過NASA風洞試驗驗證數(shù)值模擬的有效性；最后通過最小二乘法結合氣動模擬數(shù)據(jù)進行參數(shù)辨識，完成翼傘氣動模型的修正工作.

1 襟翼偏轉翼傘氣動模型

1.1 翼傘氣動外形

翼傘的氣動外形如圖1所示，其中氣動弦長c為2.133 6 m，展長b為6.400 8 m，弧面下反角β為9.55°，傘繩R為9.601 2 m，前緣切口高度h為0.242 1 m.并對翼傘做如下簡化：1)翼傘近似剛體，僅在襟翼偏轉時考慮形變；2)傘衣結構不透氣；3)忽略傘繩及負載.在保證數(shù)值模擬精度的前提下，上述簡化可以降低計算成本，縮短數(shù)值模擬周期.

圖1 翼傘氣動外形參數(shù)

1.2 翼傘氣動模型

傳統(tǒng)翼傘模型進行氣動計算時，將傘衣沿展向?qū)ΨQ等長分為8片，每片升力系數(shù)從外到里乘以0.60、1.00、1.16、1.24作為修正因子[12]，對應氣動方程為:

式中：i=1,2,…,8；Faero、Maero分別為翼傘氣動力與氣動力矩；Lo-i為翼傘質(zhì)心到各分片i質(zhì)心的矢量；TRi-o為分片i坐標系到翼傘坐標系的轉換矩陣，矩陣中的γi為分片i轉換角；FL,i、FD,i分別為分片i的升力與阻力；CL、CD分別為升、阻力系數(shù)；ρ為空氣密度；S為翼傘特征面積；v為傘體坐標系下氣流速度.

翼傘氣動模型難點在于氣動系數(shù)CL,CD的確定，傳統(tǒng)計算方法為：

上述計算方式存在較大缺陷，未考慮襟翼偏轉帶來的附加氣動力.襟翼偏轉δ僅通過增大迎角α近似模擬，忽略了傘衣形變影響，造成傳統(tǒng)氣動計算存在較大誤差.為提高翼傘的氣動計算精度，考慮襟翼偏轉的附加氣動力，本文對傳統(tǒng)氣動模型進行修正，如:

針對偏轉系數(shù)與迎角之間的非線性關系，使用二階模型即可滿足：

至此翼傘模型中的襟翼偏轉氣動計算問題，轉為對12個偏轉系數(shù)因子的辨識問題，本文通過CFD數(shù)值模擬結合最小二乘法對其進行辨識.

2 CFD求解方法與分析

2.1 CFD模型與網(wǎng)格劃分

共使用5個不同形態(tài)的翼傘CFD模型，如圖2所示.模型1中襟翼無偏轉，模型2中襟翼雙側1/3偏轉，模型3中襟翼雙側2/3偏轉，模型4中襟翼雙側全偏轉，模型5中襟翼單側2/3偏轉.其中模型1、2、3、4用于雀降階段，模型1、3、5用于轉向階段.

翼傘以弦長75%處為軸彎折[3]，即后緣25%的翼尖向下彎折，彎折角度稱為下折角.1/3偏轉對應下折角為25°，2/3偏轉對應下折角為50°，全偏轉對應下折角為75°，展向兩側偏轉量大，展向中部偏轉量小.

圖2 翼傘CFD模型

采用三維非結構化化網(wǎng)格進行空間劃分，如圖3所示，流場邊界條件應盡量遠離擾動源(翼傘)，但實際計算不允許邊界條件取得過遠，根據(jù)計算經(jīng)驗確定邊界距離翼傘分別為10c,5c,10c，c為翼傘的氣動弦長.翼傘近壁面網(wǎng)格需要加密，網(wǎng)格法向拉伸率不宜過大[13]，控制在1.12以內(nèi)，網(wǎng)格總數(shù)大致為150萬.

2.2 流場設置與求解方法

流場環(huán)境設置如下：流體為不可壓空氣，流速為12.192 m/s，入口為速度入口，出口為壓力出口，非進出口為自由邊界，翼傘為壁面邊界.翼傘氣流迎角設置范圍為0°～18°，每隔2°模擬1次，1種模型需要模擬10組數(shù)據(jù)，全部數(shù)據(jù)共70組.每次模擬時間大致為120 min，在失速迎角附近需要大致180 min.

使用的控制方程為三維笛卡爾坐標系下的時均Navier-Stokes(RANS)方程[14]，其張量指標形式表示為：

式中：i,j取值范圍為(1,2,3)，表示三維笛卡爾坐標系的3個方向;u為流體速度；u′為流體脈沖速度；ρ為流體密度；p為流體動能；μ為擴散系數(shù)；S為源項.

圖3 流場空間網(wǎng)格劃分

空間離散方式為二階迎風格式，求解器選擇基于分離式求解算法的SIMPLE求解器.湍流模型選擇RNGk-ε二方程模型[15]，對近壁面逆壓梯度變化有較好捕捉效果.

2.3 CFD有效性驗證

使用文獻[3]中風洞試驗數(shù)據(jù)對CFD的有效性進行驗證.試驗使用的ND模型為傘繩剛性連接，NASA模型為傘繩柔性連接.

翼傘滑翔階段，不存在襟翼偏轉，傘衣未發(fā)生形變.如圖4(a)所示，升力系數(shù)與風洞試驗數(shù)據(jù)接近，失速迎角偏大，是風洞試驗翼傘上翼面氣流提前分離造成的，這些偏差在翼傘數(shù)值模擬中普遍存在，不影響翼傘建模.阻力系數(shù)與風洞試驗數(shù)據(jù)重合較好.翼傘雀降階段，存在襟翼偏轉，傘衣發(fā)生形變.如圖4(b)所示，襟翼2/3偏轉與完全偏轉情況下，升阻力系數(shù)與風洞試驗數(shù)據(jù)相對接近，偏差存在是風洞試驗傘衣厚度、粗糙度、彈性形變等因素造成的.在數(shù)值模擬過程中這些偏差的存在是合理的，從而驗證了CFD數(shù)值模擬的有效性.

圖4 風洞試驗與數(shù)值模擬對比

3 襟翼偏轉影響分析

3.1 雀降階段分析

翼傘雀降階段，襟翼雙側偏轉，增大阻力實現(xiàn)減速.由于翼傘為柔性材料，襟翼偏轉帶來翼傘傘衣形變，對應氣動性能發(fā)生改變.

CFD模型選擇模型2、3、4.如圖5所示，壓力對稱分布，隨著襟翼偏轉量的增加，低壓區(qū)范圍緩慢增加.展向中部低壓區(qū)明顯，在下折角附近沒有新的低壓區(qū)產(chǎn)生，而展向邊緣低壓區(qū)不明顯，在下折角附近有新的低壓區(qū)產(chǎn)生.相對而言，高壓區(qū)范圍迅速增加，展向中部高壓范圍大于展向邊緣，集中在前緣開口與下折角部分，是阻力的主要來源.

如圖6所示，傘衣表面最低壓均出現(xiàn)在上翼面圓弧過渡區(qū)域，最高壓均出現(xiàn)在前緣切口區(qū)域.襟翼雙側無偏轉時，翼傘上翼面與下翼面的弦向壓力分布，與氣動外形即為相似，壓力差與傘衣厚度相關，壓力最小值為-101.2 Pa，最大值為59.7 Pa.襟翼雙側2/3偏轉時，翼傘上翼面與下翼面壓力分布，與氣動外形存在較大差異，壓力差相比之下有明顯增加，壓力最小值為-145.3 Pa，最大值為48.6 Pa.

對不同襟翼偏轉量分別數(shù)值模擬，如圖7所示，隨著偏轉量的增加，翼傘失速迎角逐漸減小，分別為14°、12°、10°、8°.相同迎角的升力系數(shù)迅速增大，增長速度由快變慢甚至由正變負；相同迎角的阻力系數(shù)呈單調(diào)遞增，增長速度相對平穩(wěn).因此，襟翼偏轉對翼傘氣動性能存在非線性影響，進行偏轉參數(shù)因子的辨識是必要的.

圖5 雙側襟翼偏轉下翼傘壓力分布

圖6 翼傘弦向壓力分布對比

通過最小二乘法進行參數(shù)辨識，可得到表1中修正因子的數(shù)值.如圖8所示，通過修正氣動系數(shù)的計算值與CFD數(shù)據(jù)重合很好，最大偏差與CFD數(shù)值相差5.47%，位置在無偏轉時失速迎角附近.不同的襟翼偏轉量對應不同的升阻力系數(shù)曲線，比傳統(tǒng)不考慮襟翼偏轉的翼傘氣動方程更客觀科學，為翼傘雀降階段建模提供參考.

表1 雙側襟翼偏轉翼傘模型修正

圖7 雙側襟翼偏轉翼傘氣動系數(shù)

圖8 雙側襟翼偏轉氣翼傘動系數(shù)曲面

3.2 轉向階段分析

翼傘轉向階段，襟翼單側偏轉，使翼傘向偏轉一側轉向，偏轉量越大，轉彎半徑越小.由于壓力不對稱，翼傘展向環(huán)量不再滿足橢圓分布，傳統(tǒng)翼傘氣動方程同樣不適用.

CFD模型選擇模型1、3、5.如圖9所示，單側偏轉時上翼面低壓區(qū)范圍介于無偏轉低壓區(qū)與雙側偏轉低壓區(qū)之間，均集中在圓弧過渡區(qū)域.相對而言，下翼面高壓區(qū)分布形狀存在不規(guī)則性，并且向偏轉側偏移，是翼傘下翼面兩側的迎風面積不同造成的，同時也是翼傘偏航力矩的主要來源.

圖9 單、雙側襟翼偏轉翼傘壓力分布對比

如圖10所示，襟翼單側無偏轉時，傘衣表面壓力沿展向?qū)ΨQ分布，壓力最小值為-101.2 Pa，集中在上翼面展向中部的圓弧區(qū)域，壓力最大值為59.7 Pa，集中在展向中部前緣切口區(qū)域.襟翼單側2/3偏轉時，傘衣表面壓力沿展向非對稱分布，壓力最小值為-145.3 Pa，集中在上翼面展向偏轉側的圓弧區(qū)域，壓力最大值為58.8 Pa，集中在展向兩側前緣切口區(qū)域.對于展向邊緣的低壓突峰，是由上翼面新低壓區(qū)造成的.

圖10 翼傘展向壓力分布對比

對不同襟翼偏轉量分別數(shù)值模擬，如圖11所示，隨著偏轉量的增加，翼傘失速迎角逐漸減小，分別為14°、13°、12°、10°.相同迎角的升力系數(shù)迅速增大，增長速度由快變慢甚至由正變負；相同迎角的阻力系數(shù)呈單調(diào)遞增，增長速度相對平穩(wěn).與雀降情況相似，襟翼偏轉對翼傘氣動性能同樣存在非線性影響.

圖11 單側不同偏轉量氣動系數(shù)變化

通過最小二乘法進行參數(shù)辨識，可得到表2中修正因子的數(shù)值.如圖12所示，修正氣動系數(shù)計算值與CFD數(shù)據(jù)重合很好，最大偏差與CFD數(shù)值相差4.93%，位置在無偏轉時失速迎角附近.不同的襟翼偏轉量對應不同的升阻力系數(shù)曲線，比傳統(tǒng)不考慮襟翼偏轉的翼傘模型更客觀科學，為翼傘轉彎階段建模提供參考.

表2 單側襟翼偏轉翼傘模型修正

3.3 傳統(tǒng)模型與修正模型對比

傳統(tǒng)模型無論翼傘處于滑翔、轉向、雀降階段，均通過增大迎角模擬襟翼偏轉，如圖13所示，由于為考慮襟翼偏轉帶來的傘衣形變，其準確性必然下降.表3數(shù)據(jù)表明，雙側襟翼無偏轉，傳統(tǒng)模型與修正模型氣動系數(shù)差別可忽略不計；雙側襟翼完全偏轉，傳統(tǒng)模型通過增大迎角模擬，修正模型通過設置偏轉量模擬，傳統(tǒng)模型升力系數(shù)偏大17.3%，阻力系數(shù)偏小18.5%.傳統(tǒng)模型計算誤差隨著襟翼偏轉量的增加而增加，影響翼傘模型精度.通過CFD與最小二乘法得到的修正模型，對翼傘精確建模存在重要意義.

圖12 單側襟翼偏轉氣動系數(shù)曲面

Fig.12 Aerodynamic coefficient surface of parafoil with unilateral flap deflection

圖13 翼傘襟翼偏轉對比

表3 傳統(tǒng)模型與修正模型氣動系數(shù)對比

4 結 論

1) 襟翼偏轉引起翼傘氣動性能急劇變化，隨著偏轉量的增加，翼傘壓強分布改變，失速迎角減小，氣動系數(shù)激增，相比無偏轉階段存在很大差異.

2) 結合CFD數(shù)值模擬與最小二乘法，完成翼傘氣動模型的修正工作.對比模擬數(shù)據(jù)，修正后模型較好描述了翼傘氣動性能與襟翼偏轉的變化規(guī)律.

3)針對傳統(tǒng)氣動模型在翼傘轉向與雀降階段存在的升力系數(shù)偏高、阻力系數(shù)偏低的問題，修正后氣動模型的計算精度得到較大提高，為翼傘的精確建模提供參考.

[1] TAO Jin, SUN Qinglin, ZHU Erlin, et al. Quantum genetic algorithm based homing trajectory planning of parafoil system[C]//Proceedings of the 34th Chinese Control Conference,CCC 2015. Hangzhou, China: IEEE, 2015: 2523-2528.DOI：DOI: 10.1109/ChiCC.2015.7260028.

[2] ZHU Erlin, SUN Qinglin, TAN Panlong, et al. Modeling of powered parafoil based on Kirchhoff motion equation[J]. Nonlinear Dynamics, 2015, 79(1): 617-629. DOI: 10.1007/ s11071-014-1690-9.

[3] NICOLAIDES J D. Parafoil wind tunnel tests[R]. Indiana: University of Notre Dame, 1971.

[4] BALAJI R, MITTAL S, RAI A K. Effect of leading edge cut on the aerodynamics of ram-air parachutes[J]. International Journal for Numerical Methods in Fluids, 2005, 47(1): 1-17. DOI: 10.1002/fld.779.

[5] MITTAL S, SAXENA P, SINGH A. Computation of two-dimensional flows past ram-air parachutes[J]. International Journal for Numerical Methods in Fluids, 2001, 35(6): 643-667. DOI: 10.1002/1097-0363(20010330)35:6.

[6] HAN Yihua, YANG Congxin, WANG Yuwei, et al. Aerodynamics simulation of a large multi-cells parafoil[C]// Proceedings of the 20th AIAA Aerodynamic Decelerator Systems Technology Conference. Seattle, WA: AIAA, 2009: 2978-2991.

[7] SINGH R, BAEDER J D. Direct calculation of three-dimensional indicial lift response using computational fluid dynamics[J]. Journal of Aircraft, 1997, 34(4): 465-471. DOI: 10.2514/2.2214.

[8] KALRO V, TEZDUYAR T. A parallel 3D computational method for fluid-structure interactions in parachute systems[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 190(3/4): 321-332. DOI: 10.1016/S0045-7825(00)00204-8.

[9] CAO Yihua, ZHU Xu. Effects of characteristic geometric parameters on parafoil lift and drag[J]. Aircraft Engineering & Aerospace Technology, 2013, 85(4): 280-292. DOI: 10.1108/AEAT-Jun-2011-0096.

[10]POTVIN J, HURST B, PEEK G. Semi-numerical derivation of the opening shock factor and inflation time for slider-reefed parafoils[C]//Proceedings of the 19th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar. Williamburg, VA: AIAA, 2007: 2508-2517.DOI: 10.2514/6.2007-2508.

[11]LINGARD J. The aerodynamics of gliding parachutes[C]// Proceedings of the 9th Aerodynamic Decelerator and Balloon Technology Conference. Albuquerque, NM: AIAA, 1986: 1-13.DOI: 10.2514/6.1986-2427.

[12]熊菁. 翼傘系統(tǒng)動力學與歸航方案研究[D]. 長沙: 國防科學技術大學, 2005. XIONG Jing. Research on the dynamics and homing project of parafoil system[D]. Changsha: National University of Defense Technology, 2005.

[13]ONATE E, VALLS A, GARCA J. Computation of turbulent flows using a finite calculus-finite element formulation[J]. International Journal for Numerical Methods in Fluids, 2007, 54(6/7/8): 609-637. DOI: 10.1002/fld.1476.

[14]王福軍. 計算流體動力學分析: CFD 軟件原理與應用[M]. 北京: 清華大學出版社, 2004: 115-116. WANG Fujun. Computational fluid dynamics analysis: principles and applications of CFD software[M]. Beijing: Tsinghua University Press, 2004: 115-116.

[15]MANSOUR N, KIM J, MOIN P. Near wall k-epsilon turbulence modeling[J]. AIAA Journal, 1989, 27(8): 1068-1073.DOI: 10.2514/3.10222.

(編輯 張 紅)

Numerical simulation analysis for aerodynamic performance of parafoil with flap deflection

SUN Qinglin1, LIANG Wei1, CHEN Zengqiang1, HE Yingping2

(1.College of Computer and Control Engineering, Nankai University, Tianjin 300071, China; 2.Dept. of Experiment, Aerospace Life-Support Equipment Co., Ltd., Xiangyang 441022, Hubei, China)

To investigate the effect of flap deflection on aerodynamic performance of the parafoil, CFD models are established for different flap deflections. By using the finite volume method, the spatial discretization is carried out and the RANS equation is solved. The aerodynamic performance of the parafoil in the stage of turning and flaring is simulated, and then combining with the least square procedure, parameter identification is carried out to implement the modification of the parafoil aerodynamic model. Simulation results show that the flap deflection can cause the change of the parafoil pressure distribution, the decrease of the stall angle of attack, the sudden increase of the lift drag coefficient, and the complex effect on aerodynamic performance of the parafoil. Compared with the traditional aerodynamic model, the modified parafoil aerodynamic model can better describe the regulation of the flap deflection and the aerodynamic performance of the parafoil, and the calculation precision is improved effectively. It provides a theoretical reference for precise modeling of the parafoil in the stage of turning and flaring.

parafoil; flap deflection; aerodynamic performance; numerical simulation; aerodynamic model

10.11918/j.issn.0367-6234.201508089

2015-08-27

國家自然科學基金(61273138)；天津市重點基金(14JC2DJC39300)

孫青林(1963—)，男，教授，博士生導師； 陳增強(1964—)，男，教授，博士生導師

孫青林，sunql@nankai.edu.cn

V211.3

A

0367-6234(2017)04-0048-07