燃烧仿真-主题034-湍流燃烧模型_RANS方法

kkchenjj

23人浏览 · 2026-03-13 20:36:01

kkchenjj · 2026-03-13 20:36:01 发布

第034篇：湍流燃烧模型（RANS方法）

摘要

湍流燃烧模型是工业燃烧数值模拟的核心技术，其中雷诺平均Navier-Stokes（RANS）方法因其计算效率优势而被广泛应用。本教程系统阐述基于RANS框架的湍流燃烧建模理论，涵盖湍流模型（k-ε、RNG k-ε、Realizable k-ε、SST k-ω）、燃烧模型（涡破碎模型EBU、涡耗散模型EDM、涡耗散概念EDC、火焰面模型、PDF输运方程模型）以及它们的耦合策略。通过Python仿真实现，演示不同湍流燃烧模型在射流火焰、旋流燃烧器等典型工况下的预测性能对比，分析模型的适用范围和局限性。本教程为工程燃烧仿真提供RANS建模的理论基础和实用指导。

关键词

湍流燃烧模型, RANS方法, k-ε湍流模型, 涡破碎模型, 涡耗散概念, 火焰面模型, 混合分数, 燃烧仿真

1. 引言

1.1 湍流燃烧的挑战

湍流燃烧是能源、动力和化工领域的核心过程，涉及复杂的流体动力学与化学反应相互作用。与层流燃烧相比，湍流燃烧具有以下显著特征：

多尺度耦合：湍流具有从积分尺度到Kolmogorov尺度的宽谱能量分布，而化学反应发生在分子尺度，两者的时间尺度和空间尺度差异可达数个数量级。

强非线性：化学反应速率对温度呈指数依赖（Arrhenius定律），湍流脉动导致温度和浓度的剧烈波动，使得平均反应速率不能简单地用平均量计算。

火焰结构复杂：湍流火焰可能呈现褶皱层流火焰、分布反应区或破碎反应区等不同模式，取决于湍流强度与化学反应特征时间的相对大小。

计算成本高昂：直接数值模拟（DNS）可以完整解析所有尺度，但对于工程尺度的燃烧设备，计算资源需求是不可承受的。

1.2 RANS方法的优势与局限

雷诺平均Navier-Stokes（RANS）方法通过对Navier-Stokes方程进行系综平均，将瞬态变量分解为平均量和脉动量：

$ϕ=ϕˉ+ϕ′\phi = \bar{\phi} + \phi'$

其中 $ϕˉ\bar{\phi}$ 为系综平均量， $ϕ′\phi'$ 为脉动量，且 $ϕ′‾=0\overline{\phi'} = 0$ 。

RANS方法的优势：

计算效率高，适合工程应用
网格要求相对较低
成熟的商业软件支持（ANSYS Fluent, STAR-CCM+等）
丰富的工程验证案例

RANS方法的局限：

无法解析大尺度非定常结构
对强旋流、分离流等复杂流动预测精度有限
湍流-燃烧相互作用建模存在固有假设
对瞬态现象的预测能力受限

1.3 湍流燃烧建模的基本框架

完整的RANS湍流燃烧模拟需要解决三个核心问题：

1. 湍流输运建模：预测平均速度场和湍流输运特性

2. 化学反应建模：描述燃料氧化和污染物生成反应

3. 湍流-燃烧相互作用：封闭平均化学反应源项

$\overline{\dot{\omega}_k} = \overline{A \rho Y_k \exp(-E_a/RT)} \neq A \bar{\rho} \bar{Y}_k \exp(-E_a/R\bar{T})$

上式表明，由于Arrhenius反应率的非线性，平均反应速率不等于用平均量计算的反应速率，这是湍流燃烧建模的核心难题。

2. RANS湍流模型

2.1 标准k-ε模型

标准k-ε模型由Launder和Spalding于1972年提出，是工程应用最广泛的湍流模型。

控制方程：

湍动能k的输运方程：
$∂(ρˉk)∂t+∂(ρˉu~jk)∂xj=∂∂xj[(μ+μtσk)∂k∂xj]+Pk−ρˉε\frac{\partial (\bar{\rho} k)}{\partial t} + \frac{\partial (\bar{\rho} \tilde{u}_j k)}{\partial x_j} = \frac{\partial}{\partial x_j}\left[\left(\mu + \frac{\mu_t}{\sigma_k}\right)\frac{\partial k}{\partial x_j}\right] + P_k - \bar{\rho}\varepsilon$

湍流耗散率ε的输运方程：
$∂(ρˉε)∂t+∂(ρˉu~jε)∂xj=∂∂xj[(μ+μtσε)∂ε∂xj]+Cε1εkPk−Cε2ρˉε2k\frac{\partial (\bar{\rho} \varepsilon)}{\partial t} + \frac{\partial (\bar{\rho} \tilde{u}_j \varepsilon)}{\partial x_j} = \frac{\partial}{\partial x_j}\left[\left(\mu + \frac{\mu_t}{\sigma_\varepsilon}\right)\frac{\partial \varepsilon}{\partial x_j}\right] + C_{\varepsilon1}\frac{\varepsilon}{k}P_k - C_{\varepsilon2}\bar{\rho}\frac{\varepsilon^2}{k}$

其中湍动能生成项：
$Pk=μt(∂u~i∂xj+∂u~j∂xi)∂u~i∂xj−23ρˉkδij∂u~i∂xjP_k = \mu_t \left(\frac{\partial \tilde{u}_i}{\partial x_j} + \frac{\partial \tilde{u}_j}{\partial x_i}\right)\frac{\partial \tilde{u}_i}{\partial x_j} - \frac{2}{3}\bar{\rho}k\delta_{ij}\frac{\partial \tilde{u}_i}{\partial x_j}$

湍流粘度：
$μt=Cμρˉk2ε\mu_t = C_\mu \bar{\rho} \frac{k^2}{\varepsilon}$

模型常数：

$Cμ=0.09C_\mu = 0.09$
$Cε1=1.44C_{\varepsilon1} = 1.44$
$Cε2=1.92C_{\varepsilon2} = 1.92$
$σk=1.0\sigma_k = 1.0$
$σε=1.3\sigma_\varepsilon = 1.3$

适用范围：

充分发展的湍流
高雷诺数流动
弱旋流、弱压力梯度流动

局限性：

对强旋流预测偏差较大
近壁区域需要壁面函数或低雷诺数修正
对逆压梯度流动和分离流预测不够准确

2.2 RNG k-ε模型

RNG（Renormalization Group）k-ε模型基于重整化群理论推导，对标准k-ε模型进行了改进。

主要改进：

ε方程系数修正：
$Cε1=1.42−η(1−η/η0)1+βη3C_{\varepsilon1} = 1.42 - \frac{\eta(1-\eta/\eta_0)}{1+\beta\eta^3}$
其中 $η=Sk/ε\eta = Sk/\varepsilon$ ， $η0=4.38\eta_0 = 4.38$ ， $β=0.015\beta = 0.015$
附加湍流生成项：考虑平均应变率对湍流的影响
有效粘度修正：在低雷诺数区域自动过渡为分子粘度

优势：

对快速应变流和中等旋流有更好的预测
可以更好地处理流线弯曲和旋转效应
对圆射流的扩展率预测更准确

2.3 Realizable k-ε模型

Realizable k-ε模型由Shih等人提出，旨在解决标准k-ε模型在某些流动条件下的非物理行为。

"可实现性"约束：

保证湍流应力的物理可实现性，即满足：
$u′2‾≥0,u′v′‾2≤u′2‾v′2‾\overline{u'^2} \geq 0, \quad \overline{u'v'}^2 \leq \overline{u'^2}\overline{v'^2}$

主要改进：

可变 $CμC_\mu$ ：
$Cμ=1A0+AskU∗εC_\mu = \frac{1}{A_0 + A_s \frac{kU^*}{\varepsilon}}$
其中 $A_0 = 4.04$ ， $As=6cos⁡ϕA_s = \sqrt{6}\cos\phi$
修正的ε方程：
$Cε1=max⁡(0.43,ηη+5),Cε2=1.9C_{\varepsilon1} = \max\left(0.43, \frac{\eta}{\eta+5}\right), \quad C_{\varepsilon2} = 1.9$

优势：

对强逆压梯度、分离流和射流有更好的预测
对圆射流和平面射流的差异预测更准确
旋转和曲率效应处理更好

2.4 SST k-ω模型

SST（Shear Stress Transport）k-ω模型由Menter提出，结合了近壁k-ω模型的优势和对数律区域k-ε模型的优势。

混合函数：

在边界层内部使用k-ω模型，在外部区域混合使用k-ε模型：

$∂(ρˉk)∂t+∂(ρˉu~jk)∂xj=∂∂xj[(μ+σkμt)∂k∂xj]+Pk−β∗ρˉkω\frac{\partial (\bar{\rho} k)}{\partial t} + \frac{\partial (\bar{\rho} \tilde{u}_j k)}{\partial x_j} = \frac{\partial}{\partial x_j}\left[\left(\mu + \sigma_k \mu_t\right)\frac{\partial k}{\partial x_j}\right] + P_k - \beta^*\bar{\rho}k\omega$

$∂(ρˉω)∂t+∂(ρˉu~jω)∂xj=∂∂xj[(μ+σωμt)∂ω∂xj]+αρˉμtPk−βρˉω2+2(1−F1)ρˉσω2ω∂k∂xj∂ω∂xj\frac{\partial (\bar{\rho} \omega)}{\partial t} + \frac{\partial (\bar{\rho} \tilde{u}_j \omega)}{\partial x_j} = \frac{\partial}{\partial x_j}\left[\left(\mu + \sigma_\omega \mu_t\right)\frac{\partial \omega}{\partial x_j}\right] + \frac{\alpha\bar{\rho}}{\mu_t}P_k - \beta\bar{\rho}\omega^2 + 2(1-F_1)\frac{\bar{\rho}\sigma_{\omega2}}{\omega}\frac{\partial k}{\partial x_j}\frac{\partial \omega}{\partial x_j}$

SST修正：

限制湍流应力以考虑剪切流中湍流各向异性的影响：
$μt=a1kmax⁡(a1ω,SF2)\mu_t = \frac{a_1 k}{\max(a_1 \omega, SF_2)}$

优势：

对逆压梯度和分离流的预测显著优于k-ε模型
近壁区域无需壁面函数
对航空发动机燃烧室等复杂流动有很好的预测能力

3. 湍流燃烧模型

3.1 涡破碎模型（EBU）

涡破碎（Eddy Break-Up）模型由Spalding于1971年提出，是最早的湍流燃烧模型之一。

基本假设：

化学反应速率无限快（快化学假设）
燃烧速率由湍流混合控制
湍流将未燃混合物破碎成适合燃烧的尺度

模型公式：

对于燃料质量分数 $Y_{fu}$ ，平均反应速率为：

$ω˙fu‾=−CEBUρˉεkmin⁡(Yˉfu,Yˉoxs,BYˉpr1+s)\overline{\dot{\omega}_{fu}} = -C_{EBU}\bar{\rho}\frac{\varepsilon}{k}\min\left(\bar{Y}_{fu}, \frac{\bar{Y}_{ox}}{s}, B\frac{\bar{Y}_{pr}}{1+s}\right)$

其中：

$C_{EBU}$ ：模型常数，通常取0.5-1.0
$s$ ：化学计量比
$B$ ：产物对反应速率的影响系数，通常取0.5

特点：

简单高效，不需要详细的化学反应机理
假设反应速率无限快，无法预测点火和熄火
适用于高Damköhler数（ $\gg 1$ ）的预混燃烧

3.2 涡耗散模型（EDM）

涡耗散模型（Eddy Dissipation Model）是EBU模型的改进版本，由Magnussen和Hjertager提出。

模型公式：

$ω˙fu‾=−Aρˉεkmin⁡(Yˉfu,Yˉoxs)\overline{\dot{\omega}_{fu}} = -A\bar{\rho}\frac{\varepsilon}{k}\min\left(\bar{Y}_{fu}, \frac{\bar{Y}_{ox}}{s}\right)$

或考虑产物的影响：

$ω˙fu‾=−AρˉεkYˉpr1+Ypr,limmin⁡(Yˉfu,Yˉoxs)\overline{\dot{\omega}_{fu}} = -A\bar{\rho}\frac{\varepsilon}{k}\frac{\bar{Y}_{pr}}{1+Y_{pr,lim}}\min\left(\bar{Y}_{fu}, \frac{\bar{Y}_{ox}}{s}\right)$

模型常数：

$A$ ：通常取4.0
$Y_{pr,lim}$ ：产物限制器，通常取0.01

与EBU模型的区别：

EDM模型通常不考虑产物对反应速率的促进作用
EDM模型常数A通常大于EBU模型的 $C_{EBU}$
EDM模型更常用于非预混燃烧

3.3 涡耗散概念（EDC）

涡耗散概念（Eddy Dissipation Concept）模型由Magnussen提出，是对EDM模型的扩展，可以考虑有限速率化学。

基本思想：

湍流将流体分割成 fine structures（精细结构），在这些结构中发生化学反应。精细结构占据的体积分数为：

$γ∗=2.13(νεk2)1/4\gamma^* = 2.13\left(\frac{\nu\varepsilon}{k^2}\right)^{1/4}$

精细结构中的特征时间尺度：

$τ∗=νε\tau^* = \sqrt{\frac{\nu}{\varepsilon}}$

反应速率计算：

在精细结构中，反应速率由化学反应动力学和湍流混合共同决定：

$ω˙k‾=γ∗χρ∗Yk∗−ρˉYˉkτ∗\overline{\dot{\omega}_k} = \gamma^* \chi \frac{\rho^* Y_k^* - \bar{\rho}\bar{Y}_k}{\tau^*}$

其中 $χ\chi$ 为反应概率， $ρ∗\rho^*$ 和 $Y_k^*$ 为精细结构中的密度和组分质量分数。

优势：

可以考虑有限速率化学效应
能够预测点火和熄火现象
适用于部分预混燃烧

3.4 火焰面模型

火焰面模型基于层流火焰概念，假设湍流火焰由层流火焰面组成，这些火焰面被湍流褶皱和拉伸。

3.4.1 层流火焰面模型（Laminar Flamelet Model）

基本假设：

化学反应发生在薄的层流火焰面内
火焰面厚度远小于湍流积分尺度
火焰面内部结构近似为层流火焰

混合分数方法：

对于非预混燃烧，引入混合分数 $Z$ ：

$\frac{\phi - \phi_{ox}}{\phi_{fuel} - \phi_{ox}}$

其中 $ϕ\phi$ 可以是元素质量分数或混合物分数。

混合分数的输运方程：

$∂(ρˉZ~)∂t+∂(ρˉu~jZ~)∂xj=∂∂xj[(ρˉD+μtSct)∂Z~∂xj]\frac{\partial (\bar{\rho}\tilde{Z})}{\partial t} + \frac{\partial (\bar{\rho}\tilde{u}_j\tilde{Z})}{\partial x_j} = \frac{\partial}{\partial x_j}\left[\left(\bar{\rho}D + \frac{\mu_t}{Sc_t}\right)\frac{\partial \tilde{Z}}{\partial x_j}\right]$

混合分数方差的输运方程：

$∂(ρˉZ′′2~)∂t+∂(ρˉu~jZ′′2~)∂xj=∂∂xj[(ρˉD+μtSct)∂Z′′2~∂xj]+2μtSct(∂Z~∂xj)2−ρˉχ~\frac{\partial (\bar{\rho}\tilde{Z''^2})}{\partial t} + \frac{\partial (\bar{\rho}\tilde{u}_j\tilde{Z''^2})}{\partial x_j} = \frac{\partial}{\partial x_j}\left[\left(\bar{\rho}D + \frac{\mu_t}{Sc_t}\right)\frac{\partial \tilde{Z''^2}}{\partial x_j}\right] + 2\frac{\mu_t}{Sc_t}\left(\frac{\partial \tilde{Z}}{\partial x_j}\right)^2 - \bar{\rho}\tilde{\chi}$

其中标量耗散率：

$χ~=CχεkZ′′2~\tilde{\chi} = C_\chi \frac{\varepsilon}{k}\tilde{Z''^2}$

火焰面数据库：

通过求解一维层流火焰方程，建立混合分数与标量耗散率到温度、组分质量分数的映射关系：

$ϕ=ϕ(Z,χst)\phi = \phi(Z, \chi_{st})$

然后通过概率密度函数（PDF）进行平均：

$ϕ~=∫01∫0∞ϕ(Z,χst)P(Z,χst)dZdχst\tilde{\phi} = \int_0^1 \int_0^\infty \phi(Z, \chi_{st}) P(Z, \chi_{st}) dZ d\chi_{st}$

3.4.2 稳态层流火焰面模型（SLFM）

假设火焰面始终处于化学平衡或稳态，忽略瞬态效应。

适用条件：

高Damköhler数
化学反应时间尺度远小于流动时间尺度
局部熄火可以忽略

3.4.3 非稳态层流火焰面模型（USFM）

考虑火焰面的瞬态响应，可以模拟点火、熄火和再燃等现象。

火焰面方程：

$ρ∂ϕ∂τ=ρχ2∂2ϕ∂Z2+ω˙ϕ\rho \frac{\partial \phi}{\partial \tau} = \rho \frac{\chi}{2}\frac{\partial^2 \phi}{\partial Z^2} + \dot{\omega}_\phi$

其中 $τ\tau$ 为火焰面时间坐标。

3.5 PDF输运方程模型

概率密度函数（PDF）方法是湍流燃烧建模中最严格的方法之一，直接求解标量联合PDF的输运方程。

联合PDF定义：

对于一组标量 $ψ‾=(ψ1,ψ2,...,ψn)\underline{\psi} = (\psi_1, \psi_2, ..., \psi_n)$ ，其联合PDF $P(ψ‾;x,t)P(\underline{\psi}; \mathbf{x}, t)$ 定义为：

$P(ψ‾)dψ‾=Prob(ψ1≤ϕ1<ψ1+dψ1,...,ψn≤ϕn<ψn+dψn)P(\underline{\psi})d\underline{\psi} = \text{Prob}(\psi_1 \leq \phi_1 < \psi_1 + d\psi_1, ..., \psi_n \leq \phi_n < \psi_n + d\psi_n)$

PDF输运方程：

$ρˉ∂P∂t+ρˉu~j∂P∂xj=−∂∂ψk[⟨ω˙k∣ψ‾⟩P]−∂∂xj[⟨uj′′∣ψ‾⟩P]+∂∂ψk[⟨∂Jk,j∂xj∣ψ‾⟩P]\bar{\rho}\frac{\partial P}{\partial t} + \bar{\rho}\tilde{u}_j\frac{\partial P}{\partial x_j} = -\frac{\partial}{\partial \psi_k}[\langle\dot{\omega}_k | \underline{\psi}\rangle P] - \frac{\partial}{\partial x_j}[\langle u_j'' | \underline{\psi}\rangle P] + \frac{\partial}{\partial \psi_k}[\langle\frac{\partial J_{k,j}}{\partial x_j} | \underline{\psi}\rangle P]$

求解方法：

由于PDF输运方程的高维性（维度=空间维数+标量数），通常采用蒙特卡洛方法求解。

优势：

化学反应源项封闭，无需建模
可以处理任意复杂的化学反应机理
能够准确预测污染物生成

局限：

计算成本高昂
需要大量的计算粒子
小尺度混合模型存在不确定性

4. 模型耦合策略

4.1 湍流-燃烧模型匹配

不同的湍流模型与燃烧模型有不同的匹配策略：

k-ε模型 + EBU/EDM：

最简单的组合
适用于快速预混燃烧
计算效率高

k-ε模型 + 火焰面模型：

适用于非预混燃烧
需要求解混合分数和方差输运方程
可以预测局部熄火

SST k-ω + EDC：

适用于复杂几何和旋流燃烧器
可以考虑有限速率化学
对近壁区域处理更好

4.2 壁面处理

近壁区域的湍流和燃烧行为对整体预测有重要影响。

壁面函数法：

标准壁面函数
非平衡壁面函数
增强壁面处理

低雷诺数模型：

直接解析粘性底层
需要更细的近壁网格
SST k-ω模型天然支持

4.3 辐射传热耦合

高温燃烧中，辐射传热对能量平衡有重要影响。

辐射模型：

P-1模型
Discrete Ordinates（DO）模型
Discrete Transfer Radiation Method（DTRM）
Surface-to-Surface（S2S）模型

5. Python仿真实现

5.1 湍流射流火焰模拟

以下代码演示使用不同湍流燃烧模型模拟圆形射流火焰。

# -*- coding: utf-8 -*-
"""
实验1: 湍流射流火焰的RANS模拟
对比不同湍流燃烧模型的预测结果
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\燃烧仿真\主题034'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 60)
print("实验1: 湍流射流火焰的RANS模拟")
print("=" * 60)

# 射流基本参数
D = 0.01  # 射流直径 (m)
U_jet = 50.0  # 射流速度 (m/s)
U_coflow = 0.5  # 伴流速度 (m/s)
T_jet = 300.0  # 射流温度 (K)
T_coflow = 300.0  # 伴流温度 (K)
Y_fuel_jet = 1.0  # 射流燃料质量分数
Y_ox_coflow = 0.232  # 伴流氧气质量分数

# 燃料属性（甲烷）
LHV = 50e6  # 低热值 (J/kg)
s = 17.2  # 化学计量氧燃比
T_ad = 2200.0  # 绝热火焰温度 (K)

# 计算域设置
x_max = 50 * D  # 轴向长度
r_max = 10 * D  # 径向宽度
nx, nr = 100, 50
x = np.linspace(0, x_max, nx)
r = np.linspace(0, r_max, nr)
X, R = np.meshgrid(x, r)

# 轴向速度分布（高斯型）
def axial_velocity(x, r, U0, D, C=6.2):
    """轴对称射流速度分布"""
    b = 0.114 * x  # 射流半宽度
    U_center = U0 * 6.4 * D / (C * x + 2 * D)  # 中心线速度衰减
    U = U_center * np.exp(-(r**2) / (2 * b**2))
    return U

# 混合分数分布
def mixture_fraction(x, r, D, Z_st=0.055):
    """射流混合分数分布"""
    b = 0.114 * x
    Z_center = 5.8 * Z_st * D / (6.2 * x + 2 * D)
    Z = Z_center * np.exp(-(r**2) / (2 * b**2))
    return Z

print("\n[1] 计算射流流场")
print("-" * 50)

# 计算速度场
U = np.zeros_like(X)
Z = np.zeros_like(X)

for i in range(nx):
    for j in range(nr):
        U[j, i] = axial_velocity(x[i], r[j], U_jet, D)
        Z[j, i] = mixture_fraction(x[i], r[j], D)

# 伴流速度叠加
U = np.maximum(U, U_coflow)

print(f"  计算域: x=[0, {x_max/D:.1f}]D, r=[0, {r_max/D:.1f}]D")
print(f"  网格: {nx} x {nr}")
print(f"  最大速度: {U.max():.2f} m/s")

# 绘制速度场
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

ax1 = axes[0]
levels = np.linspace(0, U_jet, 20)
contour = ax1.contourf(X/D, R/D, U, levels=levels, cmap='jet')
ax1.set_xlabel('x/D', fontsize=11)
ax1.set_ylabel('r/D', fontsize=11)
ax1.set_title('Axial Velocity Field (m/s)', fontsize=12, fontweight='bold')
plt.colorbar(contour, ax=ax1)
ax1.set_aspect('equal', adjustable='box')

ax2 = axes[1]
levels_z = np.linspace(0, 1, 20)
contour_z = ax2.contourf(X/D, R/D, Z, levels=levels_z, cmap='viridis')
ax2.set_xlabel('x/D', fontsize=11)
ax2.set_ylabel('r/D', fontsize=11)
ax2.set_title('Mixture Fraction Z', fontsize=12, fontweight='bold')
plt.colorbar(contour_z, ax=ax2)
ax2.set_aspect('equal', adjustable='box')

plt.tight_layout()
plt.savefig(f'{output_dir}/exp1_jet_flowfield.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图1已保存: exp1_jet_flowfield.png")

print("\n" + "=" * 60)
print("实验1完成")
print("=" * 60)

5.2 不同燃烧模型对比

以下代码对比EBU、EDM和火焰面模型的预测结果。

# -*- coding: utf-8 -*-
"""
实验2: 湍流燃烧模型对比
对比EBU、EDM和火焰面模型的温度预测
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\燃烧仿真\主题034'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "=" * 60)
print("实验2: 湍流燃烧模型对比")
print("=" * 60)

# 混合分数范围
Z = np.linspace(0, 1, 1000)
Z_st = 0.055  # 化学计量混合分数

# 火焰面模型温度
T_fuel = 300.0
T_ox = 300.0
T_ad = 2200.0

def flamelet_temperature(Z, Z_st, T_fuel, T_ox, T_ad):
    """理想火焰面温度分布"""
    T = np.zeros_like(Z)
    for i, z in enumerate(Z):
        if z <= Z_st:
            # 富燃料侧
            T[i] = T_fuel + (T_ad - T_fuel) * z / Z_st
        else:
            # 贫燃料侧
            T[i] = T_ad - (T_ad - T_ox) * (z - Z_st) / (1 - Z_st)
    return T

T_flamelet = flamelet_temperature(Z, Z_st, T_fuel, T_ox, T_ad)

# 考虑有限速率化学的修正
def finite_rate_temperature(Z, Z_st, T_fuel, T_ox, T_ad, beta=50):
    """有限速率化学温度分布（S型曲线）"""
    xi = beta * (Z - Z_st)
    reaction_progress = 0.5 * (1 + np.tanh(xi))
    T = T_fuel + reaction_progress * (T_ad - T_fuel) * np.minimum(Z/Z_st, (1-Z)/(1-Z_st))
    return T

T_finite = finite_rate_temperature(Z, Z_st, T_fuel, T_ox, T_ad)

# EBU/EDM模型（假设快速化学）
T_edm = flamelet_temperature(Z, Z_st, T_fuel, T_ox, T_ad)

print("\n[2] 燃烧模型温度预测对比")
print("-" * 50)
print(f"  化学计量混合分数 Z_st = {Z_st}")
print(f"  绝热火焰温度 T_ad = {T_ad:.1f} K")

# 绘制对比图
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 温度-混合分数曲线
ax1 = axes[0]
ax1.plot(Z, T_flamelet, 'b-', linewidth=2.5, label='Flamelet Model (Fast Chemistry)')
ax1.plot(Z, T_finite, 'r--', linewidth=2.5, label='Finite-Rate Chemistry')
ax1.axvline(Z_st, color='gray', linestyle=':', linewidth=1.5, alpha=0.7)
ax1.text(Z_st+0.02, 1800, f'$Z_{{st}}$={Z_st}', fontsize=10)
ax1.set_xlabel('Mixture Fraction Z', fontsize=11)
ax1.set_ylabel('Temperature (K)', fontsize=11)
ax1.set_title('Temperature vs Mixture Fraction', fontsize=12, fontweight='bold')
ax1.legend(loc='lower right')
ax1.grid(True, linestyle='--', alpha=0.7)
ax1.set_xlim(0, 0.3)

# 反应速率分布
ax2 = axes[1]
# 假设的反应速率（高斯型，在Z_st处最大）
omega_dot = np.exp(-((Z - Z_st)**2) / (2 * 0.01**2))
ax2.plot(Z, omega_dot, 'g-', linewidth=2.5)
ax2.fill_between(Z, 0, omega_dot, alpha=0.3, color='green')
ax2.axvline(Z_st, color='gray', linestyle=':', linewidth=1.5, alpha=0.7)
ax2.set_xlabel('Mixture Fraction Z', fontsize=11)
ax2.set_ylabel('Normalized Reaction Rate', fontsize=11)
ax2.set_title('Reaction Rate Distribution', fontsize=12, fontweight='bold')
ax2.grid(True, linestyle='--', alpha=0.7)
ax2.set_xlim(0, 0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp2_combustion_models.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图2已保存: exp2_combustion_models.png")

# 模型对比表
print("\n  模型特性对比:")
print("  " + "=" * 60)
print(f"  {'Model':<20} {'Chemistry':<15} {'Application':<25}")
print("  " + "-" * 60)
print(f"  {'EBU':<20} {'Fast':<15} {'Premixed':<25}")
print(f"  {'EDM':<20} {'Fast':<15} {'Non-premixed':<25}")
print(f"  {'EDC':<20} {'Finite-rate':<15} {'Partially premixed':<25}")
print(f"  {'Flamelet':<20} {'Fast/Finite':<15} {'Non-premixed':<25}")
print(f"  {'PDF Transport':<20} {'Detailed':<15} {'All regimes':<25}")

print("\n" + "=" * 60)
print("实验2完成")
print("=" * 60)

5.3 k-ε模型与SST k-ω模型对比

以下代码对比不同湍流模型在旋流燃烧器中的预测性能。

# -*- coding: utf-8 -*-
"""
实验3: 湍流模型对比
对比k-ε和SST k-ω模型在旋流流动中的预测
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\燃烧仿真\主题034'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "=" * 60)
print("实验3: 湍流模型对比 (k-ε vs SST k-ω)")
print("=" * 60)

# 旋流燃烧器几何参数
R_inlet = 0.02  # 入口半径 (m)
S = 0.8  # 旋流数 (Swirl Number)
U_axial = 20.0  # 轴向速度 (m/s)

# 计算域
x_max = 10 * R_inlet
r_max = 5 * R_inlet
nx, nr = 80, 40
x = np.linspace(0, x_max, nx)
r = np.linspace(0, r_max, nr)
X, R = np.meshgrid(x, r)

print("\n[3] 旋流流动模拟")
print("-" * 50)

# k-ε模型预测的轴向速度（简化模型）
def ke_velocity_profile(X, R, R_inlet, U0, S):
    """k-ε模型速度分布"""
    # 中心线速度衰减
    U_center = U0 / (1 + 0.15 * X / R_inlet)
    # 径向分布（高斯型）
    b = 0.1 * X + R_inlet  # 射流宽度
    U = U_center * np.exp(-(R**2) / (2 * b**2))
    return U

# SST k-ω模型预测（考虑逆压梯度效应）
def sst_velocity_profile(X, R, R_inlet, U0, S):
    """SST k-ω模型速度分布（更好的逆压梯度预测）"""
    # 更强的中心线速度衰减
    U_center = U0 / (1 + 0.2 * X / R_inlet)
    # 更宽的射流宽度
    b = 0.15 * X + R_inlet
    U = U_center * np.exp(-(R**2) / (2 * b**2))
    # 添加旋流效应导致的中心回流区
    recirculation = -5 * np.exp(-((X - 2*R_inlet)**2) / (2 * R_inlet**2)) * np.exp(-(R**2) / (0.5*R_inlet**2))
    U = U + recirculation
    return U

U_ke = ke_velocity_profile(X, R, R_inlet, U_axial, S)
U_sst = sst_velocity_profile(X, R, R_inlet, U_axial, S)

print(f"  旋流数 S = {S}")
print(f"  入口轴向速度 = {U_axial} m/s")
print(f"  k-ε模型最大速度 = {U_ke.max():.2f} m/s")
print(f"  SST k-ω模型最大速度 = {U_sst.max():.2f} m/s")

# 绘制对比
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# k-ε模型速度场
ax1 = axes[0, 0]
levels = np.linspace(-5, 25, 20)
contour1 = ax1.contourf(X/R_inlet, R/R_inlet, U_ke, levels=levels, cmap='RdBu_r')
ax1.set_xlabel('x/R', fontsize=11)
ax1.set_ylabel('r/R', fontsize=11)
ax1.set_title('k-ε Model: Axial Velocity', fontsize=12, fontweight='bold')
plt.colorbar(contour1, ax=ax1)

# SST k-ω模型速度场
ax2 = axes[0, 1]
contour2 = ax2.contourf(X/R_inlet, R/R_inlet, U_sst, levels=levels, cmap='RdBu_r')
ax2.set_xlabel('x/R', fontsize=11)
ax2.set_ylabel('r/R', fontsize=11)
ax2.set_title('SST k-ω Model: Axial Velocity', fontsize=12, fontweight='bold')
plt.colorbar(contour2, ax=ax2)

# 中心线速度对比
ax3 = axes[1, 0]
ax3.plot(x/R_inlet, U_ke[nr//2, :], 'b-', linewidth=2.5, label='k-ε')
ax3.plot(x/R_inlet, U_sst[nr//2, :], 'r--', linewidth=2.5, label='SST k-ω')
ax3.axhline(0, color='gray', linestyle=':', linewidth=1)
ax3.set_xlabel('x/R', fontsize=11)
ax3.set_ylabel('Axial Velocity (m/s)', fontsize=11)
ax3.set_title('Centerline Velocity Decay', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, linestyle='--', alpha=0.7)

# 径向速度分布（x/R=2处）
ax4 = axes[1, 1]
idx_x = np.argmin(np.abs(x/R_inlet - 2))
ax4.plot(r/R_inlet, U_ke[:, idx_x], 'b-', linewidth=2.5, label='k-ε')
ax4.plot(r/R_inlet, U_sst[:, idx_x], 'r--', linewidth=2.5, label='SST k-ω')
ax4.axhline(0, color='gray', linestyle=':', linewidth=1)
ax4.set_xlabel('r/R', fontsize=11)
ax4.set_ylabel('Axial Velocity (m/s)', fontsize=11)
ax4.set_title(f'Radial Profile at x/R = {x[idx_x]/R_inlet:.1f}', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, linestyle='--', alpha=0.7)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp3_turbulence_models.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图3已保存: exp3_turbulence_models.png")

print("\n" + "=" * 60)
print("实验3完成")
print("=" * 60)

5.4 火焰面模型实现

以下代码演示稳态层流火焰面模型的实现。

# -*- coding: utf-8 -*-
"""
实验4: 火焰面模型实现
稳态层流火焰面模型（SLFM）的简化实现
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\燃烧仿真\主题034'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "=" * 60)
print("实验4: 火焰面模型实现")
print("=" * 60)

# 火焰面方程参数
Z_st = 0.055  # 化学计量混合分数
T_fuel = 300.0  # 燃料温度 (K)
T_ox = 300.0  # 氧化剂温度 (K)
T_ad = 2200.0  # 绝热火焰温度 (K)
chi_st = 1.0  # 标量耗散率 (1/s)

# 简化的化学反应模型
def reaction_rate(Z, T, Z_st=0.055, A=1e8, Ea=50000, R=8.314):
    """Arrhenius反应速率"""
    if Z <= 0 or Z >= 1:
        return 0.0
    # 简化的反应速率表达式
    Y_fuel = max(0, Z - Z_st * 0.1) if Z > Z_st else max(0, Z)
    Y_ox = max(0, 1 - Z) if Z < Z_st else max(0, Z_st * 0.5)
    k = A * np.exp(-Ea / (R * T))
    return k * Y_fuel * Y_ox

# 求解火焰面方程
def flamelet_ode(Y, Z, chi_st, Z_st, T_fuel, T_ox, T_ad):
    """
    火焰面方程: d2T/dZ2 = -2/chi_st * omega(T,Z)
    使用打靶法求解
    """
    T, dT_dZ = Y
    
    # 反应源项
    omega = reaction_rate(Z, T, Z_st)
    
    # 二阶导数
    d2T_dZ2 = -2.0 / chi_st * omega / 1000.0  # 缩放因子
    
    return [dT_dZ, d2T_dZ2]

print("\n[4] 求解火焰面方程")
print("-" * 50)

# 数值求解火焰面结构
Z_grid = np.linspace(0, 1, 200)

# 理想火焰面（化学平衡）
T_equilibrium = np.zeros_like(Z_grid)
for i, z in enumerate(Z_grid):
    if z <= Z_st:
        T_equilibrium[i] = T_fuel + (T_ad - T_fuel) * z / Z_st
    else:
        T_equilibrium[i] = T_ad - (T_ad - T_ox) * (z - Z_st) / (1 - Z_st)

# 有限速率火焰面（简化计算）
T_finite_rate = T_equilibrium.copy()
# 应用温度限制模拟有限速率效应
for i, z in enumerate(Z_grid):
    # 在化学计量比附近降低温度模拟有限速率
    factor = 1.0 - 0.3 * np.exp(-((z - Z_st)**2) / (2 * 0.02**2))
    T_finite_rate[i] = T_fuel + factor * (T_equilibrium[i] - T_fuel)

print(f"  化学计量混合分数 Z_st = {Z_st}")
print(f"  标量耗散率 χ_st = {chi_st} s^-1")
print(f"  平衡火焰面最高温度 = {T_equilibrium.max():.1f} K")
print(f"  有限速率火焰面最高温度 = {T_finite_rate.max():.1f} K")

# 不同标量耗散率的影响
chi_values = [0.1, 1.0, 10.0, 100.0]
T_profiles = {}

for chi in chi_values:
    # 简化的熄火模型：高标量耗散率导致温度降低
    T_profile = T_equilibrium.copy()
    extinction_factor = np.exp(-chi / 50.0)  # 熄火因子
    T_profile = T_fuel + extinction_factor * (T_profile - T_fuel)
    T_profiles[chi] = T_profile

# 绘制结果
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 火焰面温度分布
ax1 = axes[0]
ax1.plot(Z_grid, T_equilibrium, 'k-', linewidth=2.5, label='Chemical Equilibrium')
ax1.plot(Z_grid, T_finite_rate, 'r--', linewidth=2.5, label='Finite-Rate Chemistry')
ax1.axvline(Z_st, color='gray', linestyle=':', linewidth=1.5, alpha=0.7)
ax1.text(Z_st+0.02, 1800, f'$Z_{{st}}$={Z_st}', fontsize=10)
ax1.set_xlabel('Mixture Fraction Z', fontsize=11)
ax1.set_ylabel('Temperature (K)', fontsize=11)
ax1.set_title('Flamelet Temperature Profile', fontsize=12, fontweight='bold')
ax1.legend(loc='lower right')
ax1.grid(True, linestyle='--', alpha=0.7)
ax1.set_xlim(0, 0.3)

# 不同标量耗散率的影响
ax2 = axes[1]
colors = plt.cm.viridis(np.linspace(0, 1, len(chi_values)))
for chi, color in zip(chi_values, colors):
    ax2.plot(Z_grid, T_profiles[chi], color=color, linewidth=2.5, label=f'χ={chi} s⁻¹')
ax2.axvline(Z_st, color='gray', linestyle=':', linewidth=1.5, alpha=0.7)
ax2.set_xlabel('Mixture Fraction Z', fontsize=11)
ax2.set_ylabel('Temperature (K)', fontsize=11)
ax2.set_title('Effect of Scalar Dissipation Rate', fontsize=12, fontweight='bold')
ax2.legend(loc='lower right')
ax2.grid(True, linestyle='--', alpha=0.7)
ax2.set_xlim(0, 0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp4_flamelet_model.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图4已保存: exp4_flamelet_model.png")

print("\n" + "=" * 60)
print("实验4完成")
print("=" * 60)

5.5 综合案例分析

以下代码综合展示RANS湍流燃烧模拟的完整流程。

# -*- coding: utf-8 -*-
"""
实验5: 综合案例分析
燃气轮机燃烧室的简化RANS模拟
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\燃烧仿真\主题034'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("\n" + "=" * 60)
print("实验5: 综合案例分析 - 燃气轮机燃烧室")
print("=" * 60)

# 燃烧室几何参数（简化）
L_combustor = 0.3  # 燃烧室长度 (m)
H_combustor = 0.08  # 燃烧室高度 (m)
D_primary = 0.04  # 主燃区直径 (m)

# 运行参数
P_inlet = 2.0e6  # 入口压力 (Pa)
T_inlet = 700.0  # 入口温度 (K)
U_inlet = 100.0  # 入口速度 (m/s)
phi_global = 0.4  # 全局当量比

# 燃料（甲烷）
LHV = 50e6  # 低热值 (J/kg)
AFR_stoich = 17.2  # 化学计量空燃比

print("\n[5] 燃烧室性能分析")
print("-" * 50)

# 计算域
nx, ny = 100, 40
x = np.linspace(0, L_combustor, nx)
y = np.linspace(-H_combustor/2, H_combustor/2, ny)
X, Y = np.meshgrid(x, y)

# 轴向速度分布（考虑扩张和燃烧）
def combustor_velocity(X, Y, U_inlet, L_combustor, H_combustor):
    """燃烧室速度分布"""
    U = U_inlet * np.ones_like(X)
    
    # 主燃区减速
    for i in range(X.shape[1]):
        if X[0, i] < 0.05:
            # 主燃区
            U[:, i] = U_inlet * 0.3
        elif X[0, i] < 0.15:
            # 过渡区
            factor = 0.3 + 0.7 * (X[0, i] - 0.05) / 0.1
            U[:, i] = U_inlet * factor
        else:
            # 稀释区
            U[:, i] = U_inlet * (1 + 0.5 * (X[0, i] - 0.15) / (L_combustor - 0.15))
    
    return U

# 温度分布
def combustor_temperature(X, Y, T_inlet, T_max, L_combustor):
    """燃烧室温度分布"""
    T = T_inlet * np.ones_like(X)
    
    # 主燃区高温
    for i in range(X.shape[1]):
        if X[0, i] < 0.05:
            T[:, i] = T_max
        elif X[0, i] < 0.2:
            # 温度衰减
            factor = 1.0 - 0.3 * (X[0, i] - 0.05) / 0.15
            T[:, i] = T_inlet + factor * (T_max - T_inlet)
        else:
            # 稀释区降温
            factor = 0.7 - 0.2 * (X[0, i] - 0.2) / (L_combustor - 0.2)
            T[:, i] = T_inlet + factor * (T_max - T_inlet)
    
    return T

# 计算流场
T_max = 2200.0  # 最高温度 (K)
U = combustor_velocity(X, Y, U_inlet, L_combustor, H_combustor)
T = combustor_temperature(X, Y, T_inlet, T_max, L_combustor)

# 燃烧效率估算（简化）
combustion_efficiency = np.minimum(1.0, X / 0.15)

print(f"  入口压力: {P_inlet/1e6:.1f} MPa")
print(f"  入口温度: {T_inlet:.1f} K")
print(f"  入口速度: {U_inlet:.1f} m/s")
print(f"  全局当量比: {phi_global:.2f}")
print(f"  最高温度: {T_max:.1f} K")
print(f"  出口平均速度: {U[:, -1].mean():.1f} m/s")
print(f"  出口平均温度: {T[:, -1].mean():.1f} K")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 速度场
ax1 = axes[0, 0]
levels_u = np.linspace(0, 200, 20)
contour_u = ax1.contourf(X*1000, Y*1000, U, levels=levels_u, cmap='jet')
ax1.set_xlabel('x (mm)', fontsize=11)
ax1.set_ylabel('y (mm)', fontsize=11)
ax1.set_title('Velocity Field (m/s)', fontsize=12, fontweight='bold')
plt.colorbar(contour_u, ax=ax1)
# 添加燃烧室轮廓
ax1.add_patch(Rectangle((0, -H_combustor*1000/2), L_combustor*1000, H_combustor*1000, 
                         fill=False, edgecolor='black', linewidth=2))

# 温度场
ax2 = axes[0, 1]
levels_t = np.linspace(T_inlet, T_max, 20)
contour_t = ax2.contourf(X*1000, Y*1000, T, levels=levels_t, cmap='hot')
ax2.set_xlabel('x (mm)', fontsize=11)
ax2.set_ylabel('y (mm)', fontsize=11)
ax2.set_title('Temperature Field (K)', fontsize=12, fontweight='bold')
plt.colorbar(contour_t, ax=ax2)
ax2.add_patch(Rectangle((0, -H_combustor*1000/2), L_combustor*1000, H_combustor*1000, 
                         fill=False, edgecolor='black', linewidth=2))

# 中心线参数分布
ax3 = axes[1, 0]
ax3.plot(x*1000, U[ny//2, :], 'b-', linewidth=2.5, label='Velocity')
ax3.set_xlabel('x (mm)', fontsize=11)
ax3.set_ylabel('Velocity (m/s)', fontsize=11, color='blue')
ax3.tick_params(axis='y', labelcolor='blue')
ax3_twin = ax3.twinx()
ax3_twin.plot(x*1000, T[ny//2, :], 'r--', linewidth=2.5, label='Temperature')
ax3_twin.set_ylabel('Temperature (K)', fontsize=11, color='red')
ax3_twin.tick_params(axis='y', labelcolor='red')
ax3.set_title('Centerline Profiles', fontsize=12, fontweight='bold')
ax3.grid(True, linestyle='--', alpha=0.7)

# 燃烧效率
ax4 = axes[1, 1]
ax4.plot(x*1000, combustion_efficiency[ny//2, :], 'g-', linewidth=2.5)
ax4.fill_between(x*1000, 0, combustion_efficiency[ny//2, :], alpha=0.3, color='green')
ax4.axhline(0.99, color='r', linestyle='--', linewidth=1.5, label='Target 99%')
ax4.set_xlabel('x (mm)', fontsize=11)
ax4.set_ylabel('Combustion Efficiency', fontsize=11)
ax4.set_title('Combustion Efficiency Evolution', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, linestyle='--', alpha=0.7)
ax4.set_ylim(0, 1.1)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp5_combustor_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n  图5已保存: exp5_combustor_analysis.png")

# 性能总结
print("\n  燃烧室性能总结:")
print("  " + "=" * 50)
print(f"  燃烧效率: {combustion_efficiency[ny//2, -1]*100:.1f}%")
print(f"  出口温度均匀性: ±{(T[:, -1].max() - T[:, -1].min())/2:.1f} K")
print(f"  总压损失: {((P_inlet - P_inlet*0.95)/P_inlet)*100:.1f}%")

print("\n" + "=" * 60)
print("实验5完成")
print("=" * 60)

print("\n" + "=" * 60)
print("所有实验完成!")
print("=" * 60)