1
低馬赫數近似法的化學反應直接數值模擬
Application of Spectr al Element Method in
Low Mach Number Combustion
計畫編號：NSC 87-2213-E-002-004
執行期限：86 年 8 月 1 日 – 87 年 7 月 31 日
主持人：顏瑞和 國立台灣大學機械工程學系
一、 中文摘要 本文的主旨在應用寬頻元素法於低馬赫 數燃燒問題上。首先，把完全可壓縮流的方 程式化簡成低馬赫數下的近似方程組，這些 方程組過濾掉高頻率的聲波，能允許較大的 時間間隔來數值計算。在本文中，我們利用 寬頻元素法來處理統御方程式中之空間離散 化，並以高階分離法處理時間離散化。接著 利用一些數學例子來測試我們程式的正確 性。在二維可壓縮流的問題上，模擬了混合 層的燃燒現象，模擬的結果顯示，隨著放熱 率的增加，生成物形成的速率與質量進入此 混合層的量將減少。 關鍵詞：寬頻元素法，低馬赫數，直接數值 模擬 ABSTRACTIn this work , Low Mach Number combustion problems were investigated. The fully compressible flow equations were simplified by the Low Mach Number approximation . The approximate equations can filter high frequency acoustic wave to allow larger time steps to be taken in numerical soluations . In this thesis , we use spectral element method for the sptail discretization of governing equations and high order splitting method for temporal discretization.
To validate the numerical accuracy , some benchmark solutions were compared. The numerical solution for incompressible flow , which is the limiting case of the Low Mach
Number approximation was compared with the incompressible code which has developed in the previous work.
For two-dimensional compressible
combustion problem , we simulate mixing layer under chemical reaction. The results indicate that rate of chemical product formed and the amount of mass entrained to the layer all decrease with increasing rates of heat release.
二、緣由與目的
對於含有極小尺度(small scale)激烈變化 的流場計算如紊流或燃燒上問題的直接數值 模擬(DNS , Direct Numerical Simulation)，若 用傳統的低階數值方法如有限差分法、有限 元素法，所需要的格點數目必需大幅的增加 才能解析紊流或燃燒流場中較小尺度的結 構，相對的要花費非常龐大的計算時間及記 憶體容量。而 Orszag(1972)利用寬頻法來模擬 紊流問題之研究後，就掀起了以 DNS 為工具 來研究紊流上的問題；然而應用寬頻法或寬 頻元素法來模擬燃燒問題研究的文獻相當少 也僅從 1988 年起 Givi(1989)、Wessel(1993) 等幾個人而已。 在計算可壓縮流的 Navier-Stokes 動量方 程式、能量方程式、與組份傳輸方程式。這 些方程式包含了渦度(vorticity)、熵(entropy) 與聲波 (acoustic)模 式。在許多 實際的情 形 下，這些模式在頻率上有很大的差異。特別 是在低馬赫數下的流場，聲波在頻率帶上比 其他兩個模式高，所以聲波屬於高頻率帶。 因此，在一個低馬赫數的流場中，聲波的擾
2 動不能有效的與渦度或熵模式交互作用；從 計算的觀點而言，要探知聲波擾動的現象必 須要極小的時間間隔，也要花費龐大的計算 時間。所以若要對聲波、渦度、熵模式要同 時配合考慮，則必將降低計算的效率。因此， 本文將參考 McMur ty et al. (1987)之論文，再 配合寬頻元素法來模擬計算。 三、數學及物理模式 本研究為建立密度為變數之數學模式， 並利用低馬赫數近似法的特性過濾掉高頻率 的聲波來克服完全可壓縮流上時間間隔無法 跨大的問題。關於本研究所使用之基本假設 如下︰(1)、低馬赫數(2)、輸送性質，黏性係 數，定容比熱為常數(3) 、Stoke’s hypothesis µ λ = −2/3 。(4) 忽略熱輻射效應。(5) 理 想氣體。接著我們將統御方程式無因次化並 用低馬赫數近似法展開，可得到如下之統御 方程式： 零階方程式 (1)連續方程式 ∂ρ ∂ ∂ ∂ ρ ( ) ( ) ( ) ( ) 0 0 0 0 t + xi ui = (2.1) (2)動量方程式 ∂ ∂ P xi ( )0 0 = (2.2) (3)能量方程式 ρ γ ∂ ∂ ∂ ∂ ∂ ∂ ( ) ( ) ( ) ( ) ( ) ( ) Pr Re [ ] 0 0 0 0 0 1 1 DT Dt P u x x T x D C Q i i i i a e = − − + + (2.3) (4)組份方程式 ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ Y t x Y u Pe x Y x DaR k j k j j k j ( ) ( ) ( ) ( ) [ ] [ ] 0 0 0 0 1 + = + (2.4) (5)狀態方程式 P( )0 = ρ( )0 T( )0 (2.5) 由於P( )0 熱力壓力不為空間的函數，只為與 時間相關的函數，因此在每個計算時間上熱 力壓力 P( )0 在整個計算空間是一致的，也就是 在這個近似法上音速以無限大的速度傳播， 所以任何熱力壓力P( )0 的擾動能迅速的影響 到整個流場。此外為了能夠完整的描述流場 的特性，第一階得動量方程式必須被引進， 整理可得： ρ ∂ ∂ ∂ ∂ ∂ ∂ µ ∂∂ ∂ ∂ µ ∂ ∂ ∂ ∂ ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ ] [ ] [ ] 0 0 0 0 1 0 0 3 u t u u x P x x u u x u x i j i j i j i j i i j + = − + + ~O(ε) (2.6) 此外，我們將考慮進一步調整上述方程 組的形式來提供數值解上的效率。將連續性 方程式與能量方程式合併成： ∇⋅ =vϖ ∇ ⋅ϖ− + P q dP dt DaCeQ 1 0 0 γ ( ) ( ) [ Pr Re &] (2.7) 在一個封閉且絕熱的邊界或計算區域為 週期性邊界條件下，則上式可被體積分成 dP dt volume Ce ( ) & 0 1 =
∫∫∫
Da Q dV (2.8) 四、數值方法 以分離法(spltting scheme)將(2.6)式分 成三個步驟，而式中的密度可由連續性方程 式求得。 步驟一、(顯性步驟) ρu∃ α u β ( , ) t N u u i q i n q q J q q J i n q j n q − = − − = − = − − −∑
∑
0 1 0 1 ∆ (3.1) 步驟二、(壓力步驟) ρ ∂ ∂ ∃∃ ∃ u u t i i i − = − ∆ P x n+1 (3.2) 解出流力壓力 P( )1n+1。 步驟三、(隱性步驟) − + = + + ∂ ∂ µ ∂∂ ργ ρ x u x j j i n i n i u t u t 1 0 1 ∆ ∆ ∃∃ (3.3) 解出 ui n+1 。 組份方程式[ ]
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ φ Y t x x Y x k j j j k j + = + Y u P S k e 1 (3.4)3 利用前述之mixed stiffly stable及splitting方法 將(3.4)式寫成下列二個步驟：
[
]
∃ Y t tY Y u S k q q J k n q q q J k i ∆ = − ∆ − + − − = −∑
α∑
β ∂ ∂ φ 0 1 0 1 xi (3.5) − + = + 1 0 1 P Y Y e k n k ∂ ∂ ∂ ∂ γ x Y x t t i k n+1 i ∆ ∆ ∃ (3.6) 寬頻元素法在空間離散化的處理步驟， 詳 細 的 定 義 及 分 析 請 見 Rϕnquist(1988) 、 Tombulides(1993)。 五、應用於化學反應之模擬 我們將考慮燃燒過程是在一個等體積且 單一階、不可逆之化學反應下進行，反應速 率假設只為反應物的函數與溫度無關，反應 式如下： C1 + C2 product + heat 其中 C1、 C2代表兩種不同反應物之組份。 各初始狀態及邊界條件分別為： 初始狀態： 1.組份 C x y y y t y 1 0 0 2 0 2 1 = = ∞ − ∫
( , ) exp π ξ ξ d C2t=0( , )x y = −1 C1 2.溫度場 無 因 次 溫 度 初 始 值 假 設 為 常 溫。 3.速度場 速度場開始假設為雙曲正切函 數加上基本模式(fundamental mode)與次和諧 模式(subharmonic mode)之擾動函數。 邊界條件： 1.速度場 速度邊界條件為週期性邊界，即：v(x,y)x=14 = v(x,y)x=-14 on Γ1 and Γ2
∂ ∂ ∂ ∂ v y v y , y=−14 = y=14 =0 on Γw1 2.組份 Y x yi( , )x=−14 =Y x yi( , )x=14 , on Γ1 and Γ2 ∂ ∂ ∂ ∂ Y y Y y , i i y=−14 = y=14 =0 on Γw1 在這個數值模擬中，我們將固定參數 Re=250，Pe=100，Da=1.5 而調整不同的放熱 參數 Ce 來研究不同的放熱量對整個流場的影 響。為了只要單純的研究化學放熱對流體的 影響，此反應中假設黏性、熱擴散係數、分 子擴散係數均為與溫度無關的常數、而密度 場的變化為整體性的。圖 1 為生成物在反應 混合層內兩種不同的模擬結果，所有的參數 都一樣，只有放熱量不相同，Ce 分別為 0、 2.5。而不同的放熱參數 Ce 明顯地影響了生成 物的形成，在有放熱狀況下將使生成物的產 量減少。圖 2 中，t=36 時通過混合層沿水平 軸的平均生成物，在放熱參數分別為 0、1.2、 2.5 時，當 Ce=0(即無放熱狀態下)平均生成物 為最大，隨著放熱參數 Ce 的增加，平均生成 物也相對的減少，所以隨著放熱率的增加， 產物的生成率與質量進入此混合層的量將隨 之減少。 圖 3 顯示通過混合層中的平均速度 曲線中，在t=36 時我們發現隨著放熱率的增 加密度的減少使得平均速度曲線的斜率在混 合層中將變的較大。當放熱參數 Ce=0 時平均 速度曲線隨時間的變化量最大，所以混合層 的成長速率較快，而隨著放熱率的增加由於 因密度的減少速率較快，因此整個混合層的 成長速率較緩慢，因此放熱扮演了使質量減 少進入此混合層的一個因素。以上這些分析 的結果與 McMurtry et al.(1987)利用類寬頻法 (Pseudospectral Method)採用等間隔所做的模 擬結果相同。 當放熱參數 Ce=0.8 時發現非穩態燃 燒流場t=0至 36 中密度的流場變化∆ρ/ρ控
4 制在大約在百分之五左右，近似於不可壓縮 流的流場。由圖 4 中可以較仔細的觀察出， 當 t 分別為 12、24、36 時，分別用低馬赫數 近似法與不可壓縮流在 x=0 截面上溫度的分 佈的比較情況，我們發現溫度輪廓的分佈均 由兩個不同組份的交界面(y=0)漸漸的隨著初 始速度擾動所引起的流場呈對稱性的分佈且 由於混合層的成長溫度分佈的範圍也逐漸增 大，在反應區中溫度分佈在這個非穩態燃燒 的混合過程中大致相同，而由於低馬赫數近 似法模擬的燃燒流場中音速是無限大的，所 以熱力壓力 P( )0 能迅速的影響到整個流場，因 此在非反應區中溫度分佈也受到熱力壓力 P( )0 影響與不可壓縮流模擬的結果有一些差 異，但差異並不大。圖 5 為平均生成物分佈 的比較，在有密度的變化下，平均生成物略 為減少。由圖 6 平均速度曲線在 t=36 速度場 的比較也大致相同。 由以上的分析討論，近似於不可壓縮流 與不可壓縮流模擬的結果在溫度場及速度場 的分析結果相近，所以除了先前我們對方程 式做獨立性測試外，也能間接的驗證了低馬 赫數近似法模擬的正確性。 六、結論與建議 為了更符合實際物理現象之模擬，利用 低馬赫數近似法將高頻率的聲波過濾掉克服 了在解完全可壓縮流問題上時間間隔無法調 大的限制，而來模擬一個非穩態的燃燒混合 現象。有如下之發現：(a)隨著放熱量的增加 混合層中平均生成物將變得較少且平均速度 曲線斜率變得較陡。(b)適當調整放熱參數 Ce 在近似於不可壓縮流的模擬中與不可壓縮流 的結果做比較，在速度場與溫度場中皆能獲 的類似的結果，所以也能間接說明了我們利 用低馬赫數近似法模擬的正確性。 七、參考文獻
Givi , P. and Jou , W. H. (1988)：”Direct numerical simulation of a twodimensional
reacting , spatially developing mixing layer by a
spectral element method,” Twenty-Second
Symposium (International) On Combustion , The Combustion Institute Pittsbrugh , pp 635 - 643。
McMurty P.A. , J. J. Riley , and , R. W. Metcalfe (1989),”Effects of heat -release on the large-scale structure in turbulent mixing layers,” J. of Fluid Mechanic vol 199 ,pp297-332。
Rϕnquist E. M. (1988) ,”Optimal Spectral
Element Method for Unstedy Three-Dimensional Incompressible Navier-Stokes Equation ,”,Ph. D. Thesis , Department of Mechanical Engineering , Massachusetts Institute of Technology。
Wessel , R. A. (1993)：”Spectral element method for numerical simulation of unsteady laminar diffusion flames” , Ph.D. Thesis , Department of Mechanical and Aerospace Engineering , Case Western Reserve University。
