传热学仿真-主题098-火灾蔓延传热模拟

kkchenjj

83人浏览 · 2026-03-21 05:15:50

kkchenjj · 2026-03-21 05:15:50 发布

第九十八篇：火灾蔓延传热模拟

摘要

火灾是威胁人类生命财产安全的主要灾害之一，其蔓延过程涉及复杂的热传递机制，包括热对流、热辐射和可燃物热解等多种物理化学过程。本文系统阐述火灾动力学基础理论、火焰蔓延机理和传热模型，重点介绍森林火灾和建筑火灾的数值模拟方法。通过Python仿真实现火焰蔓延速度计算、热释放速率分析、温度场演化和烟气扩散模拟，为火灾安全工程、消防设计和应急救援提供理论基础和数值工具。

关键词

火灾动力学，火焰蔓延，热释放速率，森林火灾，建筑火灾，CFD模拟，烟气扩散，燃烧传热

在这里插入图片描述

1. 火灾动力学基础理论

1.1 火灾的基本特征

火灾是一种复杂的非定常燃烧现象，具有以下基本特征：

链式反应机制：火灾燃烧遵循自由基链式反应机理，包括链引发、链传递和链终止三个阶段。可燃物在高温下分解产生自由基，自由基与氧气反应释放热量并产生新的自由基，形成自持燃烧。

多尺度特性：火灾现象跨越多个时间和空间尺度：

微观尺度（10⁻⁶~10⁻³ m）：化学反应动力学、分子扩散
介观尺度（10⁻³~10⁻¹ m）：火焰结构、湍流燃烧
宏观尺度（10⁻¹~10³ m）：房间火灾、森林蔓延

多物理场耦合：火灾涉及流体力学、传热学、燃烧学和材料科学的深度耦合：

流场：浮力驱动的自然对流、烟气羽流
温度场：高温火焰区与低温环境的热交换
浓度场：氧气消耗与燃烧产物生成
辐射场：火焰对周围物体的热辐射

1.2 火灾三角与四面体

传统火灾理论强调火灾三角——可燃物、氧气、点火源三要素。现代火灾科学在此基础上增加了链式反应，形成火灾四面体：

$火灾=f(可燃物,氧气,点火源,链式反应)\text{火灾} = f(\text{可燃物}, \text{氧气}, \text{点火源}, \text{链式反应})$

可燃物：包括固体（木材、塑料、纺织品）、液体（汽油、酒精）和气体（天然气、液化石油气）。可燃物的热解特性、燃烧热和火焰传播速度决定火灾发展规律。

氧气：空气中的氧气体积分数约为21%，当氧气浓度低于15%时，大多数可燃物无法持续燃烧。

点火源：提供初始活化能，包括明火、电火花、高温表面等。点火能量必须超过可燃物的最小点火能。

链式反应：燃烧过程中的自由基反应是维持火焰传播的关键。

1.3 火灾发展阶段

典型火灾经历四个发展阶段：

初期增长阶段（Incipient Stage）：

可燃物被点燃，火焰局限于点火源附近
热释放速率（HRR）较低，呈指数增长
持续时间：数秒至数分钟

全面发展阶段（Free-Burning Stage）：

火焰蔓延至整个可燃表面
热释放速率达到最大值
可能出现轰燃（Flashover）现象

衰减阶段（Decay Stage）：

可燃物逐渐消耗
热释放速率下降
温度逐渐降低

熄灭阶段（Extinction）：

可燃物耗尽或氧气不足
火焰熄灭，余热持续释放

2. 火焰蔓延传热机理

2.1 火焰蔓延的基本模式

火焰蔓延是指火焰从已燃区向未燃区传播的过程，主要有三种模式：

顺流蔓延（Wind-Driven Spread）：
火焰在风或气流作用下向下游蔓延。蔓延速度 $v$ 与风速 $u$ 的关系可表示为：

$v_0 + k \cdot u$

其中， $v_0$ 为无风时的蔓延速度， $k$ 为风速影响系数。

逆流蔓延（Against-Wind Spread）：
火焰逆着风向蔓延，主要依靠热传导和辐射传热。蔓延速度较慢，通常满足：

$v_0 - \alpha \cdot u$

当风速超过临界值时，逆流蔓延可能停止。

侧向蔓延（Lateral Spread）：
火焰垂直于主风向蔓延，蔓延速度介于顺流和逆流之间。

2.2 传热机制分析

火焰蔓延过程中的热量传递包括三种基本模式：

热传导：
通过固体可燃物的热传导预热未燃区。一维热传导方程为：

$ρcp∂T∂t=k∂2T∂x2+q˙ext′′\rho c_p \frac{\partial T}{\partial t} = k \frac{\partial^2 T}{\partial x^2} + \dot{q}_{ext}''$

其中， $q˙ext′′\dot{q}_{ext}''$ 为外部热流密度。

热对流：
高温烟气与可燃物表面的对流换热：

$q˙conv′′=h(Tg−Ts)\dot{q}_{conv}'' = h (T_g - T_s)$

对流换热系数 $h$ 在自然对流条件下约为 5-25 W/(m²·K)，在强制对流条件下可达 50-500 W/(m²·K)。

热辐射：
火焰对未燃区的辐射传热是远距离火焰蔓延的主要机制：

$q˙rad′′=εfσTf4F12\dot{q}_{rad}'' = \varepsilon_f \sigma T_f^4 F_{12}$

其中， $εf\varepsilon_f$ 为火焰发射率（通常 0.7-0.95）， $F_{12}$ 为视角因子。

2.3 火焰蔓延速度模型

Rothermel模型（森林火灾）：
Rothermel提出了著名的森林火灾蔓延速度公式：

$\frac{I_R \xi (1 + \phi_w + \phi_s)}{\rho_b \varepsilon Q_{ig}}$

其中：

$I_R$ ：反应强度（kW/m²）
$ξ\xi$ ：传播通量比
$ϕw\phi_w$ ：风速修正因子
$ϕs\phi_s$ ：坡度修正因子
$ρb\rho_b$ ：燃料床堆积密度（kg/m³）
$ε\varepsilon$ ：有效加热系数
$Q_{ig}$ ：点燃单位质量燃料所需热量（kJ/kg）

Thomas模型（建筑火灾）：
对于建筑内部火灾，Thomas提出了基于开口因子的蔓延模型：

$Q˙=1260AwHw\dot{Q} = 1260 A_w \sqrt{H_w}$

其中， $A_w$ 为开口面积（m²）， $H_w$ 为开口高度（m）。

3. 热释放速率与燃烧模型

3.1 热释放速率（HRR）

热释放速率是表征火灾强度的核心参数，定义为：

$Q˙=m˙⋅ΔHc\dot{Q} = \dot{m} \cdot \Delta H_c$

其中， $m˙\dot{m}$ 为质量燃烧速率（kg/s）， $ΔHc\Delta H_c$ 为燃烧热（kJ/kg）。

t²火灾模型：
工程上常用t²模型描述火灾增长：

$Q˙(t)=α(t−t0)2\dot{Q}(t) = \alpha (t - t_0)^2$

根据增长速率，t²火灾分为：

慢速： $α=0.00293\alpha = 0.00293$ kW/s²
中速： $α=0.01172\alpha = 0.01172$ kW/s²
快速： $α=0.0469\alpha = 0.0469$ kW/s²
超快速： $α=0.1876\alpha = 0.1876$ kW/s²

3.2 燃烧模型

单步反应模型：
最简单的燃烧模型假设燃料与氧气一步反应生成产物：

$Fuel+νO2O2→Products\text{Fuel} + \nu_{O_2} \text{O}_2 \rightarrow \text{Products}$

反应速率采用Arrhenius定律：

$ω˙=AρaYfuelbYO2cexp⁡(−EaRT)\dot{\omega} = A \rho^a Y_{fuel}^b Y_{O_2}^c \exp\left(-\frac{E_a}{RT}\right)$

涡耗散模型（EDM）：
对于湍流燃烧，Magnussen提出涡耗散模型：

$ω˙=CEDMρεkmin⁡(Yfuel,YO2s,CprodYprod1+s)\dot{\omega} = C_{EDM} \rho \frac{\varepsilon}{k} \min\left(Y_{fuel}, \frac{Y_{O_2}}{s}, C_{prod}\frac{Y_{prod}}{1+s}\right)$

其中， $k$ 为湍动能， $ε\varepsilon$ 为耗散率， $s$ 为化学计量比。

4. Python仿真实现

4.1 火灾蔓延基础仿真

"""
主题098：火灾蔓延传热模拟
包含：火焰蔓延、热释放速率、温度场、烟气扩散
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Circle, FancyBboxPatch
from matplotlib.collections import LineCollection
import matplotlib
matplotlib.use('Agg')
import warnings
warnings.filterwarnings('ignore')
import os
from scipy.ndimage import gaussian_filter
from matplotlib.colors import LinearSegmentedColormap

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\传热学仿真\主题098'
os.makedirs(output_dir, exist_ok=True)

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

print("=" * 60)
print("火灾蔓延传热模拟系统")
print("=" * 60)

# ============================================
# 1. t²火灾增长模型
# ============================================
print("\n" + "=" * 60)
print("1. t²火灾增长模型分析")
print("=" * 60)

def fire_growth_analysis():
    """
    分析不同增长速率的t²火灾
    """
    # 时间参数
    t = np.linspace(0, 600, 1000)  # 0-600秒
    
    # 不同火灾类型的增长系数
    fire_types = {
        'Slow': 0.00293,
        'Medium': 0.01172,
        'Fast': 0.0469,
        'Ultra-fast': 0.1876
    }
    
    colors = ['green', 'blue', 'orange', 'red']
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：热释放速率-时间曲线
    ax1 = axes[0, 0]
    for (ftype, alpha), color in zip(fire_types.items(), colors):
        Q = alpha * t**2 / 1000  # 转换为MW
        ax1.plot(t, Q, color=color, linewidth=2, label=f'{ftype} (α={alpha})')
    
    ax1.set_xlabel('Time (s)', fontsize=11)
    ax1.set_ylabel('Heat Release Rate (MW)', fontsize=11)
    ax1.set_title('t² Fire Growth Model', fontsize=12, fontweight='bold')
    ax1.legend(loc='upper left')
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(0, 600)
    
    # 图2：热释放速率-时间对数坐标
    ax2 = axes[0, 1]
    for (ftype, alpha), color in zip(fire_types.items(), colors):
        Q = alpha * t**2 / 1000
        ax2.semilogy(t[1:], Q[1:], color=color, linewidth=2, label=ftype)
    
    ax2.set_xlabel('Time (s)', fontsize=11)
    ax2.set_ylabel('Heat Release Rate (MW)', fontsize=11)
    ax2.set_title('HRR vs Time (Log Scale)', fontsize=12, fontweight='bold')
    ax2.legend(loc='lower right')
    ax2.grid(True, which="both", ls="-", alpha=0.3)
    ax2.set_xlim(1, 600)
    
    # 图3：达到特定HRR所需时间
    ax3 = axes[1, 0]
    target_HRR = [1, 5, 10]  # MW
    x_pos = np.arange(len(fire_types))
    width = 0.25
    
    for i, Q_target in enumerate(target_HRR):
        times = []
        for alpha in fire_types.values():
            t_target = np.sqrt(Q_target * 1000 / alpha)
            times.append(t_target)
        ax3.bar(x_pos + i*width, times, width, label=f'{Q_target} MW')
    
    ax3.set_xlabel('Fire Type', fontsize=11)
    ax3.set_ylabel('Time to Reach HRR (s)', fontsize=11)
    ax3.set_title('Time to Reach Target HRR', fontsize=12, fontweight='bold')
    ax3.set_xticks(x_pos + width)
    ax3.set_xticklabels(fire_types.keys())
    ax3.legend()
    ax3.grid(True, alpha=0.3, axis='y')
    
    # 图4：火灾功率曲线（示意）
    ax4 = axes[1, 1]
    
    # 完整的火灾曲线（增长-稳定-衰减）
    t_full = np.linspace(0, 1200, 2000)
    alpha_fast = 0.0469
    t_growth = 300  # 增长阶段结束
    t_decay_start = 600  # 衰减开始
    
    Q_full = np.zeros_like(t_full)
    Q_max = alpha_fast * t_growth**2 / 1000  # MW
    
    for i, ti in enumerate(t_full):
        if ti <= t_growth:
            Q_full[i] = alpha_fast * ti**2 / 1000
        elif ti <= t_decay_start:
            Q_full[i] = Q_max
        else:
            # 线性衰减
            Q_full[i] = Q_max * (1 - (ti - t_decay_start) / 600)
            if Q_full[i] < 0:
                Q_full[i] = 0
    
    ax4.plot(t_full, Q_full, 'r-', linewidth=2.5)
    ax4.fill_between(t_full, Q_full, alpha=0.3, color='red')
    
    # 标注阶段
    ax4.axvspan(0, t_growth, alpha=0.1, color='green', label='Growth')
    ax4.axvspan(t_growth, t_decay_start, alpha=0.1, color='yellow', label='Steady')
    ax4.axvspan(t_decay_start, 1200, alpha=0.1, color='blue', label='Decay')
    
    ax4.set_xlabel('Time (s)', fontsize=11)
    ax4.set_ylabel('Heat Release Rate (MW)', fontsize=11)
    ax4.set_title('Complete Fire Curve', fontsize=12, fontweight='bold')
    ax4.legend(loc='upper right')
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/fire_growth_analysis.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ t²火灾增长模型分析完成")
    print(f"  - 分析了4种火灾增长类型")
    print(f"  - 超快速火灾100秒可达HRR = {0.1876 * 100**2 / 1000:.1f} MW")

fire_growth_analysis()

# ============================================
# 2. 森林火灾蔓延模拟（Rothermel模型）
# ============================================
print("\n" + "=" * 60)
print("2. 森林火灾蔓延模拟")
print("=" * 60)

def forest_fire_spread():
    """
    基于Rothermel模型的森林火灾蔓延模拟
    """
    # 燃料参数（典型针叶林）
    fuel_load = 1.5  # kg/m²，燃料负荷
    fuel_depth = 0.3  # m，燃料床深度
    fuel_moisture = 0.08  # 含水率（干基）
    packing_ratio = 0.004  # 堆积比
    
    # 环境参数
    wind_speeds = np.linspace(0, 10, 50)  # m/s，风速
    slopes = np.linspace(0, 30, 50)  # 坡度（度）
    
    # 简化的Rothermel模型计算
    def rothermel_rate(wind, slope, moisture):
        """简化Rothermel蔓延速度计算（m/min）"""
        # 基础蔓延速度
        R0 = 1.5  # m/min，无风无坡度
        
        # 风速修正
        phi_w = 0.5 * wind**1.5
        
        # 坡度修正
        phi_s = 0.3 * slope**1.2 if slope > 0 else 0
        
        # 含水率修正
        moisture_factor = max(0.1, 1 - moisture * 5)
        
        R = R0 * (1 + phi_w + phi_s) * moisture_factor
        return R
    
    # 计算蔓延速度矩阵
    spread_rates = np.zeros((len(slopes), len(wind_speeds)))
    for i, slope in enumerate(slopes):
        for j, wind in enumerate(wind_speeds):
            spread_rates[i, j] = rothermel_rate(wind, slope, fuel_moisture)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：蔓延速度-风速关系
    ax1 = axes[0, 0]
    for slope in [0, 10, 20, 30]:
        rates = [rothermel_rate(w, slope, fuel_moisture) for w in wind_speeds]
        ax1.plot(wind_speeds, rates, linewidth=2, label=f'Slope = {slope}°')
    
    ax1.set_xlabel('Wind Speed (m/s)', fontsize=11)
    ax1.set_ylabel('Spread Rate (m/min)', fontsize=11)
    ax1.set_title('Fire Spread Rate vs Wind Speed', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 图2：蔓延速度-坡度关系
    ax2 = axes[0, 1]
    for wind in [0, 2, 5, 10]:
        rates = [rothermel_rate(wind, s, fuel_moisture) for s in slopes]
        ax2.plot(slopes, rates, linewidth=2, label=f'Wind = {wind} m/s')
    
    ax2.set_xlabel('Slope (degrees)', fontsize=11)
    ax2.set_ylabel('Spread Rate (m/min)', fontsize=11)
    ax2.set_title('Fire Spread Rate vs Slope', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 图3：蔓延速度等高线图
    ax3 = axes[1, 0]
    X, Y = np.meshgrid(wind_speeds, slopes)
    contour = ax3.contourf(X, Y, spread_rates, levels=20, cmap='hot')
    ax3.contour(X, Y, spread_rates, levels=10, colors='black', linewidths=0.5, alpha=0.5)
    plt.colorbar(contour, ax=ax3, label='Spread Rate (m/min)')
    ax3.set_xlabel('Wind Speed (m/s)', fontsize=11)
    ax3.set_ylabel('Slope (degrees)', fontsize=11)
    ax3.set_title('Fire Spread Rate Contours', fontsize=12, fontweight='bold')
    
    # 图4：火灾蔓延可视化
    ax4 = axes[1, 1]
    
    # 模拟火灾蔓延
    n_steps = 50
    dt = 1  # 分钟
    
    # 初始火点
    fire_x, fire_y = [0], [0]
    
    # 风向（从左向右）
    wind_direction = 0  # 度，0表示向右
    wind_speed = 5  # m/s
    
    for step in range(n_steps):
        # 计算当前蔓延速度
        R = rothermel_rate(wind_speed, 0, fuel_moisture)
        
        # 顺流蔓延距离
        dx_forward = R * dt * np.cos(np.radians(wind_direction))
        dy_forward = R * dt * np.sin(np.radians(wind_direction))
        
        # 侧向蔓延距离（较慢）
        R_lateral = R * 0.3
        
        # 更新火场边界（简化为椭圆）
        time = (step + 1) * dt
        a = R * time  # 长半轴（顺流）
        b = R_lateral * time  # 短半轴（侧向）
        
        theta = np.linspace(0, 2*np.pi, 100)
        x_ellipse = a * np.cos(theta)
        y_ellipse = b * np.sin(theta)
        
        if step % 10 == 0:
            color_intensity = step / n_steps
            ax4.plot(x_ellipse, y_ellipse, color=(1, 1-color_intensity*0.5, 0), 
                    linewidth=2, alpha=0.7, label=f't={time}min' if step == 0 or step == n_steps-1 else '')
            ax4.fill(x_ellipse, y_ellipse, alpha=0.1, color='red')
    
    # 绘制风向箭头
    ax4.arrow(-20, 15, 10, 0, head_width=2, head_length=2, fc='blue', ec='blue')
    ax4.text(-15, 18, 'Wind', fontsize=10, color='blue')
    
    ax4.set_xlabel('Distance (m)', fontsize=11)
    ax4.set_ylabel('Distance (m)', fontsize=11)
    ax4.set_title('Fire Spread Pattern', fontsize=12, fontweight='bold')
    ax4.set_aspect('equal')
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/forest_fire_spread.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 森林火灾蔓延模拟完成")
    print(f"  - 最大蔓延速度: {np.max(spread_rates):.1f} m/min")
    print(f"  - 风速影响显著，坡度加剧蔓延")

forest_fire_spread()

# ============================================
# 3. 建筑火灾温度场模拟
# ============================================
print("\n" + "=" * 60)
print("3. 建筑火灾温度场模拟")
print("=" * 60)

def building_fire_simulation():
    """
    建筑房间火灾温度场模拟
    """
    # 房间尺寸
    Lx, Ly = 10, 8  # m
    Nx, Ny = 100, 80
    dx, dy = Lx/Nx, Ly/Ny
    
    # 火源参数
    fire_x, fire_y = Nx//2, Ny//4  # 火源位置
    Q_fire = 500  # kW，热释放速率
    
    # 时间参数
    dt = 0.1  # s
    n_steps = 500
    
    # 物性参数
    rho = 1.2  # kg/m³
    cp = 1005  # J/(kg·K)
    k = 0.025  # W/(m·K)
    alpha = k / (rho * cp)  # 热扩散系数
    
    # 稳定性条件
    Fo = alpha * dt / dx**2
    print(f"  - 傅里叶数 Fo = {Fo:.4f}")
    
    # 初始化温度场
    T = np.ones((Nx, Ny)) * 20  # 初始温度20°C
    T_history = []
    
    # 模拟
    for step in range(n_steps):
        T_new = T.copy()
        
        # 内部节点（有限差分）
        for i in range(1, Nx-1):
            for j in range(1, Ny-1):
                # 热源区域
                if (i - fire_x)**2 + (j - fire_y)**2 < 25:
                    q_source = Q_fire * 1000 / (25 * dx * dy)  # W/m³
                else:
                    q_source = 0
                
                # 热传导
                T_new[i, j] = T[i, j] + Fo * (
                    T[i+1, j] + T[i-1, j] + T[i, j+1] + T[i, j-1] - 4*T[i, j]
                ) + q_source * dt / (rho * cp)
        
        # 边界条件（对流换热）
        h_conv = 10  # W/(m²·K)
        T_ambient = 20
        
        # 边界
        T_new[0, :] = T_new[1, :] - h_conv * dx / k * (T_new[1, :] - T_ambient)
        T_new[-1, :] = T_new[-2, :] - h_conv * dx / k * (T_new[-2, :] - T_ambient)
        T_new[:, 0] = T_new[:, 1] - h_conv * dy / k * (T_new[:, 1] - T_ambient)
        T_new[:, -1] = T_new[:, -2] - h_conv * dy / k * (T_new[:, -2] - T_ambient)
        
        T = T_new.copy()
        
        if step % 100 == 0:
            T_history.append(T.copy())
    
    # 绘制结果
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：最终温度场
    ax1 = axes[0, 0]
    im1 = ax1.contourf(T.T, levels=20, cmap='hot', origin='lower')
    plt.colorbar(im1, ax=ax1, label='Temperature (°C)')
    ax1.set_xlabel('X (grid)', fontsize=11)
    ax1.set_ylabel('Y (grid)', fontsize=11)
    ax1.set_title('Final Temperature Field', fontsize=12, fontweight='bold')
    
    # 图2：温度随时间演变（中心点）
    ax2 = axes[0, 1]
    times = np.arange(0, n_steps * dt, dt * 100)
    center_temps = [T_hist[fire_x, fire_y] for T_hist in T_history]
    
    ax2.plot(times, center_temps, 'r-', linewidth=2)
    ax2.set_xlabel('Time (s)', fontsize=11)
    ax2.set_ylabel('Temperature (°C)', fontsize=11)
    ax2.set_title('Temperature at Fire Center', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 图3：温度剖面
    ax3 = axes[1, 0]
    x_profile = np.linspace(0, Lx, Nx)
    y_profile = np.linspace(0, Ly, Ny)
    
    ax3.plot(x_profile, T[:, fire_y], 'r-', linewidth=2, label='Horizontal')
    ax3.plot(y_profile, T[fire_x, :], 'b-', linewidth=2, label='Vertical')
    ax3.set_xlabel('Distance (m)', fontsize=11)
    ax3.set_ylabel('Temperature (°C)', fontsize=11)
    ax3.set_title('Temperature Profiles', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 图4：温度场动画帧
    ax4 = axes[1, 1]
    im4 = ax4.imshow(T_history[2].T, cmap='hot', origin='lower', vmin=20, vmax=300)
    plt.colorbar(im4, ax=ax4, label='Temperature (°C)')
    ax4.set_xlabel('X (grid)', fontsize=11)
    ax4.set_ylabel('Y (grid)', fontsize=11)
    ax4.set_title(f'Temperature at t={times[2]:.0f}s', fontsize=12, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/building_fire_temperature.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 建筑火灾温度场模拟完成")
    print(f"  - 火源中心最高温度: {np.max(T):.1f}°C")
    print(f"  - 模拟时间: {n_steps * dt:.0f} s")

building_fire_simulation()

# ============================================
# 4. 烟气扩散模拟
# ============================================
print("\n" + "=" * 60)
print("4. 烟气扩散模拟")
print("=" * 60)

def smoke_dispersion():
    """
    火灾烟气羽流和扩散模拟
    """
    # 空间网格
    Lx, Ly = 20, 15  # m
    Nx, Ny = 200, 150
    dx, dy = Lx/Nx, Ly/Ny
    
    # 火源参数
    fire_x, fire_y = Nx//2, Ny//4
    Q = 1000  # kW，热释放速率
    
    # 时间参数
    dt = 0.05
    n_steps = 300
    
    # 烟气浓度场
    C = np.zeros((Nx, Ny))
    
    # 速度场（简化的浮力驱动流）
    U = np.zeros((Nx, Ny))
    V = np.zeros((Nx, Ny))
    
    # 生成速度场（羽流上升）
    for i in range(Nx):
        for j in range(Ny):
            r = np.sqrt((i - fire_x)**2 + (j - fire_y)**2) * dx
            z = j * dy
            
            if r < 3 and z > 1:
                # 羽流中心上升速度
                V[i, j] = 2.0 * np.exp(-r**2 / (0.5 + z*0.1))
            
            # 卷吸速度
            if r > 0.5 and r < 4:
                U[i, j] = -0.5 * (i - fire_x) * dx / r * np.exp(-r/2)
    
    # 扩散系数
    D = 0.1  # m²/s
    
    # 模拟
    C_history = []
    
    for step in range(n_steps):
        C_new = C.copy()
        
        # 对流-扩散方程
        for i in range(1, Nx-1):
            for j in range(1, Ny-1):
                # 对流项（上风差分）
                if U[i, j] > 0:
                    dCdx = (C[i, j] - C[i-1, j]) / dx
                else:
                    dCdx = (C[i+1, j] - C[i, j]) / dx
                
                if V[i, j] > 0:
                    dCdy = (C[i, j] - C[i, j-1]) / dy
                else:
                    dCdy = (C[i, j+1] - C[i, j]) / dy
                
                # 扩散项
                d2Cdx2 = (C[i+1, j] - 2*C[i, j] + C[i-1, j]) / dx**2
                d2Cdy2 = (C[i, j+1] - 2*C[i, j] + C[i, j-1]) / dy**2
                
                C_new[i, j] = C[i, j] - dt * (U[i, j] * dCdx + V[i, j] * dCdy) + \
                             D * dt * (d2Cdx2 + d2Cdy2)
        
        # 火源持续释放烟气
        if step < 200:
            C_new[fire_x-2:fire_x+3, fire_y:fire_y+3] += 0.5
        
        C = np.clip(C_new, 0, 10)
        
        if step % 60 == 0:
            C_history.append(C.copy())
    
    # 绘制结果
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：最终烟气浓度
    ax1 = axes[0, 0]
    im1 = ax1.contourf(C.T, levels=15, cmap='gray_r', origin='lower', vmin=0, vmax=5)
    plt.colorbar(im1, ax=ax1, label='Smoke Concentration')
    ax1.set_xlabel('X (grid)', fontsize=11)
    ax1.set_ylabel('Y (grid)', fontsize=11)
    ax1.set_title('Smoke Concentration Field', fontsize=12, fontweight='bold')
    
    # 图2：速度场
    ax2 = axes[0, 1]
    skip = 10
    ax2.quiver(np.arange(0, Nx, skip), np.arange(0, Ny, skip), 
              U[::skip, ::skip].T, V[::skip, ::skip].T, scale=20)
    ax2.set_xlabel('X (grid)', fontsize=11)
    ax2.set_ylabel('Y (grid)', fontsize=11)
    ax2.set_title('Velocity Field', fontsize=12, fontweight='bold')
    ax2.set_aspect('equal')
    
    # 图3：不同时刻的烟气扩散
    ax3 = axes[1, 0]
    times = [0, 1, 2, 3, 4]
    colors = plt.cm.Reds(np.linspace(0.3, 1, len(times)))
    
    for i, (t_idx, color) in enumerate(zip(times, colors)):
        if t_idx < len(C_history):
            # 取中心剖面
            ax3.plot(C_history[t_idx][:, fire_y + 10], color=color, 
                    linewidth=2, label=f't={t_idx*3}s')
    
    ax3.set_xlabel('X (grid)', fontsize=11)
    ax3.set_ylabel('Concentration', fontsize=11)
    ax3.set_title('Smoke Concentration Profiles', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 图4：烟气层高度估计
    ax4 = axes[1, 1]
    
    # 简化的烟气层高度计算
    t_height = np.linspace(0, 15, 100)
    H_room = 3.0  # 房间高度
    
    # 经验公式：烟气层下降
    z_clear = H_room * np.exp(-0.1 * t_height)
    z_clear = np.maximum(z_clear, 0.5)  # 最低保持0.5m
    
    ax4.fill_between(t_height, 0, H_room - z_clear, alpha=0.3, color='gray', label='Smoke Layer')
    ax4.fill_between(t_height, H_room - z_clear, H_room, alpha=0.3, color='cyan', label='Clear Layer')
    ax4.plot(t_height, H_room - z_clear, 'k-', linewidth=2, label='Interface')
    
    # 安全高度线
    ax4.axhline(y=1.8, color='r', linestyle='--', linewidth=2, label='Safe Height')
    
    ax4.set_xlabel('Time (min)', fontsize=11)
    ax4.set_ylabel('Height (m)', fontsize=11)
    ax4.set_title('Smoke Layer Development', fontsize=12, fontweight='bold')
    ax4.legend(loc='upper right')
    ax4.set_ylim(0, H_room)
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/smoke_dispersion.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 烟气扩散模拟完成")
    print(f"  - 烟气羽流上升速度: {np.max(V):.1f} m/s")
    print(f"  - 模拟了烟气层下降过程")

smoke_dispersion()

# ============================================
# 5. 辐射传热与火焰高度
# ============================================
print("\n" + "=" * 60)
print("5. 辐射传热与火焰高度")
print("=" * 60)

def radiation_heat_transfer():
    """
    火焰辐射传热分析
    """
    # 火焰参数
    Q_range = np.linspace(100, 5000, 100)  # kW，热释放速率
    
    # Heskestad火焰高度公式
    def flame_height(Q):
        """计算火焰高度（m）"""
        Q_kW = Q  # kW
        H_f = 0.235 * Q_kW**(2/5) - 1.02 * 1.0  # 假设直径1m
        return max(H_f, 0.5)
    
    # 计算火焰高度
    heights = [flame_height(Q) for Q in Q_range]
    
    # 辐射热流计算
    def radiation_heat_flux(Q, H_f, r):
        """
        计算距离火焰r处的辐射热流
        Q: kW
        H_f: m
        r: m
        """
        # 火焰温度
        T_f = 1000  # K
        
        # 火焰表面积（简化为圆柱体）
        D_f = 1.0  # m，火焰直径
        A_f = np.pi * D_f * H_f + np.pi * D_f**2 / 4
        
        # 视角因子（简化）
        F = min(1.0, (D_f * H_f) / (4 * r**2))
        
        # 发射率
        epsilon = 0.8
        
        # 辐射热流
        sigma = 5.67e-8  # W/(m²·K⁴)
        q_rad = epsilon * sigma * T_f**4 * F / 1000  # kW/m²
        
        return q_rad
    
    # 不同距离处的辐射热流
    distances = [1, 2, 5, 10]  # m
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：火焰高度-热释放速率
    ax1 = axes[0, 0]
    ax1.plot(Q_range, heights, 'r-', linewidth=2)
    ax1.set_xlabel('Heat Release Rate (kW)', fontsize=11)
    ax1.set_ylabel('Flame Height (m)', fontsize=11)
    ax1.set_title('Flame Height vs HRR', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 添加参考点
    ref_points = [1000, 2500, 5000]
    for Q in ref_points:
        H = flame_height(Q)
        ax1.scatter([Q], [H], s=100, color='blue', zorder=5)
        ax1.annotate(f'H={H:.1f}m', xy=(Q, H), xytext=(10, 10), 
                    textcoords='offset points', fontsize=9)
    
    # 图2：辐射热流-距离关系
    ax2 = axes[0, 1]
    r_range = np.linspace(1, 20, 100)
    Q_test = 2000  # kW
    H_test = flame_height(Q_test)
    
    q_rad_values = [radiation_heat_flux(Q_test, H_test, r) for r in r_range]
    
    ax2.semilogy(r_range, q_rad_values, 'b-', linewidth=2)
    
    # 标注临界热流
    ax2.axhline(y=12.5, color='r', linestyle='--', linewidth=2, label='Wood Ignition (12.5 kW/m²)')
    ax2.axhline(y=4.0, color='orange', linestyle='--', linewidth=2, label='Pain Threshold (4 kW/m²)')
    ax2.axhline(y=1.6, color='green', linestyle='--', linewidth=2, label='Safe (1.6 kW/m²)')
    
    ax2.set_xlabel('Distance from Fire (m)', fontsize=11)
    ax2.set_ylabel('Radiative Heat Flux (kW/m²)', fontsize=11)
    ax2.set_title('Radiation Heat Flux vs Distance', fontsize=12, fontweight='bold')
    ax2.legend(loc='upper right')
    ax2.grid(True, alpha=0.3, which='both')
    
    # 图3：不同HRR下的辐射热流
    ax3 = axes[1, 0]
    for Q in [500, 1000, 2000, 5000]:
        H = flame_height(Q)
        q_values = [radiation_heat_flux(Q, H, r) for r in r_range]
        ax3.semilogy(r_range, q_values, linewidth=2, label=f'Q={Q} kW')
    
    ax3.set_xlabel('Distance from Fire (m)', fontsize=11)
    ax3.set_ylabel('Radiative Heat Flux (kW/m²)', fontsize=11)
    ax3.set_title('Radiation from Fires of Different Sizes', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3, which='both')
    
    # 图4：安全距离计算
    ax4 = axes[1, 1]
    
    # 计算达到临界热流的距离
    q_critical = 12.5  # kW/m²，木材点燃
    safe_distances = []
    
    for Q in Q_range:
        H = flame_height(Q)
        # 数值求解达到临界热流的距离
        for r in np.linspace(1, 50, 500):
            q = radiation_heat_flux(Q, H, r)
            if q < q_critical:
                safe_distances.append(r)
                break
        else:
            safe_distances.append(50)
    
    ax4.plot(Q_range, safe_distances, 'purple', linewidth=2.5)
    ax4.fill_between(Q_range, 0, safe_distances, alpha=0.2, color='red', label='Ignition Risk Zone')
    
    ax4.set_xlabel('Heat Release Rate (kW)', fontsize=11)
    ax4.set_ylabel('Safe Distance (m)', fontsize=11)
    ax4.set_title('Safe Distance from Fire', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    
    # 添加参考
    ax4.scatter([1000, 2500, 5000], 
               [safe_distances[20], safe_distances[49], safe_distances[-1]], 
               s=100, color='red', zorder=5)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/radiation_heat_transfer.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 辐射传热分析完成")
    print(f"  - 2000kW火灾火焰高度: {flame_height(2000):.1f} m")
    print(f"  - 木材点燃安全距离(2000kW): ~{safe_distances[40]:.1f} m")

radiation_heat_transfer()

# ============================================
# 6. 综合火灾场景模拟
# ============================================
print("\n" + "=" * 60)
print("6. 综合火灾场景模拟")
print("=" * 60)

def comprehensive_fire_scenario():
    """
    综合火灾场景：多室建筑火灾模拟
    """
    fig = plt.figure(figsize=(16, 12))
    gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)
    
    # 场景参数
    room_dims = (10, 8, 3)  # m
    fire_hrr = 2500  # kW
    
    # 图1：建筑平面图
    ax1 = fig.add_subplot(gs[0, 0])
    
    # 绘制房间布局
    rooms = [
        Rectangle((0, 0), 5, 4, fill=True, facecolor='lightblue', edgecolor='black', linewidth=2),
        Rectangle((5, 0), 5, 4, fill=True, facecolor='lightgreen', edgecolor='black', linewidth=2),
        Rectangle((0, 4), 5, 4, fill=True, facecolor='lightyellow', edgecolor='black', linewidth=2),
        Rectangle((5, 4), 5, 4, fill=True, facecolor='lightcoral', edgecolor='black', linewidth=2),
    ]
    
    for room in rooms:
        ax1.add_patch(room)
    
    # 火源位置
    ax1.scatter([2.5], [2], s=300, color='red', marker='*', zorder=5, label='Fire')
    
    # 门
    ax1.plot([5, 5], [1.5, 2.5], 'brown', linewidth=4)
    ax1.plot([2, 3], [4, 4], 'brown', linewidth=4)
    
    ax1.set_xlim(-0.5, 10.5)
    ax1.set_ylim(-0.5, 8.5)
    ax1.set_aspect('equal')
    ax1.set_title('Building Floor Plan', fontsize=11, fontweight='bold')
    ax1.legend()
    ax1.set_xlabel('X (m)')
    ax1.set_ylabel('Y (m)')
    
    # 图2：HRR曲线
    ax2 = fig.add_subplot(gs[0, 1])
    t = np.linspace(0, 600, 1000)
    alpha = 0.0469  # 快速火灾
    Q_curve = alpha * t**2
    Q_curve = np.minimum(Q_curve, fire_hrr)  # 限制最大HRR
    
    ax2.plot(t, Q_curve/1000, 'r-', linewidth=2.5)
    ax2.fill_between(t, Q_curve/1000, alpha=0.3, color='red')
    ax2.axhline(y=fire_hrr/1000, color='black', linestyle='--', linewidth=1.5, label='Max HRR')
    
    ax2.set_xlabel('Time (s)', fontsize=10)
    ax2.set_ylabel('HRR (MW)', fontsize=10)
    ax2.set_title('Fire HRR Curve', fontsize=11, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 图3：温度时间曲线
    ax3 = fig.add_subplot(gs[0, 2])
    
    # 不同位置的温度
    positions = ['Fire Room', 'Adjacent', 'Corridor', 'Exit']
    delays = [0, 30, 60, 120]  # 延迟时间
    max_temps = [800, 400, 200, 100]
    
    for pos, delay, max_temp in zip(positions, delays, max_temps):
        temp = np.zeros_like(t)
        for i, ti in enumerate(t):
            if ti > delay:
                temp[i] = min(max_temp, 20 + (max_temp - 20) * (1 - np.exp(-(ti - delay) / 100)))
        ax3.plot(t, temp, linewidth=2, label=pos)
    
    ax3.axhline(y=60, color='orange', linestyle='--', linewidth=1.5, label='Tenability Limit')
    ax3.set_xlabel('Time (s)', fontsize=10)
    ax3.set_ylabel('Temperature (°C)', fontsize=10)
    ax3.set_title('Temperature at Different Locations', fontsize=11, fontweight='bold')
    ax3.legend(fontsize=8)
    ax3.grid(True, alpha=0.3)
    
    # 图4：烟气层高度
    ax4 = fig.add_subplot(gs[1, 0])
    
    H_room = room_dims[2]
    z_interface = H_room * np.exp(-0.003 * t)
    z_interface = np.maximum(z_interface, 0.5)
    
    ax4.fill_between(t, 0, H_room - z_interface, alpha=0.4, color='gray', label='Smoke')
    ax4.fill_between(t, H_room - z_interface, H_room, alpha=0.4, color='cyan', label='Clear')
    ax4.plot(t, H_room - z_interface, 'k-', linewidth=2, label='Interface')
    ax4.axhline(y=1.8, color='r', linestyle='--', linewidth=2, label='Safe Height')
    
    ax4.set_xlabel('Time (s)', fontsize=10)
    ax4.set_ylabel('Height (m)', fontsize=10)
    ax4.set_title('Smoke Layer Height', fontsize=11, fontweight='bold')
    ax4.legend(fontsize=8)
    ax4.set_ylim(0, H_room)
    ax4.grid(True, alpha=0.3)
    
    # 图5：CO浓度
    ax5 = fig.add_subplot(gs[1, 1])
    
    CO_ppm = 1000 * (1 - np.exp(-t / 200))
    CO_ppm = np.minimum(CO_ppm, 5000)
    
    ax5.plot(t, CO_ppm, 'purple', linewidth=2)
    ax5.fill_between(t, CO_ppm, alpha=0.3, color='purple')
    
    ax5.axhline(y=1200, color='red', linestyle='--', linewidth=2, label='Dangerous')
    ax5.axhline(y=250, color='orange', linestyle='--', linewidth=1.5, label='Warning')
    
    ax5.set_xlabel('Time (s)', fontsize=10)
    ax5.set_ylabel('CO Concentration (ppm)', fontsize=10)
    ax5.set_title('CO Concentration', fontsize=11, fontweight='bold')
    ax5.legend(fontsize=8)
    ax5.grid(True, alpha=0.3)
    
    # 图6：能见度
    ax6 = fig.add_subplot(gs[1, 2])
    
    visibility = 30 * np.exp(-t / 150)
    visibility = np.maximum(visibility, 0.5)
    
    ax6.plot(t, visibility, 'brown', linewidth=2)
    ax6.fill_between(t, 0, visibility, alpha=0.3, color='brown')
    
    ax6.axhline(y=10, color='red', linestyle='--', linewidth=2, label='Walking Limit')
    ax6.axhline(y=5, color='orange', linestyle='--', linewidth=1.5, label='Critical')
    
    ax6.set_xlabel('Time (s)', fontsize=10)
    ax6.set_ylabel('Visibility (m)', fontsize=10)
    ax6.set_title('Visibility in Smoke', fontsize=11, fontweight='bold')
    ax6.legend(fontsize=8)
    ax6.grid(True, alpha=0.3)
    
    # 图7：ASET/RSET分析
    ax7 = fig.add_subplot(gs[2, 0])
    
    # 可用安全疏散时间（ASET）
    ASET = 180  # s
    
    # 所需安全疏散时间（RSET）
    RSET_detection = 30
    RSET_alarm = 30
    RSET_pre_movement = 60
    RSET_travel = 120
    RSET_total = RSET_detection + RSET_alarm + RSET_pre_movement + RSET_travel
    
    categories = ['Detection', 'Alarm', 'Pre-movement', 'Travel', 'Total RSET']
    times_rset = [RSET_detection, RSET_alarm, RSET_pre_movement, RSET_travel, RSET_total]
    colors_bar = ['blue', 'green', 'orange', 'red', 'purple']
    
    bars = ax7.bar(categories, times_rset, color=colors_bar, alpha=0.7)
    ax7.axhline(y=ASET, color='red', linestyle='--', linewidth=2, label=f'ASET = {ASET}s')
    
    ax7.set_ylabel('Time (s)', fontsize=10)
    ax7.set_title('ASET/RSET Analysis', fontsize=11, fontweight='bold')
    ax7.legend()
    ax7.set_xticklabels(categories, rotation=45, ha='right')
    ax7.grid(True, alpha=0.3, axis='y')
    
    # 图8：人员疏散模拟
    ax8 = fig.add_subplot(gs[2, 1])
    
    # 疏散曲线
    n_people = 100
    t_evac = np.linspace(0, 300, 100)
    evacuated = n_people * (1 - np.exp(-t_evac / 60))
    
    ax8.plot(t_evac, evacuated, 'b-', linewidth=2.5, label='Evacuated')
    ax8.fill_between(t_evac, evacuated, alpha=0.3, color='blue')
    ax8.axhline(y=n_people, color='green', linestyle='--', linewidth=2, label='Total')
    
    ax8.set_xlabel('Time (s)', fontsize=10)
    ax8.set_ylabel('Number of People', fontsize=10)
    ax8.set_title('Evacuation Progress', fontsize=11, fontweight='bold')
    ax8.legend()
    ax8.grid(True, alpha=0.3)
    
    # 图9：风险评估矩阵
    ax9 = fig.add_subplot(gs[2, 2])
    
    risk_matrix = np.array([
        [1, 2, 3],
        [2, 3, 4],
        [3, 4, 5]
    ])
    
    im = ax9.imshow(risk_matrix, cmap='RdYlGn_r', aspect='auto', vmin=1, vmax=5)
    
    labels = [['Low', 'Low-Med', 'Medium'],
              ['Low-Med', 'Medium', 'Med-High'],
              ['Medium', 'Med-High', 'High']]
    
    for i in range(3):
        for j in range(3):
            ax9.text(j, i, labels[i][j], ha='center', va='center', 
                    fontsize=9, fontweight='bold')
    
    # 标记当前风险
    ax9.scatter([1], [2], s=300, color='blue', marker='*', 
               edgecolors='white', linewidths=2, zorder=5)
    
    ax9.set_xticks([0, 1, 2])
    ax9.set_yticks([0, 1, 2])
    ax9.set_xticklabels(['Low', 'Medium', 'High'])
    ax9.set_yticklabels(['Low', 'Medium', 'High'])
    ax9.set_xlabel('Fire Load', fontsize=10)
    ax9.set_ylabel('Occupant Density', fontsize=10)
    ax9.set_title('Fire Risk Assessment', fontsize=11, fontweight='bold')
    
    plt.savefig(f'{output_dir}/comprehensive_fire_scenario.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 综合火灾场景模拟完成")
    print(f"  - 模拟了多室建筑火灾")
    print(f"  - ASET = {ASET}s, RSET = {RSET_total}s")
    print(f"  - 安全裕度: {ASET - RSET_total} s")

comprehensive_fire_scenario()

print("\n" + "=" * 60)
print("火灾蔓延传热模拟全部完成")
print("=" * 60)
print(f"\n生成的图表文件：")
print(f"  1. fire_growth_analysis.png - t²火灾增长模型")
print(f"  2. forest_fire_spread.png - 森林火灾蔓延")
print(f"  3. building_fire_temperature.png - 建筑火灾温度场")
print(f"  4. smoke_dispersion.png - 烟气扩散模拟")
print(f"  5. radiation_heat_transfer.png - 辐射传热分析")
print(f"  6. comprehensive_fire_scenario.png - 综合火灾场景")
print(f"\n输出目录: {output_dir}")

### 4.2 仿真2：烟气层高度演化分析

```python
"""
仿真2：烟气层高度演化分析
模拟火灾中烟气层的下降过程
"""

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

# 房间参数
H_room = 3.0  # m，房间高度
A_room = 100  # m²，房间面积
V_room = H_room * A_room  # 房间体积

# 火灾参数
Q_fire = 2000  # kW，热释放速率
m_dot = 0.05 * Q_fire  # kg/s，质量流率（简化模型）

# 烟气层参数
rho_smoke = 0.8  # kg/m³，烟气密度
cp_smoke = 1005  # J/(kg·K)
T_smoke = 200 + 273.15  # K，烟气温度
T_ambient = 20 + 273.15  # K，环境温度

# 时间参数
t_max = 600  # s
t = np.linspace(0, t_max, 1000)

# 烟气层高度模型（简化）
def smoke_layer_height(t, H, A, Q, rho, cp, T_s, T_a):
    """计算烟气层高度随时间变化"""
    # 填充时间常数
    tau = rho * cp * H * A * (T_s - T_a) / (Q * 1000)
    
    # 烟气层高度（从天花板向下）
    h_smoke = H * (1 - np.exp(-t / tau))
    
    return h_smoke, tau

h_smoke, tau_fill = smoke_layer_height(t, H_room, A_room, Q_fire, 
                                       rho_smoke, cp_smoke, T_smoke, T_ambient)

# 安全高度
h_safe = 1.8  # m，人员安全高度

axes[0, 0].plot(t, h_smoke, 'r-', linewidth=2, label='烟气层高度')
axes[0, 0].axhline(y=h_safe, color='g', linestyle='--', label='安全高度')
axes[0, 0].fill_between(t, 0, h_smoke, alpha=0.3, color='gray', label='烟气区')
axes[0, 0].set_xlabel('时间 (s)', fontsize=11)
axes[0, 0].set_ylabel('烟气层高度 (m)', fontsize=11)
axes[0, 0].set_title('烟气层高度演化', fontsize=12)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 不同热释放速率的影响
Q_range_fire = [500, 1000, 2000, 5000]  # kW
colors_q = ['blue', 'green', 'red', 'purple']

for Q, color in zip(Q_range_fire, colors_q):
    h_s, _ = smoke_layer_height(t, H_room, A_room, Q, 
                                rho_smoke, cp_smoke, T_smoke, T_ambient)
    axes[0, 1].plot(t, h_s, color=color, linewidth=2, label=f'Q={Q}kW')

axes[0, 1].axhline(y=h_safe, color='k', linestyle='--', alpha=0.5)
axes[0, 1].set_xlabel('时间 (s)', fontsize=11)
axes[0, 1].set_ylabel('烟气层高度 (m)', fontsize=11)
axes[0, 1].set_title('热释放速率对烟气层的影响', fontsize=12)
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 可用安全疏散时间（ASET）
# 当烟气层降到安全高度以下时
def calculate_ASET(H, h_safe, tau):
    """计算ASET"""
    if h_safe >= H:
        return np.inf
    t_aseet = -tau * np.log(1 - h_safe / H)
    return t_aseet

ASET_values = []
for Q in Q_range_fire:
    _, tau_q = smoke_layer_height(0, H_room, A_room, Q, 
                                  rho_smoke, cp_smoke, T_smoke, T_ambient)
    aset = calculate_ASET(H_room, h_safe, tau_q)
    ASET_values.append(aset)

axes[1, 0].bar(range(len(Q_range_fire)), ASET_values, 
               color=colors_q, alpha=0.7)
axes[1, 0].set_xticks(range(len(Q_range_fire)))
axes[1, 0].set_xticklabels([f'{q}kW' for q in Q_range_fire])
axes[1, 0].set_xlabel('热释放速率', fontsize=11)
axes[1, 0].set_ylabel('ASET (s)', fontsize=11)
axes[1, 0].set_title('可用安全疏散时间', fontsize=12)
axes[1, 0].grid(True, alpha=0.3, axis='y')

# 烟气温度演化
# 假设烟气温度随时间上升
T_smoke_evolution = T_ambient + (T_smoke - T_ambient) * (1 - np.exp(-t / 200))

axes[1, 1].plot(t, T_smoke_evolution - 273.15, 'r-', linewidth=2)
axes[1, 1].axhline(y=60, color='orange', linestyle='--', label='耐受极限')
axes[1, 1].axhline(y=120, color='red', linestyle='--', label='危险温度')
axes[1, 1].set_xlabel('时间 (s)', fontsize=11)
axes[1, 1].set_ylabel('烟气温度 (°C)', fontsize=11)
axes[1, 1].set_title('烟气温度演化', fontsize=12)
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/smoke_layer_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("烟气层分析图已保存")

print("\n烟气层分析结果:")
print(f"房间高度: {H_room}m, 面积: {A_room}m²")
print(f"填充时间常数: {tau_fill:.1f}s")
print(f"ASET (Q=2000kW): {ASET_values[2]:.1f}s")

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

【继电保护】小电流接地系统故障仿真-中性点不接地与经消弧线圈接地仿真模型（Simulink仿真实现）

小电流接地系统故障仿真是电力系统中非常重要的研究领域，特别是针对中性点不接地和经消弧线圈接地的情况。这两种故障情况在电力系统中都可能发生，因此对其进行仿真模型研究具有重要意义。中性点不接地故障是指变压器或发电机中性点没有接地，这种情况下，如果出现了单相接地故障，会导致系统中产生零序电流，可能对设备和系统造成严重损坏。因此，针对中性点不接地故障，需要建立相应的仿真模型，研究其对电力系统的影响，以及采

AtomGit开源社区

BEVFormer环境搭建

AtomGit开源社区

微电网两阶段鲁棒优化经济调度方法[3]【升级优化版本】（Matlab代码实现）

针对微电网内可再生能源和负荷的不确定性，建立了min-max-min 结构的两阶段鲁棒优化模型，可得到最恶劣场景下运行成本最低的调度方案。模型中考虑了储能、需求侧负荷及可控分布式电源等的运行约束和协调控制,并引入了不确定性调节参数，可灵活调整调度方案的保守性。基于列约束生成算法和强对偶理论,可将原问题分解为具有混合整数线性特征的主问题和子问题进行交替求解,从而得到原问题的最优解。