趙曉龍，張君安，董 皓，劉 波

(西安工業(yè)大學(xué) 機(jī)電工程學(xué)院，西安 710021)

近年來(lái)，隨著航空航天事業(yè)的發(fā)展，超聲速飛行器的戰(zhàn)略價(jià)值受到了各科技大國(guó)的高度重視，成為各國(guó)近期的研究熱點(diǎn)[1].超聲速與低速最大的不同之一便是激波.在超聲速流場(chǎng)中，氣流的黏性效應(yīng)強(qiáng)，激波與物面附近的邊界層流動(dòng)會(huì)發(fā)生相互干擾，形成近壁區(qū)的復(fù)雜流場(chǎng)，導(dǎo)致干擾區(qū)內(nèi)激波的位置、形狀、激波前后的流動(dòng)參量以及邊界層特性均發(fā)生變化，甚至有可能支配整個(gè)流場(chǎng)[2-5].由于超聲速的地面實(shí)驗(yàn)存在難度大，成本高，周期長(zhǎng)和飛行條件不易控制等困難，因此，工程中常常依靠數(shù)值方法來(lái)彌補(bǔ)這些不足.

研究人員針對(duì)超聲速計(jì)算問(wèn)題做了一些研究.文獻(xiàn)[6]研究了超聲速大攻角旋轉(zhuǎn)飛行器氣動(dòng)特性的數(shù)值模擬方法，以Woodward有限基本解方法為基礎(chǔ)，計(jì)算了旋轉(zhuǎn)對(duì)彈體和彈翼渦流生成、發(fā)展的影響.文獻(xiàn)[7]對(duì)跨聲速和超聲速流中激波/邊界層干擾進(jìn)行了數(shù)值模擬，針對(duì)湍流模型中的Bradshaw 常數(shù)進(jìn)行了修正，并對(duì)跨聲速和超聲速流中激波/邊界層干擾進(jìn)行了數(shù)值模擬研究，計(jì)算得到的壁面壓力分布、分離區(qū)長(zhǎng)度和速度剖面均與實(shí)驗(yàn)值吻合較好.文獻(xiàn)[8]使用雷諾平均方法研究了平板超聲速流中的激波和邊界層問(wèn)題.文獻(xiàn)[9]研究了超聲速二維底部壓力的簡(jiǎn)化解析計(jì)算方法，應(yīng)用質(zhì)量和能量守恒原理，推導(dǎo)出了一個(gè)超聲速二維底部壓力的簡(jiǎn)化解析算法.文獻(xiàn)[10]針對(duì)超聲速氣流繞局部透氣凹角流動(dòng)問(wèn)題進(jìn)行了數(shù)值模擬，用具有壓力分裂形式的簡(jiǎn)化納維-斯托克斯方程為控制方程，數(shù)值模擬了超聲速來(lái)流條件下的激波-邊界層干擾被動(dòng)控制，數(shù)值計(jì)算時(shí)采用了多重掃描法對(duì)控制方程差分離散.文獻(xiàn)[11]在迎風(fēng)型格式和矢通量分裂技術(shù)的基礎(chǔ)之上，對(duì)捕捉激波方法進(jìn)行一種新的嘗試，該方法抑制非物理波動(dòng)，將其轉(zhuǎn)換成守恒形式，得到了基本上無(wú)振蕩的激波捕捉格式.

綜上所述，目前對(duì)于超聲速問(wèn)題計(jì)算主要從數(shù)學(xué)模型、控制方程及差分格式上進(jìn)行研究.由于超聲速問(wèn)題本身的復(fù)雜性，本文采用簡(jiǎn)化的納維-斯托克斯方程作為控制方程難以真實(shí)地模擬出氣體流動(dòng)的實(shí)際狀態(tài).本文采用直接數(shù)值模擬方法，計(jì)算連續(xù)方程、含有黏度項(xiàng)的納維-斯托克斯方程和能量方程，不對(duì)應(yīng)力分量和黏度進(jìn)行簡(jiǎn)化，同時(shí)聯(lián)立完全氣體狀態(tài)等其他基本方程，使用高精度的差分格式進(jìn)行離散，保證計(jì)算數(shù)據(jù)的準(zhǔn)確性，同時(shí)研究分析了不同雷諾數(shù)對(duì)尖前緣平板超聲速繞流狀態(tài)的影響.

1 計(jì)算模型及控制方程

本文的物理模型如圖1所示，一個(gè)帶有尖角的二維平板放置在馬赫數(shù)Ma大于1的超聲速氣流中.其中板長(zhǎng)為L(zhǎng)，初始來(lái)流的壓強(qiáng)p、密度ρ、速度V和溫度T均為標(biāo)準(zhǔn)值，超聲速氣流流過(guò)平板上表面.

圖1 物理模型

建立包含了黏性效應(yīng)和應(yīng)力分量的連續(xù)方程、動(dòng)量方程和能量方程的向量形式為

(1)

其中

式中：U，E，F(xiàn)分別為所構(gòu)建向量形式的三個(gè)分量；t為時(shí)間；x,y為位置坐標(biāo)；ρ為密度；u為x方向的速度；v為y方向的速度；V為速度；Et為單位體積流體動(dòng)能和內(nèi)能之和；e為內(nèi)能；p為流場(chǎng)內(nèi)壓力；τxx,τxy,τyx,τyy為不同方向的應(yīng)力；qx,qy為熱流分量；k為熱傳動(dòng)率；λ為第二黏性系數(shù)；μ為動(dòng)力黏度；T為溫度.

完全氣體假定的狀態(tài)方程為

由此可見(jiàn)，“301條款”是一項(xiàng)美國(guó)制定的單方面發(fā)動(dòng)調(diào)查的國(guó)內(nèi)法規(guī)則，美國(guó)總統(tǒng)憑此規(guī)定可采取單方行動(dòng)使美國(guó)從他國(guó)對(duì)其不公平貿(mào)易行為中解脫出來(lái)。

p=ρRT

(2)

其中R為理想氣體常數(shù).

進(jìn)一步假定為定壓比熱的完全氣體，則有

e=cvT

(3)

(4)

式中：|V|為速度的模；u為x方向的速度；v為y方向的速度；cv為定容熱容.

假定氣體是定壓比熱的完全氣體，則動(dòng)力黏度為

(5)

其中μ0，T0分別為標(biāo)準(zhǔn)黏度和標(biāo)準(zhǔn)溫度.

(6)

式中：Pr為普朗特常數(shù)；cp為定壓比熱.

由式(1)～(6)聯(lián)立求解，采用2階精度的數(shù)值差分，可以求得U.其差分格式為

(7)

時(shí)間步長(zhǎng)的計(jì)算公式為

Δt=min[K(Δt0)i,j]

(8)

2 數(shù)值計(jì)算分析

為使計(jì)算方便，本文選取L=0.01 mm進(jìn)行直接計(jì)算，則該模型的計(jì)算坐標(biāo)如圖2所示，x為平板表面坐標(biāo)位置；y為邊界層位置.計(jì)算網(wǎng)格既能光滑無(wú)振蕩地捕捉激波，又具有相當(dāng)高精度以保證能精確描述邊界層內(nèi)不穩(wěn)定因素的微小變化.本文選用的計(jì)算參數(shù)見(jiàn)表1.

圖2 計(jì)算坐標(biāo)

計(jì)算區(qū)域進(jìn)行初始化，同時(shí)分別計(jì)算該區(qū)域內(nèi)的邊界條件、動(dòng)力黏度、熱傳導(dǎo)系數(shù)和應(yīng)力；使用時(shí)間推進(jìn)的計(jì)算方法，單獨(dú)求解每次推進(jìn)的時(shí)間步長(zhǎng)；設(shè)定收斂準(zhǔn)則.本文數(shù)值計(jì)算的收斂準(zhǔn)則為：在每一時(shí)間步后，若所有網(wǎng)格點(diǎn)上的密度較前一時(shí)間步長(zhǎng)上的密度值的變化不超過(guò)1×10-6kg·m-3，則認(rèn)為收斂，若不收斂，則返回第二步再次計(jì)算推進(jìn)時(shí)間步長(zhǎng).比較入口質(zhì)量流量和出口質(zhì)量流量?jī)烧咂?，解耦輸出所需參?shù).

表1 計(jì)算參數(shù)Tab.1 Calculation parameters

根據(jù)控制方程，計(jì)算馬赫數(shù)等于2，4，6，8時(shí)坐標(biāo)內(nèi)氣體速度、壓力、溫度和密度的分布狀態(tài)，可得氣體狀態(tài)如圖3所示.

由圖3可知，在計(jì)算坐標(biāo)區(qū)域內(nèi)，氣體的速度、壓力、溫度和密度產(chǎn)生了明顯的非線性變化，存在明顯的速度、壓力、溫度和密度分解線，這個(gè)分界線就是通過(guò)數(shù)值計(jì)算捕捉到的激波.在激波的分界線的上方氣體的各個(gè)狀態(tài)基本不發(fā)生變化，約等于流入平板時(shí)的氣體狀態(tài)；在激波分界線下方，隨著越靠近平板的位置，來(lái)流的氣體狀態(tài)變化越強(qiáng)烈，導(dǎo)致這種變化的原因就是因?yàn)闅怏w黏度的存在，氣流的速度越高，黏度對(duì)氣體狀態(tài)的影響越大.

圖3 不同馬赫數(shù)時(shí)計(jì)算坐標(biāo)內(nèi)的氣體狀態(tài)

圖4為馬赫數(shù)分別等于2，4，6和8時(shí)平板表面的壓力分布，由圖4可知，當(dāng)超聲速的氣流流過(guò)平板表面時(shí)，平板表面的壓力先激增后急劇下降到某一值，然后趨于穩(wěn)定；且馬赫數(shù)越大，沿平板表面的壓力分布越大.當(dāng)馬赫數(shù)等于8時(shí)，平板表面的最大壓力與標(biāo)準(zhǔn)壓力的比值為5.65倍；馬赫數(shù)等于2時(shí)，平板表面的最大壓力與標(biāo)準(zhǔn)壓力的比值為2.05倍.隨著馬赫數(shù)的增大，在沿平板表面位置上，壓力急劇下降后，又略微升高，然后緩慢下降.當(dāng)馬赫數(shù)等于8時(shí)，在平板表面位置等于0.13×10-5m時(shí)，平板表面的最大壓力與標(biāo)準(zhǔn)壓力比值為3.1倍；當(dāng)平板表面位置等于0.2×10-5m時(shí)，平板表面的最大壓力與標(biāo)準(zhǔn)壓力比值增大至3.3倍.

圖5為不同馬赫數(shù)時(shí)邊界層的流體狀態(tài).隨著邊界層坐標(biāo)數(shù)值增大，流場(chǎng)內(nèi)的氣體與標(biāo)準(zhǔn)壓力和標(biāo)準(zhǔn)溫度的比值先增大后減??；馬赫數(shù)持續(xù)增大，直至等于外界來(lái)流氣體的馬赫數(shù).當(dāng)馬赫數(shù)等于2時(shí)，計(jì)算區(qū)域的邊界層為0.12×10-5m，流場(chǎng)內(nèi)最大壓力與標(biāo)準(zhǔn)壓力比值為2倍，最高溫度與標(biāo)準(zhǔn)溫度比值為1.2倍，在邊界層坐標(biāo)等于0.8×10-5m時(shí)馬赫數(shù)達(dá)到最大；馬赫數(shù)等于4時(shí)，計(jì)算區(qū)域的邊界層為0.8×10-6m，最大壓力與標(biāo)準(zhǔn)壓力比值為2.9倍，最高溫度與標(biāo)準(zhǔn)溫度比值為1.5倍，在邊界層坐標(biāo)等于4.3×10-6m時(shí)馬赫數(shù)達(dá)到最大；馬赫數(shù)等于6時(shí)，計(jì)算區(qū)域的邊界層為6.83×10-6m，最大壓力與標(biāo)準(zhǔn)壓力比值為4.1倍，最高溫度與標(biāo)準(zhǔn)溫度比值為2.4倍，在邊界層坐標(biāo)等于3.3×10-6m時(shí)馬赫數(shù)達(dá)到最大；馬赫數(shù)等于8時(shí)，計(jì)算區(qū)域的邊界層為5.7×10-6m，最大壓力與標(biāo)準(zhǔn)壓力比值為5.6倍，最高溫度與標(biāo)準(zhǔn)溫度比值為3.3倍，在邊界層坐標(biāo)等于2.9×10-6m時(shí)馬赫數(shù)達(dá)到最大.

圖4 平板表面壓力分布

圖5 不同馬赫數(shù)時(shí)邊界層的流體狀態(tài)

3 結(jié) 論

1) 應(yīng)用直接數(shù)值模擬求解含有黏度項(xiàng)的納維-斯托克斯方程，反應(yīng)了所求模型的實(shí)際流體狀態(tài)，直接數(shù)值模擬所求模型的精度優(yōu)于其他方法.

2) 氣體以超聲速的速度流入平板上時(shí)，在平板上產(chǎn)生激波.隨著馬赫數(shù)的增大，產(chǎn)生激波的范圍越來(lái)越小，靠近平板位置的速度小于流入平板的速度，但壓力整體大于遠(yuǎn)離平板位置的壓力.

3) 流場(chǎng)內(nèi)的氣體的最大壓力和最大溫度出現(xiàn)在平板表面的邊界層中，隨著馬赫數(shù)增大，出現(xiàn)的位置越來(lái)越靠近平板表面，隨著邊界層坐標(biāo)數(shù)值增大，馬赫數(shù)持續(xù)增大，直至等于外界來(lái)流氣體的馬赫數(shù).