徐林程,王 剛,武 潔,葉正寅
(西北工業(yè)大學(xué),翼型、葉柵空氣動(dòng)力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710072)
其中:
翼型模型幾何誤差對(duì)氣動(dòng)性能影響的自動(dòng)微分分析方法
徐林程,王 剛,武 潔,葉正寅
(西北工業(yè)大學(xué),翼型、葉柵空氣動(dòng)力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710072)
基于自動(dòng)微分原理和 NS方程有限體積方法建立了一套翼型敏感性導(dǎo)數(shù)計(jì)算方法和程序,可以一次性獲得翼型不同氣動(dòng)力系數(shù)、壓力分布對(duì)模型幾何外形誤差的敏感性導(dǎo)數(shù)和不確定度。計(jì)算結(jié)果表明,在跨聲速范圍內(nèi),即使翼型的外形誤差只有63μm(弦長1m),也可以給翼型壓力分布帶來0.312(以來流動(dòng)壓為參考)量級(jí)的不確定度,而激波處的流動(dòng)最為敏感。這種自動(dòng)微分方法對(duì)于分析數(shù)值模擬結(jié)果的分散性、風(fēng)洞試驗(yàn)結(jié)果的分散性或不確定性具有很好的指導(dǎo)意義。
幾何誤差;不確定性;敏感性分析;自動(dòng)微分;跨聲速;翼型
氣動(dòng)外形決定了飛行器的氣動(dòng)性能,無論是數(shù)值計(jì)算還是風(fēng)洞試驗(yàn),人們一直致力于提升對(duì)某個(gè)外形氣動(dòng)性能的預(yù)測精度。從風(fēng)洞試驗(yàn)角度,外形加工誤差、模型結(jié)構(gòu)本身變形(可以是因?yàn)檩d荷作用、也可以是因?yàn)闅夂蜃兓瘞淼奈锢碜冃危┒紩?huì)引起外形的變化;從數(shù)值計(jì)算角度講,幾何外形的連續(xù)曲線會(huì)被離散為折線段,網(wǎng)格劃分策略和網(wǎng)格數(shù)量就決定了計(jì)算外形與實(shí)際外形的差別,而且在網(wǎng)格生成過程中,插值運(yùn)算本身也會(huì)在某種程度上引起幾何誤差,這些外形的改變到底會(huì)給氣動(dòng)性能的預(yù)測帶來多大的影響?這是一個(gè)有實(shí)際意義的基本問題,雖然問題“顯而易見”,但長期以來沒有一個(gè)定量的分析方法研究該問題。解決這個(gè)問題的方案可以從國外近年來采用的不確定性分析途徑進(jìn)行分析[1-7]。在不確定性分析過程中,最重要的技術(shù)環(huán)境就是獲得不同的敏感性導(dǎo)數(shù),對(duì)于導(dǎo)數(shù)的求解方法,最直接的方式就是利用差商法獲取,隨著CFD技術(shù)的發(fā)展和成熟,利用CFD工具計(jì)算這些敏感性導(dǎo)數(shù)理論上是可以實(shí)現(xiàn)的,但是,對(duì)于眾多的敏感性導(dǎo)數(shù),如果采用簡單的差分算法意味著巨大的計(jì)算量。而且,對(duì)于幾何外形很小的變化,要準(zhǔn)確地計(jì)算出他們之間的氣動(dòng)性能差異,對(duì)計(jì)算方法和軟件是一個(gè)很大的挑戰(zhàn)。此外,CFD軟件計(jì)算的結(jié)果還受到網(wǎng)格數(shù)量、收斂精度等因素的影響,而差分算法中步長的選取也對(duì)敏感性導(dǎo)數(shù)值產(chǎn)生影響。為了更有效地獲得敏感性導(dǎo)數(shù),國外引入自動(dòng)微分 方法[8-11],這種 方 法直 接 伴隨 CFD 方 法 求解 空氣動(dòng)力學(xué)基本方程的實(shí)際流程,敏感性導(dǎo)數(shù)的計(jì)算也是數(shù)學(xué)意義上嚴(yán)格的微分概念,更有價(jià)值的是,只要開發(fā)出的計(jì)算程序設(shè)計(jì)合適,可以一次性同時(shí)獲得大量的敏感性導(dǎo)數(shù),而且敏感性導(dǎo)數(shù)的收斂精度與流場的收斂精度達(dá)到相同的量級(jí),為不確定性分析開辟了一條新的技術(shù)途徑。
本文以跨聲速翼型為對(duì)象,在課題組CFD程序基礎(chǔ)上,應(yīng)用自動(dòng)微分工具 Tapenade對(duì)程序進(jìn)行改造,對(duì)典型超臨界翼型的幾何敏感性導(dǎo)數(shù)進(jìn)行計(jì)算,給出壓力系數(shù)分布對(duì)各處幾何外形的不確定性分析。這項(xiàng)工作對(duì)于認(rèn)識(shí)數(shù)值計(jì)算和風(fēng)洞試驗(yàn)中外形參數(shù)變化引起的性能預(yù)測誤差有重要的實(shí)際意義。
自動(dòng)微分是面向計(jì)算程序,應(yīng)用求導(dǎo)的鏈?zhǔn)椒▌t,求取程序輸出量對(duì)輸入量偏導(dǎo)數(shù)的一種方法。任何計(jì)算程序是由有限個(gè)賦值序列組成,即
而編譯器會(huì)將所有的數(shù)值運(yùn)算分解為計(jì)算機(jī)硬件能夠執(zhí)行的加,減,乘,除四則基本運(yùn)算,由于四則基本運(yùn)算在數(shù)學(xué)上都是顯式析可微且可導(dǎo)的,因此所有計(jì)算機(jī)程序能表達(dá)的數(shù)學(xué)運(yùn)算,理論上是至少單側(cè)可微的,即偏導(dǎo)數(shù)
都解析存在,并且可以顯示表達(dá)。令
其中D為求導(dǎo)算子,D f為程序中賦值語句基于自身數(shù)學(xué)結(jié)構(gòu)對(duì)所有變量求導(dǎo),Dτ為變量基于程序全局?jǐn)?shù)量關(guān)系對(duì)所有變量求導(dǎo),展開后得可得:
對(duì)Dτ的每一項(xiàng)應(yīng)用鏈?zhǔn)椒▌t:
表示為矩陣形式,即為
對(duì)上式兩端同時(shí)進(jìn)行轉(zhuǎn)置運(yùn)算,可得
其中I為單位矩陣。
方程(1)給出了自動(dòng)微分前向求導(dǎo)模式的求解順序和算法
方程(2)給出了自動(dòng)微分逆向求導(dǎo)模式的求解順序和算法
由于上述推導(dǎo)過程沒有引入數(shù)學(xué)或物理假設(shè),不包含除法運(yùn)算,并且每個(gè)環(huán)節(jié)都是嚴(yán)格精確的解析過程,因此,自動(dòng)微分求導(dǎo)過程結(jié)果沒有截?cái)嗾`差,精度與原程序相同。對(duì)于線性問題,自動(dòng)微分過程的收斂速度,穩(wěn)定性與原程序一致[13]。此外,單目標(biāo)函數(shù)的逆向自動(dòng)微分過程計(jì)算量最多不超過原程序計(jì)算量的5倍已得到理論上的證明[15]。
本文以跨聲速馬赫數(shù)范圍的二維超臨界翼型問題為研究對(duì)象,以課題組非結(jié)構(gòu)網(wǎng)格求解雷諾平均N-S方程的計(jì)算程序?yàn)榛A(chǔ),應(yīng)用Tapenade[15]軟件進(jìn)行自動(dòng)微分方法的進(jìn)一步改造。原程序以雷諾平均的可壓縮、非定常N-S方程為流動(dòng)基本方程,采用目前認(rèn)為通用性較好的SA湍流模型,應(yīng)用有限體積法對(duì)空氣動(dòng)力學(xué)基本方程進(jìn)行離散求解[12],并且運(yùn)用了當(dāng)?shù)貢r(shí)間步長、隱式時(shí)間推進(jìn)等加速收斂措施。
流動(dòng)基本方程的積分形式為:
其中:
ρ、u、v、e分別為空氣密度,x、y方向的速度分量和單位體積的總內(nèi)能,n為線積分的法向單位向量,V 為積分域,?V 為積分域的邊界,F(xiàn)為通量項(xiàng),它包括無粘項(xiàng)FE和粘性項(xiàng)Fv,V 為速度矢量,i、j分別為x、y方向的單位矢量,τ為粘性應(yīng)力。
圖1 翼型附近的非結(jié)構(gòu)混合網(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 計(jì)算過程的殘值收斂歷程Fig.4 Convergence course of computing process
圖1~圖4是在迎角α=2.31°,馬赫數(shù) Ma= 0.729,雷諾數(shù)Re=650萬的狀態(tài)下,原程序計(jì)算超臨界翼型RAE2822跨聲速流場的網(wǎng)格和結(jié)果。為了模擬附面層流動(dòng),在翼型壁面附近采用了層狀的網(wǎng)格結(jié)構(gòu),第一層網(wǎng)格高度為4.0e-6??梢钥闯?,原CFD程序具有良好的收斂特性,同時(shí)計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)[16]吻合程度高,物面Y+值都小于1.2,因此網(wǎng)格是比較合理的,并且原CFD程序具有較好的計(jì)算精度。
為了驗(yàn)證在空氣動(dòng)力學(xué)方程求解程序基礎(chǔ)上改造的自動(dòng)微分程序,首先通過傳統(tǒng)的差分方法計(jì)算出翼型前緣附近某點(diǎn)的壓力系數(shù)對(duì)該點(diǎn)坐標(biāo)變量的敏感性導(dǎo)數(shù),然后采用改造后的自動(dòng)微分程序?qū)υ撁舾行詫?dǎo)數(shù)進(jìn)行計(jì)算。其結(jié)果如表1所示,可以看出,采用自動(dòng)微分的方法與采用傳統(tǒng)差分方法的計(jì)算結(jié)果接近。此外,壓力系數(shù)對(duì)x坐標(biāo)的敏感性要大于對(duì)y坐標(biāo)敏感性導(dǎo)數(shù)一個(gè)量級(jí)。
表1 不同方式計(jì)算的對(duì)翼型前緣坐標(biāo)的敏感性導(dǎo)數(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為自動(dòng)微分(Automatic Differentiation)。
在上述方法驗(yàn)證的基礎(chǔ)上,以翼型表面坐標(biāo)值為輸入變量,以翼型表面的壓力系數(shù)為輸出變量,由Tapenade軟件對(duì)源程序生成自動(dòng)微分前向求導(dǎo)程序,進(jìn)行計(jì)算。
圖5顯示了翼型表面各點(diǎn)壓力系數(shù)對(duì)物面各點(diǎn)坐標(biāo)偏導(dǎo)數(shù)的總體分布,(a)圖針對(duì)x坐標(biāo),(b)圖針對(duì)y坐標(biāo)。圖6是圖5兩圖對(duì)應(yīng)的沿x方向視圖,圖中每條曲線描述的是翼型表面某一點(diǎn)的Cp對(duì)各點(diǎn)坐標(biāo)的敏感程度,圖中的曲線族在橫軸某區(qū)間的聚集程度和幅值大小,表征著該區(qū)間坐標(biāo)對(duì)物面各處壓力系數(shù)的總體影響能力。圖7是圖5兩圖對(duì)應(yīng)的沿y方向視圖,圖中每條曲線描述的是翼型表面某一點(diǎn)的坐標(biāo)對(duì)各點(diǎn)壓力系數(shù)的影響能力,圖中曲線族在橫軸某區(qū)間的聚集程度和幅值大小,表征著該區(qū)域壓力系數(shù)對(duì)物面坐標(biāo)的總體敏感性。
圖5 Cp對(duì)物面坐標(biāo)導(dǎo)數(shù)的整體分布Fig.5 Derivatives of pressure coeffecients with respect to surface coordinates
可以看到,圖5(a、b)兩圖具有一些相似的基本特征。首先,偏導(dǎo)數(shù)總體量級(jí)很小,但在沿著XY平面內(nèi)的對(duì)角線上有明顯的突起,這說明每個(gè)點(diǎn)坐標(biāo)主要對(duì)當(dāng)?shù)鼗蛘呦嘟c(diǎn)附近的流場有明顯影響。其次,偏導(dǎo)數(shù)在激波位置存在明顯的帶狀突起,在翼型前緣和后緣有尖銳的峰值,結(jié)合圖6,可以看到,在x=0附近的前緣附近,x=1.0的后緣附近和有激波存在的0.5<x<0.6區(qū)域內(nèi),曲線高密度聚集,并且存在幅值峰,可見,這三個(gè)區(qū)域的物面坐標(biāo)對(duì)流場有很強(qiáng)的影響能力;結(jié)合圖7,可以看到,在x=0的前緣附近和x=1.0的后緣附近和有激波存在的0.5<x<0.6區(qū)域內(nèi),曲線高密度聚集,并且存在幅值峰,可見,這三個(gè)區(qū)域的流場對(duì)物面擾動(dòng)的感知能力很強(qiáng)。此外,在翼型中后段,圖5(b)圖在對(duì)角線上的突起要比圖5(a)強(qiáng)很多,這主要是因?yàn)樵摬糠值奈锩娣ㄏ蚺cy′坐標(biāo)方向接近。
為了進(jìn)一步量化幾何誤差對(duì)翼型壓力分布的影響,考慮到采用7級(jí)精度加工弦長一米的翼型模型,表面公差會(huì)達(dá)到63μm,利用已獲得的x、y方向的偏導(dǎo)數(shù)合成各物面點(diǎn)垂直物面方向的偏導(dǎo)數(shù),取物面在法向的攝動(dòng)為63μm,則壓力分布的不確定度如圖8所示,其中,圖8(a)是翼型壓力系數(shù)不確定性分布的單獨(dú)顯示,圖8(b)是翼型壓力系數(shù)不確定性附著在壓力系數(shù)上的顯示。
圖6 物面坐標(biāo)影響曲線云圖Fig.6 Effect curves of surface coordinates to pressure coeffecients
對(duì)照?qǐng)D8中的 (a)、(b)兩圖,不難發(fā)現(xiàn)翼型上表面壓力分布的不確定度在前緣和激波位置出現(xiàn)峰值,并且在激波處取得最大值0.312,顯然,在這兩個(gè)位置之間的超聲速區(qū)域,壓力系數(shù)不確定度明顯高于其余的亞聲速區(qū)域,這說明,翼型在跨聲速流中,超聲速區(qū)比亞聲速區(qū)要對(duì)物面坐標(biāo)敏感。
圖7 Cp敏感性曲線云圖Fig.7 Sensitivity curves of pressure coeffecients with respect to surface coordinates
圖8 物面坐標(biāo)不確定度誘導(dǎo)的壓力分布不確定度Fig.8 Uncertainties of surface pressure coeffecients caused by surface coordinates perturbations
圖8(a)表明,在物面63μm 的攝動(dòng)下,壓力系數(shù)的不確定度會(huì)達(dá)到0.312(以來流動(dòng)壓為參考)的量級(jí),可以從翼型模型加工和數(shù)值模擬兩個(gè)角度看待這個(gè)結(jié)果。
從翼型模型加工的角度看,要獲得更高精度的風(fēng)洞試驗(yàn)結(jié)果,必須采用更高精度的加工方法制作翼型模型。此外,由于不同區(qū)域物面坐標(biāo)對(duì)流場的影響能力相差懸殊,在加工翼型模型過程中,對(duì)敏感區(qū)域(需要預(yù)先進(jìn)行幾何外形的敏感性分析),例如:翼型前緣,激波位置(需要預(yù)先估測)等,進(jìn)行特殊處理,單獨(dú)提高其加工精度,將對(duì)提高風(fēng)洞試驗(yàn)結(jié)果的精度取得事半功倍的效果。
從數(shù)值模擬的角度看,由于物面總是被離散成有限個(gè)點(diǎn),用折線代替曲線,如圖9所示,顯然,這個(gè)過程會(huì)引起外形誤差。根據(jù)幾何關(guān)系可得:
其中,δ、ρ、Δl分別表示外形誤差,當(dāng)?shù)厍拾霃?,網(wǎng)格長度。令ε為計(jì)算結(jié)果所需要達(dá)到的精度,χ 為流場對(duì)物面幾何誤差的敏感程度(這里為物面壓力系數(shù)對(duì)物面幾何誤差敏感性的最大值,即0.312×106/63 ≈5000),則當(dāng)?shù)鼐W(wǎng)格長度應(yīng)滿足如下關(guān)系定量關(guān)系:
上式是從敏感性角度出發(fā),得到的當(dāng)?shù)鼐W(wǎng)格長度最大值的表達(dá)式,以下簡稱為物面網(wǎng)格密度的敏感性準(zhǔn)則??梢?,CFD數(shù)值模擬要獲得指定精度的結(jié)果,物面網(wǎng)格密度必須滿足式(6)才是可能的。
圖9 曲線的離散公差Fig.9 Geometry error of curves caused by dicretisation
采用有限體積法求解雷諾平均 Navier-Stokes方程的程序?yàn)榛A(chǔ),運(yùn)用自動(dòng)微分方法建立了一套翼型敏感性導(dǎo)數(shù)計(jì)算方法和程序,在此基礎(chǔ)上,研究了在跨聲速流場中模型幾何誤差對(duì)翼型表面壓力分布的影響。研究表明,在跨聲速范圍內(nèi):
物面坐標(biāo)主要影響當(dāng)?shù)氐牧鲌?,但是,激波附近的物面坐?biāo)會(huì)顯著影響物面各處的流場,同時(shí),激波位置的流場能夠明顯感受到物面各處坐標(biāo)的存在,并且,在翼型前緣區(qū)域和后緣區(qū)域,物面坐標(biāo)對(duì)流場的影響能力和流場對(duì)物面坐標(biāo)的敏感性都很強(qiáng)。
7級(jí)精度的加工誤差,在本文所取的來流狀態(tài)下,會(huì)給翼型表面壓力系數(shù)分布帶來0.312(以來流動(dòng)壓為參考)量級(jí)的不確定度。降低加工誤差對(duì)風(fēng)洞試驗(yàn)結(jié)果的影響需要更高精度的加工手段。數(shù)值模擬中,物面網(wǎng)格密度滿足敏感性準(zhǔn)則(詳見式(6)),計(jì)算結(jié)果達(dá)到預(yù)期精度才能有實(shí)際物理意義。同時(shí),翼型上表面的超聲速區(qū)域流場明顯比翼型表面亞聲速區(qū)域流場對(duì)幾何誤差更敏感,并且翼型前緣和激波位置會(huì)出現(xiàn)敏感性尖峰,這從一個(gè)方面解釋了為什么在風(fēng)洞試驗(yà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)蔣躍文,葉正寅,王剛.基于非結(jié)構(gòu)網(wǎng)格的高效求解方法研究[J].計(jì)算力學(xué)學(xué)報(bào),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
國家自然科學(xué)基金(11272264)
徐林程(1989-),男,研究生,研究方向:理論與計(jì)算流體力學(xué).E-mail:993644844@qq.com
徐林程,王 剛,武 潔,等.翼型模型幾何誤差對(duì)氣動(dòng)性能影響的自動(dòng)微分分析方法[J].空氣動(dòng)力學(xué)學(xué)報(bào),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.