指定流量分配系數的多回流出口邊界算法

許開龍，劉再剛，姜勝利，王星，張磐

1．中物院高性能數值模擬軟件中心，北京 100088

2．北京應用物理與計算數學研究所，北京 100088

在航空發動機燃燒室數值模擬中，計算域存在燃氣和冷卻氣等多個出口邊界，其出口流量分配關系通常根據工程設計要求或實驗獲得，模擬中需要引入特殊的邊界處理方法以獲得邊界物理量分布作為定解條件，為邊界處理方法帶來較大挑戰。此外，燃燒室中火焰筒等局部的燃燒組織是數值模擬的關注重點之一，截取燃燒室的部分結構作為計算域可顯著降低數值模擬的計算量，但常規出口邊界條件算法因出口存在回流或未充分發展無法保證計算穩定性和精度。因此，研發可指定流量分配系數、允許出口回流的出口邊界算法對提升發動機燃燒室數值模擬軟件的工程適用性和仿真置信度有重要意義。

出口邊界的處理方式對數值模擬的穩定性和計算準確性有顯著影響，是不可壓縮流體數值模擬中的一個挑戰性問題［1-2］。經過數十年的發展，已有諸多不可壓縮流動模擬中出口邊界處理方 法 的 研 究。Sani 和Gresho［3］綜 述 了 經 典 的 出口邊界處理方法，其中基于物理量局部單向化假設的對流邊界條件［4］及零法向梯度的充分發展邊界條件應用較為廣泛［5-6］，但該方法在邊界上存在渦團或回流時，會出現數值不穩定現象甚至直接發散，且降低時間步長或增加網格密度并不能改善這種數值不穩定性［7］。在實際應用中，常將出口邊界面置于主流區下游足夠遠處，使得流動產生的渦團在到達出口邊界前充分耗散。但主流區渦團的發展和流動雷諾數相關，計算域需要隨著雷諾數的增加而擴展，造成了很大的計算資源開 銷，如Steggel 和Rockliff［8］采 用 充 分 發 展 邊 界條件計算二維圓柱擾流問題，出口邊界位于下游125 倍圓柱直徑處，主流區僅占計算域一小部分。此外，在海洋環流、大氣流動、血液流動等問題中并不存在明確的出口，必須將計算域截斷。

針對出口流動未充分發展的情況，Papanastasiou 等［9］首次在有限元離散下提出了自由邊界條件，該方法將控制方程的離散拓展到邊界上，并將邊界上物理量和計算域內的值關聯，無需用戶額外提供流場信息。何吉歡［10］在有限元框架下進一步討論了自由邊界條件對復雜流場的處理 能 力。Li 和Tao［11］提 出 了 法 向 速 度 滿 足 局 部單元質量守恒、切向速度滿足齊次Neumann 條件的邊界處理方法，并根據全局質量守恒進一步修正法向速度大小，成功實現出口邊界上存在回流時傳熱問題的正確模擬；該方法簡單實用，僅需修改速度邊界條件而無需調整求解流程。Hasan等［1］通過改進邊界上速度插值方法，實現了計算域縮減至6～8 倍特征尺寸時鈍體繞流問題的正確計算。Xue 和Barton［12］提出了利用質量流量處理進/出口邊界條件的方法，將出口邊界處理方法拓展到了多出口的情形，并在二維T 形分叉管算例中驗證了計算域大幅縮減、出口邊界上存在回流情況時該方法的表現；該方法的優點是，以實際工程中可測量的出口質量流量作為邊界條件輸入信息，適用于多出口邊界情形，提高了工程適用性，但該方法需要修正質量通量的計算方法，對整體的求解流程有一定影響。

不可壓縮流體數值模擬中的多出口邊界和自由回流出口邊界的處理是充滿挑戰的問題，已有研究發展了適用于不同離散格式的方法，并取得了一定的應用效果，但面向復雜流動［13］、多相流［14］、多開口邊界［12］等特征的方法仍在不斷發展中。本文擬基于工程中可測量的流量分配系數確定各出口質量流率，發展了可指定出口流量分配系數、出口存在回流的出口邊界條件算法，并集成在自主研發的高性能計算流體數值模擬軟件中，驗證了該出口邊界條件的適用性和精度。

1 控制方程和數值方法

1.1 控制方程組及其離散方法

本文求解的不可壓縮Navier-Stokes（N-S）方程形式為

式中：ρ為流體密度；u為速度；p為壓力；Fl為第l個體積力；μ為流體動力黏度。

采用非結構網格有限體積法對式（1）和式（2）進行空間離散，速度和壓力的離散使用基于單元的同位網格和Rhie-Chow 插值。整個計算域被劃分為互不重復的微元控制體，取其中一個微元控制體Ω，包圍該控制體的內部面元集合為I(Ω)，邊界面元集合為B(Ω)，面元面積矢量為S，如圖1 所示。使用高斯定理將控制體上散度積分轉為通量的表面積分，可以得到積分形式的控制方程為

圖1 計算域邊界單元示意圖Fig.1 Schematic of computational domain boundary cell

式中：下標“i”表示內部面元索引號；下標“b”表示邊界面元索引號；=ρiui·Si和=ρbub·Sb分別為內部面元和邊界面元上的流體質量通量，ρi、ui、Si和ρb、ub、Sb分別為內部面元和邊界面元上的密度、速度和面積；|Ω|為微元控制體體積。針對離散后的N-S 方程，采用SIMPLEC 算法實現半隱式速度-壓力的耦合計算，采用二階Crank-Nicolson 格式離散時間項，采用二階線性迎風格式計算對流通量，采用中心差分格式計算擴散通量。

1.2 出口邊界條件

式（3）和 式（4）中、、-(μ?u)b·Sb分別表示邊界面元b∈B(Ω)上質量通量、動量對流通量和動量擴散通量，是邊界條件處理關注的重點。在本文中，壓力的出口邊界條件設置為Dirichlet 邊界條件，而速度u邊界條件是研究的主要內容。針對圖1 所示控制體，可以采用以下線性關系式描述任意邊界面元上的速度和其鄰接單元上速度的關系：

式中：uF為邊界面元中心F點的速度；quF為單位面積法向通量。為處理網格畸變嚴重的情形，需要用式（7）將單元中心P點的值重構到P'得到uP'再進行計算。邊界條件處理就是確定各面元上的值。

在常規邊界處理中，邊界上物理量的值可根據實驗測量、物理觀察等得到，如速度入口邊界、壁面邊界、對稱邊界等，進而確定4 個邊界系數的值。但對于出流邊界上，尤其是未達到充分發展的情況下，邊界上物理量的值是未知的，且很難通過實驗測量獲得邊界上的值或者通量分布。

以圖1 中F點所在邊界面元上邊界處理為例，將邊界面元速度分解到局部坐標系（n，τ，η）下，如式（8）所示。假設該出口邊界上所有邊界面元切向速度分量沿法向的梯度為0，而法向速度分量沿法向的變化率為一個常數，記為R，可以得到邊界面元上3 個速度分量的計算式為

該方法中實際是將常規速度Neumann 型邊界條件變為了Robin 型邊界條件。

法向速度變化率R的確定方法需要重點關注。與Li 和Tao［11］的思路類似，本文基于全局質量守恒確定不可壓縮流動問題中的出口邊界條件。在多出口條件下，引入出口質量流量分配系數αj進行處理：假設計算域有J個出口邊界，計算域所有入口邊界質量流量及計算域內新增質量（質量源空間積分）之和為Mtot，各出口邊界流量分配系數為αj（j=1， 2，…J；ΣJj=1αj=1.0），在不可壓縮流體模擬中，全局質量守恒是保證求解穩定的必要條件［12］，可以得到出口目標質量流量Mj為

根據各出口面目標質量約束，結合式（9）、式（12）可以得到第j個出口邊界上法向速度變化系數Rj為

結合式（9）～式（11）、式（13），可以得到邊界法向速度、切向速度的迭代更新方法為

式中：k表示當前迭代步。

上述出口邊界速度值的構造方法實際上與Papanastasiou 等［9］、Li和Tao［11］實現的自由邊界條件的思想是一致的，即雖然出口邊界面元上速度量是未知的，但當前時間步邊界單元的速度場是嚴格按照N-S 方程的求解得到的，是符合物理規律的、自然的。基于該邊界單元速度場構造邊界速度值，并以出口面質量流量偏差作為邊界速度的修正依據，也是合理的。相比于Xue 和Barton［12］提出的方法，該方法可按照常規的速度邊界條件進行實現，無需改動N-S 方程求解流程，易于實現。在充分發展情況下，整體質量守恒已經滿足，對應于Rj=0 的情形，結果是完全一致的。在應用實踐中發現，上述速度邊界條件計算方法在出口面上存在回流且網格畸變較為嚴重時，可能出現數值不穩定現象，基于P'點法向速度方向修正F點速度方向可以有效改善該不穩定性，修正方法為

式中：sgn(x)為符號函數，x＞ 0 時取值為1，否則為-1。

1.3 算法集成

考慮到實際工程問題構型復雜、計算規模大的特點，本文作者依托中物院高性能數值模擬軟件中心自主研制的非結構網格并行編程框架JAUMIN［15］研制了不可壓N-S 方程的SIMPLE，SIMPLEC 和PISO 求解算法，隱式歐拉和Crank-Nicolson 時間離散格式，一階迎風、二階線性迎風和中心差分空間離散格式，研制形成了面向工程實用和物理模型高效集成的高性能計算流體數值模擬軟件。該軟件基于非結構網格求解三維、不可壓近似、穩態/非穩態形式的N-S 方程，支持常見類型的三維凸多面體網格，集成了工程常用的雷諾時均類湍流模型、大渦模擬類湍流模型、傳熱模型、輸運方程求解功能、時空離散格式以及進口、出口、壁面、對稱邊界等典型物理邊界條件。目前，該軟件已成功應用于復雜結構湍流、航空發動機兩相湍流燃燒、流-固耦合傳熱等問題的數值模擬，并基于廣州“天河二號”超級計算機實現萬核、1.6 億四面體-六面體混合網格的流動模擬，并行計算效率超過70%，表現出良好的大規模并行擴展能力。

本文將發展的出口邊界處理算法集成到自主研制的高性能計算流體模擬軟件中，基于典型算例進行了方法的正確性驗證。

2 方法驗證

本節介紹了研究中開展的系列數值實驗及結果分析，分別從流量分配、網格類型、計算域截斷程度、雷諾數等維度驗證了上述出口邊界條件，并基于TECFLAM 旋流燃燒器流場模擬驗證了所給方法在復雜構型、流動特征中的表現。在2.1 節，針對90°T 形分叉管、60°Y 形分叉管2 種網格類型算例，設定系列流量分配系數，驗證了該方法在不同網格類型下流量分配計算的正確性；在2.2 節，針對六面體網格T 形分叉管算例，依次縮減側面出口段長度，使得出口面流動狀態在充分發展和強回流之間變化，考察了出口流動未充分發展時速度場的預測結果。Dong 和Shen［7］指出，流動雷諾數是影響出口邊界條件數值穩定性的重要因素，在2.3 節，對比研究了雷諾數在200～50 000 之間變化時該方法的適用性。在2.4 節，介紹了該方法在TECFLAM 旋流燃燒器冷態流動模擬中的模擬驗證情況。

2.1 流量分配計算驗證

2.1.1 六面體網格算例：90°T 形分叉管

采用Hayes 等［16］提出 的90°T 形分叉 管算例驗證六面體網格流量分配的正確性，算例構型如圖2 所示。該構型入口段寬度為W，長度為Lin，設定為速度入口邊界條件；主出口段長度為Lm，設定為本文中的出口邊界條件，流量分配系數為α1；側面出口段長度為Ls，設定為本文發展的出口邊界條件，流量分配系數為α2；垂直于流動平面方向的2 個面設定為對稱面，其余部分為壁面邊界條件。通過組合入口、出口段寬度、出口流量分配系數以及雷諾數等，可以產生系列復雜的流動行為，如流動分離、再附著等現象。

圖2 90°T 形分叉管算例構型及網格剖分示意圖Fig.2 Schematic of computational domain of a duct with a 90° T-branch discretized by structured meshes

基于該構型，本節數值實驗中重點關注流量分配系數和側面出口段長度的影響，固定W=1 m，Lin=3W，Lm=10W，入口處速度設定為充分發展狀態下近似拋物線分布uin=4x(1-x)，在Ls=10W和Ls=W這2 種 情 況 下，α1∶α2在0∶1.0～1.0∶0 之間變化。采用六面體網格離散計算域，網格數分別為7 464 和3 768，考慮到該設置下雷諾數較小，采用穩態求解器進行模擬。

圖3展示了不同流量分配系數下速度大小的分布。圖3（a）為Ls=10W的結果，在α1∶α2=0∶1.0 的設置下，相當于主分支出口被封閉，流體全部從右側分支出口流出，進出口的速度大小及分布非常接近。隨著α1∶α2的增加，主分支出口速度值逐漸增加，更多的流體從主分支出口流出，該規律及速度場的分布與預期是相符的。圖3（b）展示了Ls=W的結果，從圖中可以看出該組速度場計算結果定性合理，速度分布與圖3（a）中對應計算域結果非常相似。需要注意的是，在入口段與右側分支交叉拐角處存在低速回流區，說明盡管該組設置中側面分支計算域僅為圖3（a）結果的1/10，出口流動也遠未達到充分發展的要求，但本文提出的出口邊界處理方法能很好地預測流量分配結果。關于計算域縮減情況下流場的預測，將在2.2 節中詳細闡述。

圖3 90°T 形分叉管不同流量分配下的速度分布Fig.3 Velocity distribution of duct case with 90° T-branch at various mass flow ratios

表1 進一步列出了Ls=W，α1∶α2在0∶1.0～1.0∶0 這2 個極端之間變化時，計算收斂時入口流量、主分支出口1 流量、主分支出口2 流量的統計結果，以及實際計算得到的流量分配系數α1和α2。從統計結果來看，采用本文發展的出口邊界條件計算得到的流量分配結果與設定值高度吻合，出口流量的最大相對偏差在0.001%量級。

表1 90°T 形分叉管算例流量統計結果Table 1 Mass flow rate statistics results of duct case with 90° T-branch

2.1.2 四面體網格湍流算例：60°Y 形分叉管

采用如圖4 所示的三維60°Y 形分叉管算例考察本文方法在四面體網格中的表現。圖中管道直徑Φ=10 mm，入口段長度50 mm，2 個分叉管呈60°夾角，分叉段長度均為50 mm，分叉管上分支出口記為出口1，流量分配系數為α1，下分支出口記為出口2，流量分配系數為α2。

計算中采用均勻入口邊界條件，速度大小為10 m/s，上下2 個分支出口采用本文的出口邊界條件，其余邊界為無滑移壁面邊界條件。采用邊界附近棱柱網格、內部四面體網格剖分計算域，包含900 744 個四面體網格，流動處于湍流狀態，采用穩態求解器、SST（Shear Stress Transport）k-ω湍流模型進行該算例的計算。該算例可以驗證本文算法在四面體網格、復雜流動狀態下流量分配和出口回流的處理能力。

圖5 繪制了α1∶α2在0∶1.0～0.5∶0.5 范圍內變化時Y 形管軸向剖面上的速度分布，計算得到的出口流量和設定出口流量的比較結果如表2 所示。在上分支流量分配系數設置為0 時，計算得到的上分支出口接近于0。隨著設定的上分支流量分配系數增加，相應流量逐漸增加，最終上下分支的速度呈對稱分布。值得注意的是，在該算例設置下，上下分支起始段會出現一個回流區，其大小隨著流量分配系數的增加而不斷擴展。出口流量統計結果表明，本算例最大誤差出現在α1∶α2=0∶1.0 這種極端流量分配情況下，最大相對誤差在0.001%量級，說明本文發展的方法同樣適用于典型三維四面體網格、湍流流動情形。

表2 60°Y 形分叉管算例流量統計結果Table 2 Mass flow rate statistics results of duct case with 60° Y-branch

圖5 60°Y 形分叉管算例不同流量分配下速度分布Fig.5 Velocity distribution of duct case with 60° Y-branch at various mass flow ratios

2.2 計算域縮減結果一致性驗證

本文出口邊界處理方法的一個應用目標是，數值仿真中可以根據實際需要將出口邊界置于任意位置，出口流動未充分發展、甚至存在回流的情形下數值計算仍然能夠滿足正確和穩定收斂的要求。為了進一步驗證這一特性，研究中針對上文提到的90°T 形分叉管算例開展了系列數值實驗進行驗證，算例如表3所列。

表3 計算域縮減影響算例和取樣位置設置Table 3 Test cases and sampling positions for truncation of computational domain

針對圖2 所示的構型，在表3 中，算例1 中設置Lm=10W，Ls=10W。這種情形下主分支出口和右側分支出口計算域均足夠長，滿足充分發展的假設，采用常規的出流邊界進行計算，可以獲得充分發展邊界條件下流場解，并得到2 個出口上的質量流量，分別為α1和α2。在此基礎上，算例2 采用本文方法處理出口邊界條件，流量分配系數保持和算例1 一致，和算例1 的對比結果用于驗證本文方法在計算域足夠長的情況下和常規出流邊界條件的一致性；在算例2 的基礎上，進一步縮減右側分支長度，Ls/W分別取5、3、2、1、0.5，得到算例3～算例7，使得右側分支出口處流動依次從充分發展狀態向回流狀態變化，流量分配系數仍保持和算例1 一致。算例2～算例7和算例1 的對比結果，可以驗證本文方法在充分發展狀態、未充分發展狀態以及存在回流狀態下的適用性。

圖6展示了表3中x= 1.0， 1.5， 3.0， 6.0 m的采樣位置，用于定量分析出口邊界位置對計算域內部場計算結果的影響。表3 中取樣位置所在列“√”表示該算例設置下，當前取樣位置在計算域設置范圍之內，可進行取樣比對。例如算例1中，側分支長度Ls= 10W，4 個取樣位置均位于算例計算域內。但在算例6 中，側分支長度Ls=W，僅有x= 1.0， 1.5 m 這2 個有效取樣位置，而x= 3.0， 6.0 m 這2 個位置因超出計算域不再對比分析。

圖6 取樣位置示意圖（x = 1.0， 1.5， 3.0， 6.0 m）Fig.6 Schematic of sampling locations （x = 1.0， 1.5， 3.0， 6.0 m）

需要特別說明的是，在上述算例設置下，流動的物理特征由流量分配系數保證，在計算域改變時，2 個分支流量分配、物理特性等并未發生變化，重點考察算法能否保證不同計算域選擇下結果的一致性。

圖7 展示了算例1～算例7 的速度分布云圖，為了方便比對，圖中僅截選了實際計算域的一部分。圖7（a）和圖7（b）分別為Ls= 10W構型設置下充分發展邊界條件和本文方法中的結果，出口流量分配系數為α1=0.59，α2=0.41。可以看出，在出口流動達到充分發展情況下，本文方法預測的速度場和常規出流邊界預測的速度場是定性一致的。圖7（c）～圖7（g）為Ls/W=5， 3， 2， 1， 0.5 的結果，在右側分支不斷縮減的情況下，入口段、主出口分支以及右側分支部分的速度分布均和圖7（a）中的結果是一致的，即便是Ls/W= 0.5 構型設置、出口處存在強回流情形下，計算中沒有出現數值振蕩或流場形態的顯著變化。

圖8 進一步顯示了算例1～算例7 在不同取樣位置上的速度大小分布，各算例取樣設置如表3 和圖6 所示。從圖中可以看出，Ls/W= 10構型設置、不同取樣位置處，充分發展邊界條件下的模擬結果和使用本文自由出口邊界條件算法得到的速度大小分布基本一致，說明本文方法可以在更小的計算域復現出常規充分發展出流邊界條件的模擬結果，即模擬結果與計算域出口邊界位置的選擇無關。

圖8 中x= 1.5 m 和x= 3.0 m 取樣處的速度分布顯示y= 0.2～0.4 m 之間存在一個速度零點，結合圖7 中的速度云圖分布，可以看出該速度零點實際上是法向速度發生方向變化的臨界點，即該區域存在回流區。圖9 繪制了Lm/W= 10、Ls/W= 10 構型設置下2 個分支交叉點附近的流線分布，該流線圖明確顯示出主出口分支和側出口分支均存在回流區。Ls/W= 0.5， 1， 2幾種構型設置及預測結果說明了本文的出口邊界處理方法同樣適用于出口邊界存在回流的情況。

圖7 不同計算域截斷情況下速度場云圖Fig.7 Velocity contours of duct case with successively truncated computational domain

圖8 不同計算域截斷程度下典型取樣位置的速度大小分布（Lm/W = 10）Fig.8 Velocity profiles of duct case with successively truncated computational domains at various sampling locations （Lm/W = 10）

圖9 90°T形管拐角處流線分布（Lm/W = 10， Ls/W = 10）Fig.9 Streamlines distribution at corner of duct case with 90° T-branch （Lm/W = 10， Ls/W = 10）

2.3 流動雷諾數影響驗證

Dong 和Shen［7］指出，雷諾數（Re）是影響出口邊界條件數值穩定性的重要因素。為驗證本方法在不同Re下的適用情況，研究中開展了Re=200～ 50 000 范圍內的流場模擬數值實驗，其中Re通過修改入口速度uin實現，入口速度均滿足uin=Cx(1-x)的分布形式，C取值4、20、50、100、200、500、1 000，依次對應于Re為200、1 000、 2 500、5 000、10 000、25 000、50 000 的情況。構型參數選用Ls/W= 10 和Ls/W= 1 這2 種，用于表征不同出口位置選取情況。算例邊界條件設置和上述T 形管保持一致，出口流量分配系數固定為0.6∶0.4。需要指出的是，該系列算例中流動特征變化較大，覆蓋了Re= 200 的層流狀態到Re= 50 000 的旺盛湍流狀態，這里選用匹配于較寬Re范圍的SSTk-ω湍流模型。

數值實驗結果顯示，本文提出的出口邊界方法在Re= 200～50 000 范圍內均表現出良好的數值收斂性，所有算例未出現數值發散的情況。圖10 顯示了Re= 25 000 情況下2 種構型對應的速度分布云圖，相比于圖7 中的結果，高Re下的2 個分支回流區顯著擴展，這也是相同構型、采用常規出流邊界條件時高Re下更容易出現數值不穩定的主要原因。圖11 繪制了不同Re下2 種構型x= 1.0 m 處的速度剖面，其中實線為Ls=10W構型結果，實心方塊為Ls=W構型結果，圖中為了方便展示不同Re結果，縱軸采用對數坐標系。從圖中可以看出，無論是低Re的層流狀態還是高Re的旺盛湍流狀態，仿真設置的計算域大小均不會影響速度場結果。

圖10 Re = 25 000 時2 種構型下速度云圖Fig.10 Velocity contours for two configurations with Re = 25 000

圖11 不同Re 下2 種構型x = 1.0 m 處速度分布Fig.11 Velocity profiles at x = 1.0 m for two configurations with various Re

2.4 TECFLAM 復雜旋流算例驗證

將1.2 節提出的方法應用于德國TECFLAM旋流燃燒器冷態流場［17-18］的模擬，檢驗該方法在復雜構型流動特征中的適用性。該燃燒室屬于強旋流動，有豐富的冷態流動、甲烷-空氣預混燃燒實驗數據，被廣泛用于數值仿真的結果驗證［17-18］，其幾何構型和網格剖分結果如圖12 所示。TECFLAM 旋流燃燒器主體結構是環形通道包圍著一柱狀鈍體，鈍體直徑為30 mm，環形通道（內徑30 mm，外徑60 mm）通空氣或空氣-甲烷混合物，旋流由環形通道上游一可移動塊產生，旋流數可在 0～2.0 范圍內連續調節，默認旋流數為0.75。整個燃燒器置于同軸空氣伴流中，伴流速度為0.5 m/s。

圖12 TECFLAM 旋流燃燒器構型和網格剖分Fig.12 Configuration and meshing of TECFLAM swirl burner

計算域選定為燃燒器出口上部空間，出口處高度記為z=0 mm。燃燒室內產生的旋流等效為燃燒室出口處的速度邊界條件，即u=-5.23y/r，v=5.23x/r，w=5.0 m/s，其中r為指定點距離中心軸線的徑向距離。出口邊界位于z=H處，H分別取值120 mm 和360 mm，計算域剖分為45 萬網格和110 萬網格。前者計算域覆蓋了實驗和數值研究的重點關注區域（z= 1～90 mm），但出口處存在較強的徑向和周向運動，無法直接應用常規出口邊界條件；后者出口位于主流區下游較遠處，適用于傳統邊界處理方法，但網格量是前者的2.4 倍。

采用大渦模擬方法、動態Smagorinsky 亞格子應力封閉模型，基于上述2 種計算域模擬TECFLAM 冷 態 流 場，將z= 1， 10， 20， 30， 60 mm 處時均速度3 個分量的預測值與實驗結果比對，結果如圖13 所示。其中，實驗值取自文獻數據［17-18］，設置1 對應于H= 120 mm 的計算域，設置2 對應于H= 360 mm 的計算域。在H= 360 mm 時，速度分量的分布說明了計算域中的流動特征：氣體從噴口噴出后，切向分量和軸向分量相當（～ 5 m/s），主導整體運動。沿著徑向方向，切向速度分量先增加后減小，并在某個徑向位置處出現最大切向速度，構成旋渦區。隨著流動的發展，旋渦區域逐步沿著徑向方向發展（～ 1 m/s）。圖13 中設置2 和實驗值結果表明，不同高度處速度3 個分量的預測值與實驗值符合得較好。當H從360 mm 減小到120 mm 時，計算域大幅縮減，本文方法能夠保證穩定收斂。從對比結果來看，在z= 1～60 mm 范圍內，除了部分流場變化劇烈區域速度預測值出現一定的偏差，速度3 個分量的預測值仍然和實驗數據保持較高的一致性。

圖13 TECFLAM 旋流燃燒器模擬的時均速度分布與實驗結果［17-18］對比Fig.13 Comparison of time-averaged velocity distribution in TECFLAM swirl burner simulation with experimental［17-18］ data

3 結 論

本文發展了一種基于質量流量分配系數的多出口計算域邊界條件處理方法，實現了不可壓縮流動模擬中多出口邊界、存在回流情形的高效、穩定模擬，并基于自主研發的高性能計算流體數值模擬軟件驗證了該方法的有效性，主要結論如下：

1）針對典型六面體網格、四面體網格算例，該方法實現了出口質量流量分配系數在0～1.0范圍內的穩定模擬，計算得到的質量流量分配系數和設定值最大相對偏差在0.001%量級。

2）在保持出口流量分配系數不變的情況下，該方法正確預測了出口流動未達到充分發展、存在回流或旋流下的流場，出口邊界位置對計算結果的影響很小；在計算域出口流動達到充分發展狀態時，其預測結果和常規充分發展邊界條件的結果一致。

3）通過2 種典型構型設置、系列Re數值實驗，驗證了該方法適用于Re處于200～50 000 范圍、出口流場未充分發展或存在回流的情形，覆蓋了低Re層流到高Re旺盛湍流。

4）使用TECFLAM 旋流燃燒器冷態流場大渦模擬驗證自由出口邊界條件算法，結果表明計算區域縮減至1/3 后模擬仍能保證穩定收斂，預測值和實驗數據保持較高的一致性。相比于使用常規出口邊界的方法，本文的自由出口邊界條件算法支持仿真人員根據需要選擇計算域，可顯著減小數值計算開銷。