鄭忠華，范周琴，王子昂，余彬，張斌,*

1. 中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所 高超聲速沖壓發(fā)動機技術(shù)重點實驗室，綿陽 621000 2. 上海交通大學 航空航天學院，上海 200240

作為一類基礎(chǔ)流動模型，兩個同向旋轉(zhuǎn)渦旋之間的相互作用首先由飛機尾跡渦系抽象提出[1]。研究表明，合理控制外部自由來流經(jīng)過機翼產(chǎn)生的尾跡相互作用不僅有助于減小誘導阻力，更能從渦系不穩(wěn)定性的角度將尾跡對周圍氣流擾動的影響降至最低，提高飛行安全性[2]。除廣泛實際工程應(yīng)用外，深入認識渦相互作用流動機理對理解剪切層、湍流等流體前沿難題同樣具有重要科學意義，至今已吸引眾多頂尖學者通過實驗[3]及數(shù)值手段[4]進行研究。

首先，針對渦旋發(fā)生相互作用的條件，Saffman和Szeto[5]提出兩個渦旋距離與渦旋半徑之比(定義為渦對展弦比)是決定渦對內(nèi)部應(yīng)變強弱的重要參數(shù)。Melander等[6]對兩個渦量均勻分布的渦旋進行數(shù)值模擬，發(fā)現(xiàn)渦相互作用發(fā)生的臨界條件為渦對展弦比在0.3左右。Meunier等[7]通過渦旋旋轉(zhuǎn)角動量構(gòu)建渦旋有效半徑的概念，理論推導指出渦對臨界展弦比范圍為0.22～0.26且與渦旋雷諾數(shù)以及具體渦旋類型無關(guān)。關(guān)于渦相互作用的物理機制，Melander和Mcwilliams[8]在與渦對同速旋轉(zhuǎn)的參考坐標系中通過流線觀察渦相互作用過程指出虛擬渦旋(Ghost Vortex)與懸臂(Filament)形成是推動渦心相互靠近的無黏機制，并將渦對融合歸納為對稱-反對稱-再對稱的過程。在此基礎(chǔ)上為定量描述渦對靠近的快慢程度，Williamson[3]將渦量場分解為對稱分量與反對稱分量并證明渦對靠近的速度僅由反對稱分量誘導決定。渦量分解理論描述的渦對間距演化規(guī)律與模擬結(jié)果吻合良好。關(guān)于渦相互作用流動特征的詳細研究進展可參考Leweke等[9]的綜述文獻。

區(qū)別于外部尾跡流動速度遠低于聲速的不可壓縮流場特征，以沖壓發(fā)動機為代表的內(nèi)部流動軸向速度接近甚至遠遠超過當?shù)芈曀賉10]，隨之而來的強可壓縮性對流場演化過程的影響不可忽略。針對如何將燃料與氧化物有效混合燃燒這一長期制約沖壓發(fā)動機性能的核心問題[11]，近幾年，利用渦相互作用過程中引入的應(yīng)變作用提出了一種新型塔式燃料噴射可以有效將燃料通過渦流發(fā)生器卷入渦中，從而達到與周圍空氣混合增強的效果[12-13]。在此工程背景與實際需求下，研究流場可壓縮性對渦相互作用過程的影響這一基礎(chǔ)理論問題對優(yōu)化設(shè)計塔式燃料噴射裝置具有一定借鑒意義。本文基于可壓縮Navier-Stokes方程求解算法，結(jié)合可壓縮壓力泊松方程設(shè)置初始條件模擬不同馬赫數(shù)渦對的相互作用過程，通過比較流動結(jié)構(gòu)及動力學特征的異同，重點探索可壓縮性對渦相互作用流場演化過程的影響，并在此基礎(chǔ)上初步提出服務(wù)于工程設(shè)計的可壓縮時間尺度修正關(guān)系。

1 數(shù)值模擬方法

1.1 流動控制方程與數(shù)值方法

渦相互作用問題流場演化的控制方程是可壓縮黏性Navier-Stokes方程，其在二維流動中表達形式為

(1)

式中：U代表由守恒參數(shù)組成的向量；向量場F、G和Fv、Gv分別代表對流通量和耗散通量，且

U=[ρ,ρu,ρv,ρE,ρY1,…，ρYNS-1]T

F=[ρ,ρu2+p,ρuv,(ρE+p)u,ρ1u,…,

ρNS-1u]T

G=[ρ,ρv2+p,ρuv,(ρE+p)v,ρ1v,…,

ρNS-1v]T

(2)

式中：ρ、p、E分別代表混合物密度、壓力、單位質(zhì)量總能量；u和v為混合物速度矢量的兩個分量；ρi和Yi為組分i的密度和質(zhì)量分數(shù)；NS為總組分數(shù)。黏性力τ和熱流量q為

(3)

其中：T代表溫度；H代表比焓；μ和λ分別代表氣動黏性和混合物熱傳導率，由Wilke 半經(jīng)驗公式計算得到；對于高速流動問題，質(zhì)量擴散系數(shù)Di為

(4)

二元擴散系數(shù)Dij為

(5)

其中：Xj為組分j的摩爾分數(shù)；文中假設(shè)施密特數(shù)Sc取0.5[14]。引入理想氣體狀態(tài)方程封閉方程組：

(6)

式中：Ru為普適氣體常量；Mi為組分i對應(yīng)的分子量?？偰芰縀包含氣體動能及氣體內(nèi)能，即

(7)

cpi=Ru(a1i+a2iT+a3iT2+a4iT3+a5iT4)

(8)

其中：系數(shù)a1i,a2i,…,a5i可以由權(quán)威的NASA 熱力學參數(shù)數(shù)據(jù)庫查得。上述數(shù)學模型在無量綱后通過有限體積法進行離散。時間推進采用三階精度的TVD Runge-Kutta方法，對流項采用五階WENO格式離散[15]，黏性項采用中心差分法離散。

1.2 基于可壓縮壓力泊松方程設(shè)置的初始條件

在基于Navier-Stokes方程的非定常計算中，初始壓力場設(shè)置對模擬結(jié)果至關(guān)重要。不合理的初始條件設(shè)置會在流場演化初期額外引入非物理的噪聲轉(zhuǎn)捩，由此可能導致完全錯誤的后續(xù)流場發(fā)展過程[16]。本文考慮兩個完全相同(強度、有效半徑均相等且旋轉(zhuǎn)方向相同)Lamb-Oseen渦之間的相互作用過程。初始速度場由每個渦旋各自依據(jù)Biot-Savart定律的誘導速度場疊加而成：

(9)

式中：r1與r2為流場中某一點到兩個渦心的距離；(x1,y1) 與(x2,y2)為兩個渦心的空間位置；Γ為單個渦旋強度；渦旋有效半徑[7]a為

(10)

其中:ω為渦量；r為到渦心的空間距離；A為積分區(qū)域。在給定初始速度場基礎(chǔ)上，初始熱力學參數(shù)分布基于壓力泊松方程結(jié)合等熵流動假設(shè)構(gòu)建。該方法首先由Harlow和Welch基于不可壓縮Navier-Stokes方程組推導提出并應(yīng)用于不可壓縮流動模擬中[17]，具有不可壓縮條件內(nèi)在相容、數(shù)值計算穩(wěn)定、方程形式簡潔易于構(gòu)造高精度格式等優(yōu)點。在忽略黏性項與非定常項后壓力泊松方程可直接應(yīng)用于初始條件設(shè)置中，具體形式為

P,ii=-ρ0ui,juj,i

(11)

式中：求和記號為愛因斯坦標記法。在不可壓縮流場中，密度變化不大，近似為一常數(shù)ρ0。

在可壓縮流動中，為盡可能消除流場演化初始階段由于壓力場與速度場不相容導致的噪聲轉(zhuǎn)捩現(xiàn)象，通常假設(shè)初始流場滿足等熵假設(shè)[18]。另外聯(lián)立可壓縮壓力泊松方程求解獲得所有初始熱力學參數(shù)分布：

(12)

根據(jù)等熵條件，壓力與密度之間存在如下聯(lián)系：

(13)

關(guān)于泊松方程數(shù)值算法，采用文獻[19]中的緊致差分格式進行離散求解。在二維矩形網(wǎng)格(x及y方向上的網(wǎng)格精度分別為h和k)上考慮含源項f的泊松方程，形如：

φxx+φyy=f(x,y)

(14)

基于Keriss緊致差分公式[20]，可以得到九結(jié)點矩形網(wǎng)絡(luò)下的高精度優(yōu)化差分格式為

(15)

式中：r=h/k依賴于x與y方向上網(wǎng)格精度；上標n代表第n次迭代的值。泊松方程求解離散網(wǎng)格編號示意圖如圖1所示，算法的具體構(gòu)造細節(jié)詳見文獻[19]。在每次迭代中均需調(diào)用一次壓力邊界條件完成全場計算，最后以相鄰兩次迭代中全場最大壓力變化Δpmax<10-6為判據(jù)認為迭代收斂，終止計算。

圖1 泊松方程求解離散網(wǎng)格編號示意圖

1.3 數(shù)值驗證

為驗證上述Navier-Stokes非定常模擬算法與初始條件設(shè)置對模擬渦相互作用問題的有效性，采用文獻[21]中的算例進行數(shù)值驗證。初始渦對幾何參數(shù)、渦量分布形式、計算域大小及網(wǎng)格精度均與文獻[21]相同。初始熱力學參數(shù)分布通過求解壓力泊松方程設(shè)置。本文模擬流場演化過程與文獻中大渦模擬結(jié)果對比如圖2所示。在兩個渦旋相互連接靠近的過程中渦對發(fā)生扭曲變形并向外傳輸渦量形成懸臂結(jié)構(gòu)。在黏性擴散作用的推動下兩個渦旋中心逐漸合并最終融合為單個渦旋。

另外，本文根據(jù)渦對間距與渦旋有效半徑(根據(jù)式(10)計算)隨時間的演化規(guī)律定量反映計算結(jié)果與文獻之間的差異，如圖3所示?？梢园l(fā)現(xiàn)在演化前中期本文渦對間距與渦旋有效半徑計算結(jié)果均與文獻吻合良好(最大誤差不超過5%)，在演化后期(無量綱渦對間距b/b0<0.3)由于渦對充分靠近難以準確識別渦心位置。此時根據(jù)式(10)計算的數(shù)據(jù)無法真實反映渦旋有效半徑，不再具有參考比較意義。總體而言無論是定性的渦量場發(fā)展過程還是定量的渦對幾何參數(shù)演化規(guī)律都基本證明本文采用的數(shù)值方法可用于研究渦相互作用物理問題。

圖2 本文數(shù)值模擬渦量云圖與文獻結(jié)果對比

圖3 無量綱渦對間距、渦旋有效半徑數(shù)值驗證與文獻結(jié)果對比

2 算例設(shè)置條件

控制渦相互作用進程的3個重要無量綱參數(shù)分別為：渦對展弦比、渦旋雷諾數(shù)與渦旋馬赫數(shù)[4]。其中，展弦比決定渦對的幾何構(gòu)型，而雷諾數(shù)與馬赫數(shù)反映流場演化的動力學機制。本文重點聚焦的問題是可壓縮效應(yīng)對相同強度渦相互作用過程的影響，因而在保證渦對展弦比與雷諾數(shù)不變的情況下，通過改變渦旋馬赫數(shù)以調(diào)整流場可壓縮性。為在清晰呈現(xiàn)結(jié)果的前提下節(jié)省計算資源，所有算例選取初始渦對展弦比為0.177，渦旋雷諾數(shù)為5 000，以避免相互作用進程過于冗長。另一方面在該雷諾數(shù)下橢圓不穩(wěn)定性主導的湍流三維效應(yīng)尚不顯著[22]，只適用于二維數(shù)值計算方法。本文模擬算例的計算域及計算條件如圖4所示，其中渦旋模型為Lamb-Oseen渦，渦量分布為

(16)

式中：U0為渦內(nèi)最大速度，在渦核半徑r0處取得。渦旋強度以環(huán)量Γ描述；渦旋馬赫數(shù)及雷諾數(shù)可分別定義為

圖4 計算域及模擬條件設(shè)置

Ma=U0/c∞

(17)

式中：c∞為無窮遠處聲速。

ReΓ=Γ/υ

(18)

式中：υ為流體動力黏度。

計算域大小及網(wǎng)格解析度對于模擬渦相互作用同樣具有顯著影響。采用均勻網(wǎng)格進行數(shù)值計算(網(wǎng)格精度Δx/b0=Δy/b0=0.01)，計算域大小為20b0×20b0且計算域邊界均設(shè)置為外推條件以消除壁面效應(yīng)對結(jié)果的潛在作用[17]。時間步長為Δt=10-9s。經(jīng)過驗證，該網(wǎng)格精度與時間步長能夠保證充分的空間與時間離散，基本能夠消除數(shù)值黏性對結(jié)果的影響。

3 結(jié)果與討論

3.1 不同馬赫數(shù)渦對融合演化過程

不同馬赫數(shù)渦對的渦量場形態(tài)演化過程如圖5所示。為保證算例結(jié)果之間的對比性與普適性，物理時間采用無量綱化處理:

(19)

式中：特征時間t*的物理意義是兩個具有相同強度Γ、相距b0的點渦在無黏環(huán)境中相互旋轉(zhuǎn)運動的周期；T*為無量綱演化時間。可以發(fā)現(xiàn)，不同馬赫數(shù)渦對相互作用進程總體相似：首先在黏性擴散的作用下渦旋半徑逐漸增加，此時渦對僅呈現(xiàn)細微的變形且兩個渦旋在空間中沒有直接接觸(T*=0.41)；當渦旋半徑不斷變大與渦對間距之比達到臨界閾值后，渦對系統(tǒng)開始融合。此時兩個渦旋已相互連接并受相互誘導的應(yīng)變作用發(fā)生扭曲變形，部分渦量在旋轉(zhuǎn)速度場下向外輸運形成典型的懸臂結(jié)構(gòu)。在整個發(fā)展過程中流場均保持中心對稱特征。另一方面，在相同無量綱時間下不同馬赫數(shù)渦對的旋轉(zhuǎn)角度存在顯著的差異，這意味著以馬赫數(shù)表征的流場可壓縮性能夠改變渦對旋轉(zhuǎn)速度從而對融合過程產(chǎn)生一定的影響。此外，在渦對形態(tài)方面高馬赫數(shù)渦旋在渦心處拉伸效應(yīng)更明顯，而低馬赫數(shù)渦對渦量分布相對集中。

圖5 不同馬赫數(shù)渦對相互作用流場演化規(guī)律

3.2 不同馬赫數(shù)渦對間距演化規(guī)律

為了更加深入地刻畫可壓縮性對渦相互作用流場的影響，計算了不同算例渦對間距隨時間演化規(guī)律，如圖6所示。在第1階段中，渦對間距基本上保持不變，可以近似地認為渦對維持在相互圍繞旋轉(zhuǎn)的穩(wěn)定狀態(tài)，這也與渦量云圖中直觀獲得的結(jié)論一致。在第2階段中，渦對間距急劇減小，其相互靠近的速率可通過Biot-Savart誘導定律確定。隨著兩個渦旋之間距離的減小，自誘導產(chǎn)生對流效應(yīng)對流場發(fā)展的主導地位逐漸削弱。隨后渦對受無黏對流效應(yīng)與黏性擴散效應(yīng)競爭機制的影響，渦對間距出現(xiàn)振蕩并且在黏性擴散的推動下完成最終的融合。比較不同馬赫數(shù)渦對間距演化過程的差異發(fā)現(xiàn)：在相同無量綱時間下，低馬赫數(shù)渦對間距較小，融合較為充分；換而言之，若將渦對從初始狀態(tài)減小到相同無量綱間距，高馬赫數(shù)渦旋所需的時間更長?？梢哉J為可壓縮性具有抑制渦相互作用發(fā)展的影響。類似現(xiàn)象在可壓縮剪切層中引起了廣泛的研究，并指出可壓縮性起到維持流場穩(wěn)定的作用[23]。值得注意的是，圖中圓圈標注的數(shù)據(jù)為不可壓縮渦相互作用的融合過程，可視作渦旋馬赫數(shù)Ma趨向于0的極限理想情況，與低馬赫數(shù)算例(可壓縮性影響不明顯)的模擬結(jié)果基本吻合。

圖6 不同馬赫數(shù)渦對間距隨時間演化規(guī)律

3.3 不同馬赫數(shù)渦對融合臨界展弦比

在關(guān)于渦相互作用模型的理論研究中，確定系統(tǒng)開始融合的臨界條件對于進一步認識渦相互作用物理過程及其規(guī)律奠定了基礎(chǔ)。本文基于數(shù)值模擬結(jié)果，討論可壓縮性對于渦對臨界展弦比的影響。不同馬赫數(shù)渦對有效半徑以及渦對間距隨時間變化規(guī)律如圖7所示：以無量綱渦對間距b/b0=0.9為臨界標志[3](即相互作用第1、第2階段的分界點)可以發(fā)現(xiàn)，隨著馬赫數(shù)的增加，渦對系統(tǒng)的臨界展弦比變化不大，基本保持在理論推導獲得的(a/b0)cr=0.24左右[7]。因而可以認為，以渦對展弦比描述的系統(tǒng)融合準則對于渦旋類型、渦旋雷諾數(shù)以及渦旋馬赫數(shù)均具有較好的魯棒性。然而，另一方面，相互作用第二階段開始的臨界無量綱時間隨著馬赫數(shù)的增加顯著提高，這也從側(cè)面印證了可壓縮性具有推遲相互作用進程的作用。進一步地，為更深入地認識并預測可壓縮渦相互作用發(fā)展的物理規(guī)律，本文下一節(jié)將針對不可壓縮渦相互作用過程中的特征時間尺度，引入可壓縮性修正，以統(tǒng)一不同馬赫數(shù)渦對融合的進程。

圖7 不同馬赫數(shù)渦對有效半徑及間距時間演化規(guī)律

4 相互作用時間尺度的可壓縮性修正

在不可壓縮渦相互作用模型中(Ma→0的極限情況)，流動特征時間定義為

(20)

根據(jù)上文結(jié)果可以發(fā)現(xiàn)，該特征時間難以對不同馬赫數(shù)渦對相互作用過程進行歸一，其原因是可壓縮流動與低速流動的內(nèi)在機制存在差異，而該特征時間無法體現(xiàn)出可壓縮性強弱對于流場演化的影響?？紤]到可壓縮性具有推遲相互作用進程的影響，本節(jié)在不可壓縮特征時間的定義中引入依賴于渦旋馬赫數(shù)的可壓縮時間尺度修正，為進一步分析可壓縮渦相互作用的動力學行為提供一定的理論基礎(chǔ)。首先，本文通過不同馬赫數(shù)渦心處密度的時間歷程從物理角度說明流場中可壓縮性相對強弱的變化規(guī)律，如圖8所示。其中渦心密度的時間變化率采用中心差分格式近似，即

(21)

式中：F為本文考慮隨時間變化渦中心處的密度；Δt為時間步長。

圖8 不同馬赫數(shù)渦對渦心處無量綱密度和渦心處無量綱密度時間變化率的演化過程

渦心密度的時間變化率關(guān)于時間步長具有二階精度。經(jīng)過與高階數(shù)值離散公式結(jié)果比較可以認為計算精度基本滿足要求。

可以發(fā)現(xiàn)，不同馬赫數(shù)渦對渦心處密度演化規(guī)律關(guān)于時間具有較好的統(tǒng)一性，其中不同馬赫數(shù)渦心密度出現(xiàn)拐點的時間相互一致，如圖8(b)所示。該結(jié)果表明，僅根據(jù)馬赫數(shù)即可反映可壓縮性的強弱而不依賴于流場具體演化過程的影響。因此在構(gòu)造可壓縮時間尺度修正時首先將馬赫數(shù)視作給定常數(shù)而不考慮其與時間尺度之間的耦合作用。

關(guān)于可壓縮修正項的具體表達形式，待定參數(shù)少且滿足同類物理量綱相加原則的冪律模化(Power Law)在剪切層[24]、湍流[25]、稀薄氣體流動[26]建模中應(yīng)用廣泛。其中渦相互作用模型大量存在于剪切層流動中可以認為內(nèi)在符合相似的模化規(guī)律，然目前在可壓縮問題中尚未推廣。通過渦旋馬赫數(shù)的冪律描述可壓縮性強弱與流動特征時間的關(guān)系，可壓縮時間尺度T表示為

(22)

圖9 通過引入可壓縮時間尺度修正后統(tǒng)一的無量綱渦對間距以及渦對有效半徑演化規(guī)律

5 結(jié) 論

針對流動可壓縮性對渦相互作用動力學特征的影響這一科學問題，基于經(jīng)過驗證的Navier-Stokes算法以及壓力泊松方程初始條件設(shè)置方法數(shù)值模擬了不同馬赫數(shù)渦對相互作用流場演化過程，主要結(jié)論如下：

1) 可壓縮與低速渦對相互作用流動結(jié)構(gòu)與整體過程相似。另一方面，以高馬赫數(shù)表征的可壓縮性在提高渦量分布拉伸程度的同時具有抑制渦對系統(tǒng)融合，延遲相互作用進程的作用。

2) 以渦對展弦比描述的系統(tǒng)融合準則對流場可壓縮性具有較強魯棒性，然而臨界無量綱時間隨馬赫數(shù)增加明顯提高。

為在無量綱時間意義下統(tǒng)一不同馬赫數(shù)渦對相互作用的進程，初步提出了可壓縮時間尺度修正理論，其在不同渦旋雷諾數(shù)條件的適用性及其物理機制值得深入研究。

致 謝

感謝上海交通大學超算“π”為本研究提供的計算資源。