超聲速半球形多孔探針大氣脫體激波仿真研究

張陽春，周樹道,2，姚 韜

(1.國防科技大學(xué) 氣象海洋學(xué)院，江蘇 南京 211101；2.南京信息工程大學(xué) 氣象災(zāi)害預(yù)警與評估協(xié)同創(chuàng)新中心，江蘇 南京 210044)

1 引 言

半球形多空探針是一種適用于大氣三維流場探測的傳感器，具有結(jié)構(gòu)簡單、靈敏度高等優(yōu)點，其外形呈球頭-圓柱體結(jié)構(gòu)。多孔探針的傳統(tǒng)標(biāo)定方法是基于探針孔壓的多項式函數(shù)擬合法，這種方法是在不考慮探針頭部繞流模型情況下建立起來一種風(fēng)速、風(fēng)向與孔壓間的多項式函數(shù)關(guān)系，其標(biāo)定精度主要取決于探針孔數(shù)和多項式階數(shù)，標(biāo)定過程也依賴于大量的風(fēng)洞實驗，成本較大[1～4]。目前針對多孔探針的研究大多是在亞聲速下，在超聲速和可壓縮條件下的研究較少。為改進(jìn)多孔探針的標(biāo)定方法，提高測量精度和范圍，降低標(biāo)定成本，需要對多超聲速條件下探針頭部流體繞流流場進(jìn)行研究，對超聲速大氣條件下探針頭部脫體激波的研究必不可少。

針對超聲速下鈍頭體脫體激波形狀和脫體距離的研究，國內(nèi)外學(xué)者以大量實驗數(shù)據(jù)為基礎(chǔ)得到了各種工程算法，并發(fā)現(xiàn)不同的工程算法只在自身實驗條件下適用[5～9]。本文采用CFD方法，針對半球形探針在超聲速大氣條件下產(chǎn)生的脫體激波形狀及位置開展研究，在1.2～1.7馬赫數(shù)條件下進(jìn)行數(shù)值仿真，使用最小二乘法進(jìn)行數(shù)據(jù)處理，發(fā)現(xiàn)了脫體激波形狀與脫體距離受馬赫數(shù)的影響規(guī)律，建立了半球形探針在超聲速下脫體激波曲線的參數(shù)化方程。

2 探針外形和網(wǎng)格

2.1 多孔探針結(jié)構(gòu)

半球形多孔探針尺寸較多，以王偲臣等提出的一種可用于葉輪機(jī)內(nèi)部流場測量的半球形高頻氣動五孔探針為例[10]，探針正面及剖面圖如圖1所示，其機(jī)構(gòu)由頭部的半球體和尾部的圓柱體構(gòu)成，內(nèi)部有5個平行于軸線的測壓孔，測壓孔直徑是探針直徑的1/10。在研究探針頭部脫體激波時，考慮到測壓孔的孔徑較探針直徑而言很小、孔內(nèi)部與外界間氣體通量為零，且脫體激波與探針表面有一定距離，因此可以忽略測壓孔的影響，將探針看作由半球和圓柱體組成的實心幾何體，在此基礎(chǔ)上對流場網(wǎng)格進(jìn)行劃分。

圖1 半球形多孔探針結(jié)構(gòu)圖Fig.1 Hexical porous probe structure diagram

2.2 網(wǎng)格處理

當(dāng)迎角為0°時流場在空間呈軸對稱分布,僅對四分之一流場仿真即可。圖2(a)所示為四分之一流體仿真域，邊界包括出口、入口、壁面、對稱面和探針表面，仿真域內(nèi)存在比較均勻的自由流和速度、壓強(qiáng)、密度、溫度等物理量梯度較大的激波，使用自適應(yīng)和誤差估計功能，基于內(nèi)置誤差估計自動調(diào)整網(wǎng)格尺寸的自適應(yīng)網(wǎng)格劃分方法可以很好的預(yù)測激波的位置并解析激波[11]。在常規(guī)尺寸網(wǎng)格基礎(chǔ)上對于激波位置處建立細(xì)化的網(wǎng)格，得到的自適應(yīng)網(wǎng)格，計算網(wǎng)格數(shù)約50萬個，網(wǎng)格頂點約9萬個。圖2(b)所示為網(wǎng)格劃分后的探針表面和對稱面，圖中網(wǎng)格的疏密程度變化反映了流場中不同位置處物理量的梯度大小。

圖2 流體仿真域和自適應(yīng)網(wǎng)格Fig.2 Fluid simulation domain and adaptive grid

3 數(shù)值方法及驗證

3.1 控制方程組

在馬赫數(shù)大于0.3的情況下，須考慮氣體的可壓縮性，研究脫體激波現(xiàn)象可認(rèn)為氣體是無粘可壓縮的[12]。N-S方程(Navier-Stokes方程)是描述粘性流體動量守恒的運動方程，忽略方程中的粘性項后可用來描述無粘流體。流體的可壓縮性可以用連續(xù)性方程描述。激波的產(chǎn)生會使氣體溫度發(fā)生劇烈變化，建立控制方程組時必須考慮氣體動能與內(nèi)能之間的轉(zhuǎn)化。鑒于此，得到適用于氣體激波無粘可壓縮流體控制方程組如下：

(1)

(2)

(3)

p=ρRT

(4)

e=γT

(5)

式(1)為N-S方程忽略粘性項時的無粘流體動量方程，式(2)為可壓縮流體的連續(xù)性方程，式(3)為流體能量方程，式(4)為氣體狀態(tài)方程，式(5)為氣體內(nèi)能公式。其中：ρ是流體密度；t是時間；u是流體速度矢量；e是氣體內(nèi)能；p是氣體壓強(qiáng)；R是氣體常數(shù)；γ是氣體比熱比；T是氣體溫度。

同時假設(shè)流體在壁面和探針表面滿足絕熱和滑移(無摩擦)條件：

nu=0

(6)

nq=0

(7)

式中：n是壁面和探針表面法相單位矢量；q是熱量。

3.2 數(shù)值方法和算例驗證

基于目前的全球大氣觀測數(shù)據(jù)，對流層頂氣壓場介于100～320 hPa之間[13],以300 hPa、273 K大氣環(huán)境值為數(shù)值模擬的計算工況，采用有限元方法對式(1)～式(5)所示的控制方程組進(jìn)行離散,為提高計算效率采用三次迭代求解，邊界條件取遠(yuǎn)場邊界無反射邊界條件，探針表面和壁面為無摩擦滑移、絕熱，全場取自由來流值為初始條件。自由來流速度在1.2到1.7馬赫之間，每間隔0.1馬赫做一組數(shù)值采樣。

圖3是來流馬赫數(shù)為1.5、攻角為0°時探針表面和對稱面的速度、壓強(qiáng)云圖，在探針頭部流體的速度和壓強(qiáng)出現(xiàn)明顯的弧形間斷面，脫體激波位置和形狀清晰可見，探針尾部形成尾流和斜激波，仿真流場與實驗現(xiàn)象相符。

圖3 速度和壓強(qiáng)云圖Fig.3 Velocity and pressure cloud

3.3 數(shù)據(jù)提取與處理

為得到激波曲線的具體位置，需要對仿真數(shù)據(jù)進(jìn)行篩選。流體經(jīng)過激波時壓強(qiáng)突增，壓強(qiáng)梯度云圖可以很好地識別激波，在數(shù)值仿真得到的流場壓強(qiáng)云圖基礎(chǔ)上對壓強(qiáng)求梯度得到如圖4(a)所示的壓強(qiáng)梯度云圖。從圖中可以看出脫體激波在球頭處呈弧線分布，沿直線向無窮遠(yuǎn)處延伸，與球面間的距離為脫體距離。

將網(wǎng)格節(jié)點數(shù)據(jù)導(dǎo)出并重新繪制得到如圖4(b)所示的壓強(qiáng)梯度節(jié)點離散圖，各節(jié)點的值代表節(jié)點處壓強(qiáng)梯度。按照一定的閾值將壓強(qiáng)梯度值較小的節(jié)點，即與激波曲線相關(guān)性較弱的節(jié)點剔除，得到如圖4(c)所示的與激波位置相關(guān)性較高的節(jié)點離散圖，作為激波曲線采樣點。

4 激波曲線擬合

在超聲速且迎角為零時，球-圓柱體形成的激波形狀為一軸對稱曲面，子午面如圖5所示，研究這一曲面只需考慮子午面的母線即可。來流馬赫數(shù)Ma=1時球-圓柱體頭部脫體激波是一無限遠(yuǎn)處的直線，隨著Ma的增大，激波的脫體距離D和曲率半徑R逐漸縮小，在球面處形成一條間距極窄的圓弧段，后沿以直線向無窮遠(yuǎn)處延伸，這一特征與雙曲線類似，激波脫體距離和曲率半徑可看作馬赫數(shù)的單值函數(shù)。

圖4 數(shù)據(jù)處理過程 Fig.4 Data processing

圖5 脫體激波直角坐標(biāo)關(guān)系Fig.5 Schematic diagram of detached shock wave

鄭之初等使用了曲率中心位于球心的一組圓錐曲線來表示激波形狀[14]。考慮到曲率中心與球心并不必然重合，將曲率半徑和脫體距離分別看作Ma的單值函數(shù)，分別擬合激波的位置和形狀，使用如圖(5)所示的以球心為原點的雙曲線直角坐標(biāo)參數(shù)化方程表示脫體激波，其表達(dá)式為：

(8)

(9)

(10)

D=|a+d|-1

(11)

式(8)為激波曲線的直角坐標(biāo)參數(shù)化方程。其中:x、y表示激波曲線上的點到坐標(biāo)軸的距離;r為球體半徑;x/r、y/r表示將激波曲線以r為單位長度等比例變化;a、b、d為待求解的未知量。

式(9)、式(10)、式(11)分別表示激波曲線的離心率e、曲率半徑R和脫體距離D與未知量a、b、d之間的關(guān)系。將e、R、D看作馬赫數(shù)Ma的單值函數(shù)，基于仿真數(shù)據(jù)擬合可得到e、R、D與Ma間的函數(shù)關(guān)系。對任意給定的Ma，聯(lián)立式(9)、式(10)、式(11)可得到式(8)中未知量a、b、d的值，代入式(8)中，任意Ma對應(yīng)的激波曲線都得以確定。

4.1 參數(shù)擬合

在1.2到1.7馬赫范圍內(nèi)進(jìn)行6組數(shù)值仿真，以不同馬赫數(shù)時激波曲線的采樣點為樣本數(shù)據(jù)，如圖4(c)所示，式(8)為回歸模型，采用最小二乘法擬合不同馬赫數(shù)時的激波曲線，得到公式中參數(shù)的估計值，擬合結(jié)果如圖6所示，樣本點大于200個，參數(shù)a、b、d的擬合方差如表1給出。

表1 參數(shù)a,b,d擬合方差Tab.1 Fitting variance of parameters a,b and d

由圖6擬合曲線中的參數(shù)a、b、d的值及式(9)、式(10)、式(11)計算得到不同馬赫數(shù)時參數(shù)e、R、D的值。

圖7給出了離心率e、曲率半徑R、脫體距離D在不同馬赫數(shù)下的值和各參數(shù)與馬赫數(shù)Ma的函數(shù)關(guān)系擬合曲線。

圖6 激波曲線擬合Fig.6 Shock curve fitting

擬合結(jié)果顯示隨著馬赫數(shù)的增大激波離心率、曲率半徑和脫體距離呈下降趨勢，且下降速度逐漸減緩，這與實驗現(xiàn)象相符。參數(shù)e、R、D的具體擬合公式為：

(12)

(13)

(14)

圖7 參數(shù)與馬赫數(shù)擬合曲線Fig.7 Fitting curve of parameters and Mach number

4.2 結(jié)果比較與分析

對不同的馬赫數(shù)Ma和球體半徑r，聯(lián)立式(9)、式(10)、式(11)、式(12)、式(13)、式(14)可得參數(shù)a、b、d的值，代入式(8)中得到由參數(shù)化方程預(yù)測的激波曲線。

圖8是本文得到的激波參數(shù)化方程曲線與仿真得到的對稱面壓強(qiáng)梯度云圖的對比，在1.2～1.7馬赫內(nèi),由參數(shù)化方程給出的曲線與仿真結(jié)果極為相符。

圖8 對稱面壓強(qiáng)梯度云圖與參數(shù)化方程曲線對比Fig.8 Symmetrical face pressure gradient cloud chart is compared with the paramethy equation curve

Billig基于試驗數(shù)據(jù)提出一種適用于球頭-圓柱脫體激波形狀的經(jīng)驗公式[15]，該公式適用于較低溫度下的理想氣體狀態(tài)。

圖9將經(jīng)驗公式與本文得到的參數(shù)化方程進(jìn)行了比較，從圖9中可以看出，當(dāng)流速介于1.2～1.4馬赫之間時2條曲線重合度較低，在1.5～1.7馬赫之間時兩曲線重合度較高。在同一馬赫數(shù)下，距離球面越近曲線重合度越高，隨著曲線遠(yuǎn)離球面重合度降低。兩條曲線對激波脫體距離的預(yù)測基本一致，預(yù)測結(jié)果的差值小于球體半徑的5%。

圖9 脫體激波隨馬赫數(shù)變化Fig.9 Changes of oobody shock wave with Mach number

距離頭部較近的曲線重合度高的原因，是在處理數(shù)據(jù)時，按照壓強(qiáng)梯度值來選取節(jié)點的，由于距離較遠(yuǎn)處的激波強(qiáng)度小，尤其當(dāng)馬赫數(shù)較低時激波強(qiáng)度衰減更快。節(jié)點處壓強(qiáng)梯度達(dá)不到閾值而被剔除，最終的擬合曲線更加接近球頭-圓柱體頭部強(qiáng)度較大的激波的形狀。說明該參數(shù)曲線對于近球面激波情況更加適用。

5 結(jié) 論

本文采用數(shù)值計算仿真了半球形多孔探針在超聲速大氣條件下的脫體激波，并研究了脫體激波形狀和脫體距離受大氣流體馬赫數(shù)影響規(guī)律。得到了在1.2～1.7馬赫范圍內(nèi)與數(shù)值仿真結(jié)果一致性較好的激波曲線參數(shù)化方程，方程對激波脫體距離預(yù)測準(zhǔn)確，在近球面處精度更高。