Мобильная версияhttps://www.rivm.nl/bibliotheek/rapporten/2017-0197.pdf
Выпустил эту версию модели и это лучшее, что мне досталось.
Код: Выделить всё
# Import necessary libraries
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# -------------------------------
# Constants
# -------------------------------
R_J_per_mol_K = 8.314 # Universal gas constant [J/(mol·K)]
# -------------------------------
# Scenario Selection
# -------------------------------
selected_scenario = 1 # Change this to select scenarios (1, 2, or 3)
# -------------------------------
# Input Parameters
# -------------------------------
# General Parameters (Common across all scenarios)
body_weight_kg = 60 # Body weight [kg]
frequency_per_year = 1 # Frequency of exposure [per year]
# Scenario Specific Parameters
if selected_scenario == 1:
# Scenario 1 (Base Scenario)
product_surface_area_cm2 = 10 # Product surface area [cm²]
product_thickness_mm = 5 # Product thickness [mm]
product_density_g_per_cm3 = 0.1 # Product density [g/cm³]
diffusion_coefficient_m2_per_s = 1E-08 # Diffusion coefficient [m²/s]
weight_fraction_substance = 0.6 # Weight fraction of the substance [fraction]
product_air_partition_coefficient = 0.5 # Product/air partition coefficient [-]
room_volume_m3 = 0.5 # Room volume [m³]
ventilation_rate_per_hour = 58 # Ventilation rate [per hour]
inhalation_rate_m3_per_hour = 10 # Inhalation rate [m³/hr]
start_exposure_h = 1 # Start exposure time [hours]
exposure_duration_h = 2 # Exposure duration [hours]
mass_transfer_coefficient_m_per_h = 10 # Mass transfer coefficient [m/h]
#Absorption
absorption_fraction = 1 # Absorption fraction [fraction]
# -------------------------------
# Unit Conversions
# -------------------------------
product_surface_area_m2 = product_surface_area_cm2 / 10000 # Convert cm² to m²
product_thickness_m = product_thickness_mm / 1000 # Convert mm to m
product_density_kg_per_m3 = product_density_g_per_cm3 * 1000 # Convert g/cm³ to kg/m³
mass_transfer_coefficient_m_per_s = mass_transfer_coefficient_m_per_h / 3600 # Convert m/h to m/s
ventilation_rate_per_s = ventilation_rate_per_hour / 3600 # Convert per hour to per second
start_exposure_s = start_exposure_h * 3600 # Convert hours to seconds
exposure_duration_s = exposure_duration_h * 3600 # Convert hours to seconds
total_emission_time_s = start_exposure_s + exposure_duration_s # Total time to simulate [s]
frequency_per_day = frequency_per_year / 365 # Frequency per day
inhalation_rate_m3_per_s = inhalation_rate_m3_per_hour / 3600 # Convert inhalation rate from m³/h to m³/s
# -------------------------------
# Initial Concentrations
# -------------------------------
volume_material_m3 = product_surface_area_m2 * product_thickness_m # Volume of material [m³]
mass_material_kg = volume_material_m3 * product_density_kg_per_m3 # Mass of material [kg]
mass_substance_initial_kg = mass_material_kg * weight_fraction_substance # Initial mass of substance [kg]
concentration_initial_kg_per_m3 = mass_substance_initial_kg / volume_material_m3 # Initial concentration [kg/m³]
N = 1000 # Reduced number of layers for faster computation
dx = product_thickness_m / (N - 1) # Spatial step size [m]
t_eval = np.linspace(0, total_emission_time_s, 2000) # Time points for evaluation
# -------------------------------
# Emission Model Definition
# -------------------------------
# Precompute constants
D = diffusion_coefficient_m2_per_s
hm = mass_transfer_coefficient_m_per_s
K = product_air_partition_coefficient
S = product_surface_area_m2
V = room_volume_m3
q = ventilation_rate_per_s
def emission_model(t, y):
C = y[:-1] # Concentration in the material layers [kg/m³]
y_air = y[-1] # Air concentration [kg/m³]
dCdt = np.zeros_like(C)
# No-flux boundary condition at x=0 (back surface)
dCdt[0] = D * (C[1] - C[0]) / dx**2 # Diffusion at back surface [kg/(m²·s²)]
# Diffusion in the material
for i in range(1, N - 1):
dCdt[i] = D * (C[i + 1] - 2 * C[i] + C[i - 1]) / dx**2 # Diffusion in layers [kg/(m²·s²)]
# Flux boundary condition at x=L (front surface)
J_diff = -D * (C[N - 1] - C[N - 2]) / dx # Diffusion flux at surface [kg/(m²·s)]
J_air = hm * (C[N - 1] / K - y_air) # Mass transfer to air [kg/(m²·s)]
dCdt[N - 1] = D * (C[N - 2] - C[N - 1]) / dx**2 - J_air / dx # Change in surface layer concentration [kg/m³/s]
dy_air_dt = (S / V) * hm * (C[N - 1] / K - y_air) - q * y_air # Change in air concentration [kg/m³/s]
return np.concatenate((dCdt, [dy_air_dt]))
# -------------------------------
# Initial Conditions
# -------------------------------
C0 = np.full(N, concentration_initial_kg_per_m3) # Initial concentration in material layers [kg/m³]
y_air_0 = 0.0 # Initial air concentration [kg/m³]
y0 = np.concatenate((C0, [y_air_0])) # Initial condition vector
# -------------------------------
# Solve the ODE system with BDF solver for speed and accuracy
# -------------------------------
solution = solve_ivp(emission_model, [0, total_emission_time_s], y0, t_eval=t_eval, method='BDF', vectorized=False)
# Extract results
C_material = solution.y[:-1, :] # Material concentrations over time [kg/m³]
y_air = solution.y[-1, :] # Air concentrations over time [kg/m³]
time = solution.t # Time vector [s]
# -------------------------------
# Extract Air Concentration during Exposure
# -------------------------------
exposure_indices = np.where((time >= start_exposure_s) & (time = start_exposure_s # Mask for exposure period
incremental_external_dose_mg_per_kg_bw = np.zeros_like(t_s_plot) # Initialize dose array
incremental_external_dose_mg_per_kg_bw[exposure_mask] = (C_t_mg_per_m3[exposure_mask] * delta_inhaled_volume_m3) / body_weight_kg # Incremental dose [mg/kg bw]
cumulative_external_dose_mg_per_kg_bw = np.cumsum(incremental_external_dose_mg_per_kg_bw) # Cumulative external dose [mg/kg bw]
cumulative_internal_dose_mg_per_kg_bw = cumulative_external_dose_mg_per_kg_bw * absorption_fraction # Cumulative internal dose [mg/kg bw]
# -------------------------------
# Visualization
# -------------------------------
fig, axs = plt.subplots(1, 3, figsize=(24, 6)) # Create a figure with three subplots
# 1. Air Concentration Over Time
axs[0].plot(t_s_plot / 3600, C_t_mg_per_m3, label='Air Concentration (mg/m³)', color='blue')
axs[0].set_title('Air Concentration Over Time')
axs[0].set_xlabel('Time (hours)')
axs[0].set_ylabel('Concentration (mg/m³)')
axs[0].grid(True)
axs[0].legend()
axs[0].axvline(x=start_exposure_s / 3600, color='red', linestyle='--', label='Exposure Start')
axs[0].axvline(x=(start_exposure_s + 15 * 60) / 3600, color='green', linestyle='--', label='15 min after Exposure Start')
axs[0].legend()
# 2. Cumulative External Dose Over Time
axs[1].plot(t_s_plot / 3600, cumulative_external_dose_mg_per_kg_bw, label='Cumulative External Dose (mg/kg bw)', color='orange')
axs[1].set_title('External Dose Over Time')
axs[1].set_xlabel('Time (hours)')
axs[1].set_ylabel('Dose (mg/kg bw)')
axs[1].grid(True)
axs[1].legend()
axs[1].axvline(x=start_exposure_s / 3600, color='red', linestyle='--', label='Exposure Start')
# 3. Cumulative Internal Dose Over Time (Adjusted for absorption)
axs[2].plot(t_s_plot / 3600, cumulative_internal_dose_mg_per_kg_bw, label='Cumulative Internal Dose (mg/kg bw)', color='green')
axs[2].set_title('Internal Dose Over Time')
axs[2].set_xlabel('Time (hours)')
axs[2].set_ylabel('Dose (mg/kg bw)')
axs[2].grid(True)
axs[2].legend()
axs[2].axvline(x=start_exposure_s / 3600, color='red', linestyle='--', label='Exposure Start')
plt.tight_layout() # Adjust layout for clarity
plt.show() # Show the plots
# -------------------------------
# End of Script
# -------------------------------
Код: Выделить всё
Mean Event Concentration: 1.3e-01 mg/m³
Mean Concentration on Day of Exposure: 1.1e-02 mg/m³
Year Average Concentration: 2.9e-05 mg/m³
External Event Dose: 4.3e-02 mg/kg bw
External Dose on Day of Exposure: 4.3e-02 mg/kg bw
Internal Event Dose: 4.3e-02 mg/kg bw
Internal Dose on Day of Exposure: 4.3e-02 mg/kg bw/day
Internal Year Average Dose: 1.2e-04 mg/kg bw/day
ОЖИДАЕТСЯ:
Код: Выделить всё
Mean event concentration
1.2 × 10⁻¹ mg/m³
average air concentration on exposure event. Note: depends strongly on chosen exposure duration
Mean concentration on day of exposure
9.6 × 10⁻³ mg/m³
average air concentration over the day (accounts for the number of events on one day)
Year average concentration
2.6 × 10⁻⁵ mg/m³
mean daily air concentration averaged over a year
External event dose
3.8 × 10⁻² mg/kg bw
the amount that can potentially be absorbed per kg body weight during one event
External dose on day of exposure
3.8 × 10⁻² mg/kg bw
the amount that can potentially be absorbed per kg body weight during one day
Internal event dose
3.8 × 10⁻² mg/kg bw
absorbed dose per kg body weight during one exposure event
Internal dose on day of exposure
3.8 × 10⁻² mg/kg bw/day
absorbed dose per kg body weight during one day. Note: these can be higher than the ‘event dose’ for exposure frequencies larger than 1 per day.
Internal year average dose
1.1 × 10⁻⁴ mg/kg bw/day
daily absorbed dose per kg body weight averaged over a year.
- Результаты не совпадают настолько, насколько хотелось бы.
- Он может загружаться очень медленно. Я адаптировал его так, чтобы он работал быстрее, но все же. Если я увеличу количество слоев, чтобы получить более точные результаты, это будет очень медленно и не намного точнее.
Так что на самом деле у меня всего несколько вопросов.
Любой, кто знаком с PDE/ODE и коэффициентами массопередачи. Как лучше всего решить эту проблему?
Я считаю, что этот сценарий, возможно, очень неэффективный способ воспроизвести модель.
Подробнее здесь: https://stackoverflow.com/questions/790 ... but-quickl
