細長旋成體大攻角非對稱渦模擬的擾動引入方式研究

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

（清華大學航天航空學院，北京 100084）

0 引 言

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

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

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

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

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

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

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

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

圖1 導彈模型尺寸示意圖Fig.1 Schematic diagram of missile model size

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

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

2 數(shù)值方法

2.1 原始的SST-IDDES方法

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

式中：為密度；為時間；為湍動能；為層流黏性；μ為渦黏性；P為湍動能的生成項；IDDES 的長度尺度表達為

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

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

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

2.2 剪切層自適應長度尺度

剪切層自適應長度尺度公式如式（6）所示。

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

3 結(jié)果與分析

圖4展示了計算得到的無量綱準則瞬時等值面，等值面采用馬赫數(shù)進行著色。準則定量反映了流場中的相干結(jié)構(gòu)，可以用應變率和渦量的函數(shù)來表示，即

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

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

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

圖5 有/無擾動導彈壓力分布Fig.5 Pressure distribution of missiles with/without perturbation

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

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

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

4 結(jié)束語

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