徐林程，王 剛，武 潔，葉正寅

（西北工業(yè)大學，翼型、葉柵空氣動力學國防科技重點實驗室，陜西 西安 710072）

其中：

翼型模型幾何誤差對氣動性能影響的自動微分分析方法

徐林程，王 剛，武 潔，葉正寅

（西北工業(yè)大學，翼型、葉柵空氣動力學國防科技重點實驗室，陜西 西安 710072）

基于自動微分原理和 NS方程有限體積方法建立了一套翼型敏感性導數(shù)計算方法和程序，可以一次性獲得翼型不同氣動力系數(shù)、壓力分布對模型幾何外形誤差的敏感性導數(shù)和不確定度。計算結果表明，在跨聲速范圍內(nèi)，即使翼型的外形誤差只有63μm（弦長1m），也可以給翼型壓力分布帶來0.312（以來流動壓為參考）量級的不確定度，而激波處的流動最為敏感。這種自動微分方法對于分析數(shù)值模擬結果的分散性、風洞試驗結果的分散性或不確定性具有很好的指導意義。

幾何誤差；不確定性；敏感性分析；自動微分；跨聲速；翼型

0 引 言

氣動外形決定了飛行器的氣動性能，無論是數(shù)值計算還是風洞試驗，人們一直致力于提升對某個外形氣動性能的預測精度。從風洞試驗角度，外形加工誤差、模型結構本身變形（可以是因為載荷作用、也可以是因為氣候變化帶來的物理變形）都會引起外形的變化；從數(shù)值計算角度講，幾何外形的連續(xù)曲線會被離散為折線段，網(wǎng)格劃分策略和網(wǎng)格數(shù)量就決定了計算外形與實際外形的差別，而且在網(wǎng)格生成過程中，插值運算本身也會在某種程度上引起幾何誤差，這些外形的改變到底會給氣動性能的預測帶來多大的影響？這是一個有實際意義的基本問題，雖然問題“顯而易見”，但長期以來沒有一個定量的分析方法研究該問題。解決這個問題的方案可以從國外近年來采用的不確定性分析途徑進行分析［1-7］。在不確定性分析過程中，最重要的技術環(huán)境就是獲得不同的敏感性導數(shù)，對于導數(shù)的求解方法，最直接的方式就是利用差商法獲取，隨著CFD技術的發(fā)展和成熟，利用CFD工具計算這些敏感性導數(shù)理論上是可以實現(xiàn)的，但是，對于眾多的敏感性導數(shù)，如果采用簡單的差分算法意味著巨大的計算量。而且，對于幾何外形很小的變化，要準確地計算出他們之間的氣動性能差異，對計算方法和軟件是一個很大的挑戰(zhàn)。此外，CFD軟件計算的結果還受到網(wǎng)格數(shù)量、收斂精度等因素的影響，而差分算法中步長的選取也對敏感性導數(shù)值產(chǎn)生影響。為了更有效地獲得敏感性導數(shù)，國外引入自動微分 方法［8-11］，這種 方 法直 接 伴隨 CFD 方 法 求解 空氣動力學基本方程的實際流程，敏感性導數(shù)的計算也是數(shù)學意義上嚴格的微分概念，更有價值的是，只要開發(fā)出的計算程序設計合適，可以一次性同時獲得大量的敏感性導數(shù)，而且敏感性導數(shù)的收斂精度與流場的收斂精度達到相同的量級，為不確定性分析開辟了一條新的技術途徑。

本文以跨聲速翼型為對象，在課題組CFD程序基礎上，應用自動微分工具 Tapenade對程序進行改造，對典型超臨界翼型的幾何敏感性導數(shù)進行計算，給出壓力系數(shù)分布對各處幾何外形的不確定性分析。這項工作對于認識數(shù)值計算和風洞試驗中外形參數(shù)變化引起的性能預測誤差有重要的實際意義。

1 自動微分方法

自動微分是面向計算程序，應用求導的鏈式法則，求取程序輸出量對輸入量偏導數(shù)的一種方法。任何計算程序是由有限個賦值序列組成，即

而編譯器會將所有的數(shù)值運算分解為計算機硬件能夠執(zhí)行的加，減，乘，除四則基本運算，由于四則基本運算在數(shù)學上都是顯式析可微且可導的，因此所有計算機程序能表達的數(shù)學運算，理論上是至少單側可微的，即偏導數(shù)

都解析存在，并且可以顯示表達。令

其中D為求導算子，D f為程序中賦值語句基于自身數(shù)學結構對所有變量求導，Dτ為變量基于程序全局數(shù)量關系對所有變量求導，展開后得可得：

對Dτ的每一項應用鏈式法則：

表示為矩陣形式，即為

對上式兩端同時進行轉置運算，可得

其中I為單位矩陣。

方程（1）給出了自動微分前向求導模式的求解順序和算法

方程（2）給出了自動微分逆向求導模式的求解順序和算法

由于上述推導過程沒有引入數(shù)學或物理假設，不包含除法運算，并且每個環(huán)節(jié)都是嚴格精確的解析過程，因此，自動微分求導過程結果沒有截斷誤差，精度與原程序相同。對于線性問題，自動微分過程的收斂速度，穩(wěn)定性與原程序一致［13］。此外，單目標函數(shù)的逆向自動微分過程計算量最多不超過原程序計算量的5倍已得到理論上的證明［15］。

2 空氣動力學計算方法

本文以跨聲速馬赫數(shù)范圍的二維超臨界翼型問題為研究對象，以課題組非結構網(wǎng)格求解雷諾平均N-S方程的計算程序為基礎，應用Tapenade［15］軟件進行自動微分方法的進一步改造。原程序以雷諾平均的可壓縮、非定常N-S方程為流動基本方程，采用目前認為通用性較好的SA湍流模型，應用有限體積法對空氣動力學基本方程進行離散求解［12］，并且運用了當?shù)貢r間步長、隱式時間推進等加速收斂措施。

流動基本方程的積分形式為：

其中：

ρ、u、v、e分別為空氣密度，x、y方向的速度分量和單位體積的總內(nèi)能，n為線積分的法向單位向量，V 為積分域，?V 為積分域的邊界，F(xiàn)為通量項，它包括無粘項FE和粘性項Fv，V 為速度矢量，i、j分別為x、y方向的單位矢量，τ為粘性應力。

圖1 翼型附近的非結構混合網(wǎng)格Fig.1 Unstructured hybrid mesh attached to the airfoil

圖2 第一層網(wǎng)格格心的y＋分布Fig.2 Distribution of y＋value of the first grid layer

圖3 翼型的壓力系數(shù)分布Fig.3 Pressure coeffecients distribution of the airfoil

圖4 計算過程的殘值收斂歷程Fig.4 Convergence course of computing process

圖1～圖4是在迎角α=2.31°，馬赫數(shù) Ma= 0.729，雷諾數(shù)Re=650萬的狀態(tài)下，原程序計算超臨界翼型RAE2822跨聲速流場的網(wǎng)格和結果。為了模擬附面層流動，在翼型壁面附近采用了層狀的網(wǎng)格結構，第一層網(wǎng)格高度為4.0e-6。可以看出，原CFD程序具有良好的收斂特性，同時計算結果與實驗數(shù)據(jù)［16］吻合程度高，物面Y+值都小于1.2，因此網(wǎng)格是比較合理的，并且原CFD程序具有較好的計算精度。

3 計算結果分析

為了驗證在空氣動力學方程求解程序基礎上改造的自動微分程序，首先通過傳統(tǒng)的差分方法計算出翼型前緣附近某點的壓力系數(shù)對該點坐標變量的敏感性導數(shù)，然后采用改造后的自動微分程序對該敏感性導數(shù)進行計算。其結果如表1所示，可以看出，采用自動微分的方法與采用傳統(tǒng)差分方法的計算結果接近。此外，壓力系數(shù)對x坐標的敏感性要大于對y坐標敏感性導數(shù)一個量級。

表1 不同方式計算的對翼型前緣坐標的敏感性導數(shù)Table 1 Derivatives of surface pressure coefficient with respect to local coordinates from different differentiation methods

其中，In V 為自變量（Independent Variables），Var為變量（Variables），DS為差分步長（Difference Step），F(xiàn)D為前向差分（Forward Difference），AD為自動微分（Automatic Differentiation）。

在上述方法驗證的基礎上，以翼型表面坐標值為輸入變量，以翼型表面的壓力系數(shù)為輸出變量，由Tapenade軟件對源程序生成自動微分前向求導程序，進行計算。

圖5顯示了翼型表面各點壓力系數(shù)對物面各點坐標偏導數(shù)的總體分布，（a）圖針對x坐標，（b）圖針對y坐標。圖6是圖5兩圖對應的沿x方向視圖，圖中每條曲線描述的是翼型表面某一點的Cp對各點坐標的敏感程度，圖中的曲線族在橫軸某區(qū)間的聚集程度和幅值大小，表征著該區(qū)間坐標對物面各處壓力系數(shù)的總體影響能力。圖7是圖5兩圖對應的沿y方向視圖，圖中每條曲線描述的是翼型表面某一點的坐標對各點壓力系數(shù)的影響能力，圖中曲線族在橫軸某區(qū)間的聚集程度和幅值大小，表征著該區(qū)域壓力系數(shù)對物面坐標的總體敏感性。

圖5 Cp對物面坐標導數(shù)的整體分布Fig.5 Derivatives of pressure coeffecients with respect to surface coordinates

可以看到，圖5（a、b）兩圖具有一些相似的基本特征。首先，偏導數(shù)總體量級很小，但在沿著XY平面內(nèi)的對角線上有明顯的突起，這說明每個點坐標主要對當?shù)鼗蛘呦嘟c附近的流場有明顯影響。其次，偏導數(shù)在激波位置存在明顯的帶狀突起，在翼型前緣和后緣有尖銳的峰值，結合圖6，可以看到，在x=0附近的前緣附近，x=1.0的后緣附近和有激波存在的0.5＜x＜0.6區(qū)域內(nèi)，曲線高密度聚集，并且存在幅值峰，可見，這三個區(qū)域的物面坐標對流場有很強的影響能力；結合圖7，可以看到，在x=0的前緣附近和x=1.0的后緣附近和有激波存在的0.5＜x＜0.6區(qū)域內(nèi)，曲線高密度聚集，并且存在幅值峰，可見，這三個區(qū)域的流場對物面擾動的感知能力很強。此外，在翼型中后段，圖5（b）圖在對角線上的突起要比圖5（a）強很多，這主要是因為該部分的物面法向與y′坐標方向接近。

為了進一步量化幾何誤差對翼型壓力分布的影響，考慮到采用7級精度加工弦長一米的翼型模型，表面公差會達到63μm，利用已獲得的x、y方向的偏導數(shù)合成各物面點垂直物面方向的偏導數(shù)，取物面在法向的攝動為63μm，則壓力分布的不確定度如圖8所示，其中，圖8（a）是翼型壓力系數(shù)不確定性分布的單獨顯示，圖8（b）是翼型壓力系數(shù)不確定性附著在壓力系數(shù)上的顯示。

圖6 物面坐標影響曲線云圖Fig.6 Effect curves of surface coordinates to pressure coeffecients

對照圖8中的 （a）、（b）兩圖，不難發(fā)現(xiàn)翼型上表面壓力分布的不確定度在前緣和激波位置出現(xiàn)峰值，并且在激波處取得最大值0.312，顯然，在這兩個位置之間的超聲速區(qū)域，壓力系數(shù)不確定度明顯高于其余的亞聲速區(qū)域，這說明，翼型在跨聲速流中，超聲速區(qū)比亞聲速區(qū)要對物面坐標敏感。

圖7 Cp敏感性曲線云圖Fig.7 Sensitivity curves of pressure coeffecients with respect to surface coordinates

圖8 物面坐標不確定度誘導的壓力分布不確定度Fig.8 Uncertainties of surface pressure coeffecients caused by surface coordinates perturbations

圖8（a）表明，在物面63μm 的攝動下，壓力系數(shù)的不確定度會達到0.312（以來流動壓為參考）的量級，可以從翼型模型加工和數(shù)值模擬兩個角度看待這個結果。

從翼型模型加工的角度看，要獲得更高精度的風洞試驗結果，必須采用更高精度的加工方法制作翼型模型。此外，由于不同區(qū)域物面坐標對流場的影響能力相差懸殊，在加工翼型模型過程中，對敏感區(qū)域（需要預先進行幾何外形的敏感性分析），例如：翼型前緣，激波位置（需要預先估測）等，進行特殊處理，單獨提高其加工精度，將對提高風洞試驗結果的精度取得事半功倍的效果。

從數(shù)值模擬的角度看，由于物面總是被離散成有限個點，用折線代替曲線，如圖9所示，顯然，這個過程會引起外形誤差。根據(jù)幾何關系可得：

其中，δ、ρ、Δl分別表示外形誤差，當?shù)厍拾霃剑W(wǎng)格長度。令ε為計算結果所需要達到的精度，χ 為流場對物面幾何誤差的敏感程度（這里為物面壓力系數(shù)對物面幾何誤差敏感性的最大值，即0.312×106／63 ≈5000），則當?shù)鼐W(wǎng)格長度應滿足如下關系定量關系：

上式是從敏感性角度出發(fā)，得到的當?shù)鼐W(wǎng)格長度最大值的表達式，以下簡稱為物面網(wǎng)格密度的敏感性準則?？梢?，CFD數(shù)值模擬要獲得指定精度的結果，物面網(wǎng)格密度必須滿足式（6）才是可能的。

圖9 曲線的離散公差Fig.9 Geometry error of curves caused by dicretisation

4 結 論

采用有限體積法求解雷諾平均 Navier-Stokes方程的程序為基礎，運用自動微分方法建立了一套翼型敏感性導數(shù)計算方法和程序，在此基礎上，研究了在跨聲速流場中模型幾何誤差對翼型表面壓力分布的影響。研究表明，在跨聲速范圍內(nèi)：

物面坐標主要影響當?shù)氐牧鲌?，但是，激波附近的物面坐標會顯著影響物面各處的流場，同時，激波位置的流場能夠明顯感受到物面各處坐標的存在，并且，在翼型前緣區(qū)域和后緣區(qū)域，物面坐標對流場的影響能力和流場對物面坐標的敏感性都很強。

7級精度的加工誤差，在本文所取的來流狀態(tài)下，會給翼型表面壓力系數(shù)分布帶來0.312（以來流動壓為參考）量級的不確定度。降低加工誤差對風洞試驗結果的影響需要更高精度的加工手段。數(shù)值模擬中，物面網(wǎng)格密度滿足敏感性準則（詳見式（6）），計算結果達到預期精度才能有實際物理意義。同時，翼型上表面的超聲速區(qū)域流場明顯比翼型表面亞聲速區(qū)域流場對幾何誤差更敏感，并且翼型前緣和激波位置會出現(xiàn)敏感性尖峰，這從一個方面解釋了為什么在風洞試驗和數(shù)值模擬中，上表面壓力分布分散性比下表面大；翼型吸力峰值附近和激波附近的壓力系數(shù)分散性大。

［1］ WEAVER A B，ALEXEENKO A A，et al.Flow field uncertainty analysis for hypersonic CFD simulations［R］.AIAA Paper 2010-1180.

［2］ BETTIS B R，HOSDER S.Uncertainty quantification in hypersonic reentry flows due to aleatory and epistemic uncertainties ［R］.AIAA Paper 2011-0252.

［3］ DEEPAK BOSE D，BROWN J L，et al.Uncertainty assessment of hypersonic aerothermodynamics prediction capability ［J］.Journal of Spacecraft and Rockets，2013，50（1）：12-18.

［4］ KUCHI-ISHI S.Uncertainty evaluation of thermocouple aeroheating measurements for hypersonic wind-tunnel tests［J］.Journal of Spacecraft and Rockets，2006，43（3）：698-700.

［5］ KAMMEYER M E，RUEGER M L.On the classification of errors：systematic，random，and replication level［R］.AIAA Paper 2000-2203

［6］ ULBRICH N.Test data uncertainty analysis algorithm of NASA Ames wind tunnels［J］.AIAA Journal，2003，41（8）：1616-1619.

［7］ KAMMEYER M E，WOZNIAK R W.An uncertainty analysis for low-speed wind tunnel pressure measurements［R］.AIAA Paper 2004-2196.

［8］ MADER C A，MARTINS J R R A.Computation of aircraft stability derivatives using an automatic differentiation adjoint approach［J］.AIAA Journal，2011，49（12）：2737-2750.

［9］ JONES D，MUELLER J D，BAYYUK S.CFD development with automatic differentiation［R］.AIAA Paper 2012-0573.

［10］ZHOUJIE LYU Z，GAETAN K W，KENWAY G K W，et al.Automatic differentiation adjoint of the reynolds-averaged Navier-Stokes equations with a turbulence model［R］.AIAA Paper 2013-2581.

［11］DALLE D J，DRISCOLL J F.Continuous differentiation of complex systems applied to a hypersonic vehicle［R］.AIAA Paper 2012-4958.

［12］JIANG Y W，YE Z Y，WANG G.Efficient solution of Euler／N-S equations on unstructrured grids［J］.Chinese Journal of Computional Mechanics，2012，29（2）：217-233.（in Chinese）蔣躍文，葉正寅，王剛.基于非結構網(wǎng)格的高效求解方法研究［J］.計算力學學報，2012，29（2）：217-223.

［13］GILES M B，GHATE D，DUTA M C.Using automatic differentiation for adjoint CFD code development［R］.Oxford：Oxford University Computing Laboratory.2005.

［14］GRIEWANK A.On automatic differentiation.in：mathematical programming’88［M］.Kluwer Academic Publishers，Dordrecht，1989.

［15］HASCOЁT L，PASCUAL V.TAPENADE 2.1 user’s guide ［R］.INRIA，Sophia-Anpipolis，2004

［16］THIBERT J J，GRANJACQUES M，OHMAN L H.RAE2822 airfoil agard advisory report No.138［R］.Experimental Data Base for Computer Program Assessment，1979，A182.

An automatic differentiation method for uncertainty analysis due to airfoil configuration variation

XU Lincheng，WANG Gang，WU Jie，YE Zhengyin

（National Key Laboratory of Aerodynamic Design and Research，NWPU，Xi’an 710072，China）

Focused on the quantification of the uncertainties of areodynamics performance of airfoils with respect to geometry error，with a set of CFD program based on finite volume algorithm solving the Reynolds-Averaged Navier-Stokes equations with S-A turbulent model，adopting automatic differentiation method to reform the program simultaniously，all kinds of sensitive derivatives，uncertainties of all kinds of aerodynamic coefficients and pressure coefficients distribution resulting from geometry error could be obtained in one course of computation.As the computational results show，even if the geometry error is only 63 microns （while the length of chord is 1 meter），the pressure distribution of the walls could be influenced obviously with uncertainty quantity reaching 0.312（taking dynamic pressure of the flow as reference）for an airfois in transonic flow，moreover，pressure attached to the place where shock wave stationed bears peak uncertainty.the results of method of automatic differentiation account for the dispersity of results of numeric simulations and wind tunnel experiments well.

geometry error；uncertainty；sensitivity analysis；automatic differentiation；transonic flow；airfoils

V211.3

Adoi：10.7638／kqdlxxb-2013.0097

0258-1825（2014）04-0551-06

2013-10-15；

2014-01-02

國家自然科學基金（11272264）

徐林程（1989-），男，研究生，研究方向：理論與計算流體力學.E-mail：993644844＠qq.com

徐林程，王 剛，武 潔，等.翼型模型幾何誤差對氣動性能影響的自動微分分析方法［J］.空氣動力學學報，2014，32（4）：551-556.

10.7638／kqdlxxb-2013.0097. XU L C，WANG G，WU J，et al.An automatic differentiation method for uncertainty analysis due to airfoil configuration variation［J］.ACTA Aerodynamica Sinica，2014，32（4）：551-556.