蒸发暴露数学模型复现失败求助：Python实现与预期不符

模型来源 #

基于RIVM报告中4.1.1.4章节《Exposure to Vapour: Evaporation》的数学模型复现。

核心公式 #

文档中的核心控制方程分为两个阶段：

1. 释放阶段：物质从表面蒸发进入空气，同时通风排出
   $$\frac{dA_{air}}{dt} = K \cdot A_{release} \cdot \frac{M}{R \cdot T} \cdot (P_{eq} - P_{air}) - \frac{Q}{3600} \cdot V \cdot \frac{A_{air}}{V}$$
   $$\frac{dA_{prod}}{dt} = -\frac{dA_{air}}{dt}$$
2. 释放结束后：仅通风排出空气中的物质
   $$\frac{dA_{air}}{dt} = -\frac{Q}{3600} \cdot V \cdot \frac{A_{air}}{V}$$
   $$\frac{dA_{prod}}{dt} = 0$$

其中参数定义：

• $A_{air}$：空气中物质质量（kg）
• $A_{prod}$：表面残留物质质量（kg）
• $K$：传质系数（m/s）
• $A_{release}$：释放面积（m²）
• $M$：物质摩尔质量（kg/mol）
• $R$：气体常数（8.314 J/mol·K）
• $T$：环境温度（K）
• $P_{eq}$：平衡蒸气压（Pa）
• $P_{air}$：空气中物质分压（Pa）
• $Q$：通风速率（h⁻¹）
• $V$：房间体积（m³）

已实现的Python代码 #

import math
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

class ExposureToVapourEvaporation:
    def __init__(self):
        # Input parameters
        self.frequency = 197  # per year
        self.exposure_duration = 0.75  # minute
        self.product_amount = 500  # g (amount of diluted product applied on a surface)
        self.weight_fraction_substance = 0.2  # fraction in the product
        self.room_volume = 1  # m³
        self.ventilation_rate = 0.5  # per hour
        self.inhalation_rate = 22.9 / 1000  # m³/min (converted from L/min)
        self.vapour_pressure = 0.0106  # Pa
        self.molecular_weight = 46.1 / 1000  # kg/mol (converted from g/mol)
        self.release_area = 0.002  # m²
        self.release_duration = 0.3  # minute
        self.application_temperature = 20 + 273.15  # K (converted from °C)
        self.mass_transfer_coefficient = 10 / 3600  # m/s (converted from m/h)
        self.body_weight = 9.8  # kg
        
        # Optional parameters
        self.is_product_used_in_dilution = True
        self.dilution = 1  # times (only used if is_product_used_in_dilution is True)
        self.molecular_weight_matrix = 22 / 1000  # kg/mol (converted from g/mol)
        self.is_pure_substance = True
        
        # Derived parameters
        self.is_constant_surface_area = True
        self.weight_fraction_solution = self.calculate_weight_fraction_solution()

    def calculate_weight_fraction_solution(self):
        if self.is_product_used_in_dilution:
            return self.weight_fraction_substance / self.dilution
        return self.weight_fraction_substance

    def calculate_equilibrium_vapour_pressure(self):
        return self.vapour_pressure

    def evaporation_ode(self, y, t):
        A_air, A_prod = y
        K = self.mass_transfer_coefficient
        P_eq = self.calculate_equilibrium_vapour_pressure()
        P_air = A_air * 8.314 * self.application_temperature / (self.molecular_weight * self.room_volume)
        
        if t < self.release_duration * 60:  # During release
            dA_air_dt = K * self.release_area * (self.molecular_weight / (8.314 * self.application_temperature)) * (P_eq - P_air) - (self.ventilation_rate / 3600) * self.room_volume * A_air
            dA_prod_dt = -dA_air_dt

        elif t == self.release_duration * 60:  # At the exact end of release
            dA_air_dt = 100 * K * self.release_area * (self.molecular_weight / (8.314 * self.application_temperature)) * (P_eq - P_air) - (self.ventilation_rate / 3600) * self.room_volume * A_air
        
        else:  # After release
            dA_air_dt = -(self.ventilation_rate / 3600) * self.room_volume * A_air
            dA_prod_dt = 0
        
        if A_prod <= 0:
            dA_prod_dt = 0
            dA_air_dt = -(self.ventilation_rate / 3600) * self.room_volume * A_air
        
        return [dA_air_dt, dA_prod_dt]

    def solve_evaporation(self):
        t = np.linspace(0, self.exposure_duration * 60, 1000)  # seconds
        t = np.sort(np.unique(np.append(t, self.release_duration * 60)))  # Ensure we have a point at the end of release
        y0 = [0, (self.product_amount / 1000) * self.weight_fraction_solution]  # kg
        solution = odeint(self.evaporation_ode, y0, t)
        return t, solution

    def calculate_metrics(self):
        t, solution = self.solve_evaporation()
        A_air = solution[:, 0]
        
        concentrations = A_air / self.room_volume * 1e6  # Convert to mg/m³
        mean_concentration = np.mean(concentrations)
        peak_concentration = np.max(concentrations)
        
        # Calculate TWA 15 min
        twa_15_min = mean_concentration if self.exposure_duration <= 15 else np.mean(concentrations[-int(15*60/self.exposure_duration):])
        
        # Calculate daily and yearly averages
        daily_average = mean_concentration * (self.exposure_duration / 1440)
        yearly_average = daily_average * (self.frequency / 365)
        
        # Calculate doses
        event_dose = (mean_concentration * self.inhalation_rate * self.exposure_duration) / self.body_weight
        day_dose = event_dose  # Assuming one event per day
        
        metrics = {
            "mean_event_concentration": mean_concentration,
            "peak_concentration_twa_15_min": twa_15_min,
            "mean_concentration_day": daily_average,
            "year_average_concentration": yearly_average,
            "external_event_dose": event_dose,
            "external_day_dose": day_dose
        }
        
        return metrics, t, concentrations

    def plot_air_concentration(self, t, concentrations):
        plt.figure(figsize=(10, 6))
        plt.plot(t / 60, concentrations)  # Convert time to minutes
        plt.axvline(x=self.release_duration, color='r', linestyle='--', label='End of Release')
        plt.axvline(x=self.exposure_duration, color='g', linestyle='--', label='End of Exposure')
        plt.title('Air Concentration Over Time')
        plt.xlabel('Time (minutes)')
        plt.ylabel('Concentration (mg/m³)')
        plt.legend()
        plt.grid(True)
        plt.show()

# Create an instance of the model and calculate metrics
model = ExposureToVapourEvaporation()
metrics, t, concentrations = model.calculate_metrics()

print("Calculated Results:")
for key, value in metrics.items():
    print(f"{key}: {value:.2e} mg/m³" if "concentration" in key else f"{key}: {value:.2e} mg/kg bw")

# Plot the air concentration over time
model.plot_air_concentration(t, concentrations)

结果对比 #

• 预期结果：浓度在0-0.3分钟释放阶段快速上升至约1.2×10⁻³ mg/m³，0.3-0.75分钟阶段缓慢下降
• 当前结果：浓度始终维持在约1×10⁻⁶ mg/m³的极低水平，无明显上升趋势

问题说明 #

已尝试调整代码中的数学逻辑，但始终无法匹配预期结果，怀疑原文档公式可能存在表述错误，寻求技术排查与修正帮助。

内容的提问来源于stack exchange，提问作者Nick