流場可壓縮性對渦相互作用影響的數值研究

2020-03-25 11:02:04鄭忠華范周琴王子昂余彬張斌

航空學報 2020年2期

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

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

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

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

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

1 數值模擬方法

1.1 流動控制方程與數值方法

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

(1)

式中：U代表由守恒參數組成的向量；向量場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分別代表混合物密度、壓力、單位質量總能量；u和v為混合物速度矢量的兩個分量；ρi和Yi為組分i的密度和質量分數；NS為總組分數。黏性力τ和熱流量q為

(3)

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

(4)

二元擴散系數Dij為

(5)

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

(6)

式中：Ru為普適氣體常量；Mi為組分i對應的分子量。總能量E包含氣體動能及氣體內能，即

(7)

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

(8)

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

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

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

(9)

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

(10)

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

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

(11)

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

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

(12)

根據等熵條件，壓力與密度之間存在如下聯系：

(13)

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

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

(14)

基于Keriss緊致差分公式[20]，可以得到九結點矩形網絡下的高精度優化差分格式為

(15)

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

圖1 泊松方程求解離散網格編號示意圖

1.3 數值驗證

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

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

圖2 本文數值模擬渦量云圖與文獻結果對比

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

2 算例設置條件

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

(16)

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

圖4 計算域及模擬條件設置

Ma=U0/c∞

(17)

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

ReΓ=Γ/υ

(18)

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

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

3 結果與討論

3.1 不同馬赫數渦對融合演化過程

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

(19)

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

圖5 不同馬赫數渦對相互作用流場演化規律

3.2 不同馬赫數渦對間距演化規律

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

圖6 不同馬赫數渦對間距隨時間演化規律

3.3 不同馬赫數渦對融合臨界展弦比

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

圖7 不同馬赫數渦對有效半徑及間距時間演化規律

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

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

(20)

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

(21)

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

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

渦心密度的時間變化率關于時間步長具有二階精度。經過與高階數值離散公式結果比較可以認為計算精度基本滿足要求。

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

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

(22)