龔安龍, 劉曉文, 劉 周, 周偉江, 紀(jì)楚群

(中國航天空氣動力技術(shù)研究院, 北京 100074)

高超聲速壁面黏性力快速計算方法

龔安龍*, 劉曉文, 劉 周, 周偉江, 紀(jì)楚群

(中國航天空氣動力技術(shù)研究院, 北京 100074)

當(dāng)前高超聲速飛行器壁面黏性力的預(yù)測方法主要有精確的Novier-Stokes方程CFD方法(計算流體力學(xué))和快速工程計算方法。前者計算精度高，但效率很低；后者效率極高，但精度難以滿足工程應(yīng)用要求。如何將兩種方法有效結(jié)合，發(fā)展具有一定精度且效率也較高的壁面黏性力計算方法，一直是計算空氣動力學(xué)研究的重要方向之一。為此，本文基于無黏Euler方程CFD方法獲得的壁面流場參數(shù)，發(fā)展了一種更加精確的高超聲速飛行器壁面黏性力快速工程計算方法。該方法以經(jīng)典參考溫度法為基礎(chǔ)，以壁面當(dāng)?shù)亓鲌鰠?shù)替代自由來流參數(shù)，更加充分地考慮了高超聲速復(fù)雜流動的特性，從理論上提高了壁面黏性力的預(yù)測精度。為驗(yàn)證該方法的準(zhǔn)確性，以典型高超聲速高升阻比飛行器外形為研究對象，分別采用精確Navier-Stokes方程CFD方法和本文發(fā)展的快速工程計算方法，預(yù)測了高空高馬赫數(shù)飛行環(huán)境下飛行器壁面的黏性力；通過對比分析表明，本文所建立的快速計算方法預(yù)測偏差為10～20%，能夠滿足工程初步設(shè)計的要求。

高超聲速；壁面黏性力；工程計算方法；參考溫度法；CFD方法

0 引 言

隨著計算機(jī)運(yùn)算速度的大幅提高和計算流體力學(xué)(Computational Fluid Dynamics，CFD)方法的日趨完善[1]，快速求解Navier-Stokes方程已成為可能但快速的工程方法在飛行器設(shè)計中仍有很大的應(yīng)用前景。無黏Euler方程CFD方法與壁面黏性力快速工程計算方法相結(jié)合能夠快速預(yù)測流場特性，并保證良好的計算精度[2]。因此目前工程應(yīng)用中，特別是針對復(fù)雜飛行器外形繞流情況，常常采用該方法快速獲得滿足工程設(shè)計精度要求的氣動性能[3-4]。

壁面黏性力工程計算方法通?；谶吔鐚永碚摚宰杂蓙砹鲄?shù)作為邊界層外緣參數(shù)進(jìn)行計算，并需要針對飛行器幾何形狀進(jìn)行簡化和近似。比如已經(jīng)獲得廣泛應(yīng)用的等效平板層流Blassius近似方法和湍流Van Driest II方法[5]，工程應(yīng)用于一定馬赫數(shù)和雷諾數(shù)范圍內(nèi)均具有較好的預(yù)測精度[6]。而對于高超聲速高溫氣體作用下的壁面黏性力預(yù)測，參考溫度法則提供了一種基本滿足工程要求的預(yù)測能力[7-8]。然而，對于復(fù)雜的飛行器外形，特別是在高超聲速情況下，自由來流受到飛行器各部件的干擾非常嚴(yán)重，使得壁面邊界層外緣的流場參數(shù)與自由來流參數(shù)相差很大，因而上述壁面黏性力工程計算方法從理論上就會產(chǎn)生較大的偏差。而事實(shí)上，基于無黏Euler方程的CFD方法在進(jìn)行無黏流場計算時，就已經(jīng)比較準(zhǔn)確的獲得了邊界層外緣的流場參數(shù)。鑒于此，本文從無黏Euler方程CFD方法獲得的壁面流場參數(shù)出發(fā)，建立了基于當(dāng)?shù)亓骶€長度和流場參數(shù)的參考溫度法，應(yīng)用于高超聲速飛行器壁面黏性力的預(yù)測。將本文計算方法獲取的壁面黏性力計算結(jié)果與完全Navier-Stokes方程CFD方法計算結(jié)果進(jìn)行了比較，分析表明該方法相比于基于自由來流參數(shù)的經(jīng)典參考溫度計算方法具有更高的精度，能夠更好的滿足工程應(yīng)用需求。

1 基于無黏流場的壁面流線計算方法

無黏Euler方程計算中，壁面采用無穿透邊界條件，即速度方向平行于壁面。通過壁面的三個速度分量對時間的積分得到壁面流線上的坐標(biāo)值。假設(shè)P=(x,y,z)T為流線上的點(diǎn)，該點(diǎn)上的速度為：

根據(jù)流線[9]的定義：

初值條件為：x(t0)=x0,y(t0)=y0,z(t0)=z0

采用4階顯示Runge-Kutta方法求解該初值問題，具體公式如下：

h為步長，通?？梢匀∫粋€較小的常數(shù)，反映時間推進(jìn)的間隔。

由于壁面流線的特殊性，流線計算推進(jìn)時需要解決一些特殊問題[10]，主要包括：時間步長的合理選取；流線推進(jìn)點(diǎn)在壁面的投影；壁面速度場的插值。下面分別進(jìn)行討論。

1.1 變時間步長處理方法

時間步長h的選取決定了流線計算的精度和效率。如果時間步長取值過大，必然造成流線計算的偏差很大；如果時間步長取值過小，流線計算精度雖然很高，但會極大的增加計算時間。通常，Runge-Kutta方法中時間步長取某一合理的常數(shù)，但其合理性也只是從整個流場大概來說的。另一方面，無黏流場計算獲得的壁面流動參數(shù)與來流狀態(tài)和氣動外形密切相關(guān)，因此合理的時間步長還需要根據(jù)具體的流動特性來確定。鑒于此，本文提出了流場自適應(yīng)的變時間步長處理方法。

流場自適應(yīng)的變時間步長處理方法基本思路是：從流線點(diǎn)Pn推進(jìn)到下一個點(diǎn)Pn+1時，4階Runge-Kutta方法中的時間步長h取點(diǎn)Pn所在當(dāng)?shù)貑卧鲃犹卣鲿r間Tc。該特征時間定義為單元四條邊中以各自平均速度通過的最小時間(見圖1)，即:

eij為點(diǎn)i到點(diǎn)j的單位矢量。

該時間基本反映了流動穿越該當(dāng)?shù)貑卧奶卣鲿r間，因此能夠保證每一個壁面單元都至少存在一個流線點(diǎn)，從而保證了流線精度；同時也保證同一個單元內(nèi)不會存在過多的流線點(diǎn)使得流線計算時間太長。流場自適應(yīng)的變時間步長處理方法有效解決了流線計算精度和計算效率的矛盾。

1.2 流線推進(jìn)點(diǎn)在壁面的投影

1.3 壁面流場參數(shù)的插值

流線推進(jìn)獲得壁面投影點(diǎn)的坐標(biāo)后，其該點(diǎn)的流場參數(shù)是未知的。而已知的僅僅是單元節(jié)點(diǎn)上的值。因此需要采用插值方法獲得單元內(nèi)投影點(diǎn)的流場參數(shù)。本文采用了反距離權(quán)重的插值方法。如圖4所示，三角形單元由節(jié)點(diǎn)1、2、3構(gòu)成，其上的流場參數(shù)V1、V2、V3為已知量。Pn為單元內(nèi)的投影點(diǎn)，d1、d2、d3分別為Pn到三個節(jié)點(diǎn)的距離。Pn點(diǎn)流場參數(shù)Vn采用反距離權(quán)重的插值方法表示為：

其中權(quán)重因子表示為：

V為任意流場參數(shù)，包括速度、壓力和密度等，均采用該方法進(jìn)行插值。

1.4 流線計算的驗(yàn)證算例

采用流場后處理商業(yè)軟件Tecplot對本文的流線計算方法進(jìn)行了考核驗(yàn)證。圖5給出了典型升力體外形壁面流場所確定的流線，圖中紫色圓點(diǎn)是本文方法計算得到的壁面流線，藍(lán)線為Tecplot軟件繪出的壁面流線。圖6為典型高升阻比外形在高超聲速流動情況下壁面的流線，同樣圖中用紫色圓點(diǎn)代表本文方法計算得到的壁面流線，藍(lán)線代表Tecplot軟件繪出的壁面流線。通過比較發(fā)現(xiàn)：本文計算方法獲得的流線與Tecplot軟件繪制結(jié)果完全一致。

2 基于當(dāng)?shù)亓鲌鰠?shù)的參考溫度法

2.1 參考溫度法

參考溫度法被認(rèn)為是高超聲速高溫氣體流動下黏性氣動力預(yù)測的一種較成熟的近似工程算法[7]。在高超聲速飛行器選型設(shè)計中，該方法得到了廣泛應(yīng)用[11-13]。該方法基于不可壓平板邊界層流動理論，考慮了可壓縮修正和溫度影響。

首先討論層流黏性摩擦力系數(shù)的參考溫度計算方法。對于不可壓平板層流，Blasius解得到當(dāng)?shù)啬Σ亮ο禂?shù)為：

參考溫度法引入可壓縮修正，可以將當(dāng)?shù)啬Σ亮ο禂?shù)寫成：

假設(shè)邊界層內(nèi)氣體參數(shù)仍然滿足完全氣體狀態(tài)方程，參考溫度T*下的密度ρ*可表示為：

根據(jù)邊界層內(nèi)壓力的分布特點(diǎn)，通常取壓力p*≈pe，于是

黏性系數(shù)μ*由溫度T*唯一確定，采用粘溫指數(shù)函數(shù)關(guān)系式進(jìn)行計算[11]，即：

通常取ω=0.75。

對于湍流情況，也采用與層流參考溫度法相同的處理。根據(jù)平板不可壓湍流常用的摩擦力系數(shù)關(guān)系式：

考慮可壓縮修正，可以得到：

對于層、湍流同時存在的情況，引入轉(zhuǎn)捩雷諾數(shù)來判定流動是層流還是湍流，并分段計算黏性力。轉(zhuǎn)捩雷諾數(shù)可用如下的近似公式計算[11]，

其中Rext為轉(zhuǎn)捩雷諾數(shù)。當(dāng)邊界層外緣雷諾數(shù)Reex>Rext時認(rèn)為流動為湍流，摩擦力系數(shù)采用公式(15)計算；而當(dāng)Reex≤Rext時，認(rèn)為流動為層流，摩擦力系數(shù)采用公式(9)進(jìn)行計算。

2.2 黏性力計算

通常平板流動邊界層外緣的物理量近似認(rèn)為等于自由來流的物理量。但是對于復(fù)雜外形的高超聲速流動，由于可壓縮氣流的復(fù)雜作用(激波壓縮和氣流膨脹等)，邊界層外緣流場物理量與自由來流的值差別很大，因此如果用自由來流的參數(shù)作為邊界層外緣的流動參數(shù)將使當(dāng)?shù)啬Σ亮ο禂?shù)的計算結(jié)果產(chǎn)生很大的偏差。采用無黏Euler方程CFD方法可以快速得到無黏流場的特性，將其壁面流動參數(shù)作為對應(yīng)當(dāng)?shù)剡吔鐚油饩壍牧鲃訁?shù)則能夠比較準(zhǔn)確的反映流動真實(shí)特性。

全體表面的黏性作用力需要通過對所有壁面單元的摩擦力進(jìn)行求和獲得，即:

其中：Fv為總的黏性作用力；

fi為壁面單元i的黏性摩擦力；

cfi為壁面單元i的黏性摩擦力系數(shù)；

QAi為壁面單元i的當(dāng)?shù)貏訅海?/p>

Ai為壁面單元i的面積。

單元i的黏性摩擦力系數(shù)cfi采用基于當(dāng)?shù)亓鲌鰠?shù)的參考溫度法進(jìn)行計算，見公式(9)和(15)，即:

采用當(dāng)?shù)匚恢昧骶€長度計算雷諾數(shù)，即:

式中s為從頭部附近發(fā)出的到達(dá)單元i的流線總長度(圖7)，由第1節(jié)中的流線計算方法獲得。

單元i的當(dāng)?shù)貏訅篞Ai為參考溫度T*下的動壓：

于是單元i的黏性摩擦力可以表示為：

該摩擦力fi的方向是流線的切向，也就是當(dāng)?shù)乇诿嫠俣确较?。令切向單位矢量為τi=(τx,τy,τz)，見圖7。于是單元i的摩擦力fi在直角坐標(biāo)系(x,y,z)中的三個分量為：

fix=fiτx

fiy=fiτy

于是飛行器表面總的黏性作用力在直角坐標(biāo)系的三個分量可以表示為：

無量綱化的黏性作用力系數(shù)表示為：

其中ρ∞、u∞分別為自由來流的密度和速度，AR為參考面積。

3 典型算例驗(yàn)證

為考察上述黏性力快速計算方法的準(zhǔn)確性，針對高升阻比外形(圖6)的高空高超聲速流動情況進(jìn)行了研究。其中以Navier-Stokes方程CFD數(shù)值計算方法[14-15]獲得的完全氣體模型下的黏性力作為準(zhǔn)確值，并采用基于自由來流條件的經(jīng)典平板邊界層參考溫度法和本文快速計算方法分別進(jìn)行了黏性力的計算。本文的研究選用了馬赫數(shù)10和20兩個速度情況，高度也選取了50 km和60 km兩組情況，迎角最大值達(dá)到20°。通過對Navier-Stokes方程CFD數(shù)值計算結(jié)果的分析表明，各計算狀態(tài)下飛行器表面未發(fā)生明顯流動分離現(xiàn)象。

圖8給出了三種方法所獲得的黏性軸向力分量結(jié)果的比較?？梢钥吹?，在馬赫數(shù)10情況下，經(jīng)典參考溫度法和本文計算方法所獲得的結(jié)果相差不大；馬赫數(shù)20情況下，經(jīng)典參考溫度法計算結(jié)果與準(zhǔn)確值存在較大偏差，而本文計算方法獲得的結(jié)果與準(zhǔn)確值更接近；在50 km和60 km兩種不同高度情況下，上述規(guī)律完全一致。

同時我們發(fā)現(xiàn)在所研究的馬赫數(shù)和迎角范圍內(nèi)，馬赫數(shù)越高，迎角越大，經(jīng)典參考溫度法計算結(jié)果偏差越大，而本文計算方法與準(zhǔn)確值更接近。這是因?yàn)轳R赫數(shù)越高，迎角越大，自由來流受到飛行器影響就越大，邊界層外緣實(shí)際流動參數(shù)與自由來流之間的差異則更大；采用本文方法以Euler方程CFD方法獲得的準(zhǔn)確流動參數(shù)作為邊界層外緣參數(shù)，相比于經(jīng)典參考溫度法以自由來流作為邊界層外緣參數(shù)，能夠更精確的反映控制黏性底層流動的外圍流動特性，因此本文方法計算結(jié)果相比于經(jīng)典參考溫度法會更加準(zhǔn)確。

總之，相比于經(jīng)典參考溫度法，本文計算方法的結(jié)果與Navier-Stokes方程CFD數(shù)值計算結(jié)果的一致性更好，存在約10～20%的偏差，能夠滿足工程初步設(shè)計的要求。

另一方面，從計算效率來看，對于單個計算狀態(tài)，采用Navier-Stokes方程CFD方法所需要的計算時間高出Euler方程CFD方法計算時間一個量級，而本文基于Euler方程CFD方法計算獲得的無黏流場所建立的黏性力快速計算方法在個人臺式計算機(jī)單個CPU上運(yùn)行所需的計算時間通常只有幾秒或幾十秒，相比于CFD方法的計算時間可以忽略不計。因此本文計算方法與Navier-Stokes方程CFD方法的效率之比，相當(dāng)于Euler方程CFD方法與Navier-Stokes方程CFD方法的效率之比，即達(dá)到10倍的量級?？梢姳疚姆椒ㄔ谟嬎阈史矫婢哂忻黠@的優(yōu)勢。

4 結(jié) 論

本文從工程應(yīng)用的角度出發(fā)，基于無黏Euler方程CFD方法獲得的壁面流場參數(shù)，建立了一種高超聲速飛行器壁面黏性力的快速工程計算方法。所提出的流場自適應(yīng)變時間步長流線計算方法大大提高了流線計算效率；建立的基于當(dāng)?shù)亓鲌鰠?shù)的參考溫度法，對壁面黏性力的計算相比于經(jīng)典參考溫度法具有更高的精確性，并得到了Navier-Stokes方程CFD方法的驗(yàn)證，能夠滿足工程設(shè)計的要求。

由于壁面流線建立的前提是流動在壁面沒有發(fā)生分離，因此本文的計算方法在理論上僅僅適用于沒有發(fā)生流動分離的情況。但飛行器高超聲速飛行環(huán)境下，在一定迎角范圍內(nèi)通常不會發(fā)生比較明顯的分離流動；而且在工程設(shè)計中也需要盡量避免大迎角高超聲速飛行情況。因而，對于高超聲速飛行器工程設(shè)計來講，本文建立的高超聲速壁面黏性力快速計算方法仍然具有良好的準(zhǔn)確性和實(shí)用性。

[1]Yan C, Yu J, et al. On the achievements and prospects for the methods of computational fluid dynamics[J]. Advances In Mechanics, 2011, 41(5): 562-589. (in Chinese)閻超, 于劍, 等. CFD 模擬方法的發(fā)展成就與展望[J]. 力學(xué)進(jìn)展, 2011, 41(5): 562-589.

[2]Zhu Z Q, et al. Application computational fluid dynamics[M]. Beijing: Press of Beijing University of Aeronautics and Astronautics, 1998. (in Chinese)朱自強(qiáng)等編著, 應(yīng)用計算流體力學(xué)[M]. 北京: 北京航空航天大學(xué)出版社, 1998.

[3]Wang R H, Yin G L, et al. Fast CFD tools for civil aircraft conceptual design and optimization use[J]. Aircraft Design, 2012, 32(5): 31-35. (in Chinese)王如華, 尹貴魯, 等. 快速CFD計算工具在民機(jī)概念優(yōu)化設(shè)計中的應(yīng)用[J]. 飛機(jī)設(shè)計, 2012, 32(5): 31-35.

[4]Xu R F, Deng Y J, et al. Aerodynamic optimization design and its requirement to CFD[J]. Aeronautical Science & Technology, 2011(2): 50-52. (in Chinese)許瑞飛, 鄧一菊, 等. 氣動優(yōu)化設(shè)計及其對CFD的需求[J]. 航空科學(xué)技術(shù), 2011(2): 50-52.

[5]Van Driest E R. On turbulent flow near a wall[J]. Journal of the Aeronautical Sciences, 1956, 23(11): 1007-1011.

[6]Fleeman E L. Tactical missile design, second edition[M]. Virginia: AIAA Inc., 2006.

[7]Anderson J D. Hypersonic and high temperature gas dynamics, second edition[M]. Virginia: AIAA Inc., 2006.

[8]Gnoffo P A. Computational fluid dynamics technology for hypersonic applications[R]. NASA 20040013407, 2004.

[9]Legendre R. Streamline in a steady flow[R]. NASA ADA395523, 1966.

[10]Diachin D P, Herzog J. Analytic streamline calculations on linear tetrahedral[R]. AIAA-97-1975:733-742, 1997.

[11]Corda C, Anderson J D. Viscous optimized hypersonic waveriders designed from axisymmetric flow fields[R]. AIAA-88-0369, 1988.

[12]Frederick F, Terry C, et al. The development of Waveriders from axisymmetric flowfield[R]. AIAA 2007-847.

[13]Zhang D, Tang S, et al. Engineering calculation method of viscous force for air-breathing hypersonic vehicle[J]. Journal of Solid Rocket Technology, 2013, 36(3): 291-295. (in Chinese)張棟, 唐碩, 等. 吸氣式高超聲速飛行器黏性力工程計算方法[J]. 固體火箭技術(shù), 2013, 36(3): 291-295.

[14]Gong A L, Liu Z, et al. Numerical study on hypersonic double-cone separated flow[J]. Acta Aerodynamica Sinica, 2014, 32(6): 767-771. (in Chinese)龔安龍, 劉周, 等. 高超聲速激波/邊界層干擾流動數(shù)值模擬研究[J]. 空氣動力學(xué)學(xué)報, 2014, 32(6): 767-771.

[15]Liu Z. Delayed detached eddy simulation for static and forced oscillating airfoil at high angle of attack[R]. IAC 2013, 2013.

A rapid method for hypersonic skin viscous force calculation

Gong Anlong*, Liu Xiaowen, Liu Zhou, Zhou Weijiang, Ji Chuqun

(ChinaAcademyofAerospaceAerodynamics,Beijing100074,China)

Accurate computational fluid dynamic (CFD) methods with Novier-Stokes equations and rapid engineering methods are currently two major approaches to obtain skin viscous force for hypersonic aircraft. The former is relatively accurate with poor efficiency, and the later is efficient with poor accuracy. How to combine the two methods and develop a new way with enough efficiency and accuracy is always an important aspect for computational aerodynamics. Therefor , a rapid engineering method to predict skin viscous force of hypersonic aircraft more accurately was developed with near-wall flow characteristics derived from invicid Euler equations. The method was based on the classical reference temperature method and constructed with local flow parameters instead of free flow parameters, so that could count the hypersonic flow characteristics more adequately and achieve a more accurate prediction of skin viscous forces in theory. In order to examine the accuracy of the rapid method, the skin viscous forces of a typical hypersonic aircraft with high Lift-to-Drag ratio at hypersonic velocity and high altitude was calculated by a Navier-Stokes equations CFD method and the presented rapid method. The comparison analysis finally showed a relative deviation of 10～20% from the rapid method, which could gave a good satisfaction for engineering applications.

hypersonic; skin viscous force; engineering method; reference temperature method; CFD methods

0258-1825(2017)01-0033-06

2015-04-02;

2015-05-24

國家自然科學(xué)基金(11372040)；國家自然科學(xué)基金(11472258)

龔安龍*(1977-)，男，湖北人，高級工程師，研究方向：氣動外形設(shè)計與CFD應(yīng)用.E-mail：gonganlong@tsinghua.org.cn

龔安龍, 劉曉文, 劉周, 等. 高超聲速壁面黏性力快速計算方法[J]. 空氣動力學(xué)學(xué)報, 2017, 35(1): 33-38.

10.7638/kqdlxxb-2015.0040 Gong A L, Liu X W, Liu Z, et al. A rapid method for hypersonic skin viscous force calculation[J]. Acta Aerodynamica Sinica, 2017, 35(1): 33-38.