張紹廣，肖茂超，張宇飛，陳海昕

（清華大學(xué)航天航空學(xué)院，北京 100084）

0 引 言

導(dǎo)彈和飛行器需要在大攻角工況下飛行，因此必須研究大攻角下旋成體的流動(dòng)結(jié)構(gòu)，才能保證在大攻角工況對(duì)導(dǎo)彈和飛行器的精確控制。小攻角時(shí)（0°≤≤15°），旋成體生成附著對(duì)稱的渦結(jié)構(gòu)；在中等攻角時(shí)（15°＜≤30°），旋渦脫體但仍然維持對(duì)稱結(jié)構(gòu)，此時(shí)產(chǎn)生的側(cè)向力可以忽略不計(jì)；在大攻角時(shí)（30°＜≤65°），受對(duì)流不穩(wěn)定性影響，旋渦變得不對(duì)稱和非定常，側(cè)向力也變得顯著。在大攻角時(shí)（30°＜≤65°），側(cè)向力甚至超過(guò)法向力，從而將導(dǎo)致較大的不良偏航力矩。

側(cè)向力的產(chǎn)生與很多因素有關(guān)，如雷諾數(shù)、攻角、馬赫數(shù)和導(dǎo)彈前體形狀等。Lamont等研究了不同雷諾數(shù)對(duì)側(cè)向力的影響：基于導(dǎo)彈直徑的雷諾數(shù)Re為0.2×10和4.0×10時(shí)，分別為層流分離和全湍流的分離；Re在0.4×10至0.8×10的范圍時(shí)，分離剪切層中存在轉(zhuǎn)捩；在攻角小于70°時(shí)，層流分離和全湍流的分離非常相似，相比存在轉(zhuǎn)捩的分離產(chǎn)生更加顯著的側(cè)向力。Kumar等發(fā)現(xiàn)導(dǎo)彈側(cè)向力的產(chǎn)生主要是由模型的不完美導(dǎo)致的，而且在實(shí)際生產(chǎn)制造過(guò)程中很難制造出完美的模型。Keener等在研究中指出，導(dǎo)彈在不同滾轉(zhuǎn)角時(shí)會(huì)產(chǎn)生不同的側(cè)向力，側(cè)向力的不同主要取決于導(dǎo)彈前體20%的長(zhǎng)度位置，單獨(dú)旋轉(zhuǎn)前體和旋轉(zhuǎn)整個(gè)導(dǎo)彈產(chǎn)生的側(cè)向力基本一致。有學(xué)者在數(shù)值計(jì)算中發(fā)現(xiàn)，在導(dǎo)彈的頭部設(shè)置擾動(dòng)可以計(jì)算出不對(duì)稱的渦結(jié)構(gòu)。Lim等在導(dǎo)彈前體采用不同的層流-湍流轉(zhuǎn)捩條件，捕捉到了旋渦的不對(duì)稱結(jié)構(gòu)。

主動(dòng)和被動(dòng)的流動(dòng)控制方法（如：擾流片、鈍頭、圓形轉(zhuǎn)捩帶、頭部凹坑、吹吸氣等）被用來(lái)控制導(dǎo)彈大攻角時(shí)側(cè)向力的產(chǎn)生。對(duì)這些流動(dòng)控制方法的研究，大多采用試驗(yàn)測(cè)試的方法進(jìn)行，效率往往較低。但是，目前采用數(shù)值計(jì)算的方法發(fā)展減小側(cè)向力的流動(dòng)控制方法依然很難實(shí)現(xiàn)。主要的原因在于，數(shù)值計(jì)算中需要引入微小的擾動(dòng)才可以捕捉到不對(duì)稱的流動(dòng)結(jié)構(gòu)，擾動(dòng)的位置和尺寸都需要根據(jù)試驗(yàn)值來(lái)確定。

研究導(dǎo)彈大攻角非對(duì)稱分離形態(tài)和側(cè)向力的變化規(guī)律需要準(zhǔn)確的預(yù)測(cè)模型。迄今為止，雷諾平均方程（Reynolds-averaged Navier-Stokes，RANS）仍然是計(jì)算工程湍流問(wèn)題最常用的方法。RANS 方法采取系綜平均的思路，將湍流信息分解為平均量和脈動(dòng)量，脈動(dòng)量在Navier-Stokes 方程中通過(guò)雷諾平均轉(zhuǎn)化為雷諾應(yīng)力的形式，然后利用湍流模式對(duì)其進(jìn)行?；５荝ANS方法在計(jì)算分離流動(dòng)時(shí)會(huì)產(chǎn)生較大誤差，尤其不適用于預(yù)測(cè)大分離流動(dòng)。與RANS 方法對(duì)比，大渦模擬（large-eddy simulation，LES）方法采取空間過(guò)濾的思路，將流場(chǎng)信息分解成大尺度低頻脈動(dòng)和小尺度高頻脈動(dòng)。但是，采用LES方法解析邊界層的計(jì)算量很大。Spalart指出，使用壁面解析的LES（wall-resolved LES，WRLES）百萬(wàn)雷諾數(shù)的機(jī)翼邊界層流動(dòng)，需要網(wǎng)格數(shù)量為10量級(jí)，時(shí)間推進(jìn)步數(shù)為10量級(jí)，計(jì)算量遠(yuǎn)遠(yuǎn)超過(guò)了現(xiàn)代計(jì)算機(jī)的計(jì)算能力。為了克服LES在模擬高雷諾數(shù)湍流邊界層時(shí)的計(jì)算量過(guò)大的問(wèn)題，可以采用RANS 方法?；麄€(gè)邊界層，即脫體渦模擬（detachededdy simulation，DES）方法。近年來(lái)，雖然DES 方法仍然存在一些問(wèn)題，但是該方法還是得到了廣泛的應(yīng)用。目前，關(guān)于DES 方法的研究主要集中于如何緩解DES方法的“灰區(qū)”問(wèn)題。

改進(jìn)的延遲脫體渦模擬（improved delayed detached-eddy simulation，IDDES）方法是DES方法的一種改進(jìn)型，但其依然存在“灰區(qū)”問(wèn)題。為了解決IDDES 方法存在的“灰區(qū)”問(wèn)題，本文提出了一種剪切層自適應(yīng)的IDDES 方法。選取了具有可靠試驗(yàn)數(shù)據(jù)的導(dǎo)彈標(biāo)準(zhǔn)算例，通過(guò)引入適當(dāng)?shù)臄_動(dòng)，對(duì)導(dǎo)彈大攻角非對(duì)稱的分離形態(tài)進(jìn)行了詳細(xì)分析；介紹了數(shù)值模擬中，準(zhǔn)確預(yù)測(cè)細(xì)長(zhǎng)旋成體大攻角側(cè)向力的方法。

1 幾何模型與網(wǎng)格

圖1展示了本文中數(shù)值計(jì)算采用的導(dǎo)彈模型尺寸示意圖。其前體為Tangent曲線成型，彈體的直徑為152.4 mm，軸向長(zhǎng)2；后體為軸對(duì)稱圓柱體，軸向長(zhǎng)度為13。試驗(yàn)測(cè)試了不同雷諾數(shù)、不同攻角時(shí)導(dǎo)彈側(cè)向力的變化，本文的計(jì)算中基于導(dǎo)彈直徑的雷諾數(shù)Re為3.0×10，來(lái)流馬赫數(shù)為0.2。

計(jì)算中采用結(jié)構(gòu)網(wǎng)格對(duì)計(jì)算域進(jìn)行離散，圖2展示了彈體表面網(wǎng)格、攻角平面-的網(wǎng)格和垂直于攻角平面-平面的網(wǎng)格。彈體的軸向共布置了287個(gè)節(jié)點(diǎn)，垂直于彈體轉(zhuǎn)軸的截面上共布置328 個(gè)點(diǎn)，全模的總網(wǎng)格量為2 500萬(wàn)。為了很好地捕捉導(dǎo)彈大攻角下的渦結(jié)構(gòu)，在彈體的背風(fēng)面進(jìn)行網(wǎng)格加密，背風(fēng)面的網(wǎng)格占總網(wǎng)格的65%左右。為了排除數(shù)值計(jì)算帶來(lái)的不對(duì)稱干擾，網(wǎng)格生成過(guò)程中首先生成導(dǎo)彈半模的網(wǎng)格，在將網(wǎng)格對(duì)稱生成全模的網(wǎng)格，從而保證計(jì)算網(wǎng)格的完全對(duì)稱，排除網(wǎng)格的干擾。Kumar等在研究中發(fā)現(xiàn)，在數(shù)值計(jì)算中想準(zhǔn)確地預(yù)測(cè)試驗(yàn)中測(cè)得的側(cè)向力，必須添加適當(dāng)?shù)臄_動(dòng)。本文在計(jì)算中采用了文獻(xiàn)［9］中提出的添加擾動(dòng)的方式，在軸向距離導(dǎo)彈的頭部0.08的位置，-平面內(nèi)，安裝最大高度為0.004，長(zhǎng)度為0.04，最大寬度為0.004的擾動(dòng)，幾何尺寸如圖1所示，頭部擾動(dòng)位置網(wǎng)格如圖3所示。

圖1 導(dǎo)彈模型尺寸示意圖Fig.1 Schematic diagram of missile model size

圖2 導(dǎo)彈全局網(wǎng)格和頭部網(wǎng)格Fig.2 Global mesh and mesh around the missile head

圖3 擾動(dòng)及其網(wǎng)格Fig.3 Perturbation and the grids

2 數(shù)值方法

2.1 原始的SST-IDDES方法

為了解決IDDES 方法存在的“灰區(qū)”問(wèn)題，發(fā)展了剪切層自適應(yīng)的IDDES 方法。本部分對(duì)原始的剪切應(yīng)力傳輸IDDESS（shear-stress transport-IDDES，SST-IDDES）方法的框架進(jìn)行簡(jiǎn)單的介紹，從而方便對(duì)剪切層自適應(yīng)尺度的描述。IDDES 方法通過(guò)將RANS 的長(zhǎng)度尺度替換為IDDES 長(zhǎng)度尺度得以實(shí)現(xiàn)，IDDES 的湍動(dòng)能的輸運(yùn)方程如式（1）所示。

式中：為密度；為時(shí)間；為湍動(dòng)能；為層流黏性；μ為渦黏性；P為湍動(dòng)能的生成項(xiàng)；IDDES 的長(zhǎng)度尺度表達(dá)為

式中：＞0 時(shí)為提升函數(shù)，增加RANS/LES 切換之前RANS 區(qū)的渦黏性，避免此區(qū)域的渦黏性過(guò)度減??；當(dāng)來(lái)流不存在湍流度時(shí)，=0，IDDES 表現(xiàn)為DDES；為切換函數(shù)，詳細(xì)公式見文獻(xiàn)［19］。

式中：常數(shù)=0.15；為網(wǎng)格單元中心到壁面的距離；和分別為網(wǎng)格單元最大長(zhǎng)邊和網(wǎng)格在壁面垂直方向的網(wǎng)格尺度；切換函數(shù)=max{(1-)，}，其中與流動(dòng)相關(guān)，用于控制IDDES 分支中RANS/LES 的切換，僅與網(wǎng)格相關(guān)，用于控制壁面?；拇鬁u模擬（WMLES）分支中RANS/LES 的切換。在來(lái)流或者初場(chǎng)存在湍流脈動(dòng)時(shí)，=，IDDES變現(xiàn)為WMLES，對(duì)應(yīng)的長(zhǎng)度尺度為

式中：＞0 為提升函數(shù)，增加RANS/LES 切換之前RANS 區(qū)的渦黏性，避免此區(qū)域的渦黏性過(guò)度減小。當(dāng)來(lái)流不存在湍流度時(shí)，=0，IDDES 表現(xiàn)為DDES，=1-，對(duì)應(yīng)的長(zhǎng)度尺度為

2.2 剪切層自適應(yīng)長(zhǎng)度尺度

剪切層自適應(yīng)長(zhǎng)度尺度公式如式（6）所示。

式中：Δ為渦量矢量方向的網(wǎng)格尺度；n為單位化的渦量矢量；r(，=1，2，…，8)為網(wǎng)格單元任意兩頂點(diǎn)構(gòu)成的矢量。

3 結(jié)果與分析

圖4展示了計(jì)算得到的無(wú)量綱準(zhǔn)則瞬時(shí)等值面，等值面采用馬赫數(shù)進(jìn)行著色。準(zhǔn)則定量反映了流場(chǎng)中的相干結(jié)構(gòu)，可以用應(yīng)變率和渦量的函數(shù)來(lái)表示，即

由圖4可以看出：沒有添加擾動(dòng)時(shí)，導(dǎo)彈背風(fēng)面生成穩(wěn)定的一對(duì)旋渦，旋渦由導(dǎo)彈頭部脫出，逐漸向下游發(fā)展，在整個(gè)彈體的長(zhǎng)度上都沒有出現(xiàn)渦破碎；添加擾動(dòng)之后，旋渦在三分之一彈體長(zhǎng)度的位置處出現(xiàn)渦破碎，發(fā)展成為完全不對(duì)稱的渦形態(tài)。

圖4 有/無(wú)擾動(dòng)導(dǎo)彈渦結(jié)構(gòu)示意圖Fig.4 Vortex structure diagram of missile with/without perturbation

圖5展示了導(dǎo)彈不同軸向位置截面的壓力分布，圖中黑色圓圈表示試驗(yàn)測(cè)量值，藍(lán)色點(diǎn)表示不添加擾動(dòng)時(shí)的計(jì)算值，紅色實(shí)線為添加擾動(dòng)的計(jì)算值。試驗(yàn)測(cè)得的壓力分布明顯為非對(duì)稱結(jié)構(gòu)，左側(cè)的壓力分布吸力峰值更高，右側(cè)更低。不添加擾動(dòng)時(shí)計(jì)算得到的壓力分布在不同軸向位置基本為對(duì)稱的分布，添加擾動(dòng)后計(jì)算得到的壓力分布為非對(duì)稱結(jié)構(gòu)，且與試驗(yàn)結(jié)果吻合很好。證明了本文添加的擾動(dòng)是合理的，發(fā)展的數(shù)值方法也可以準(zhǔn)確地預(yù)測(cè)導(dǎo)彈大攻角下的側(cè)向力。

圖5 有/無(wú)擾動(dòng)導(dǎo)彈壓力分布Fig.5 Pressure distribution of missiles with/without perturbation

圖6和圖7分別展示了無(wú)/有擾情況下不同展向位置截面的流線與馬赫數(shù)分布云圖：無(wú)擾動(dòng)時(shí)不同軸向位置截面幾乎都為對(duì)稱的渦結(jié)構(gòu)，隨著向下游發(fā)展渦核逐漸增大；添加擾動(dòng)后，右側(cè)分離明顯提前，速度降低，壓力升高，所以產(chǎn)生方向?yàn)榈呢?fù)方向的側(cè)向力（在不同軸向位置截面都為非對(duì)稱的渦結(jié)構(gòu)，如圖7所示）。

圖6 無(wú)擾動(dòng)時(shí)導(dǎo)彈不同軸向位置截面流線與馬赫數(shù)分布云圖Fig.6 Streamline and Mach number distribution contour at different axial sections of missile，without perturbation

圖7 有擾動(dòng)時(shí)導(dǎo)彈不同軸向位置截面流線與馬赫數(shù)分布云圖Fig.7 Streamline and Mach number distribution contour at different axial sections of missile，with perturbation

4 結(jié)束語(yǔ)

采用數(shù)值計(jì)算的方法研究了導(dǎo)彈側(cè)向力產(chǎn)生的原因，分析了擾動(dòng)對(duì)導(dǎo)彈側(cè)向力的影響。數(shù)值方法描述試驗(yàn)中產(chǎn)生側(cè)向力的現(xiàn)象，需要人為引入合適的擾動(dòng)。采用帶剪切層自適應(yīng)的IDDES 方法，對(duì)有/無(wú)擾動(dòng)的導(dǎo)彈構(gòu)型進(jìn)行了計(jì)算。采用完全軸對(duì)稱的幾何模型和對(duì)稱的網(wǎng)格，在整個(gè)彈體的背風(fēng)面生成對(duì)稱的渦結(jié)構(gòu)，不會(huì)發(fā)生渦破碎，幾乎沒有產(chǎn)生側(cè)向力。在導(dǎo)彈頭部設(shè)置擾流片引入擾動(dòng)，會(huì)在導(dǎo)彈的背風(fēng)面誘導(dǎo)出完全不對(duì)稱的渦結(jié)構(gòu)，在添加擾動(dòng)一側(cè)，分離會(huì)明顯提前。分離提前使添加擾動(dòng)一側(cè)的速度降低、壓力升高，從而產(chǎn)生側(cè)向力。本文通過(guò)添加擾流片引入擾動(dòng)，計(jì)算得到的壓力分布與試驗(yàn)結(jié)果非常吻合。但是，較大的側(cè)向力不利于導(dǎo)彈的飛行控制。需要設(shè)計(jì)合理的流動(dòng)控制手段，抑制側(cè)向力的產(chǎn)生，實(shí)現(xiàn)導(dǎo)彈大攻角下的精確控制。