關清帝，梁劍寒，張林，陳文武，陳玉俏

國防科技大學 空天科學學院 臨近空間技術研究所，長沙 410073

超燃沖壓發(fā)動機燃燒室中燃燒過程十分復雜，湍流、激波、火焰相互耦合，火焰結(jié)構一般高度褶皺、破碎，通常還伴隨局部熄火、再燃、閃回等復雜過程［1］。在超聲速條件下，由于流動速度極快，湍流流動與化學反應的特征時間尺度處于相當?shù)牧考?，湍流與化學反應緊密耦合在一起，相互作用十分劇烈。超聲速、強非定常效應等極端條件帶來了一系列值得深入研究的可壓、非平衡復雜流動燃燒行為［2］，給超聲速湍流燃燒模擬帶來了巨大的挑戰(zhàn)。又由于化學反應源項的強非線性，對大渦模擬（Large Eddy Simulation， LES）方法而言，很難對濾波后的化學反應源項進行準確建模，湍流燃燒模型是LES 準確模擬超聲速湍流燃燒的關鍵。

目前，超聲速流中比較常用的湍流燃燒模型包括火焰面進度變量模型［3］、部分攪拌反應器模型［4］、線性渦模型和一維湍流模型［5-6］以及輸運型概率密度函數(shù)（Probability Density Function， PDF）模型［7］等。在這些模型中，輸運型PDF 方法可以精確處理湍流與化學反應的相互作用，適用于任何復雜的湍流燃燒模式，因此被認為是最可信的湍流燃燒模型之一［8］。輸運型PDF 方法從湍流的隨機特性出發(fā)，基于Navier-Stokes 方程建立并求解湍流物理量的單點聯(lián)合概率密度輸運方程。湍流物理量最基本的可預測性是它的概率及其概率密度分布。已知隨機變量的概率密度分布，則可以求解包括化學反應源項在內(nèi)的任意復雜和任意高階的單點矩，因此PDF 方法理論上可以精確地處理湍流化學反應相互作用。

PDF 方法體系的建立，主要來源于Pope［7］的開創(chuàng)性工作，其對PDF 方法的理論基礎、?；椒ㄒ约扒蠼膺^程都進行了完整的論述。Givi［9］首先將PDF 方法的思想應用于LES 濾波化學反應源 項的 封 閉，Pope［10］隨 后 提出了濾 波PDF的概念，Gao 和O'brien［11］依此推導并建立了濾波PDF的輸運方程。此后，國內(nèi)外學者在不可壓流中開展了大量的LES-PDF數(shù)值模擬。Jaishree［12］采用LES-PDF 方法對Sandia Flame D、E 和F 進行了數(shù)值模擬，與實驗結(jié)果對比，不論是平均量還是脈動量，均與實驗結(jié)果符合較好。楊越等［13］采用LES-PDF 方法對非預混CO/H2射流火焰進行了模擬，與DNS（Direct Numerical Simulation）結(jié)果進行了對比驗證，結(jié)果顯示溫度以及主要組分的質(zhì)量分數(shù)與DNS 結(jié)果符合很好，LESPDF 方法較好地預測了火焰的局部熄火與再燃過程?？傮w上，PDF 方法在不可壓湍流燃燒中取得了很好的效果。

由于在模擬湍流燃燒相互作用方面的潛在優(yōu)勢，PDF 方法也開始逐步向超聲速流發(fā)展。Jaberi 等［14］首次將可壓縮LES 求解器與PDF 相結(jié)合，發(fā)展了可壓縮LES-PDF 方法，在二維無反應時間混合層算例中初步驗證了PDF 方法在可壓縮流中的可行性。Banaeizadeh 等［15］進一步發(fā)展了這種可壓縮LES-PDF 方法，在激波管和無反應超聲速同軸射流中開展了相應的數(shù)值測試。最近，Komperda 等［16］將可壓縮PDF 與間斷譜元方法相結(jié)合，對球形域中的爆炸過程進行了數(shù)值模擬。在作者團隊先前的工作［17-18］中，發(fā)展了一種守恒可壓縮LES-PDF 方法。該方法在激波管和時間反應混合層中得到了初步驗證，結(jié)果表明守恒可壓縮PDF 方法顯著提高了能量解的精度。這些工作在可壓縮LES-PDF 方法建模等方面取得了重要的進展，但是目前的大多數(shù)模型僅在無反應流或者簡單反應流中進行了驗證。超燃沖壓發(fā)動機燃燒室中燃燒過程十分復雜，很有必要在實際構型中開展模型驗證工作。

隨著可壓縮PDF 模型的不斷完善和數(shù)值方法的發(fā)展，最近的一些研究也開始將PDF 方法應用于比較實際的可壓縮湍流燃燒問題。Validi等［19］對快速壓縮機中的湍流射流火焰進行了LES-PDF 模擬，定性比較了溫度場計算結(jié)果。Ranadive［20］模 擬 了Cheng 等［21］的超聲速射流火焰，與實驗和層流有限速率模型結(jié)果進行了比較詳細的定量對比。最近，Abdulrahman 等［22］使用LES-PDF 方法模擬了非預混和預混的分布式燃燒過程，PDF 預測的軸向平均速度與實驗符合較好。雖然上述部分工作開始應用于比較實際的超聲速燃燒問題，但是采用的網(wǎng)格仍為簡單笛卡爾網(wǎng)格，該網(wǎng)格及其數(shù)值求解框架難以應用到實際的超燃沖壓發(fā)動機構型，極大地限制了LESPDF 方法在實際發(fā)動機燃燒室中的應用。

總體來說，目前PDF 方法在超聲速流中的發(fā)展還十分初步，大多數(shù)可壓縮LES-PDF 方法只在簡單笛卡爾網(wǎng)格下進行了方法驗證，可壓縮LES-PDF 方法在實際構型中的研究和應用還十分缺乏。為了將該方法應用于發(fā)動機燃燒室模擬，需要發(fā)展適用于復雜構型的LES-PDF 方法。同時復雜構型下LES-PDF 耦合方法，粒子追蹤方法以及大規(guī)模并行計算等方面，都還有待深入的研究。本文在前期可壓縮LES-PDF 工作的基礎上，基于多塊結(jié)構網(wǎng)格，實現(xiàn)了基于曲線坐標的粒子追蹤方法，建立了曲線坐標LES-PDF 方法，使可壓縮LES-PDF 方法適用于復雜構型計算。

1 基本控制方程

1.1 可壓縮LES 控制方程

LES 方法通過對湍流運動的濾波，將湍流脈動分解為可解的大尺度運動和需要封閉的小尺度脈動。可解尺度湍流運動直接求解，小尺度湍流脈動的質(zhì)量、動量和能量輸運對大尺度運動的作用采用亞格子模型進行封閉。對Navier-Stokes 方程進行Favre 濾波，得到大渦模擬控制方程為

式中：t為物理時間；x為物理坐標；i為空間維度；ui為i方向的速度分量；ρ為密度；p為流場壓力；δij為Kronecker 函數(shù)；τij為黏性應力張量；E為總能；qi為i方向熱通量；Yα、Vi，α和ω?α分別為組分α的質(zhì)量分數(shù)、i方向擴散速度以及質(zhì)量生成率；“-”表示空間濾波；“～”表示Favre 濾波；上標“sgs”表示未封閉的亞格子項。其中τsgsij為亞格子黏性應力，采用一方程Yoshizawa 模型［23］進行封閉：

其中：ksgs為亞格子湍動能；νt為渦黏系數(shù)，由渦黏輸運方程直接求解；Sij為應變率張量。

Hsgsi為亞格子焓通量，Ysgsi，α為亞格子組分對流通量，這2 項均采用梯度擴散假設進行封閉：

式中：H為總焓；Prt和Sct分別為湍流普朗特數(shù)和湍流施密特數(shù)，在本文的模擬中，均為0.72。

1.2 曲線坐標PDF 輸運方程

本文僅考慮標量PDF，在標量PDF 方法中，直接求解熱力學標量?的聯(lián)合概率密度函數(shù)輸運方程，其中?包括Ns個組分的質(zhì)量分數(shù)以及顯焓h。可壓縮流中標量PDF 定義為

式中：?和ψ分別表示標量向量和標量向量的樣本空間為笛卡爾物理空間坐標；y為x對應的積分變量；δ（·）為細粒密度函數(shù)；G（·）為濾波核函數(shù)。由標量守恒方程，可以直接推導建立笛卡爾坐標下的PDF 輸運方程：

式中：為可解速度；Γ為擴散系數(shù)，為分子擴散系數(shù)和湍流擴散系數(shù)之和；?α為組分α的質(zhì)量分數(shù)；Sα為組分α的化學反應源項；Ωm為標量混合頻率；ψα為?α對應的樣本空間。式（9）等號左邊的2 項為PDF 瞬時項和對流項。式（9）等號右側(cè)第1 項代表湍流脈動和分子擴散導致PDF 在物理空間的輸運；第2 項為分子擴散項引起的PDF在組分空間的輸運，采用最常用的IEM（Interaction by Exchange with the Mean）小尺度混合模型進行封閉；第3 項表示化學反應以及可壓縮項引起的PDF 在組分空間的輸運。

式（9）給出的輸運方程中化學反應源項Sα可以展開為

式中：ω?α為 組分α的質(zhì)量生成率；α=Ns+1 時，Sα表示顯焓方程的源項，為反應釋熱ω?T和可壓縮高速源項Sh之和；Sh的表達式為

由于坐標變換只涉及到物理空間，所以本文重點關注跟物理空間坐標相關聯(lián)的項，忽略對可壓縮項和反應源項等單點源項過程的處理。對于超聲速流中反應源項和可壓縮項的建?？梢娮髡邎F隊前期超聲速PDF 方法的建模工作［18］。

對式（9）兩端同時除以J，J的表達式為

式中：J為坐標變換系數(shù)矩陣對應的行列式值；ξx、ξy、ξz、ηx、ηy、ηz、ζx、ζy、ζz為度量系數(shù)。利用一般曲線 坐 標 系(ξ，η，ζ，τ)和笛卡爾直角坐標系(x，y，z，t)之間的轉(zhuǎn)換關系，進行坐標變換。

為簡潔起見，這里直接給出推導整理后的一般曲線坐標系下的PDF 輸運方程：

式中：ξ=[ξ，η，ζ]為一般曲線坐標系下的空間坐標；Fξ?=Fx?/J表示一般曲線坐標系下的標量概率密度函數(shù)。為逆變速度，其表達式為

其中：u、v、w為物理空間的速度分量。

對流項系數(shù)Ai為

可以看出，坐標變換只涉及物理空間，與標量空間無關。由此得到了完全封閉的曲線坐標下PDF 輸運方程，可以進行數(shù)值求解。

2 數(shù)值求解方法

2.1 曲線坐標系下LES 離散方法

采用有限差分方法對曲線坐標下LES 控制方程進行數(shù)值離散，曲線坐標下控制方程的具體形式為

式中：Q為守恒變量矢量，表達式為

Ec、Fc、Gc為對流通量矢量；Ev、Fv、Gv為黏性通量矢量。式（17）中各項形式較為復雜，這里不再詳細展開，具體表達可以參見流體力學教材。

本文采用的可壓縮LES 求解器在超聲速湍流流動、燃燒模擬中已經(jīng)進行了大量的數(shù)值驗證［24-25］。時間方向的離散采用具有TVD（Total Variation Distance）性 質(zhì) 的 三 階Runge-Kutta 方法?？臻g上對流項和黏性項的離散分開進行，對流項具有很強的非線性，采用五階WENO［26］格式離散，黏性項體現(xiàn)為物理過程中的耗散，沒有明顯的方向性，采用簡單的二階中心差分格式進行離散。

2.2 曲線坐標系下Monte Carlo 粒子方法

PDF 輸運方程的維數(shù)很高，用傳統(tǒng)的有限差分等方法求解計算量過大，一般采用Monte Carlo 粒子方法求解。在粒子方法中，流體由大量的拉格朗日粒子來描述，每個粒子根據(jù)隨機微分方程隨時間演化。為了得到對應的曲線坐標下的隨機微分方程，引入Fokker-Planck 的矢量形式。對多維擴散過程，其轉(zhuǎn)移密度存在，且滿足Fokker-Planck 方程：

式中：f為轉(zhuǎn)移密度函數(shù)；d'為轉(zhuǎn)移密度函數(shù)的維度；ys為狀態(tài)空間。式（19）右端第1 項Bij為擴散項，第2 項a為漂移項。其對應的擴散過程為

其中：Ui為隨機變量；U為隨機變量的樣本空間；ai[U(t)，t]稱為漂移速度，而bij[U(t)，t]稱為擴散速度，兩者均為隨機變量的函數(shù)；W表示維納隨機過程。

此外，民立中學歷來也有演劇祝圣的活動，直到1908年8月，其演劇祝旦活動受到王安民的指責，從此不再舉辦演劇?！缎聞∈贰吩敿氂涊d了這一事件：

結(jié)合式（13）所示的曲線坐標下PDF 輸運方程，可以得到對應的粒子隨機微分方程。這里直接給出曲線坐標下ξ方向粒子的隨機微分方程為

式中：上標“*”表示粒子所在位置的物理參數(shù)；U為ξ方向粒子的逆變速度分量；dτ為曲線坐標下的時間增量；dWi代表維納隨機過程［27］的增量。粒子在逆變速度、分子輸運以及亞格子湍流擴散作用下在計算空間移動，而粒子標量在小尺度混合、化學反應以及可壓縮高速源項的作用下發(fā)生變化。粒子隨機微分方程所對應的Fokker-Planck 方程與?；蟮臉肆縋DF 輸運方程形式一致，通過概率相似體系建立起了Monte Carlo粒子和真實流體運動的關聯(lián)。

粒子位置以及標量的演化通過王海峰等［28］發(fā)展的弱二階算子分裂格式進行。粒子隨機微分方程可以分裂為3 個具體過程，分別是計算空間的輸運（T），標量空間的混合（C）以及標量空間的單點化學反應（R）。該二階精度的算子分裂格式稱為“TCRCT”，依次包含以下5個過程：物理空間輸運半時間步，標量空間半時間步的亞格子混合，整個時間步長的化學反應，另外半個時間步的亞格子混合以及位置空間的輸運。上述分裂格式的優(yōu)勢是可以保證化學反應步不分裂，因為在所有過程中，化學反應計算占用時間最長，尤其在采用詳細化學反應機理時，化學反應計算開銷更大。

2.3 曲線坐標系下LES-PDF 耦合方法

LES 有限差分求解器和PDF 粒子求解器需要耦合求解，耦合過程如圖1 所示。LES 有限差分求解器求解大渦模擬流動和組分方程，計算得到LES“可解場”的分布。PDF 粒子求解器計算粒子隨機微分方程，并通過系綜平均得到PDF“粒子場”平均量。LES 向PDF 提供濾波速度、擴散系數(shù)、可壓縮源項等物理量用以PDF 粒子輸運，PDF 則向LES 反饋化學反應源項以完全封閉大渦模擬組分方程。

圖1 LES-PDF 方法耦合示意圖Fig. 1 Schematic diagram of LES-PDF coupling

上述耦合過程通過混合拉格朗日粒子/歐拉網(wǎng)格方法實現(xiàn)，該方法如圖2 所示，每個網(wǎng)格中分布若干PDF 粒子。LES 向PDF 粒子傳遞的濾波速度等物理量，由網(wǎng)格節(jié)點向粒子位置處進行插值，得到該粒子處的具體值。在本文中，為簡便起見，采用了線性插值方法。從PDF 粒子傳給LES 求解器的值，則由選定區(qū)域中粒子的系綜平均得到：

圖2 混合拉格朗日粒子/歐拉網(wǎng)格方法示意圖Fig. 2 Schematic diagram of hybrid Lagrangian particle/Eulerian mesh method

式中：Np為網(wǎng)格中粒子總數(shù)；ω()n為第n個粒子的權重；?(n)為第n個粒子的標量值；Q為粒子標量的單點統(tǒng)計。為了與插值過程相匹配，本文選定的系綜區(qū)域為單個網(wǎng)格，即采用PIC（Particles In Cell）系綜方法。粒子在計算空間輸運過程中，由于PDF方法的隨機效應，可能出現(xiàn)部分網(wǎng)格內(nèi)粒子數(shù)過少而部分網(wǎng)格粒子數(shù)過多的情況，這會明顯增加粒子系綜統(tǒng)計誤差，且導致計算負載的不平衡。本文通過設置平均粒子數(shù)目對網(wǎng)格中粒子的數(shù)量進行限制，保證粒子數(shù)密度在倍網(wǎng)格平均粒子數(shù)之間。本文的流動混合算例中，網(wǎng)格平均粒子數(shù)目設置為50，燃燒算例中，由于粒子組分較多且化學反應計算量較大，網(wǎng)格平均粒子數(shù)設置為30。需要說明的是，在曲線坐標框架下，基于計算坐標，便于構造高階的插值和系綜方法。在超聲速流中，流場速度等物理量變化大，高階的速度插值和系綜統(tǒng)計方法可能有助于提高計算精度。

綜上，LES-PDF 方法的基本算法流程可以概括如下：首先對LES 流場進行初始化，PDF 粒子的初始質(zhì)量、能量（顯焓）以及組分等參數(shù)由LES 流場初始化參數(shù)計算得到。然后開始計算LES 質(zhì)量、動量和能量方程，求出LES 密度、速度、溫度和壓力，在此基礎上計算得到高速源項和湍流擴散系數(shù)等需要傳遞給PDF 的流場信息。下一步將濾波速度、高速源項以及湍流擴散系數(shù)等插值到PDF粒子所在位置處。之后，PDF 粒子根據(jù)式（21）和式（22），進行位置空間的粒子追蹤和標量空間的混合和化學反應，粒子的位置和標量參數(shù)發(fā)生變化。最后，根據(jù)網(wǎng)格內(nèi)所有粒子的系綜平均求出PDF粒子場統(tǒng)計結(jié)果，并將粒子場計算的化學反應源項傳遞給LES 求解器，求解得出新的LES 可解場信息，如此完成一個時間步長的耦合求解過程。本文提出的曲線坐標LES-PDF 方法，主要對粒子位置空間的追蹤過程進行了改進，提出了基于一般曲線坐標系的粒子追蹤方法。

2.4 曲線坐標系下粒子追蹤方法

本文提出的一般曲線坐標下LES-PDF 方法中，采用多塊結(jié)構網(wǎng)格。多塊網(wǎng)格由局部結(jié)構網(wǎng)格組成，局部塊結(jié)構網(wǎng)格相互拼接形成一個全局網(wǎng)格。曲線坐標下的粒子追蹤方法也是基于多塊結(jié)構網(wǎng)格進行實現(xiàn)。

傳統(tǒng)方法基于粒子物理速度(u，v，w)在位置空間進行移動，如圖3（a）所示。物理坐標下粒子追蹤過程需要判斷粒子和網(wǎng)格之間的相對位置關系，即粒子與網(wǎng)格的4 條邊進行矢量計算才能判斷是否位于當前網(wǎng)格中。同時從A到B的移動過程中，還要進行搜索以確定粒子的移動路徑。PDF 方法中，網(wǎng)格中包含大量粒子，每一時間步均需要進行搜索和定位，在物理坐標系下進行追蹤相對繁瑣和耗時。

曲線坐標下的追蹤方法，即通過坐標變換，將傳統(tǒng)在物理坐標下的粒子追蹤過程，轉(zhuǎn)換到曲線坐標系下進行。如圖3（b）所示，粒子在計算空間根據(jù)逆變速度(U，V，W)移動。曲線坐標下粒子追蹤的一個主要特點是，在單塊網(wǎng)格中，粒子的計算坐標(ξP，ηP，ζP)與網(wǎng)格單元索引(Ic，Jc，Kc)直接相關：

圖3 單塊網(wǎng)格內(nèi)曲線坐標粒子追蹤方法Fig. 3 Curvilinear particle tracking method within a single block grid

式中：floor 表示向下取整。粒子所在位置的曲線坐標向下取整即為粒子所在網(wǎng)格的網(wǎng)格索引，從A向B的移動中，粒子不需要進行額外的搜索和定位。

如何實現(xiàn)粒子在多個網(wǎng)格分區(qū)間的連續(xù)追蹤，是多塊結(jié)構網(wǎng)格下曲線坐標粒子追蹤方法需要處理的一個關鍵問題。本文提出了一種處理多塊網(wǎng)格分區(qū)間粒子曲線坐標轉(zhuǎn)換的方法，實現(xiàn)了粒子在多塊分區(qū)間曲線坐標的連續(xù)追蹤，其基本過程如圖4 所示。

圖4 曲線坐標粒子追蹤方法流程圖Fig. 4 Flow chart of curvilinear particle tracking method

首先計算粒子到網(wǎng)格分區(qū)邊界的最小時間Δtmin，判斷粒子是否移動出當前網(wǎng)格分區(qū)。若Δtmin≥Δt，則粒子仍在當前分區(qū)中，在當前網(wǎng)格分區(qū)中的粒子不需要追蹤，直接更新粒子的曲線坐標，若當前的迭代步數(shù)n小于給定的最大迭代計算步nmax，進入下一個迭代，否則結(jié)束計算。若Δtmin＜Δt則表明粒子將要移出當前分區(qū)。對跨越網(wǎng)格分區(qū)的粒子需要進行額外的處理。

圖5給出了跨越多塊網(wǎng)格分區(qū)粒子的追蹤過程。對要移動出當前網(wǎng)格分區(qū)的粒子，首先確定粒子穿越網(wǎng)格的位置C，將粒子先輸運至分區(qū)邊界處，并在邊界處更新輸運時間Δt=Δt-Δtmin。然后在邊界C處，根據(jù)邊界兩側(cè)分區(qū)的拓撲關系，進行曲線坐標間的坐標轉(zhuǎn)換，將當前分區(qū)1 下的曲線坐標轉(zhuǎn)換為目標分區(qū)2 中的曲線坐標。除了位置坐標需要轉(zhuǎn)換外，粒子的逆變速度也要進行相應的轉(zhuǎn)換。位置矢量和逆變速度矢量的坐標轉(zhuǎn)換在數(shù)值上需要乘以坐標轉(zhuǎn)換矩陣完成，其中的坐標轉(zhuǎn)換矩陣即表征了2 個分區(qū)間的拓撲關系。最后將該粒子進行非阻塞的點對點通信，從當前分區(qū)1 發(fā)送給目標分區(qū)2，并在網(wǎng)格分區(qū)2 中從C處開始繼續(xù)對該粒子進行追蹤，重復前面的過程，直到剩余的時間步長Δt變?yōu)?。

圖5 多塊網(wǎng)格間曲線坐標粒子追蹤方法Fig. 5 Curvilinear particle tracking method between multiblock meshes

基于上述粒子追蹤方法，可以實現(xiàn)全局范圍內(nèi)粒子的連續(xù)追蹤。實際計算中，只有少部分邊界附近的粒子會穿越網(wǎng)格分區(qū)，本文提出的曲線坐標粒子追蹤方法，單塊網(wǎng)格中不需要額外追蹤，減少了單塊網(wǎng)格內(nèi)粒子追蹤的計算量，同時也很好地處理了多個網(wǎng)格分區(qū)結(jié)合位置處的粒子追蹤過程，十分適合拉格朗日粒子方法的大規(guī)模并行計算。此外，為了保證計算空間中粒子追蹤的精度，當時間步長過大，粒子穿越多個網(wǎng)格，即時，將時間步長減半，保留剩余的輸運時間，重新計算粒子位置，直到粒子剩余的時間步長為0。

3 算例驗證

為了驗證曲線坐標下LES-PDF 方法的準確性，針對曲線坐標下的關鍵問題，曲線坐標形式的PDF 輸運方程及其粒子隨機微分方程，多塊分區(qū)間的粒子坐標轉(zhuǎn)換，通信和連續(xù)追蹤，以及曲線坐標程序復雜構型的處理能力等問題，從簡單到復雜進行逐步的數(shù)值驗證。

3.1 拉伸網(wǎng)格下Sod 激波管

首先針對一維Sod 激波管問題進行驗證。無量綱計算域取為［0，1］，x方向共布置了101 個節(jié)點，第1 個和最后1 個網(wǎng)格節(jié)點為壁面邊界條件。采用的網(wǎng)格為拉伸網(wǎng)格，在激波管中間位置進行了指數(shù)加密，激波管中間位置網(wǎng)格無量綱寬度約為0.005，激波管第1 個和最后1 個網(wǎng)格的無量綱寬度為0.025。計算采用的CFL（Courant-Friedrichs-Lewy）數(shù)為0.5，無量綱初始條件設為

在該算例中只考慮無黏流，標量PDF 只求解了顯焓，主要通過溫度和密度等熱力學結(jié)果來驗證方法的準確性。初始的顯焓流場由式（25）給定的初始壓力、速度以及密度計算得到。在簡單拉伸網(wǎng)格下重點驗證曲線坐標下可壓縮模型以及顯焓方程求解的正確性，排除了其他模型和方程的干擾。圖6 分別給出了曲線坐標LES-PDF方法計算的無量綱溫度和密度結(jié)果，從圖中可以看出，PDF 計算結(jié)果與精確解符合得很好，初步驗證了曲線坐標下PDF 輸運方程和粒子隨機微分方程的準確性。

圖6 無量綱溫度和密度計算結(jié)果Fig. 6 Calculation results of dimensionless temperature and density

3.2 非正交網(wǎng)格下時間混合層

隨時間發(fā)展的平面混合層流動示意圖如圖7所示，由2 股相反的平行流發(fā)展而成。平均流速度的初值為雙曲正切分布uˉ=tanhy，擾動流函數(shù)φ取為

圖7 時間混合層示意圖Fig. 7 Schematic diagram of temporal mixing layer

對應的最不穩(wěn)定擾動波數(shù)β=0.444 6，擾動波長λ=2πβ=14.132 2。由此得到初始速度場為

式中：擾動量級ε取為0.1；無量綱初始壓力為p=1/(γMa2c)，Mac為對流馬赫數(shù)，本文中為0.4，γ=1.4為比熱比。初始溫度T由Busemann-Crocco 關系確定，T=1+(γ-1)Mac2(1-uˉ2)/2。該算例中除顯焓外，PDF 變量中還包含一個被動標量?，主要研究該被動標量的輸運和混合，對曲線坐標LES-PDF 方法進行定性和定量的驗證。初始時刻上層標量?1=1，下層標量為?2=0，標量剖面由雙曲正切函數(shù)給定：

式中：δω為初始渦量厚度。選取初始渦量厚度δω為特征長度，基于此計算的雷諾數(shù)約為 2 800。

本文考慮的為雙周期混合問題，所以x方向計算域包含2 個最不穩(wěn)定波長，y方向計算域長度與x方向相同，即Lx=Ly=2λ。x方向為周期邊界條件，y方向為出口邊界。模擬中網(wǎng)格不分區(qū)，排除多塊網(wǎng)格分區(qū)間粒子追蹤的干擾，只驗證曲線坐標方程以及單塊網(wǎng)格內(nèi)粒子追蹤方法的準確性。

分別采用傳統(tǒng)LES 方法、曲線坐標LESPDF方法和DNS（Direct Numerical Simulation）進行模擬，其中DNS 網(wǎng)格為1 000×1 000 的正交網(wǎng)格，而LES 和曲線坐標LES-PDF 方法均采用100×100 的非正交曲線網(wǎng)格，網(wǎng)格分布如圖8 所示，圖中只給出了30 個網(wǎng)格，具體的網(wǎng)格坐標為

圖8 非正交曲線網(wǎng)格Fig. 8 Non-orthogonal curvilinear grid

式中：ξ和η為計算坐標；ε為網(wǎng)格擾動量級，取為0.5。

圖9給出了t=40時刻被動標量?的分布云圖。從圖中可以看出，曲線坐標下LES-PDF方法很好地捕捉了上下2股流體的標量界面，粒子分布與DNS計算的標量結(jié)果十分接近，表明了曲線坐標LESPDF方法以及單塊網(wǎng)格內(nèi)粒子追蹤方法的準確性。

圖9 標量? 的分布云圖Fig. 9 Contours of scalar ?

從圖10 中標量一階矩?ˉ的計算結(jié)果可以看出，曲線坐標PDF 方法與DNS 結(jié)果符合得很好，準確預測了大部分的標量細節(jié)分布，而LES 的結(jié)果則較為光滑。標量二階矩主要反映了可解尺度的相互作用，可以看出，曲線坐標PDF方法計算的結(jié)果不論是特征還是幅值上，與DNS符合得很好。標量二階矩表征了標量的亞格子脈動效應，該部分傳統(tǒng)的LES 由于缺乏亞格子信息，無法模擬，PDF 方法中粒子的隨機脈動模擬了亞格子脈動效果，從結(jié)果來看，曲線坐標PDF 方法定性地捕捉了混合層亞格子脈動規(guī)律，但是在混合層中心位置，預估的亞格子脈動要高于DNS，這可能是由于混合層中心網(wǎng)格加密不足，導致其高估了亞格子脈動效應。

圖10 標量統(tǒng)計結(jié)果Fig. 10 Scalar statistics results

總體來看，曲線坐標LES-PDF 方法較好地捕捉到了亞格子脈動效應，驗證了單塊網(wǎng)格內(nèi)曲線坐標下粒子定位、追蹤以及LES-PDF 耦合等方法的準確性。

3.3 多塊分區(qū)下的圓柱繞流

二維圓柱繞流主要用以驗證網(wǎng)格分區(qū)間粒子追蹤方法的準確性。該算例中PDF 主要求解了顯焓方程，通過粒子的溫度分布來驗證粒子是否正確穿越了分區(qū)以及是否正確處理了與激波的相互作用。入口來流馬赫數(shù)選擇為2.5，來流壓力為101 325.0 Pa，來流溫度為288.15 K。在該條件下，激波正好能穿越多個拓撲分區(qū)，能更好地驗證曲線坐標下粒子是否被正確輸運。入口給定超聲速入口條件，上下以及出口邊界均采用超聲速出口條件，圓柱表面為壁面邊界條件。

計算網(wǎng)格拓撲如圖11（a）所示，多個分區(qū)結(jié)合的角點，給粒子追蹤方法帶來了更大的挑戰(zhàn)。從圖11（a）中全場粒子分布結(jié)果可以看出，在角點位置后，存在一條明顯的空粒子帶。圖11（b）給出了局部放大的粒子、流線以及網(wǎng)格分布。從流場流線可以看出，粒子主要由分區(qū)1 向分區(qū)5進行移動，存在上下2 條路徑。粒子可以正確穿越從分區(qū)1、分區(qū)4 到分區(qū)5 的路徑。而粒子由分區(qū)1，經(jīng)過分區(qū)2 和分區(qū)3，無法正確通過最后的分區(qū)邊界。數(shù)值測試發(fā)現(xiàn)，主要是因為分區(qū)2 和分區(qū)3 在角點存在較大的網(wǎng)格畸變，在進行曲線坐標轉(zhuǎn)換時，粒子逆變速度出現(xiàn)明顯偏差，改變了粒子的正確運動方向。為解決上述問題，對角點附近的粒子追蹤過程進行修正：① 當前分區(qū)和目標分區(qū)邊界處使用相同的度量系數(shù)以及雅可比行列式，對分區(qū)邊界的度量系數(shù)進行光滑；② 粒子進入新的分區(qū)后，在新分區(qū)重新計算粒子逆變速度。修正后的粒子分布如圖11（c）所示，結(jié)果顯示粒子正確穿越了分區(qū)邊界，成功解決了空粒子帶問題。

圖11 網(wǎng)格畸變導致的粒子無法正確穿越邊界的問題Fig. 11 Problem of particles not crossing boundary correctly due to mesh distortion

圖12進一步給出了傳統(tǒng)LES 方法和曲線坐標PDF 方法計算的溫度云圖對比?？梢钥闯?，LES 網(wǎng)格顯示結(jié)果和PDF 粒子計算結(jié)果基本相同，激波位置以及尾跡的分布具有很好的一致性，定性驗證了多塊分區(qū)網(wǎng)格下曲線坐標粒子追蹤方法和粒子邊界條件的準確性。

圖12 LES 和PDF 溫度結(jié)果對比Fig. 12 Comparison of LES and PDF temperature results

3.4 復雜拓撲下的球頭繞流

進一步在三維球頭繞流算例中驗證曲線坐標LES-PDF 方法在三維復雜拓撲中的計算穩(wěn)定性。圖13（a）給出了三維球頭繞流算例的網(wǎng)格拓撲，此次計算中總共分了48 個分區(qū)，網(wǎng)格拓撲相較于二維圓柱繞流進一步復雜化。同時，分區(qū)的頂點也是典型的多塊接合點問題，而且這4 個頂點位置，網(wǎng)格的畸變均比較大。球頭繞流的計算條件和設置與二維圓柱繞流基本一致，PDF 主要求解了顯焓能量方程，通過能量場（溫度）的分布定性地顯示曲線坐標LES-PDF 方法的計算結(jié)果。圖13（b）給出了PDF 粒子的分布云圖，其中粒子由溫度著色，圖13（c）給出了中心截面上PDF 溫度分布。上述結(jié)果表明當前曲線坐標下的LES-PDF 方法具備處理三維復雜拓撲下的粒子追蹤和相應計算的能力。

圖13 球頭繞流網(wǎng)格拓撲及PDF 溫度計算結(jié)果Fig. 13 Mesh topology of spherical blunt and temperature results of PDF method

4 超聲速同軸射流火焰數(shù)值模擬

4.1 計算構型和邊界條件

基于Evans等［29］的超聲速同軸射流實驗，在實際超聲速湍流燃燒問題中對曲線坐標LES-PDF方法及其數(shù)值求解框架進行全面的驗證。Evans超聲速同軸射流實驗基本構型如圖14 所示，內(nèi)外噴管形成2 股平行的超聲速射流，內(nèi)部為低溫氫氣燃料，外部為高溫的污染空氣氧化劑流。內(nèi)部噴管外徑d為 9.525 mm，唇口厚度為 1.5 mm。為簡化計算構型，本文未考慮噴管內(nèi)部的流場結(jié)構。計算域軸向長度（Lx）為30 倍內(nèi)射流直徑，徑向長度（Lr）為2 倍射流直徑，本文中計算域用內(nèi)射流直徑進行無量綱化。軸向均勻分布480 個網(wǎng)格，周向的網(wǎng)格數(shù)目均為120 個，徑向網(wǎng)格在噴管內(nèi)外壁面均進行加密，入口截面上的網(wǎng)格分布如圖15 所示，網(wǎng)格拓撲為典型的“O-Block”分布。

圖14 同軸射流火焰示意圖Fig. 14 Schematic diagram of coaxial jet flame

圖15 入口截面網(wǎng)格和無量綱流向脈動速度分布Fig. 15 Mesh and dimensionless streamwise pulsation velocity distribution in entrance section

內(nèi)外射流入口均采用超聲速入口條件，根據(jù)Evans 等［29］的實驗參數(shù)直接給定。內(nèi)外射流詳細的流動參數(shù)如表1 所示。

表1 同軸射流火焰實驗參數(shù)Table 1 Parameters of coaxial jet flame experiment

本文的計算域以及邊界條件設置與劉朝陽等［30］的工作基本一致。劉朝陽等［30］采用的網(wǎng)格總數(shù)約為3 500 萬，燃燒模型為設定型PDF 模型。本文的網(wǎng)格總數(shù)約560 萬，采用數(shù)字濾波方法給定了如圖15 所示的湍流邊界層，邊界處厚度約為 2.5 mm?；瘜W反應計算采用Conaire 等［31］的9 組 分（H2、O2、H2O、H、O、OH、HO2、H2O2、N2） 21 步氫/空氣反應機理?；谶x擇的反應機理，PDF 求解顯焓的能量方程以及9 個組分的輸運方程，在燃燒條件下全面驗證曲線坐標LESPDF 方法的準確性。

4.2 燃燒流場結(jié)構

氫氣射流與高溫燃氣射流相互作用，在下游發(fā)展形成圖16 所示的超聲速抬舉射流火焰，圖中為射流中心截面上流場瞬態(tài)分布。從圖16（a）的軸向速度分布可以看出，由于唇口的存在，內(nèi)外2股射流在入口處形成2 道膨脹波，在射流中心處相交匯聚。膨脹波反射形成2 道激波，并在正中心出現(xiàn)一個馬赫盤。內(nèi)外射流之間具有很強的速度剪切，Kelvin-Helmholtz（K-H）不穩(wěn)定性導致剪切層失穩(wěn)。圖16（b）給出了混合分數(shù)結(jié)果，射流標量界面在K-H 渦結(jié)構作用下逐漸變得褶皺，向下游的發(fā)展過程中逐漸摻混。在該實驗中，高溫燃氣的靜溫為1 495 K，遠超氫/空氣的自著火溫度。如圖16（c）所示，射流經(jīng)過一小段混合距離之后，H2O 開始大量出現(xiàn)，并且H2O 的分布主要集中在K-H 渦結(jié)構中，表明燃燒主要發(fā)生在這些混合充分的K-H 渦結(jié)構中。圖16（d）進一步給出了燃燒釋熱率分布云圖，可以明顯看出在射流入口附近存在一個明顯的釋熱率峰值區(qū)域，其表征了火焰基底位置?？傮w上，曲線坐標LES-PDF 比較準確地反應了射流火焰流場的基本結(jié)構。

圖16 軸向截面流場云圖Fig. 16 Contours of flow field in axial plane

4.3 統(tǒng)計結(jié)果對比

圖17進一步給出了曲線坐標LES-PDF 方法計算的4個軸向位置上H2、O2、H2O 和N2沿 徑向的平均組分分布，并與實驗結(jié)果和劉朝陽等［30］的設定型PDF 數(shù)值模擬結(jié)果進行對比。

圖17（a）和圖17（b）給出了近場x/d=8.26和x/d=15.5 截面上的組分分布。近場的2 個軸向位置處的組分分布規(guī)律基本一致。由于剪切層中渦擴散和分子擴散作用，H2燃料向高溫氧化劑流中擴散，同時K-H 渦結(jié)構中混合較好的位置，開始出現(xiàn)著火，進一步消耗燃料，H2質(zhì)量分數(shù)沿徑向逐漸下降。外部高溫燃氣中的O2和N2則與之相反。近場燃燒僅在薄剪切層中進行，產(chǎn)物H2O 的峰值位于剪切層內(nèi)。整體上曲線坐標LES-PDF 方法計算的結(jié)果與設定型PDF 模型的結(jié)果相似，均低估了近場剪切層中的燃燒強度。

圖17 不同軸向位置處徑向組分分布Fig. 17 Radial component distribution at different axial positions

圖17（c）和圖17（d）給出了下游x/d=21.7和x/d=27.9 火焰區(qū)內(nèi)的平均組分分布。在下游位置，從產(chǎn)物H2O 的分布可以看出，火焰區(qū)域明顯擴大。而射流中心的H2濃度比實驗值偏高，說明燃料消耗不足，預測的燃燒強度仍低于實際情況。曲線坐標LES-PDF 方法計算的產(chǎn)物H2O的分布范圍與實驗更加接近，而設定型PDF 模型的計算結(jié)果則相對低估了下游火焰區(qū)內(nèi)的H2O的生成和H2的消耗。此外，本文計算的N2結(jié)果存在一定的不足，明顯高估了N2向射流中心的擴散。總體來看，本文的網(wǎng)格數(shù)量相比劉朝陽等［30］小接近一個量級，盡管預測的被動標量N2的分布稍 差，但H2和O2等 反應標量 與劉朝陽 等［30］的 設定型PDF 模型計算的結(jié)果基本相當，尤其是下游火焰區(qū)內(nèi)反應生成物H2O 的分布明顯優(yōu)于設定型PDF 模型，這些定量結(jié)果表明了曲線坐標LES-PDF 方法的準確性。

5 結(jié) 論

1）針對超燃沖壓發(fā)動機燃燒室中超聲速湍流燃燒數(shù)值模擬的實際需求，面向發(fā)動機實際構型，建立了一般曲線坐標系下的LES-PDF 方法?；诙鄩K結(jié)構網(wǎng)格，對PDF 輸運方程和粒子隨機微分方程進行了坐標變換，實現(xiàn)了基于曲線坐標的粒子追蹤和粒子邊界條件，解決了多塊網(wǎng)格分區(qū)間粒子連續(xù)追蹤問題，實現(xiàn)了復雜構型下的曲線坐標LES-PDF 計算。

2）由簡單到復雜逐步對曲線坐標LESPDF 方法中的關鍵問題進行了充分的驗證。一維Sod 激波管和二維時間混合層算例中，曲線坐標LES-PDF 方法計算的結(jié)果與精確解和直接數(shù)值模擬結(jié)果符合得很好，表明了曲線坐標形式的PDF 輸運方程及其粒子隨機微分方程的正確性。二維圓柱繞流算例中發(fā)現(xiàn)，在多個分區(qū)接合的角點位置，由于網(wǎng)格畸變較大，導致粒子速度計算出現(xiàn)較大誤差，導致其無法正確穿越分區(qū)邊界，在網(wǎng)格分區(qū)角點后出現(xiàn)了一條明顯的空粒子帶。通過對分區(qū)間度量系數(shù)和雅可比行列式的光滑，并在新分區(qū)重新計算粒子速度，成功解決了多塊分區(qū)間粒子丟失問題。三維復雜拓撲球頭繞流結(jié)果進一步證明了當前曲線坐標LES-PDF 方法對復雜拓撲的處理能力。

3）超聲速同軸射流火焰模擬中，曲線坐標LES-PDF 方法給出了較為準確的射流火焰流場基本結(jié)構，同時不論是近場還是遠場，該方法預測的組分分布整體上與實驗均吻合較好，尤其在遠場火焰區(qū)內(nèi)，其給出了更準確的反應生成物H2O 的分布，在實際超聲速燃燒問題中定量驗證了曲線坐標LES-PDF 方法的準確性。