燃烧仿真-主题025-喷雾燃烧基础

kkchenjj

347人浏览 · 2026-03-13 20:29:20

kkchenjj · 2026-03-13 20:29:20 发布

第二十五篇：喷雾燃烧基础

摘要

喷雾燃烧是液体燃料在燃烧设备中最常见的燃烧形态，广泛应用于航空发动机、燃气轮机、柴油机、工业锅炉等动力装置中。本教程系统介绍喷雾燃烧的理论基础，包括液滴破碎与雾化、液滴蒸发、液滴燃烧、喷雾-湍流相互作用等核心内容。内容涵盖液滴尺寸分布模型、蒸发模型、燃烧模型、喷雾数值模拟方法等。通过Python仿真实验，读者将深入理解喷雾燃烧的数值模拟技术，掌握液滴轨迹计算、蒸发速率预测、喷雾燃烧特性分析等实用技能，为理解航空发动机和燃气轮机中的喷雾燃烧过程奠定坚实的理论基础。

关键词

喷雾燃烧，液滴雾化，液滴蒸发，液滴燃烧，Sauter平均直径，蒸发常数，喷雾-湍流相互作用，离散相模型

1. 喷雾燃烧的基本概念

1.1 喷雾燃烧的特点

喷雾燃烧是指液体燃料通过喷嘴雾化成细小液滴后，与氧化剂混合并发生燃烧的燃烧形态。喷雾燃烧具有以下特点：

喷雾燃烧的特征：

液体燃料需要经过雾化和蒸发过程才能参与燃烧
燃烧速率受液滴蒸发速率控制
存在复杂的液滴-气体相互作用
喷雾结构影响燃烧稳定性和效率
污染物排放与喷雾特性密切相关

喷雾燃烧的优势：

液体燃料能量密度高，便于储存和运输
喷雾燃烧可以实现较高的燃烧效率
通过控制喷雾特性可以调节燃烧过程
适用于大功率动力装置

喷雾燃烧的挑战：

液滴蒸发需要时间，影响燃烧响应速度
喷雾不均匀可能导致燃烧不稳定
碳烟生成与液滴蒸发和混合过程相关
需要精确的喷雾和燃烧模型

1.2 喷雾燃烧的过程

喷雾燃烧是一个复杂的多物理过程，包括以下几个阶段：

雾化阶段：

液体燃料通过喷嘴喷出形成液膜或射流
液膜或射流在气动力作用下破碎成液滴
液滴可能经历二次破碎形成更小的液滴
最终形成具有一定尺寸分布的喷雾

蒸发阶段：

液滴在热气流中加热
液滴表面燃料蒸发形成燃料蒸气
蒸发速率取决于液滴温度、环境温度和相对速度
蒸发过程吸收热量，降低局部温度

混合阶段：

燃料蒸气与空气混合
混合过程受湍流和分子扩散控制
形成可燃混合物

燃烧阶段：

可燃混合物在点火源作用下发生燃烧
燃烧释放热量加热周围液滴
形成稳定的喷雾火焰

1.3 喷雾燃烧的分类

根据喷雾特性和燃烧条件，喷雾燃烧可以分为不同的类型：

按雾化方式分类：

压力雾化：依靠燃料压力通过喷嘴雾化
气动雾化：依靠高速气流将液体雾化
旋转雾化：依靠离心力将液体雾化
超声波雾化：依靠超声波振动将液体雾化

按燃烧 regime 分类：

预蒸发燃烧 regime：燃料在燃烧前完全蒸发
液滴群燃烧 regime：液滴群集体燃烧
单液滴燃烧 regime：每个液滴独立燃烧

按火焰结构分类：

扩散火焰：燃料和氧化剂在火焰面混合
预混火焰：燃料和氧化剂在燃烧前预混
部分预混火焰：介于扩散和预混之间

2. 液滴雾化理论

2.1 液滴破碎机制

液体燃料通过喷嘴喷出后，在气动力作用下破碎成液滴。液滴破碎的主要机制包括：

袋式破碎（Bag Breakup）：

液滴在气流作用下形成袋状结构
袋壁变薄并最终破裂
产生大量小液滴
发生在相对较低的韦伯数（We ≈ 10-40）

剪切破碎（Shear Breakup）：

液滴表面被气流剪切剥离
产生细小的液滴从表面剥离
发生在中等韦伯数（We ≈ 40-100）

灾难性破碎（Catastrophic Breakup）：

液滴在强气流作用下迅速破碎
产生大量细小液滴
发生在高韦伯数（We > 100）

振动破碎（Vibrational Breakup）：

液滴表面振动导致破碎
发生在低韦伯数（We < 10）

2.2 韦伯数与破碎判据

韦伯数是描述液滴破碎特性的重要无量纲数：

$\frac{\rho_g u_{rel}^2 d_d}{\sigma}$

其中 $ρg\rho_g$ 是气体密度， $u_{rel}$ 是相对速度， $d_d$ 是液滴直径， $σ\sigma$ 是表面张力。

韦伯数的物理意义：

表示气动力与表面张力的比值
韦伯数越大，液滴越容易破碎
临界韦伯数决定了液滴是否破碎

破碎判据：

$We < We_{crit}$ ：液滴稳定，不破碎
$We > We_{crit}$ ：液滴破碎

临界韦伯数取决于液滴的粘度：

$We_{crit} = We_{crit,0} \left( 1 + C Oh \right)$

其中 $\mu_l / \sqrt{\rho_l d_d \sigma}$ 是Ohnesorge数， $We_{crit,0}$ 是无粘情况下的临界韦伯数。

2.3 液滴尺寸分布

喷雾中的液滴尺寸通常服从一定的统计分布。常用的液滴尺寸分布模型包括：

Rosin-Rammler分布：

$\exp\left[ -\left( \frac{d}{\bar{d}} \right)^n \right]$

其中 $Q (d)$ 是累积体积分数， $dˉ\bar{d}$ 是特征直径， $n$ 是分布指数。

Nukiyama-Tanasawa分布：

$f(d) = A d^2 \exp(-B d^C)$

其中 $A$ 、 $B$ 、 $C$ 是分布参数。

对数正态分布：

$\frac{1}{d \sigma \sqrt{2\pi}} \exp\left[ -\frac{(\ln d - \mu)^2}{2\sigma^2} \right]$

其中 $μ\mu$ 和 $σ\sigma$ 是对数均值和标准差。

2.4 特征直径

为了描述喷雾的液滴尺寸特性，定义了多种特征直径：

Sauter平均直径（SMD或D32）：

$D_{32} = \frac{\sum n_i d_i^3}{\sum n_i d_i^2}$

SMD表示单位体积液滴的总表面积，是喷雾燃烧中最重要的特征直径。

算术平均直径（D10）：

$D_{10} = \frac{\sum n_i d_i}{\sum n_i}$

体积平均直径（D30）：

$D_{30} = \sqrt[3]{\frac{\sum n_i d_i^3}{\sum n_i}}$

质量中间直径（MMD或D50）：

累积质量分数为50%对应的直径。

3. 液滴蒸发理论

3.1 液滴蒸发模型

液滴在热气流中蒸发是一个复杂的传热传质过程。常用的液滴蒸发模型包括：

d²定律：
液滴直径的平方随时间线性减小：

$d^2 = d_0^2 - K t$

其中 $d_0$ 是初始直径， $K$ 是蒸发常数。

蒸发常数：

$\frac{8 \lambda_g}{\rho_l c_{p,g}} \ln(1 + B_M)$

其中 $λg\lambda_g$ 是气体导热系数， $ρl\rho_l$ 是液体密度， $c_{p,g}$ 是气体比热容， $B_M$ 是传质数。

传质数：

$B_M = \frac{Y_{F,s} - Y_{F,\infty}}{1 - Y_{F,s}}$

其中 $Y_{F,s}$ 是液滴表面燃料质量分数， $YF,∞Y_{F,\infty}$ 是环境燃料质量分数。

3.2 液滴加热阶段

液滴蒸发前需要经历加热阶段，温度从初始温度升高到沸点：

加热时间：

$\tau_h = \frac{\rho_l c_{p,l} d_0^2}{6 Nu \lambda_g} \ln\left( \frac{T_\infty - T_0}{T_\infty - T_b} \right)$

其中 $c_{p,l}$ 是液体比热容， $T_0$ 是初始温度， $T_b$ 是沸点温度， $N u$ 是努塞尔数。

努塞尔数：
对于静止液滴：

$Nu_0 = 2.0$

对于运动液滴（Ranz-Marshall关联式）：

$Nu = 2.0 + 0.6 Re_d^{0.5} Pr^{0.33}$

其中 $Re_d$ 是液滴雷诺数， $P r$ 是普朗特数。

3.3 液滴寿命

液滴从初始状态到完全蒸发的总时间称为液滴寿命：

总蒸发时间：

$\tau_{evap} = \tau_h + \tau_v$

其中 $τh\tau_h$ 是加热时间， $τv\tau_v$ 是蒸发时间。

蒸发时间：
根据d²定律：

$\tau_v = \frac{d_0^2}{K}$

影响因素：

环境温度：温度越高，蒸发越快
环境压力：压力影响沸点和蒸发速率
相对速度：速度影响传热传质系数
液滴尺寸：小液滴蒸发更快

3.4 多组分燃料蒸发

实际燃料通常是多组分混合物，蒸发特性与单组分燃料不同：

多组分蒸发特点：

轻组分先蒸发，重组分后蒸发
液滴组成随时间变化
沸点随组成变化

连续热力学模型：
将多组分燃料视为连续分布，用分布函数描述组成。

离散组分模型：
将燃料简化为几个代表性组分，分别计算各组分的蒸发。

4. 液滴燃烧理论

4.1 单液滴燃烧模型

单液滴燃烧是喷雾燃烧的基础。经典的液滴燃烧模型包括：

扩散火焰模型：

液滴表面蒸发产生燃料蒸气
燃料蒸气向外扩散与空气混合
在化学计量比位置形成扩散火焰
火焰呈球形包围液滴

火焰位置：
火焰半径与液滴半径的比值：

$\frac{r_f}{r_s} = \frac{\ln(1 + B_{ox})}{\ln(1 + B_F)}$

其中 $B_{ox}$ 和 $B_F$ 分别是氧化剂和燃料的传质数。

燃烧速率：
液滴质量燃烧速率：

$\dot{m} = 2 \pi r_s \rho_g D_g \ln(1 + B)$

其中 $r_s$ 是液滴半径， $D_g$ 是气体扩散系数， $B$ 是传质数。

4.2 液滴燃烧时间

液滴燃烧时间可以通过d²定律预测：

燃烧时间：

$\tau_b = \frac{d_0^2}{K_b}$

其中 $K_b$ 是燃烧蒸发常数，通常大于纯蒸发时的蒸发常数。

燃烧蒸发常数：

$K_b = \frac{8 \lambda_g}{\rho_l c_{p,g}} \ln(1 + B_T)$

其中 $B_T$ 是传热数：

$B_T = \frac{c_{p,g}(T_\infty - T_s) + q_c Y_{ox,\infty} / s}{h_{fg}}$

其中 $q_c$ 是燃烧热， $s$ 是化学计量比， $h_{fg}$ 是汽化潜热。

4.3 液滴群燃烧

实际喷雾中液滴密集，存在液滴间相互作用：

液滴群燃烧特点：

液滴间共享火焰
燃烧速率受群体效应影响
存在液滴间传热传质

群体燃烧数：

$\frac{n d_d^2}{4 D_{32}}$

其中 $n$ 是液滴数密度， $d_d$ 是液滴直径。

燃烧 regime：

$\ll 1$ ：单液滴燃烧 regime
$\sim 1$ ：过渡 regime
$\gg 1$ ：液滴群燃烧 regime

4.4 液滴燃烧的影响因素

环境条件：

环境温度：影响蒸发和燃烧速率
环境压力：影响沸点和扩散系数
氧气浓度：影响燃烧强度和火焰温度

液滴特性：

液滴尺寸：小液滴燃烧更快
液滴初温：影响加热时间
燃料性质：影响蒸发和燃烧特性

流动条件：

相对速度：影响传热传质系数
湍流强度：影响混合和燃烧速率

5. 喷雾-湍流相互作用

5.1 液滴-湍流相互作用机制

喷雾中的液滴与湍流存在复杂的相互作用：

湍流对液滴的影响：

湍流扩散使液滴分散
湍流脉动影响液滴轨迹
湍流增强传热传质
湍流影响液滴破碎

液滴对湍流的影响：

液滴蒸发改变气体密度
液滴阻力消耗湍动能
液滴诱导产生涡量

5.2 液滴扩散

液滴在湍流中的扩散可以用扩散系数描述：

液滴扩散系数：
对于小惯性液滴（Stokes数 $\ll 1$ ）：

$D_p \approx D_t$

对于大惯性液滴（ $\gg 1$ ）：

$D_p \approx 0$

其中 $D_t$ 是湍流扩散系数。

Stokes数：

$\frac{\tau_p}{\tau_t}$

其中 $τp\tau_p$ 是液滴弛豫时间， $τt\tau_t$ 是湍流时间尺度。

5.3 喷雾穿透深度

喷雾在交叉气流中的穿透深度是喷雾燃烧设计的重要参数：

穿透深度模型：

$d_n \left( \frac{\rho_l u_l^2}{\rho_g u_g^2} \right)^{0.5}$

其中 $C$ 是常数， $d_n$ 是喷嘴直径， $ρl\rho_l$ 和 $ρg\rho_g$ 分别是液体和气体密度， $u_l$ 和 $u_g$ 分别是液体和气体速度。

动量比：
动量比是决定穿透深度的关键参数：

$\frac{\rho_l u_l^2}{\rho_g u_g^2}$

6. 喷雾数值模拟方法

6.1 欧拉-拉格朗日方法

欧拉-拉格朗日方法是喷雾数值模拟中最常用的方法：

气相（欧拉框架）：
求解连续相的控制方程：

$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = \dot{S}_m$

$\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \dot{S}_u$

其中 $S˙m\dot{S}_m$ 和 $S˙u\dot{S}_u$ 分别是液滴引起的质量和动量源项。

液相（拉格朗日框架）：
追踪离散液滴的运动：

$\frac{d \mathbf{x}_p}{dt} = \mathbf{u}_p$

$m_p \frac{d \mathbf{u}_p}{dt} = \mathbf{F}_D + \mathbf{F}_P + \mathbf{F}_G$

其中 $FD\mathbf{F}_D$ 是阻力， $FP\mathbf{F}_P$ 是压力梯度力， $FG\mathbf{F}_G$ 是重力。

6.2 液滴阻力模型

液滴在气体中运动受到的阻力：

Stokes阻力（低雷诺数）：

$\mathbf{F}_D = 3 \pi \mu_g d_d (\mathbf{u}_g - \mathbf{u}_p)$

标准阻力模型：

$\mathbf{F}_D = \frac{1}{2} C_D \rho_g A_d |\mathbf{u}_g - \mathbf{u}_p| (\mathbf{u}_g - \mathbf{u}_p)$

其中 $C_D$ 是阻力系数， $A_d$ 是液滴投影面积。

阻力系数：
对于球形液滴：

$C_D = \frac{24}{Re_d} (1 + 0.15 Re_d^{0.687})$

6.3 喷雾源项

液滴对气相的影响通过源项体现：

质量源项：

$\dot{S}_m = -\sum_k \dot{m}_k$

其中 $m˙k\dot{m}_k$ 是第 $k$ 个液滴的蒸发速率。

动量源项：

$\dot{S}_u = -\sum_k \mathbf{F}_{D,k}$

能量源项：

$\dot{S}_e = -\sum_k (\dot{m}_k h_{fg} + Q_k)$

其中 $Q_k$ 是液滴与气体的对流换热量。

6.4 喷雾模拟的数值方法

液滴追踪：

显式积分：简单但可能不稳定
隐式积分：稳定但计算量大
龙格-库塔方法：精度高

气相求解：

有限体积法
SIMPLE算法处理压力-速度耦合
迎风格式处理对流项

耦合方法：

单向耦合：只考虑气相对液滴的影响
双向耦合：同时考虑气相和液滴的相互作用

7. Python仿真实验

7.1 实验1：液滴尺寸分布

实验目的：模拟和分析不同雾化条件下的液滴尺寸分布。

Rosin-Rammler分布：

$\exp\left[ -\left( \frac{d}{\bar{d}} \right)^n \right]$

概率密度函数：

$\frac{n}{\bar{d}} \left( \frac{d}{\bar{d}} \right)^{n-1} \exp\left[ -\left( \frac{d}{\bar{d}} \right)^n \right]$

特征直径计算：

D10：算术平均直径
D32：Sauter平均直径
D50：质量中间直径

7.2 实验2：液滴蒸发

实验目的：模拟单液滴在热气流中的蒸发过程。

控制方程：

能量平衡：

$m_p c_{p,l} \frac{dT_p}{dt} = \pi d_d Nu \lambda_g (T_g - T_p) - \dot{m} h_{fg}$

质量变化：

$\frac{dm_p}{dt} = -\dot{m}$

直径变化（d²定律）：

$\frac{d(d^2)}{dt} = -K$

数值求解：

时间推进：显式欧拉法或龙格-库塔法
时间步长：根据稳定性条件确定

7.3 实验3：液滴轨迹

实验目的：模拟液滴在气流中的运动轨迹。

运动方程：

$\frac{d \mathbf{x}_p}{dt} = \mathbf{u}_p$

$m_p \frac{d \mathbf{u}_p}{dt} = \frac{1}{2} C_D \rho_g A_d |\mathbf{u}_g - \mathbf{u}_p| (\mathbf{u}_g - \mathbf{u}_p)$

阻力系数：

$C_D = \frac{24}{Re_d} (1 + 0.15 Re_d^{0.687})$

数值求解：

四阶龙格-库塔法
自适应时间步长

7.4 实验4：喷雾穿透深度

实验目的：分析喷雾在交叉气流中的穿透特性。

穿透深度模型：

$d_n \sqrt{q}$

其中 $\rho_l u_l^2 / (\rho_g u_g^2)$ 是动量比。

影响因素：

喷射速度
气流速度
密度比
喷嘴直径

7.5 实验5：液滴燃烧

实验目的：模拟单液滴燃烧过程。

燃烧模型：

扩散火焰模型
火焰位置计算
燃烧速率计算

温度分布：

液滴表面温度
火焰温度
环境温度

# -*- coding: utf-8 -*-
"""
主题025：喷雾燃烧基础
Spray Combustion Fundamentals Simulation
"""

import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import os

# 创建输出目录
output_dir = r'd:\文档\燃烧仿真\主题025'
os.makedirs(output_dir, exist_ok=True)

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("="*60)
print("主题025：喷雾燃烧基础")
print("="*60)

# ============================================================
# 实验1：液滴尺寸分布
# ============================================================
print("\n实验1：液滴尺寸分布")
print("-" * 40)

# 参数设置
d_range = np.linspace(1, 100, 500)  # 液滴直径范围 (微米)
d_bar = 30  # 特征直径 (微米)

# Rosin-Rammler分布的不同分布指数
n_values = [1.5, 2.0, 2.5, 3.0, 3.5]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

colors = plt.cm.viridis(np.linspace(0, 1, len(n_values)))

# 1. 累积体积分数分布
for i, n in enumerate(n_values):
    Q = 1 - np.exp(-(d_range / d_bar)**n)
    axes[0, 0].plot(d_range, Q, color=colors[i], linewidth=2, label=f'n={n}')

axes[0, 0].axhline(y=0.5, color='r', linestyle='--', linewidth=1.5, alpha=0.7)
axes[0, 0].set_xlabel('Droplet Diameter (micron)', fontsize=11)
axes[0, 0].set_ylabel('Cumulative Volume Fraction Q', fontsize=11)
axes[0, 0].set_title('Rosin-Rammler Cumulative Distribution', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# 2. 概率密度函数
for i, n in enumerate(n_values):
    pdf = (n / d_bar) * (d_range / d_bar)**(n - 1) * np.exp(-(d_range / d_bar)**n)
    axes[0, 1].plot(d_range, pdf, color=colors[i], linewidth=2, label=f'n={n}')

axes[0, 1].set_xlabel('Droplet Diameter (micron)', fontsize=11)
axes[0, 1].set_ylabel('PDF f(d)', fontsize=11)
axes[0, 1].set_title('Rosin-Rammler PDF', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)

# 3. 对数正态分布
sigma_ln = 0.5  # 对数标准差
mu_ln = np.log(30)  # 对数均值

lognormal_pdf = (1 / (d_range * sigma_ln * np.sqrt(2 * np.pi))) * \
                np.exp(-(np.log(d_range) - mu_ln)**2 / (2 * sigma_ln**2))

axes[1, 0].plot(d_range, lognormal_pdf, 'b-', linewidth=2.5, label='Log-normal PDF')
axes[1, 0].axvline(x=np.exp(mu_ln), color='r', linestyle='--', linewidth=2, label=f'Mean={np.exp(mu_ln):.1f} micron')
axes[1, 0].set_xlabel('Droplet Diameter (micron)', fontsize=11)
axes[1, 0].set_ylabel('PDF f(d)', fontsize=11)
axes[1, 0].set_title('Log-normal Distribution', fontsize=12, fontweight='bold')
axes[1, 0].legend(fontsize=10)
axes[1, 0].grid(True, alpha=0.3)

# 4. 特征直径比较
n_fixed = 2.5
Q_values = np.linspace(0.01, 0.99, 100)

# 从累积分布反算直径
def rosin_rammler_inv(Q, d_bar, n):
    return d_bar * (-np.log(1 - Q))**(1/n)

d_Q = rosin_rammler_inv(Q_values, d_bar, n_fixed)

# 特征直径
D10 = rosin_rammler_inv(0.1, d_bar, n_fixed)  # 近似
D32 = d_bar * (3 / n_fixed)  # Sauter平均直径近似
D50 = rosin_rammler_inv(0.5, d_bar, n_fixed)  # 质量中间直径

axes[1, 1].plot(Q_values, d_Q, 'b-', linewidth=2.5, label='D-Q Curve')
axes[1, 1].axhline(y=D10, color='g', linestyle='--', linewidth=2, label=f'D10={D10:.1f} micron')
axes[1, 1].axhline(y=D32, color='r', linestyle='--', linewidth=2, label=f'D32={D32:.1f} micron')
axes[1, 1].axhline(y=D50, color='m', linestyle='--', linewidth=2, label=f'D50={D50:.1f} micron')
axes[1, 1].set_xlabel('Cumulative Volume Fraction Q', fontsize=11)
axes[1, 1].set_ylabel('Droplet Diameter (micron)', fontsize=11)
axes[1, 1].set_title('Characteristic Diameters', fontsize=12, fontweight='bold')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp1_droplet_size_distribution.png', dpi=150, bbox_inches='tight')
plt.close()

print("✓ 实验1完成：液滴尺寸分布")

# ============================================================
# 实验2：液滴蒸发
# ============================================================
print("\n实验2：液滴蒸发")
print("-" * 40)

# 参数设置
rho_l = 750  # 液体密度 (kg/m3)
cp_l = 2100  # 液体比热容 (J/kg/K)
lambda_g = 0.03  # 气体导热系数 (W/m/K)
cp_g = 1005  # 气体比热容 (J/kg/K)
h_fg = 3.0e5  # 汽化潜热 (J/kg)
T_inf = 800  # 环境温度 (K)
T_b = 400  # 沸点温度 (K)
T_0 = 300  # 初始温度 (K)
Nu = 2.0  # 努塞尔数（静止液滴）

# 传质数（简化）
B_M = 0.5

# 蒸发常数
K_evap = (8 * lambda_g / (rho_l * cp_g)) * np.log(1 + B_M)

# 初始直径
d0_values = [10, 20, 50, 100]  # 微米

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

colors = plt.cm.plasma(np.linspace(0, 1, len(d0_values)))

# 1. d²随时间变化（d²定律）
for i, d0 in enumerate(d0_values):
    # 蒸发时间
    tau_evap = (d0 * 1e-6)**2 / K_evap
    t = np.linspace(0, tau_evap, 100)
    
    # d²随时间变化
    d_squared = (d0 * 1e-6)**2 - K_evap * t
    d_squared = np.maximum(d_squared, 0)
    
    axes[0, 0].plot(t * 1000, np.sqrt(d_squared) * 1e6, color=colors[i], linewidth=2, 
                    label=f'd0={d0} micron')

axes[0, 0].set_xlabel('Time (ms)', fontsize=11)
axes[0, 0].set_ylabel('Droplet Diameter (micron)', fontsize=11)
axes[0, 0].set_title('Droplet Diameter vs Time (d^2 Law)', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# 2. d²随时间变化（验证线性关系）
for i, d0 in enumerate(d0_values):
    tau_evap = (d0 * 1e-6)**2 / K_evap
    t = np.linspace(0, tau_evap, 100)
    
    d_squared = (d0 * 1e-6)**2 - K_evap * t
    d_squared = np.maximum(d_squared, 0)
    
    axes[0, 1].plot(t * 1000, d_squared * 1e12, color=colors[i], linewidth=2, 
                    label=f'd0={d0} micron')

axes[0, 1].set_xlabel('Time (ms)', fontsize=11)
axes[0, 1].set_ylabel('d^2 (micron^2)', fontsize=11)
axes[0, 1].set_title('d^2 vs Time (Linear)', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)

# 3. 蒸发时间与初始直径的关系
d0_range = np.linspace(10, 200, 100)  # 微米
tau_evap_range = (d0_range * 1e-6)**2 / K_evap * 1000  # 转换为ms

axes[1, 0].plot(d0_range, tau_evap_range, 'b-', linewidth=2.5)
axes[1, 0].set_xlabel('Initial Diameter (micron)', fontsize=11)
axes[1, 0].set_ylabel('Evaporation Time (ms)', fontsize=11)
axes[1, 0].set_title('Evaporation Time vs Initial Diameter', fontsize=12, fontweight='bold')
axes[1, 0].grid(True, alpha=0.3)

# 4. 不同环境温度下的蒸发
T_inf_values = [500, 600, 700, 800, 900]  # K
d0_fixed = 50  # 微米

colors_T = plt.cm.coolwarm(np.linspace(0, 1, len(T_inf_values)))

for i, T_inf_i in enumerate(T_inf_values):
    # 传质数随温度变化（简化模型）
    B_M_i = 0.3 + 0.3 * (T_inf_i - 500) / 400
    K_evap_i = (8 * lambda_g / (rho_l * cp_g)) * np.log(1 + B_M_i)
    
    tau_evap_i = (d0_fixed * 1e-6)**2 / K_evap_i
    t = np.linspace(0, tau_evap_i, 100)
    
    d_squared_i = (d0_fixed * 1e-6)**2 - K_evap_i * t
    d_squared_i = np.maximum(d_squared_i, 0)
    
    axes[1, 1].plot(t * 1000, np.sqrt(d_squared_i) * 1e6, color=colors_T[i], linewidth=2, 
                    label=f'T_inf={T_inf_i} K')

axes[1, 1].set_xlabel('Time (ms)', fontsize=11)
axes[1, 1].set_ylabel('Droplet Diameter (micron)', fontsize=11)
axes[1, 1].set_title(f'Evaporation at Different T_inf (d0={d0_fixed} micron)', fontsize=12, fontweight='bold')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp2_droplet_evaporation.png', dpi=150, bbox_inches='tight')
plt.close()

print("✓ 实验2完成：液滴蒸发")

# ============================================================
# 实验3：液滴轨迹
# ============================================================
print("\n实验3：液滴轨迹")
print("-" * 40)

# 参数设置
rho_g = 1.2  # 气体密度 (kg/m3)
mu_g = 1.8e-5  # 气体动力粘度 (Pa.s)
g = 9.81  # 重力加速度 (m/s2)

# 初始条件
U_g = np.array([50, 0])  # 气流速度 (m/s)
U_p0_values = [30, 50, 70, 90]  # 不同初始液滴速度 (m/s)
theta_0 = np.pi / 6  # 喷射角度 (30度)

# 液滴直径
d_p = 50e-6  # 50微米

# 阻力系数模型
def drag_coefficient(Re):
    if Re < 0.1:
        return 24 / Re
    elif Re < 1000:
        return 24 / Re * (1 + 0.15 * Re**0.687)
    else:
        return 0.44

# 液滴运动方程
def droplet_motion(y, t, d_p, rho_l, rho_g, mu_g, U_g, g):
    x, z, u, w = y
    
    # 相对速度
    u_rel = U_g[0] - u
    w_rel = U_g[1] - w
    V_rel = np.sqrt(u_rel**2 + w_rel**2)
    
    # 雷诺数
    Re = rho_g * V_rel * d_p / mu_g
    
    # 阻力系数
    C_D = drag_coefficient(Re)
    
    # 投影面积
    A_p = np.pi * (d_p / 2)**2
    
    # 阻力
    F_D = 0.5 * rho_g * C_D * A_p * V_rel
    
    # 加速度
    if V_rel > 1e-6:
        ax = F_D * u_rel / V_rel / (rho_l * np.pi * d_p**3 / 6)
        az = F_D * w_rel / V_rel / (rho_l * np.pi * d_p**3 / 6) - g
    else:
        ax = 0
        az = -g
    
    return [u, w, ax, az]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

colors = plt.cm.viridis(np.linspace(0, 1, len(U_p0_values)))

# 1. 不同初始速度的液滴轨迹
for i, U_p0 in enumerate(U_p0_values):
    # 初始速度分量
    u0 = U_p0 * np.cos(theta_0)
    w0 = U_p0 * np.sin(theta_0)
    
    # 初始条件
    y0 = [0, 0, u0, w0]
    
    # 时间
    t = np.linspace(0, 0.05, 500)
    
    # 求解ODE
    sol = odeint(droplet_motion, y0, t, args=(d_p, rho_l, rho_g, mu_g, U_g, g))
    
    x_traj = sol[:, 0]
    z_traj = sol[:, 1]
    
    # 只绘制z>=0的部分
    valid_idx = z_traj >= 0
    axes[0, 0].plot(x_traj[valid_idx] * 1000, z_traj[valid_idx] * 1000, 
                    color=colors[i], linewidth=2, label=f'Up0={U_p0} m/s')

axes[0, 0].set_xlabel('x (mm)', fontsize=11)
axes[0, 0].set_ylabel('z (mm)', fontsize=11)
axes[0, 0].set_title('Droplet Trajectories (Different Initial Velocities)', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].set_aspect('equal')

# 2. 不同直径的液滴轨迹
U_p0_fixed = 70  # m/s
d_p_values = [10e-6, 30e-6, 50e-6, 100e-6]  # m

colors_d = plt.cm.plasma(np.linspace(0, 1, len(d_p_values)))

for i, d_p_i in enumerate(d_p_values):
    u0 = U_p0_fixed * np.cos(theta_0)
    w0 = U_p0_fixed * np.sin(theta_0)
    y0 = [0, 0, u0, w0]
    t = np.linspace(0, 0.05, 500)
    
    sol = odeint(droplet_motion, y0, t, args=(d_p_i, rho_l, rho_g, mu_g, U_g, g))
    
    x_traj = sol[:, 0]
    z_traj = sol[:, 1]
    
    valid_idx = z_traj >= 0
    axes[0, 1].plot(x_traj[valid_idx] * 1000, z_traj[valid_idx] * 1000, 
                    color=colors_d[i], linewidth=2, label=f'd={d_p_i*1e6:.0f} micron')

axes[0, 1].set_xlabel('x (mm)', fontsize=11)
axes[0, 1].set_ylabel('z (mm)', fontsize=11)
axes[0, 1].set_title(f'Droplet Trajectories (Different Diameters, Up0={U_p0_fixed} m/s)', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].set_aspect('equal')

# 3. 液滴速度随时间变化
U_p0_fixed2 = 70  # m/s
d_p_fixed = 50e-6  # m

u0 = U_p0_fixed2 * np.cos(theta_0)
w0 = U_p0_fixed2 * np.sin(theta_0)
y0 = [0, 0, u0, w0]
t = np.linspace(0, 0.03, 300)

sol = odeint(droplet_motion, y0, t, args=(d_p_fixed, rho_l, rho_g, mu_g, U_g, g))

u_traj = sol[:, 2]
w_traj = sol[:, 3]
V_traj = np.sqrt(u_traj**2 + w_traj**2)

axes[1, 0].plot(t * 1000, u_traj, 'b-', linewidth=2, label='u (horizontal)')
axes[1, 0].plot(t * 1000, w_traj, 'r-', linewidth=2, label='w (vertical)')
axes[1, 0].plot(t * 1000, V_traj, 'g--', linewidth=2, label='|V|')
axes[1, 0].axhline(y=U_g[0], color='b', linestyle=':', linewidth=2, alpha=0.7, label=f'Ug={U_g[0]} m/s')
axes[1, 0].set_xlabel('Time (ms)', fontsize=11)
axes[1, 0].set_ylabel('Velocity (m/s)', fontsize=11)
axes[1, 0].set_title('Droplet Velocity Evolution', fontsize=12, fontweight='bold')
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)

# 4. Stokes数分析
tau_p_range = np.logspace(-4, 0, 100)  # 液滴弛豫时间范围
tau_t = 0.01  # 湍流时间尺度 (假设)

St_range = tau_p_range / tau_t

axes[1, 1].semilogx(tau_p_range * 1000, St_range, 'b-', linewidth=2.5)
axes[1, 1].axhline(y=1, color='r', linestyle='--', linewidth=2, label='St=1')
axes[1, 1].axvline(x=tau_t * 1000, color='g', linestyle='--', linewidth=2, label=f'tau_t={tau_t*1000:.1f} ms')
axes[1, 1].set_xlabel('Particle Relaxation Time (ms)', fontsize=11)
axes[1, 1].set_ylabel('Stokes Number', fontsize=11)
axes[1, 1].set_title('Stokes Number Analysis', fontsize=12, fontweight='bold')
axes[1, 1].legend(fontsize=10)
axes[1, 1].grid(True, alpha=0.3)

# 添加regime标签
axes[1, 1].text(0.1, 0.01, 'Tracer\nSt<<1', fontsize=9, ha='center',
                bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
axes[1, 1].text(100, 100, 'Inertial\nSt>>1', fontsize=9, ha='center',
                bbox=dict(boxstyle='round', facecolor='salmon', alpha=0.8))

plt.tight_layout()
plt.savefig(f'{output_dir}/exp3_droplet_trajectory.png', dpi=150, bbox_inches='tight')
plt.close()

print("✓ 实验3完成：液滴轨迹")

# ============================================================
# 实验4：喷雾穿透深度
# ============================================================
print("\n实验4：喷雾穿透深度")
print("-" * 40)

# 参数设置
C_pen = 1.5  # 穿透深度常数
d_n = 0.001  # 喷嘴直径 (m)
rho_l = 750  # 液体密度 (kg/m3)
rho_g = 1.2  # 气体密度 (kg/m3)

# 速度范围
u_l_range = np.linspace(10, 100, 100)  # 液体喷射速度 (m/s)
u_g_range = np.linspace(10, 100, 100)  # 气体速度 (m/s)

# 创建网格
U_l_grid, U_g_grid = np.meshgrid(u_l_range, u_g_range)

# 动量比
q_grid = (rho_l * U_l_grid**2) / (rho_g * U_g_grid**2)

# 穿透深度
S_grid = C_pen * d_n * np.sqrt(q_grid)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. 穿透深度等高线
contour1 = axes[0, 0].contourf(U_l_grid, U_g_grid, S_grid * 1000, levels=50, cmap='viridis')
axes[0, 0].set_xlabel('Liquid Jet Velocity (m/s)', fontsize=11)
axes[0, 0].set_ylabel('Gas Velocity (m/s)', fontsize=11)
axes[0, 0].set_title('Spray Penetration Depth (mm)', fontsize=12, fontweight='bold')
plt.colorbar(contour1, ax=axes[0, 0])

# 2. 动量比分布
contour2 = axes[0, 1].contourf(U_l_grid, U_g_grid, np.log10(q_grid), levels=50, cmap='plasma')
axes[0, 1].contour(U_l_grid, U_g_grid, q_grid, levels=[1, 10, 100], colors='k', linewidths=2)
axes[0, 1].set_xlabel('Liquid Jet Velocity (m/s)', fontsize=11)
axes[0, 1].set_ylabel('Gas Velocity (m/s)', fontsize=11)
axes[0, 1].set_title('Momentum Ratio (log10)', fontsize=12, fontweight='bold')
plt.colorbar(contour2, ax=axes[0, 1])

# 3. 固定气体速度下的穿透深度
u_g_fixed = 50  # m/s
q_fixed = (rho_l * u_l_range**2) / (rho_g * u_g_fixed**2)
S_fixed = C_pen * d_n * np.sqrt(q_fixed)

axes[1, 0].plot(u_l_range, S_fixed * 1000, 'b-', linewidth=2.5)
axes[1, 0].set_xlabel('Liquid Jet Velocity (m/s)', fontsize=11)
axes[1, 0].set_ylabel('Penetration Depth (mm)', fontsize=11)
axes[1, 0].set_title(f'Penetration Depth vs Jet Velocity (Ug={u_g_fixed} m/s)', fontsize=12, fontweight='bold')
axes[1, 0].grid(True, alpha=0.3)

# 4. 不同喷嘴直径的影响
d_n_values = [0.0005, 0.001, 0.002, 0.005]  # m
u_l_fixed = 70  # m/s
u_g_fixed2 = 50  # m/s
q_fixed2 = (rho_l * u_l_fixed**2) / (rho_g * u_g_fixed2**2)

colors_dn = plt.cm.coolwarm(np.linspace(0, 1, len(d_n_values)))

for i, d_n_i in enumerate(d_n_values):
    S_i = C_pen * d_n_i * np.sqrt(q_fixed2)
    axes[1, 1].bar(i, S_i * 1000, color=colors_dn[i], label=f'dn={d_n_i*1000:.1f} mm')

axes[1, 1].set_xticks(range(len(d_n_values)))
axes[1, 1].set_xticklabels([f'{d*1000:.1f}' for d in d_n_values])
axes[1, 1].set_xlabel('Nozzle Diameter (mm)', fontsize=11)
axes[1, 1].set_ylabel('Penetration Depth (mm)', fontsize=11)
axes[1, 1].set_title(f'Penetration Depth vs Nozzle Diameter', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig(f'{output_dir}/exp4_spray_penetration.png', dpi=150, bbox_inches='tight')
plt.close()

print("✓ 实验4完成：喷雾穿透深度")

# ============================================================
# 实验5：液滴燃烧
# ============================================================
print("\n实验5：液滴燃烧")
print("-" * 40)

# 参数设置
T_inf = 1500  # 环境温度 (K)
T_s = 400  # 液滴表面温度 (K，近似沸点)
q_c = 4.5e7  # 燃烧热 (J/kg)
s_stoich = 15.0  # 化学计量比
Y_ox_inf = 0.23  # 环境氧气质量分数

# 传热数
B_T = (cp_g * (T_inf - T_s) + q_c * Y_ox_inf / s_stoich) / h_fg

# 燃烧蒸发常数
K_burn = (8 * lambda_g / (rho_l * cp_g)) * np.log(1 + B_T)

# 初始直径
d0_burn_values = [50, 100, 200]  # 微米

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

colors = plt.cm.hot(np.linspace(0.3, 0.9, len(d0_burn_values)))

# 1. 燃烧过程中液滴直径变化
for i, d0 in enumerate(d0_burn_values):
    tau_burn = (d0 * 1e-6)**2 / K_burn
    t = np.linspace(0, tau_burn, 100)
    
    d_squared = (d0 * 1e-6)**2 - K_burn * t
    d_squared = np.maximum(d_squared, 0)
    
    axes[0, 0].plot(t * 1000, np.sqrt(d_squared) * 1e6, color=colors[i], linewidth=2, 
                    label=f'd0={d0} micron')

axes[0, 0].set_xlabel('Time (ms)', fontsize=11)
axes[0, 0].set_ylabel('Droplet Diameter (micron)', fontsize=11)
axes[0, 0].set_title('Droplet Diameter during Combustion', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# 2. 火焰半径与液滴半径的比值
# 简化模型：火焰半径与传质数相关
B_F = B_T  # 简化假设
B_ox = Y_ox_inf / (1 - Y_ox_inf)  # 氧化剂传质数

rf_rs = np.log(1 + B_ox) / np.log(1 + B_F)

# 不同液滴直径下的火焰半径
d_drop_range = np.linspace(10, 200, 100)  # 微米
r_flame = rf_rs * d_drop_range / 2  # 火焰半径 (微米)

axes[0, 1].plot(d_drop_range, r_flame, 'r-', linewidth=2.5, label='Flame Radius')
axes[0, 1].plot(d_drop_range, d_drop_range / 2, 'b--', linewidth=2, label='Droplet Radius')
axes[0, 1].fill_between(d_drop_range, d_drop_range / 2, r_flame, alpha=0.3, color='orange', label='Flame Zone')
axes[0, 1].set_xlabel('Droplet Diameter (micron)', fontsize=11)
axes[0, 1].set_ylabel('Radius (micron)', fontsize=11)
axes[0, 1].set_title(f'Flame Stand-off Distance (rf/rs={rf_rs:.2f})', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)

# 3. 燃烧时间与初始直径的关系（d²定律验证）
d0_burn_range = np.linspace(10, 300, 100)  # 微米
tau_burn_range = (d0_burn_range * 1e-6)**2 / K_burn * 1000  # ms

axes[1, 0].plot(d0_burn_range, tau_burn_range, 'b-', linewidth=2.5)
axes[1, 0].set_xlabel('Initial Diameter (micron)', fontsize=11)
axes[1, 0].set_ylabel('Burning Time (ms)', fontsize=11)
axes[1, 0].set_title('Burning Time vs Initial Diameter', fontsize=12, fontweight='bold')
axes[1, 0].grid(True, alpha=0.3)

# 4. 不同环境温度下的燃烧
T_inf_burn_values = [1000, 1200, 1500, 1800]  # K
d0_burn_fixed = 100  # 微米

colors_Tburn = plt.cm.coolwarm(np.linspace(0, 1, len(T_inf_burn_values)))

for i, T_inf_i in enumerate(T_inf_burn_values):
    B_T_i = (cp_g * (T_inf_i - T_s) + q_c * Y_ox_inf / s_stoich) / h_fg
    K_burn_i = (8 * lambda_g / (rho_l * cp_g)) * np.log(1 + B_T_i)
    
    tau_burn_i = (d0_burn_fixed * 1e-6)**2 / K_burn_i
    t = np.linspace(0, tau_burn_i, 100)
    
    d_squared_i = (d0_burn_fixed * 1e-6)**2 - K_burn_i * t
    d_squared_i = np.maximum(d_squared_i, 0)
    
    axes[1, 1].plot(t * 1000, np.sqrt(d_squared_i) * 1e6, color=colors_Tburn[i], linewidth=2, 
                    label=f'T_inf={T_inf_i} K')

axes[1, 1].set_xlabel('Time (ms)', fontsize=11)
axes[1, 1].set_ylabel('Droplet Diameter (micron)', fontsize=11)
axes[1, 1].set_title(f'Combustion at Different T_inf (d0={d0_burn_fixed} micron)', fontsize=12, fontweight='bold')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/exp5_droplet_combustion.png', dpi=150, bbox_inches='tight')
plt.close()

print("✓ 实验5完成：液滴燃烧")

# ============================================================
# 总结
# ============================================================
print("\n" + "="*60)
print("所有实验完成！")
print("="*60)
print(f"\n生成的图像文件：")
print(f"1. {output_dir}/exp1_droplet_size_distribution.png")
print(f"2. {output_dir}/exp2_droplet_evaporation.png")
print(f"3. {output_dir}/exp3_droplet_trajectory.png")
print(f"4. {output_dir}/exp4_spray_penetration.png")
print(f"5. {output_dir}/exp5_droplet_combustion.png")
print("\n主题025仿真完成！")