自然層流機(jī)翼氣動(dòng)外形優(yōu)化研究

馬曉永，張彥軍，段卓毅，郭洪濤，李 權(quán)

（1.中國(guó)空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川綿陽(yáng) 621000；2.中國(guó)空氣動(dòng)力研究與發(fā)展中心高速空氣動(dòng)力研究所，四川綿陽(yáng) 622762；3.中國(guó)航空工業(yè)集團(tuán)公司第一飛機(jī)設(shè)計(jì)研究院，陜西西安 710089）

自然層流機(jī)翼氣動(dòng)外形優(yōu)化研究

馬曉永1，2，＊，張彥軍3，段卓毅3，郭洪濤2，李 權(quán)3

（1.中國(guó)空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川綿陽(yáng) 621000；2.中國(guó)空氣動(dòng)力研究與發(fā)展中心高速空氣動(dòng)力研究所，四川綿陽(yáng) 622762；3.中國(guó)航空工業(yè)集團(tuán)公司第一飛機(jī)設(shè)計(jì)研究院，陜西西安 710089）

層流技術(shù)在飛行器尤其是民用客機(jī)減阻方面具有較好的應(yīng)用潛力。采用全速勢(shì)方程（TRANAIR）、序列二次規(guī)劃（Sequential Quadratic Programming，SQP）優(yōu)化算法及自適應(yīng)笛卡爾網(wǎng)格技術(shù)，用兩種不同的參數(shù)化建模方法，類函數(shù)／型函數(shù)變換（Class function／Shape function Transformation，CST）和非均勻有理B樣條（Non-Uniform Rational B-Spline，NURBS）方法，對(duì)某自然層流機(jī)翼進(jìn)行了多點(diǎn)優(yōu)化設(shè)計(jì)，設(shè)計(jì)點(diǎn)為：CL＝0.45（M：0.735，0.755，0.765，0.775），CL＝0.50（M：0.755），約束條件為：升力、俯仰力矩和翼型20%、80%及最大厚度不得低于初值。優(yōu)化結(jié)果表明：雖然CST方法在翼型描述上具有先天優(yōu)勢(shì)，NURBS方法更多用于三維曲面，而對(duì)于文中自然層流機(jī)翼算例，NURBS的優(yōu)化效果更好；在CL＝0.45（M：0.755）狀態(tài)下，優(yōu)化后上翼面轉(zhuǎn)捩位置明顯后移，層流區(qū)域變廣，總阻力減小了4.5%；采用的優(yōu)化方法在進(jìn)行五點(diǎn)優(yōu)化時(shí)仍然具有較高的計(jì)算效率，具有較強(qiáng)的工程實(shí)用能力。

自然層流機(jī)翼；氣動(dòng)外形優(yōu)化；TRANAIR；參數(shù)化方法；數(shù)值模擬

0 引 言

基于CFD（Computational Fluid Dynamics）的數(shù)值優(yōu)化是自然層流機(jī)翼設(shè)計(jì)較為有效的方法，它集成了數(shù)值模擬、優(yōu)化算法和參數(shù)化建模方法等內(nèi)容。在優(yōu)化算法方面，基于梯度的尋優(yōu)算法、進(jìn)化類算法和響應(yīng)面方法是在工程優(yōu)化研究領(lǐng)域內(nèi)較為常用的方法。由于飛行器優(yōu)化設(shè)計(jì)參數(shù)通常較多，對(duì)進(jìn)化類算法（如遺傳算法）要求的種群規(guī)模較大，而響應(yīng)面也較難準(zhǔn)確構(gòu)建。基于梯度的優(yōu)化算法是當(dāng)前較為有效的方法，如在機(jī)翼氣動(dòng)外形優(yōu)化迭代過(guò)程中，正是由這些梯度信息告知下一步尋優(yōu)方向和步長(zhǎng)選擇。已有的研究表明，整個(gè)優(yōu)化過(guò)程95%以上的時(shí)間都耗用在流場(chǎng)求解方面，包括梯度求解時(shí)的流場(chǎng)計(jì)算，因此流場(chǎng)的計(jì)算效率直接決定著整個(gè)優(yōu)化過(guò)程的時(shí)間消耗，也是優(yōu)化方法具備工程實(shí)用價(jià)值的重要考核指標(biāo)［12］。波音公司TRANAIR全速勢(shì)方法在流場(chǎng)求解方面效率較高，被廣泛應(yīng)用于波音系列飛機(jī)的優(yōu)化設(shè)計(jì)方面［13］。

本文以TRANAIR軟件為優(yōu)化平臺(tái)，采用兩種參數(shù)化建模及動(dòng)網(wǎng)格方法：非均勻有理B樣條（Non-Uniform Rational B-Spline，NURBS）和類函數(shù)／型函數(shù)變換（Class function／Shape function Transformation，CST）方法，對(duì)某自然層流機(jī)翼進(jìn)行了亞聲速條件下氣動(dòng)外形優(yōu)化設(shè)計(jì)，并對(duì)優(yōu)化效果進(jìn)行了分析和討論。

1 優(yōu)化模型

對(duì)于飛行器氣動(dòng)外形優(yōu)化設(shè)計(jì)問(wèn)題，通常定義目標(biāo)函數(shù)為阻力系數(shù)CD，約束條件為升力系數(shù)CL、俯仰力矩系數(shù)Cm和機(jī)翼幾何約束（厚度、容積等）等分別不小于它們各自初值，同時(shí)還要滿足自變量設(shè)計(jì)范圍、邊界層分離因子和物面變形曲率等其它約束條件的限制。多點(diǎn)優(yōu)化時(shí)，通過(guò)不同設(shè)計(jì)點(diǎn)下的目標(biāo)函數(shù)加權(quán)平均構(gòu)造總的目標(biāo)函數(shù)來(lái)實(shí)現(xiàn)，并滿足不同設(shè)計(jì)點(diǎn)下的氣動(dòng)、幾何約束條件。

優(yōu)化算例為某民用通用飛機(jī)設(shè)計(jì)使用的自然層流機(jī)翼，其幾何外形如圖1所示。機(jī)翼參考面積80 m2，平均氣動(dòng)弦長(zhǎng)3.05m，展弦比10.5，前緣后掠角17.5°，25%弦線后掠角15.6°，后緣內(nèi)側(cè)后掠角0°，外側(cè)后掠角9.6°，翼根、轉(zhuǎn)折和翼尖處相對(duì)厚度分布約13.5%、12%和10%。主設(shè)計(jì)點(diǎn)為M＝0.755，CL＝0.45，Re＝18×106，多點(diǎn)優(yōu)化時(shí)為五個(gè)設(shè)計(jì)點(diǎn)，其中四個(gè)為不同M 數(shù)同升力設(shè)計(jì)點(diǎn)，另外一個(gè)為主設(shè)計(jì)M數(shù)時(shí)CL＝0.50點(diǎn)，如圖2所示。

圖1 某自然層流機(jī)翼Fig.1 Natural laminar flow wing

五點(diǎn)優(yōu)化時(shí)，鑒于自然層流機(jī)翼的波阻較小，目標(biāo)函數(shù)設(shè)置時(shí)賦予較小的加權(quán)比例，突出總阻力和型阻的比重；約束條件中各設(shè)計(jì)點(diǎn)下升力和俯仰力矩不得減小，機(jī)翼各剖面20%、80%和最大厚度不得減小。其數(shù)學(xué)模型：

圖2 多點(diǎn)優(yōu)化圖示Fig.2 Multi-points optimization

其中，CD1～CD5分別為各設(shè)計(jì)點(diǎn)下的阻力：

式中，CDFar為總阻力，CDi為誘導(dǎo)阻力、CDp為外型阻力、CDw為激波阻力。設(shè)計(jì)變量為物面變形控制參數(shù)和來(lái)流迎角。在本文算例中，采用CST方法時(shí)，72個(gè)參數(shù)控制上下翼面，5個(gè)迎角，一共77個(gè)設(shè)計(jì)變量。NURBS方法中，僅采用30個(gè)參數(shù)控制上翼面（不優(yōu)化下翼面），5個(gè)迎角，一共35個(gè)設(shè)計(jì)變量。

（3）積極探索“五方聯(lián)動(dòng)”和“社工+”的資源整合的創(chuàng)新模式，初步形成系統(tǒng)化解決婚姻家庭困境的標(biāo)準(zhǔn)模式。重慶市民政局、重慶市婚管中心、各區(qū)縣民政局、婚姻登記處及社工機(jī)構(gòu)等多方參與，以“社工+”為服務(wù)平臺(tái)，通過(guò)開(kāi)展五大服務(wù)內(nèi)容，整合資源，創(chuàng)新重慶市婚姻家庭社會(huì)工作服務(wù)標(biāo)準(zhǔn)化模式，不斷提升社區(qū)家庭成員的能力建設(shè)、社區(qū)參與及社會(huì)融入。

2 參數(shù)化方法

2.1 NURBS方法

對(duì)于一個(gè)NURBS曲面，若u方向?yàn)閜階，控制點(diǎn)數(shù)目為n＋1，v方向?yàn)閝階，控制點(diǎn)數(shù)目為m＋1，則：

其中：

式中Pi，j為控制點(diǎn)矢量，總數(shù)目為（n＋1）×（m＋1），wi，j為權(quán)重系數(shù)，Ni，p（u）和Nj，q（v）為分別定義在非減節(jié)點(diǎn)矢量U（u）、V（v）上的基函數(shù)，即：

同理可以給出Nj，q（v）和V（v）的表達(dá)式。方程式的分段形式，保證了NURBS函數(shù)的局域性質(zhì)，即控制點(diǎn)的變化僅會(huì)影響到其附近區(qū)域的變化，不會(huì)影響到其它區(qū)域［12，14－15］。

2.2 CST方法

在CST方法中，首先將物面的物理坐標(biāo)用歸一化參數(shù)表示，然后把參數(shù)坐標(biāo)解析表達(dá)成型函數(shù)和類函數(shù)的乘積，其中類函數(shù)用有明確物理意義的典型翼型參數(shù)（前緣半徑，彎度等）表示，型函數(shù)則在類函數(shù)的基礎(chǔ)上定義特定的幾何形體［16－20］。其具體表達(dá)式為：

其中：

式中下標(biāo)U表示機(jī)翼上表面，LOC表示機(jī)翼當(dāng)?shù)仄拭?；c（η）為當(dāng)?shù)叵议L(zhǎng)；xLE為當(dāng)?shù)厍熬夵c(diǎn)。圖3給出了部分參數(shù)的示意圖。下翼面也可相應(yīng)給出。

在本文算例中，CST方法使用了四個(gè)主控制剖面（翼型），每個(gè)翼型上下各使用了9個(gè)控制參數(shù)，共有72個(gè)控制參數(shù)實(shí)現(xiàn)整個(gè)機(jī)翼上下翼面的參數(shù)化描述。由于NURBS方法更適合空間三維曲面表述，且對(duì)于自然層流機(jī)翼的精細(xì)優(yōu)化，上翼面的氣動(dòng)外形對(duì)阻力影響往往較大，故在NURBS方法中，僅使用30個(gè)參數(shù)優(yōu)化上翼面。

圖3 CST參數(shù)化方法［20］Fig.3 CST parameterization［20］

2.3 優(yōu)化平臺(tái)

波音公司TRANAIR軟件在工程領(lǐng)域經(jīng)過(guò)大量的測(cè)試，具有較好的可信度［13］。流場(chǎng)求解采用無(wú)粘、無(wú)旋的全速勢(shì)方程，層流轉(zhuǎn)捩采用基于線性穩(wěn)定性理論的en法［21］，并用人工粘性方法進(jìn)行湍流模擬。優(yōu)化計(jì)算過(guò)程中，使用具有自動(dòng)加密功能的笛卡爾網(wǎng)格，易于物面／空間網(wǎng)格變形及加速流場(chǎng)收斂。

優(yōu)化算例是在中航工業(yè)一飛院曙光服務(wù)器上進(jìn)行，處理器為Dual-Core AMD Opteron（tm）Processor 2220（主頻2.80GHz），內(nèi)存8G，操作系統(tǒng)為Red Hat Enterprise Linux。TRANAIR軟件的核心代碼及參數(shù)化接口均用計(jì)算能力較強(qiáng)的Fortran語(yǔ)言程序開(kāi)發(fā)。由于采用了網(wǎng)格逐漸加密和自適應(yīng)網(wǎng)格收斂加速技術(shù)，雖然兩種參數(shù)化方法設(shè)計(jì)變量數(shù)目相差一倍，但是時(shí)間消耗基本一樣。五點(diǎn)優(yōu)化采用單CPU串行計(jì)算，耗時(shí)約25～30h。

3 結(jié)果分析

圖4給出了兩種參數(shù)化方法的優(yōu)化結(jié)果，從中可以看出兩種方法優(yōu)化后阻力都有所減小。主設(shè)計(jì)點(diǎn)（CL＝0.45，M：0.755），CST方法優(yōu)化后阻力系數(shù)減小了0.00022，NURBS方法優(yōu)化后阻力減小了0.00053，同比分別減小了1.8%和4.5%。并且隨著馬赫數(shù)的增加，優(yōu)化后阻力減小更加明顯，在（CL＝0.50，M：0.775）設(shè)計(jì)點(diǎn)，NURBS方法優(yōu)化后阻力減小了0.00083之多。單純從阻力減小量值上看，NURBS方法的優(yōu)化結(jié)果明顯優(yōu)于CST方法。

圖4 不同參數(shù)化方法的優(yōu)化結(jié)果Fig.4 Optimization results of different parameterizations

圖5給出了上翼面優(yōu)化前后轉(zhuǎn)捩位置對(duì)比（綠色表示層流區(qū)域，紅色表示湍流區(qū)域）。圖6為優(yōu)化前后不同展向壓力和摩擦阻力對(duì)比。從中可以看出，原始的機(jī)翼已具備較好的阻力特性，上、下翼面分別有約25%和55%的層流區(qū)。在內(nèi)翼段，優(yōu)化后轉(zhuǎn)捩位置后移，層流區(qū)域變廣，圖6（a）10%展長(zhǎng)剖面位置對(duì)比顯示，CST方法優(yōu)化后轉(zhuǎn)捩位置后移了2.5%當(dāng)?shù)叵议L(zhǎng)，NURBS方法優(yōu)化后則后移達(dá)到了18.1%當(dāng)?shù)叵议L(zhǎng)；在48%展長(zhǎng)剖面位置（圖6（b）），CST方法優(yōu)化后轉(zhuǎn)捩位置后移了約3.5%當(dāng)?shù)叵议L(zhǎng)，NURBS方法優(yōu)化后則后移了約4.4%當(dāng)?shù)叵议L(zhǎng)；在下翼面處，CST方法優(yōu)化后轉(zhuǎn)捩位置后移了約1.1%，由于NURBS方法并沒(méi)有優(yōu)化下翼面，轉(zhuǎn)捩位置幾乎不變。

從10%剖面壓力分布上看，翼面上較好的順壓梯度是保證層流區(qū)的必要條件。對(duì)于上翼面，原始機(jī)翼在前緣附近壓力減小過(guò)快，而后順壓梯度變化較緩，導(dǎo)致轉(zhuǎn)捩過(guò)早發(fā)生。優(yōu)化后前緣順壓變化較為平緩，進(jìn)而延遲了層流轉(zhuǎn)捩，帶來(lái)較廣的層流區(qū)域。此外，優(yōu)化結(jié)果顯示，上翼面氣動(dòng)外形優(yōu)化是自然層流機(jī)翼精細(xì)優(yōu)化設(shè)計(jì)的主要方向。

圖5 優(yōu)化前后上翼面轉(zhuǎn)捩位置比較Fig.5 Laminar flow transition of up wing surface

圖6 不同展向壓力和摩阻對(duì)比Fig.6 Pressure and friction coefficient in different spanwise directions

4 結(jié) 論

通過(guò)對(duì)某自然層流機(jī)翼氣動(dòng)外形優(yōu)化設(shè)計(jì)研究，得出了以下主要結(jié)論：

（1）在本文三維層流機(jī)翼優(yōu)化算例中，與傳統(tǒng)基于剖面翼型的參數(shù)化方法相比較，面向三維曲面的NURBS方法的表現(xiàn)出了更好的優(yōu)化效果，主設(shè)計(jì)點(diǎn)優(yōu)化后阻力減小了4.5%。

（2）機(jī)翼翼面上較緩的順壓梯度分布是層流保持的必要條件。

（3）基于全速勢(shì)（TRANAIR）的自然層流機(jī)翼優(yōu)化設(shè)計(jì)方法，多點(diǎn)優(yōu)化時(shí)具有較高的計(jì)算效率，工程實(shí)用性較強(qiáng)。

［1］ Arnal D，Archambaud J P.Laminar-turbulent transition control：NLF，LFC，HLFC［R］.RTO-EN-AVT-151，2008.

［2］ Somers D M.Design and experimental results for a natural-laminar-flow airfoil for general aviation applications［R］.NASA TP-1861，1981.

［3］ Holmes B J，Obara C J，Yip L P.Natural laminar flow experiments on modern airplane surface［R］.NASA TP-2256，1984.

［4］Joslin R D.Overview of laminar control［R］.NASA TP-208705，1998.

［5］ Anderson B T，Meyer R R.Effects of wing sweep on boundarylayer transition for a smooth F-14Awing at Mach number from 0.7000to 0.825［R］.NASA TM-101712，1990.

［6］ Braslow A L.A history of suction-type laminar flow control with emphasis on flight research［R］／／Monographs in Aerospace History，1999.

［7］ Li Q，Duan Z Y，Zhang Y J，et al.Progress in research on natural laminar wing for civil aircraft［J］.Advances in Aeronautical Science and Engineering，2013，4（4）：399-406.（in Chinese）李權(quán)，段卓毅，張彥軍，等.民用飛機(jī)自然層流機(jī)翼研究進(jìn)展［J］.航空工程進(jìn)展，2013，4（4）：399-406.

［8］ Liu P Q，Ma L C，Qu Q L，et al.Numerical investigation of the laminar separation bubble control by blowing／suction on an airfoil at lowRe number［J］.Acta Aerodynamica Sinica，2013，31（4）：518-524.（in Chinese）劉沛清，馬利川，屈秋林，等.低雷諾數(shù)下翼型層流分離泡及吹吸氣控制數(shù)值研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2013，31（4）：518-524.

［9］ Xu P，Jiang C S.Optimization design of UAV laminar aerofoil based on gengtic algorithm［J］.Aero Weaonry，2009，1：7-10.（in Chinese）許平，姜長(zhǎng)生.基于遺傳算法的無(wú)人機(jī)層流翼型優(yōu)化設(shè)計(jì)［J］.航空兵器，2009，1：7-10.

［10］He J K，Bai J Q，Zhu J.Optimization design of laminar flow airfoil baesd on genetic algorithm and transition model［J］.Aeronautical Computing Technique，2010，40（3）：44-47.（in Chinese）何俊奎，白俊強(qiáng)，朱軍.基于遺傳算法及轉(zhuǎn)捩模型的層流翼型優(yōu)化設(shè)計(jì)研究［J］.航空計(jì)算技術(shù)，2010，40（3）：44-47.

［11］Collier F.NASA’s integrated system research program and environmentally responsible aviation（ERA）project［R］.Ecoaerovision，Uc Davis，2010.

［12］Ma X Y，F(xiàn)an Z L，Wu W H，et al.Aerodynamic shape optimization for wing based on NURBS［J］.Acta Aeronautica et Astronautica Sinica，2011，32（9）：1616-1621.（in Chinese）doi：CNKI：11-1929／V.20110419.1702.004.馬曉永，范召林，吳文華，等.基于非均勻有理B樣條方法的機(jī)翼氣動(dòng)外形優(yōu)化［J］.航空學(xué)報(bào)，2011，32（9）：1616-1621.

［13］Forrester T J，Edward N T，Yu N J.Thirty years of development and application of CFD at Boeing commercial airplanes，seattle［R］.AIAA 2003-3439，2003.

［14］Liu R Z，Zhang X D，An B T，et al.Application of non-uniform rational B-spline curve and knot insertion algorithm to turbine blade optimization［J］.Journal of Aerospace Power，2010，25（2）：451-458.

［15］Jason E H，David W Z.Integrated geometry parametrization and grid movement using B-spline meshes［R］.AIAA 2008-6079，2008.

［16］Sripavadkul V，Padulo M.A comparison of airfoil shape parameterization techniques for early design optimization［R］.AIAA 2010-9050，2010.

［17］Kulfan B M.A universal parametric geometry representation method-CST［R］.AIAA 2007-62，2007.

［18］Guan X H，Li Z K，Song B F.A study on CST aerodynamic shape parameterization method［J］.Acta Aeronautica et Astronautica Sinica，2012，33（4）：625-633.（in Chinese）doi：CNKI：11-1929／V.20111011.1411.005.關(guān)曉輝，李占科，宋筆鋒.CST氣動(dòng)外形參數(shù)化方法研究［J］.航空學(xué)報(bào)，2012，33（4）：625-633.

［19］Zhang L，Chen H Q.Research on genetic algorithm for aerodynamic shape optimization based on CST［J］.Aeronautical Computing Technique，2011，41（6）：53-57.（in Chinese）張磊，陳紅全.基于CST參數(shù)化的翼型優(yōu)化遺傳算法研究［J］.航空計(jì)算技術(shù)，2011，41（6）：53-57.

［20］Kulfan B M.Recent extensions and applications of the“CST”universal parametric geometry representation method［R］.AIAA 2007-7709，2007.

［21］TRANAIR User′s Manual［M］.The Boeing Company，2007.

Study of aerodynamic shape optimization for natural laminar wing

Ma Xiaoyong1，2，＊，Zhang Yanjun3，Duan Zhuoyi3，Guo Hongtao2，Li Quan3
（1.State Key Laboratory of Aerodynamics，China Aerodynamics Research and Development Center，Mianyang 621000，China；2.High Speed Aerodynamic Institute of China Aerodynamics Research and Development Center，Mianyang 622762，China；3.The First Aircraft Institute，Aviation Industry Corporation of China，Xi′an 710089，China）

Natural laminar wing is a potential technology to minimize the drag coefficient of an aircraft.The code of TRANAIR is used to optimize a natural laminar wing with subject to the constraints on lift coefficient，pitching moment coefficient and geometries height（20%，80%and max height），and the optimization are carried out with multi-points conditions，CL＝0.45（M：0.735，0.755，0.765，0.775）and CL＝0.50（M：0.755）.The nonlinear full potential equation is used to solve flow field，sequential quadratic programming and adaptive refinement Cartesian grids are applied to optimize the aerodynamic shape.Two parameter methods，non-uniform relational B-spline（NURBS）and class function／shape function transformation（CST），are used to represent wing shape and update the new geometry，although the CST method is usually used to describe aerofoil and the NURBS is more applied on three-dimensional curved surface.The optimization results shown that the NURBS method is better than the CST one，and the transition position of laminar flow to turbulent flow are moved backward obviously，the laminar region are more wide than before，and the CDof the optimized wing decreased 4.5%at CL＝0.45（M：0.755）.The optimization method is usable and effective for aerodynamic shape design such as natural laminar wing.

natural laminar wing；aerodynamic shape optimization；TRANAIR；parameterization；numerical simulation

V211.41；V224

：Adoi：10.7638／kqdlxxb-2014.0054

2014-06-10；

2014-12-02

馬曉永＊（1979-），男，河南鄭州人，博士，副研究員，主要從事飛行器氣動(dòng)布局評(píng)估及優(yōu)化設(shè)計(jì)研究.E-mail：mxycardc＠m(xù)ail.ustc.edu.cn

馬曉永，張彥軍，段卓毅，等.自然層流機(jī)翼氣動(dòng)外形優(yōu)化研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2015，33（6）：812-817.

10.7638／kqdlxxb-2014.0054 Ma X Y，Zhang Y J，Duan Z Y，et al.Study of aerodynamic shape optimization for natural laminar wing［J］.Acta Aerodynamica Sinica，2015，33（6）：812-817.

0258-1825（2015）06-0812-06