Cough Droplets Partitioning¶
Welcome to the Cough Droplets Partitioning notebookâa step-by-step walkthrough to help you set up and run multi-component aerosol simulations using the Particula Python package.
This example focuses on modeling cough droplets, a mixture of water, salts, and biological materials, as they evaporate under different humidity conditions.
đ What This Notebook Covers¶
In this notebook, youâll learn how to:
- Model realistic aerosol particles (like cough droplets)
- Set up species properties (molar masses, densities, hygroscopicity)
- Construct gas-phase and particle-phase objects
- Run isothermal evaporation simulations
- Visualize how particle size and composition change over time
đĄ Why It Matters¶
Cough droplets are central to human exposure studies:
- They can carry infectious agents
- Influence indoor air quality
- Affect aerosol dynamics in confined environments like classrooms, hospitals, or aircraft cabins
đ Getting Started¶
This notebook is beginner-friendly. No prior experience with Particula or aerosol modeling is needed. Each step is clearly explained and includes code snippets to help you build intuition and experiment with your own setups.
In [23]:
Copied!
# đ§ Install Particula if you're running this in Google Colab
# Remove the comment below to enable installation
# !pip install particula[extra] --quiet
# đŠ Import necessary libraries
import copy # to safely duplicate Python objects
import matplotlib.pyplot as plt # for plotting
import numpy as np # numerical tools
import particula as par # the main Particula package
from tqdm import tqdm # optional: for progress bars during simulation
# đš Set default plot styling using Tailwind-inspired colors from Particula
TAILWIND = par.util.colors.TAILWIND
base_color = TAILWIND["gray"]["600"]
plt.rcParams.update({
"text.color": base_color,
"axes.labelcolor": base_color,
"figure.figsize": (5, 4),
"font.size": 14,
"axes.edgecolor": base_color,
"axes.labelcolor": base_color,
"xtick.color": base_color,
"ytick.color": base_color,
"pdf.fonttype": 42,
"ps.fonttype": 42,
})
# đ§ Install Particula if you're running this in Google Colab
# Remove the comment below to enable installation
# !pip install particula[extra] --quiet
# đŠ Import necessary libraries
import copy # to safely duplicate Python objects
import matplotlib.pyplot as plt # for plotting
import numpy as np # numerical tools
import particula as par # the main Particula package
from tqdm import tqdm # optional: for progress bars during simulation
# đš Set default plot styling using Tailwind-inspired colors from Particula
TAILWIND = par.util.colors.TAILWIND
base_color = TAILWIND["gray"]["600"]
plt.rcParams.update({
"text.color": base_color,
"axes.labelcolor": base_color,
"figure.figsize": (5, 4),
"font.size": 14,
"axes.edgecolor": base_color,
"axes.labelcolor": base_color,
"xtick.color": base_color,
"ytick.color": base_color,
"pdf.fonttype": 42,
"ps.fonttype": 42,
})
đ§Ș Defining Chemical Species in Cough Droplets¶
In this section, we define the chemical makeup of the simulated cough droplets. These are multi-component aerosols made up of water, salts, buffers, proteins, and other organics commonly found in human respiratory fluids.
đŹ Chemical Components¶
We create a list of 12 representative chemicals that capture the diversity of real cough droplet composition:
- Water (main solvent)
- Inorganic salts like NaCl, KHCOâ, KHâPOâ, Ca(HCOâ)â, Mg(HCOâ)â
- Organic compounds like urea, glucose, and lactate
- Biological macromolecules like lysozyme, IgA, and cholesterol
đ Initial Mass Fractions¶
For each chemical, we define a range of initial mass fractions that represent its likely concentration in fresh respiratory droplets. These values are approximate and intended to capture real-world variability.
âïž Automatic Property Retrieval¶
Using Particulaâs get_chemical_search() and get_chemical_stp_properties() utilities, we fetch:
- Molar mass (kg/mol)
- Density (kg/mÂł)
- Surface tension (N/m)
- Vapor pressure (Pa)
- Molecular identity (e.g. SMILES, CAS)
These properties are needed to accurately simulate evaporation, condensation, and phase partitioning.
đ Manual Adjustment¶
We manually correct the surface tension for NaCl, since tabulated value was for the solid phase. Here, we use a typical value of 0.091 N/m.
đ§Ÿ Output¶
For each species, the notebook prints:
- Name, CAS identifier, SMILES string
- Molar mass, density, surface tension, and vapor pressure
This lets you verify the chemical setup before moving forward with building the aerosol system.
In [24]:
Copied!
list_of_chemicals = [
"Water",
"NaCl",
"potassium bicarbonate (KHCO3)",
"potassium dihydrogen phosphate (KH2PO4)", # K+ + H2PO4â
"calcium bicarbonate [Ca(HCO3)2]",
"magnesium bicarbonate [Mg(HCO3)2]",
"Urea",
"Lactate",
"Glucose",
"Lysozyme",
"IgA",
"cholesterol",
]
initial_mass_fraction_bounds = [
(0.985, 0.995), # Water
(0.0006, 0.0023), # NaCl
(0.0007, 0.0030), # KHCO3
(0.0002, 0.0027), # KH2PO4
(0.00002, 0.00045), # Ca(HCO3)2
(0.000006, 0.00003), # Mg(HCO3)2
(0.0002, 0.0015), # Urea
(0.00005, 0.0004), # Lactate
(0.00001, 0.0002), # Glucose
(0.000005, 0.00002), # Lysozyme
(0.00002, 0.0002), # IgA
(0.00005, 0.0005), # Cholesterol
]
density_array = np.array([])
molar_mass_array = np.array([])
surface_tension_array = np.array([])
short_name = []
chemical_dict = {}
# Get the CAS numbers for each chemical
for chem in list_of_chemicals:
cas = par.util.get_chemical_search(chem)
print(f"{chem}: {cas}")
chemical_dict[chem] = par.util.get_chemical_stp_properties(cas)
# Store the molar mass and density
molar_mass_array = np.append(
molar_mass_array, chemical_dict[chem]["molar_mass"]
)
density_array = np.append(density_array, chemical_dict[chem]["density"])
surface_tension_array = np.append(
surface_tension_array, chemical_dict[chem]["surface_tension"]
)
short_name.append(chemical_dict[chem]["name"])
# correct NaCl surface tension
surface_tension_array[1] = 0.091 # N/m, typical value
# Print the chemical properties
for chem, props in chemical_dict.items():
print(f"{chem}:")
print(f" Name: {props['name']}")
print(f" CAS: {props['cas_name']}")
print(f" Smiles: {props['smiles']}")
print(f" Molar Mass: {props['molar_mass']:.4f} kg/mol")
print(f" Density: {props['density']:.2f} kg/mÂł")
print(f" Surface Tension: {props['surface_tension']} N/m")
print(f" Vapor Pressure: {props['pure_vapor_pressure']} Pa\n")
list_of_chemicals = [
"Water",
"NaCl",
"potassium bicarbonate (KHCO3)",
"potassium dihydrogen phosphate (KH2PO4)", # K+ + H2PO4â
"calcium bicarbonate [Ca(HCO3)2]",
"magnesium bicarbonate [Mg(HCO3)2]",
"Urea",
"Lactate",
"Glucose",
"Lysozyme",
"IgA",
"cholesterol",
]
initial_mass_fraction_bounds = [
(0.985, 0.995), # Water
(0.0006, 0.0023), # NaCl
(0.0007, 0.0030), # KHCO3
(0.0002, 0.0027), # KH2PO4
(0.00002, 0.00045), # Ca(HCO3)2
(0.000006, 0.00003), # Mg(HCO3)2
(0.0002, 0.0015), # Urea
(0.00005, 0.0004), # Lactate
(0.00001, 0.0002), # Glucose
(0.000005, 0.00002), # Lysozyme
(0.00002, 0.0002), # IgA
(0.00005, 0.0005), # Cholesterol
]
density_array = np.array([])
molar_mass_array = np.array([])
surface_tension_array = np.array([])
short_name = []
chemical_dict = {}
# Get the CAS numbers for each chemical
for chem in list_of_chemicals:
cas = par.util.get_chemical_search(chem)
print(f"{chem}: {cas}")
chemical_dict[chem] = par.util.get_chemical_stp_properties(cas)
# Store the molar mass and density
molar_mass_array = np.append(
molar_mass_array, chemical_dict[chem]["molar_mass"]
)
density_array = np.append(density_array, chemical_dict[chem]["density"])
surface_tension_array = np.append(
surface_tension_array, chemical_dict[chem]["surface_tension"]
)
short_name.append(chemical_dict[chem]["name"])
# correct NaCl surface tension
surface_tension_array[1] = 0.091 # N/m, typical value
# Print the chemical properties
for chem, props in chemical_dict.items():
print(f"{chem}:")
print(f" Name: {props['name']}")
print(f" CAS: {props['cas_name']}")
print(f" Smiles: {props['smiles']}")
print(f" Molar Mass: {props['molar_mass']:.4f} kg/mol")
print(f" Density: {props['density']:.2f} kg/mÂł")
print(f" Surface Tension: {props['surface_tension']} N/m")
print(f" Vapor Pressure: {props['pure_vapor_pressure']} Pa\n")
Water: Water NaCl: NaCl potassium bicarbonate (KHCO3): potassium bicarbonate (KHCO3) potassium dihydrogen phosphate (KH2PO4): potassium dihydrogen phosphate (KH2PO4) calcium bicarbonate [Ca(HCO3)2]: calcium carbonate [usan] magnesium bicarbonate [Mg(HCO3)2]: magnesium carbonate (mgco3) Urea: Urea Lactate: Lactate Glucose: Glucose Lysozyme: soy-dome IgA: IgA cholesterol: cholesterol Water: Name: water CAS: 7732-18-5 Smiles: O Molar Mass: 0.0180 kg/mol Density: 997.06 kg/mÂł Surface Tension: 0.07197220523022962 N/m Vapor Pressure: 3169.9293388738784 Pa NaCl: Name: sodium chloride CAS: 7647-14-5 Smiles: [Na+].[Cl-] Molar Mass: 0.0584 kg/mol Density: 2169.99 kg/mÂł Surface Tension: 0.35290821232462 N/m Vapor Pressure: 2.245741871175006e-30 Pa potassium bicarbonate (KHCO3): Name: potassium bicarbonate CAS: 298-14-6 Smiles: C(=O)(O)[O-].[K+] Molar Mass: 0.1001 kg/mol Density: 2170.00 kg/mÂł Surface Tension: None N/m Vapor Pressure: None Pa potassium dihydrogen phosphate (KH2PO4): Name: potassium dihydrogen phosphate CAS: 7778-77-0 Smiles: OP(=O)(O)[O-].[K+] Molar Mass: 0.1361 kg/mol Density: 2340.01 kg/mÂł Surface Tension: None N/m Vapor Pressure: None Pa calcium bicarbonate [Ca(HCO3)2]: Name: calcium carbonate CAS: 471-34-1 Smiles: C(=O)([O-])[O-].[Ca+2] Molar Mass: 0.1001 kg/mol Density: 2710.00 kg/mÂł Surface Tension: None N/m Vapor Pressure: None Pa magnesium bicarbonate [Mg(HCO3)2]: Name: magnesium carbonate CAS: 546-93-0 Smiles: C(=O)([O-])[O-].[Mg+2] Molar Mass: 0.0843 kg/mol Density: 3010.00 kg/mÂł Surface Tension: None N/m Vapor Pressure: None Pa Urea: Name: urea CAS: 57-13-6 Smiles: C(=O)(N)N Molar Mass: 0.0601 kg/mol Density: 932.80 kg/mÂł Surface Tension: 0.06630868366869905 N/m Vapor Pressure: 9.647479358454168 Pa Lactate: Name: lactic acid CAS: 50-21-5 Smiles: CC(C(=O)O)O Molar Mass: 0.0901 kg/mol Density: 1146.11 kg/mÂł Surface Tension: 0.06892887098077884 N/m Vapor Pressure: 0.42214021514470007 Pa Glucose: Name: glucose CAS: 50-99-7 Smiles: O=C[C@H](O)[C@@H](O)[C@H](O)[C@H](O)CO Molar Mass: 0.1802 kg/mol Density: 1618.57 kg/mÂł Surface Tension: 0.10078381855423499 N/m Vapor Pressure: 3.472643152420593e-11 Pa Lysozyme: Name: hexachlorophene CAS: 70-30-4 Smiles: C1=C(C(=C(C(=C1Cl)Cl)CC2=C(C(=CC(=C2Cl)Cl)Cl)O)O)Cl Molar Mass: 0.4069 kg/mol Density: 2009.99 kg/mÂł Surface Tension: 0.07787535981784431 N/m Vapor Pressure: 9.468715811988936e-14 Pa IgA: Name: isoguanine CAS: 3373-53-3 Smiles: C1=NC2=NC(=O)NC(=C2N1)N Molar Mass: 0.1511 kg/mol Density: 815.71 kg/mÂł Surface Tension: 0.1373263034398279 N/m Vapor Pressure: 0.00012951501534944908 Pa cholesterol: Name: cholesterol CAS: 57-88-5 Smiles: CC(C)CCC[C@@H](C)[C@H]1CC[C@H]2[C@@H]3CC=C4C[C@@H](O)CC[C@]4(C)[C@H]3CC[C@]12C Molar Mass: 0.3867 kg/mol Density: 1064.37 kg/mÂł Surface Tension: 0.03942801914463943 N/m Vapor Pressure: 8.548725890318368e-07 Pa
𧟠Step 1: Set Up Particle and Species Properties¶
đČ Reproducibility¶
To ensure that simulations can be reproduced exactly (e.g. for comparison or debugging), we seed NumPyâs random number generator:
np.random.seed(100)
𧫠1a. Define Particle Sampling and Environmental Conditions¶
We now define the simulation scale and temperature conditions:
number_of_samples: number of particles weâll track individuallytotal_number_per_cm3: assumed number concentration in the airsimulation_volume: computed simulation volume based on the abovetemperature: constant system temperature (isothermal)
This configuration controls the scale and realism of the simulation.
đĄ 1b. Assign Vapor Pressure Strategies¶
Each species needs a vapor pressure strategy, which defines how readily it transitions between the gas and particle phases.
We retrieve each speciesâ vapor pressure (at room temperature), and if missing, assign a small fallback value. Then we use Particulaâs ConstantVaporPressureBuilder to build the vapor behavior for each compound.
â ïž Surface tensions are also cleaned: if missing, a default of
0.072 N/mis applied (typical for water at 25°C).
In [25]:
Copied!
# %% Step 1: Define system size, temperature, and vapor pressure strategies
# đČ Ensure reproducibility
np.random.seed(100)
# đą Simulation settings
number_of_samples = 10_000 # Number of particles to simulate
total_number_per_cm3 = 1e-4 # Particle number concentration [#/cmÂł]
simulation_volume = (
number_of_samples / total_number_per_cm3 * 1e-6
) # Simulation volume in mÂł (converted from cmÂł)
temperature = 298.15 # System temperature in Kelvin (25°C)
# đĄ Build vapor pressure strategies for each chemical
vapor_pressure_strategy_list = []
for chem_i, props in chemical_dict.items():
vapor_pressure_temp = props["pure_vapor_pressure"]
if vapor_pressure_temp is None:
print(
f"â ïž Warning: No vapor pressure data for {chem_i}, using default."
)
vapor_pressure_temp = 1e-30 # fallback for non-volatile species
vapor_pressure_i = (
par.gas.ConstantVaporPressureBuilder()
.set_vapor_pressure(vapor_pressure_temp, "Pa")
.build()
)
vapor_pressure_strategy_list.append(vapor_pressure_i)
# đ§ Ensure surface tension values are set
surface_tension_array = np.array(
[0.072 if x is None else x for x in surface_tension_array]
)
# %% Step 1: Define system size, temperature, and vapor pressure strategies
# đČ Ensure reproducibility
np.random.seed(100)
# đą Simulation settings
number_of_samples = 10_000 # Number of particles to simulate
total_number_per_cm3 = 1e-4 # Particle number concentration [#/cmÂł]
simulation_volume = (
number_of_samples / total_number_per_cm3 * 1e-6
) # Simulation volume in mÂł (converted from cmÂł)
temperature = 298.15 # System temperature in Kelvin (25°C)
# đĄ Build vapor pressure strategies for each chemical
vapor_pressure_strategy_list = []
for chem_i, props in chemical_dict.items():
vapor_pressure_temp = props["pure_vapor_pressure"]
if vapor_pressure_temp is None:
print(
f"â ïž Warning: No vapor pressure data for {chem_i}, using default."
)
vapor_pressure_temp = 1e-30 # fallback for non-volatile species
vapor_pressure_i = (
par.gas.ConstantVaporPressureBuilder()
.set_vapor_pressure(vapor_pressure_temp, "Pa")
.build()
)
vapor_pressure_strategy_list.append(vapor_pressure_i)
# đ§ Ensure surface tension values are set
surface_tension_array = np.array(
[0.072 if x is None else x for x in surface_tension_array]
)
â ïž Warning: No vapor pressure data for potassium bicarbonate (KHCO3), using default. â ïž Warning: No vapor pressure data for potassium dihydrogen phosphate (KH2PO4), using default. â ïž Warning: No vapor pressure data for calcium bicarbonate [Ca(HCO3)2], using default. â ïž Warning: No vapor pressure data for magnesium bicarbonate [Mg(HCO3)2], using default.
đŹ Step 2: Set Up the Gas-Phase Environment¶
To simulate how particles exchange mass with the air (e.g. water evaporating), we need to define the surrounding gas-phase composition and thermodynamic conditions.
đ§ Relative Humidity for Water¶
We begin by setting near-zero concentrations for all species in the gas phase to simulate dry air. Then we adjust the water vapor concentration to simulate a 30% relative humidity environment:
- The saturation concentration for water is calculated using its vapor pressure strategy.
- The RH is applied by multiplying by 0.3.
đ§Ș Building Gas-Phase Species¶
We use GasSpeciesBuilder() to define:
- Species names and molar masses
- Their ability to partition
- Vapor pressure strategies
- Initial gas-phase concentrations
đĄ Creating the Atmosphere¶
Finally, we construct the atmosphere with temperature, pressure, and link the gas-phase species. This environment will be used to govern evaporation/condensation dynamics during the simulation.
In [26]:
Copied!
# %% Step 2: Define gas-phase concentrations and thermodynamic environment
# đš Initialize gas-phase concentrations near zero
concentration_gas = np.ones(len(list_of_chemicals)) * 1e-30 # kg/mÂł
# đ§ Set water vapor to 30% RH
isat_water = vapor_pressure_strategy_list[0].saturation_concentration(
molar_mass_array[0], temperature
)
concentration_gas[0] = isat_water * 0.3 # 30% RH
# đ§Ș Build gas-phase species with vapor strategies and concentrations
gas_species = (
par.gas.GasSpeciesBuilder()
.set_name(np.array(list_of_chemicals))
.set_molar_mass(molar_mass_array, "kg/mol")
.set_partitioning(True)
.set_vapor_pressure_strategy(vapor_pressure_strategy_list)
.set_concentration(concentration_gas, "kg/m^3")
.build()
)
# đĄ Create atmosphere: includes temperature, pressure, and gas species
atmosphere = (
par.gas.AtmosphereBuilder()
.set_more_partitioning_species(gas_species)
.set_temperature(temperature, "K")
.set_pressure(1, "atm")
.build()
)
# %% Step 2: Define gas-phase concentrations and thermodynamic environment
# đš Initialize gas-phase concentrations near zero
concentration_gas = np.ones(len(list_of_chemicals)) * 1e-30 # kg/mÂł
# đ§ Set water vapor to 30% RH
isat_water = vapor_pressure_strategy_list[0].saturation_concentration(
molar_mass_array[0], temperature
)
concentration_gas[0] = isat_water * 0.3 # 30% RH
# đ§Ș Build gas-phase species with vapor strategies and concentrations
gas_species = (
par.gas.GasSpeciesBuilder()
.set_name(np.array(list_of_chemicals))
.set_molar_mass(molar_mass_array, "kg/mol")
.set_partitioning(True)
.set_vapor_pressure_strategy(vapor_pressure_strategy_list)
.set_concentration(concentration_gas, "kg/m^3")
.build()
)
# đĄ Create atmosphere: includes temperature, pressure, and gas species
atmosphere = (
par.gas.AtmosphereBuilder()
.set_more_partitioning_species(gas_species)
.set_temperature(temperature, "K")
.set_pressure(1, "atm")
.build()
)
đ Step 3: Define Particle Composition and Build the Aerosol Phase¶
Now that the gas-phase environment is ready, we define the particle-phase systemârealistic cough droplets containing a mix of water, salts, organics, and proteins.
đĄ 3a. Size Distribution of Cough Droplets¶
Human coughs produce droplets from different regions of the respiratory tract, which can be modeled using lognormal size distributions:
- Bronchiolar: smaller mode (~0.3 ÎŒm radius)
- Laryngeal: medium mode (~4 ÎŒm radius)
- Oropharyngeal: large mode (~40 ÎŒm radius)
We use par.particles.get_lognormal_sample_distribution() to sample droplet radii based on these modes, weighted by their relative occurrence.
đ§ 3b. Assign Mass Composition¶
We assume all particles start as water droplets, then apply random mass fractions to simulate different droplet compositions. These fractions are sampled from the ranges defined earlier for each chemical species.
âïž 3c. Define Activity and Surface Strategies¶
We use Îș-Köhler theory to estimate hygroscopicity, setting kappa values for each species. This is fed into ActivityKappaParameterBuilder.
We also define a surface tension strategy that accounts for the mass of each component using SurfaceStrategyMassBuilder.
𧱠Create the Resolved Particle and Aerosol System¶
We wrap everything into a ResolvedParticleMassRepresentation and then build the full Aerosol object by combining:
- Particle-phase definition (
resolved_masses) - Gas-phase environment (
atmosphere)
This sets up the full simulation-ready aerosol system.
In [27]:
Copied!
# %% Step 3: Define particle composition and create aerosol object
# đ 3a. Lognormal size distributions for cough droplet modes
# Modes: bronchiolar, laryngeal, oropharyngeal
modes = np.array([0.6, 8.0, 80.0]) / 2 * 1e-6 # radii in meters
geometric_standard_deviations = np.array([1.65, 2.4, 1.6])
relative_number = np.array([15, 2, 0.2]) # relative weights
# Generate sampled radii
radii_cough = par.particles.get_lognormal_sample_distribution(
mode=modes,
geometric_standard_deviation=geometric_standard_deviations,
number_of_particles=relative_number,
number_of_samples=number_of_samples,
)
# đ§ 3b. Assign initial water mass and scale by chemical mass fractions
mass_water = 4 / 3 * np.pi * (radii_cough**3) * density_array[0]
mass_speciation = np.repeat(
mass_water[:, np.newaxis], len(list_of_chemicals), axis=1
)
# Apply random mass fractions from defined bounds
for i, (chem, bounds) in enumerate(
zip(list_of_chemicals, initial_mass_fraction_bounds)
):
mass_fraction = np.random.uniform(bounds[0], bounds[1], number_of_samples)
mass_speciation[:, i] *= mass_fraction
# âïž 3c. Define activity strategy using kappa parameters
kappa_parameter = np.array(
[0.0, 1.1, 0.6, 0.6, 0.6, 0.6, 0.3, 0.2, 0.1, 0.1, 0.1, 0.01]
)
activity_strategy = (
par.particles.ActivityKappaParameterBuilder()
.set_molar_mass(molar_mass_array, "kg/mol")
.set_density(density_array, "kg/m^3")
.set_kappa(kappa_parameter)
.set_water_index(0)
.build()
)
# Define surface tension strategy
surface_strategy = (
par.particles.SurfaceStrategyMassBuilder()
.set_surface_tension(surface_tension_array, "N/m")
.set_density(density_array, "kg/m^3")
.set_phase_index(np.arange(len(list_of_chemicals)))
.build()
)
# 𧱠3d. Build resolved particle object
resolved_masses = (
par.particles.ResolvedParticleMassRepresentationBuilder()
.set_distribution_strategy(par.particles.ParticleResolvedSpeciatedMass())
.set_activity_strategy(activity_strategy)
.set_surface_strategy(surface_strategy)
.set_mass(mass_speciation, "kg")
.set_density(density_array, "kg/m^3")
.set_charge(0)
.set_volume(simulation_volume, "m^3")
.build()
)
# Combine gas and particle phases into full aerosol system
aerosol = par.Aerosol(atmosphere=atmosphere, particles=resolved_masses)
# Optional: print summary of aerosol system
print(aerosol)
# %% Step 3: Define particle composition and create aerosol object
# đ 3a. Lognormal size distributions for cough droplet modes
# Modes: bronchiolar, laryngeal, oropharyngeal
modes = np.array([0.6, 8.0, 80.0]) / 2 * 1e-6 # radii in meters
geometric_standard_deviations = np.array([1.65, 2.4, 1.6])
relative_number = np.array([15, 2, 0.2]) # relative weights
# Generate sampled radii
radii_cough = par.particles.get_lognormal_sample_distribution(
mode=modes,
geometric_standard_deviation=geometric_standard_deviations,
number_of_particles=relative_number,
number_of_samples=number_of_samples,
)
# đ§ 3b. Assign initial water mass and scale by chemical mass fractions
mass_water = 4 / 3 * np.pi * (radii_cough**3) * density_array[0]
mass_speciation = np.repeat(
mass_water[:, np.newaxis], len(list_of_chemicals), axis=1
)
# Apply random mass fractions from defined bounds
for i, (chem, bounds) in enumerate(
zip(list_of_chemicals, initial_mass_fraction_bounds)
):
mass_fraction = np.random.uniform(bounds[0], bounds[1], number_of_samples)
mass_speciation[:, i] *= mass_fraction
# âïž 3c. Define activity strategy using kappa parameters
kappa_parameter = np.array(
[0.0, 1.1, 0.6, 0.6, 0.6, 0.6, 0.3, 0.2, 0.1, 0.1, 0.1, 0.01]
)
activity_strategy = (
par.particles.ActivityKappaParameterBuilder()
.set_molar_mass(molar_mass_array, "kg/mol")
.set_density(density_array, "kg/m^3")
.set_kappa(kappa_parameter)
.set_water_index(0)
.build()
)
# Define surface tension strategy
surface_strategy = (
par.particles.SurfaceStrategyMassBuilder()
.set_surface_tension(surface_tension_array, "N/m")
.set_density(density_array, "kg/m^3")
.set_phase_index(np.arange(len(list_of_chemicals)))
.build()
)
# 𧱠3d. Build resolved particle object
resolved_masses = (
par.particles.ResolvedParticleMassRepresentationBuilder()
.set_distribution_strategy(par.particles.ParticleResolvedSpeciatedMass())
.set_activity_strategy(activity_strategy)
.set_surface_strategy(surface_strategy)
.set_mass(mass_speciation, "kg")
.set_density(density_array, "kg/m^3")
.set_charge(0)
.set_volume(simulation_volume, "m^3")
.build()
)
# Combine gas and particle phases into full aerosol system
aerosol = par.Aerosol(atmosphere=atmosphere, particles=resolved_masses)
# Optional: print summary of aerosol system
print(aerosol)
Gas mixture at 298.15 K, 101325.0 Pa, partitioning=['Water' 'NaCl' 'potassium bicarbonate (KHCO3)' 'potassium dihydrogen phosphate (KH2PO4)' 'calcium bicarbonate [Ca(HCO3)2]' 'magnesium bicarbonate [Mg(HCO3)2]' 'Urea' 'Lactate' 'Glucose' 'Lysozyme' 'IgA' 'cholesterol'], gas_only_species=None Particle Representation: Strategy: ParticleResolvedSpeciatedMass Activity: ActivityKappaParameter Surface: SurfaceStrategyMass Mass Concentration: 1.592e-09 [kg/m^3] Number Concentration: 1.000e+02 [#/m^3]
đ Step 4: Visualize the Initial Particle Size Distribution¶
Before simulating evaporation or chemical evolution, letâs examine the initial state of the aerosol system weâve just created.
đ§Ÿ Summary of the Initialized System¶
At this point, the system contains:
Gas Mixture at:
- Temperature:
298.15 K - Pressure:
101325 Pa(1 atm) - Partitioning Species: Water, salts, organics, proteins
- Temperature:
Particle Representation:
- Strategy:
ParticleResolvedSpeciatedMass - Activity Model:
ActivityKappaParameter - Surface Model:
SurfaceStrategyMass - Mass Concentration: ~
1.6 Ă 10â»âč kg/mÂł - Number Concentration:
1.0 Ă 10ÂČ #/mÂł
- Strategy:
đ Plotting Initial Size Distribution¶
To visualize the initial aerosol size distribution:
- We extract particle radii using
aerosol.particles.get_radius() - Convert to log10 scale (standard for aerosol data)
- Plot a histogram to show how particle sizes are distributed
In [28]:
Copied!
# %% Step 4: Plot initial particle size distribution
fig, ax = plt.subplots()
# Log-scale histogram of particle radii
ax.hist(
np.log10(aerosol.particles.get_radius()), # convert to log10
bins=50,
density=False,
alpha=0.5,
)
ax.set_xlabel("logââ(Radius [m])")
ax.set_ylabel("Bin Counts")
ax.set_title("Initial Particle Size Distribution")
plt.tight_layout()
plt.show()
# %% Step 4: Plot initial particle size distribution
fig, ax = plt.subplots()
# Log-scale histogram of particle radii
ax.hist(
np.log10(aerosol.particles.get_radius()), # convert to log10
bins=50,
density=False,
alpha=0.5,
)
ax.set_xlabel("logââ(Radius [m])")
ax.set_ylabel("Bin Counts")
ax.set_title("Initial Particle Size Distribution")
plt.tight_layout()
plt.show()
đ Step 5: Define Aerosol Dynamics â Condensation & Coagulation¶
To simulate how aerosol particles evolve in time, we need to define physical processes that change their size, mass, and composition.
In this example, we model:
đ§ 5a. Isothermal Condensation¶
We simulate mass transfer from the gas phase to particles via isothermal condensation. Particulaâs CondensationIsothermal model accounts for:
- Molar mass of each species
- Diffusion coefficient (typically ~2e-5 mÂČ/s for air)
- Accommodation coefficient (set to 1 for full sticking)
This sets up a MassCondensation process that can be used during time evolution.
âȘ 5b. Brownian Coagulation¶
We also define a Brownian coagulation process, where particles randomly collide and merge based on their sizes and thermal motion. This uses Particulaâs:
BrownianCoagulationBuilder()withparticle_resolveddistribution- Wrapped into a
Coagulationobject
In [29]:
Copied!
# %% Step 5: Set up physical aerosol processes
# đ§ Isothermal Condensation
condensation_strategy = par.dynamics.CondensationIsothermal(
molar_mass=np.array(molar_mass_array),
diffusion_coefficient=2e-5, # typical air value [mÂČ/s]
accommodation_coefficient=1, # full sticking efficiency
# skip_partitioning_indices=[] # optional: exclude species
)
condensation_process = par.dynamics.MassCondensation(condensation_strategy)
# âȘ Brownian Coagulation
coagulation_strategy = (
par.dynamics.BrownianCoagulationBuilder()
.set_distribution_type("particle_resolved") # use individual particles
.build()
)
coagulation_process = par.dynamics.Coagulation(coagulation_strategy)
# %% Step 5: Set up physical aerosol processes
# đ§ Isothermal Condensation
condensation_strategy = par.dynamics.CondensationIsothermal(
molar_mass=np.array(molar_mass_array),
diffusion_coefficient=2e-5, # typical air value [mÂČ/s]
accommodation_coefficient=1, # full sticking efficiency
# skip_partitioning_indices=[] # optional: exclude species
)
condensation_process = par.dynamics.MassCondensation(condensation_strategy)
# âȘ Brownian Coagulation
coagulation_strategy = (
par.dynamics.BrownianCoagulationBuilder()
.set_distribution_type("particle_resolved") # use individual particles
.build()
)
coagulation_process = par.dynamics.Coagulation(coagulation_strategy)
ⱠStep 6: Run the Simulation Loop¶
Now that the aerosol system and physical processes are set up, we simulate how particles change over time due to condensation and optionally coagulation.
âïž Simulation Settings¶
- Duration: 120 seconds
- Steps: 200
- Condensation sub-steps: 100 (to improve resolution for fast vapor interactions)
- Coagulation sub-steps: 1 (low priority here due to low concentration)
time_array = np.linspace(0, total_time, total_steps)
time_step = total_time / total_steps
đŸ Data Tracking¶
At each time step, we record:
- Particle size distribution (histogram of radii)
- Species mass across all particles
- Gas-phase saturation ratio for each chemical
The size distribution is stored in log-spaced bins ranging from 10â»âž to 10â»âŽ m â a typical range for atmospheric aerosol particles.
đ§Ș Optional: Toggle Physical Processes¶
Users can choose to:
- Fix the water vapor concentration (e.g. simulate a constant RH environment)
- Disable condensation or coagulation individually to isolate effects
These are included as commented lines in the code and easy to turn on/off.
In [30]:
Copied!
# %% Step 6: Run the simulation loop
print(aerosol) # print system state before simulation
# đ Deep copy for simulation
aerosol_process = copy.deepcopy(aerosol)
# Simulation parameters
total_time = 120 # seconds
total_steps = 200
condensation_sub_step = 100
coagulation_sub_step = 1
time_step = total_time / total_steps
time_array = np.linspace(0, total_time, total_steps)
time_len = len(time_array)
# đ Initialize arrays for tracking outputs
bins_lognormal = np.logspace(-8, -4, 100)
distribution_counts = np.zeros((time_len, len(bins_lognormal) - 1))
species_mass = np.zeros((time_len, len(list_of_chemicals)))
vapor_saturation = np.zeros((time_len, len(list_of_chemicals)))
# đ Time integration loop
for step, t in enumerate(
tqdm(time_array, desc="Running Sim", mininterval=0.5)
):
if step > 0:
# Optional: fix water vapor concentration (e.g. controlled RH)
# aerosol_process.atmosphere.partitioning_species.concentration[0] = (
# aerosol.atmosphere.partitioning_species.concentration[0]
# )
# đ§ Execute condensation
aerosol_process = condensation_process.execute(
aerosol=aerosol_process,
time_step=time_step,
sub_steps=condensation_sub_step,
)
# âȘ Optional: execute coagulation
# aerosol_process = coagulation_process.execute(
# aerosol=aerosol_process,
# time_step=time_step,
# sub_steps=coagulation_sub_step,
# )
# Record size distribution
distribution_counts[step, :], edges = np.histogram(
aerosol_process.particles.get_radius(clone=True),
bins=bins_lognormal,
density=False,
)
# Record species mass in particles
species_temp = aerosol_process.particles.get_species_mass(clone=True)
species_mass[step, :] = np.sum(species_temp, axis=0)
# Record vapor saturation ratios
vapor_saturation[step, :] = (
aerosol_process.atmosphere.partitioning_species.get_saturation_ratio(
temperature=temperature
)
)
# Normalize counts to number concentration (#/mÂł)
concentrations = distribution_counts / simulation_volume
# Final system state
print(aerosol_process)
# %% Step 6: Run the simulation loop
print(aerosol) # print system state before simulation
# đ Deep copy for simulation
aerosol_process = copy.deepcopy(aerosol)
# Simulation parameters
total_time = 120 # seconds
total_steps = 200
condensation_sub_step = 100
coagulation_sub_step = 1
time_step = total_time / total_steps
time_array = np.linspace(0, total_time, total_steps)
time_len = len(time_array)
# đ Initialize arrays for tracking outputs
bins_lognormal = np.logspace(-8, -4, 100)
distribution_counts = np.zeros((time_len, len(bins_lognormal) - 1))
species_mass = np.zeros((time_len, len(list_of_chemicals)))
vapor_saturation = np.zeros((time_len, len(list_of_chemicals)))
# đ Time integration loop
for step, t in enumerate(
tqdm(time_array, desc="Running Sim", mininterval=0.5)
):
if step > 0:
# Optional: fix water vapor concentration (e.g. controlled RH)
# aerosol_process.atmosphere.partitioning_species.concentration[0] = (
# aerosol.atmosphere.partitioning_species.concentration[0]
# )
# đ§ Execute condensation
aerosol_process = condensation_process.execute(
aerosol=aerosol_process,
time_step=time_step,
sub_steps=condensation_sub_step,
)
# âȘ Optional: execute coagulation
# aerosol_process = coagulation_process.execute(
# aerosol=aerosol_process,
# time_step=time_step,
# sub_steps=coagulation_sub_step,
# )
# Record size distribution
distribution_counts[step, :], edges = np.histogram(
aerosol_process.particles.get_radius(clone=True),
bins=bins_lognormal,
density=False,
)
# Record species mass in particles
species_temp = aerosol_process.particles.get_species_mass(clone=True)
species_mass[step, :] = np.sum(species_temp, axis=0)
# Record vapor saturation ratios
vapor_saturation[step, :] = (
aerosol_process.atmosphere.partitioning_species.get_saturation_ratio(
temperature=temperature
)
)
# Normalize counts to number concentration (#/mÂł)
concentrations = distribution_counts / simulation_volume
# Final system state
print(aerosol_process)
Gas mixture at 298.15 K, 101325.0 Pa, partitioning=['Water' 'NaCl' 'potassium bicarbonate (KHCO3)' 'potassium dihydrogen phosphate (KH2PO4)' 'calcium bicarbonate [Ca(HCO3)2]' 'magnesium bicarbonate [Mg(HCO3)2]' 'Urea' 'Lactate' 'Glucose' 'Lysozyme' 'IgA' 'cholesterol'], gas_only_species=None Particle Representation: Strategy: ParticleResolvedSpeciatedMass Activity: ActivityKappaParameter Surface: SurfaceStrategyMass Mass Concentration: 1.592e-09 [kg/m^3] Number Concentration: 1.000e+02 [#/m^3]
Running Sim: 100%|ââââââââââ| 200/200 [01:55<00:00, 1.73it/s]
Gas mixture at 298.15 K, 101325.0 Pa, partitioning=['Water' 'NaCl' 'potassium bicarbonate (KHCO3)' 'potassium dihydrogen phosphate (KH2PO4)' 'calcium bicarbonate [Ca(HCO3)2]' 'magnesium bicarbonate [Mg(HCO3)2]' 'Urea' 'Lactate' 'Glucose' 'Lysozyme' 'IgA' 'cholesterol'], gas_only_species=None Particle Representation: Strategy: ParticleResolvedSpeciatedMass Activity: ActivityKappaParameter Surface: SurfaceStrategyMass Mass Concentration: 1.068e-11 [kg/m^3] Number Concentration: 1.000e+02 [#/m^3]
đ Step 7: Compare Initial and Final Size Distributions¶
This plot shows how the particle size distribution evolved during the simulation.
đ What It Shows:¶
- Dashed line: Initial size distribution before simulation
- Solid line: Final distribution after 120 seconds
- X-axis (log scale): Particle radius in meters
- Y-axis: Particle number concentration (#/mÂł)
đ§Ș What to Look For:¶
- If the final peak shifts left, particles are shrinking due to evaporation
- If the peak moves right or broadens, it may indicate condensation or growth
- In this case, we observe a leftward shift for larger particles, consistent with net evaporation
In [31]:
Copied!
# %% Step 7: Plot initial vs final size distributions
fig, ax = plt.subplots(figsize=(7, 5))
# Get radii before and after simulation
initial_distribution = aerosol.particles.get_radius(clone=True)
final_distribution = aerosol_process.particles.get_radius(clone=True)
# Compute binned counts
counts_initial, _ = np.histogram(
initial_distribution, bins=bins_lognormal, density=False
)
counts_final, _ = np.histogram(
final_distribution, bins=bins_lognormal, density=False
)
# Plot initial distribution
ax.plot(
bins_lognormal[:-1],
counts_initial / simulation_volume,
label="Initial Size Distribution",
color=TAILWIND["sky"]["800"],
linestyle="--",
)
# Plot final distribution
ax.plot(
bins_lognormal[:-1],
counts_final / simulation_volume,
label="Final Size Distribution",
color=TAILWIND["rose"]["400"],
linestyle="-",
)
# Set axis scales and labels
ax.set_xscale("log")
ax.set_xlabel("Particle Radius (m)")
ax.set_ylabel("Concentration (#/mâ»Âł)")
ax.legend()
plt.tight_layout()
plt.show()
# %% Step 7: Plot initial vs final size distributions
fig, ax = plt.subplots(figsize=(7, 5))
# Get radii before and after simulation
initial_distribution = aerosol.particles.get_radius(clone=True)
final_distribution = aerosol_process.particles.get_radius(clone=True)
# Compute binned counts
counts_initial, _ = np.histogram(
initial_distribution, bins=bins_lognormal, density=False
)
counts_final, _ = np.histogram(
final_distribution, bins=bins_lognormal, density=False
)
# Plot initial distribution
ax.plot(
bins_lognormal[:-1],
counts_initial / simulation_volume,
label="Initial Size Distribution",
color=TAILWIND["sky"]["800"],
linestyle="--",
)
# Plot final distribution
ax.plot(
bins_lognormal[:-1],
counts_final / simulation_volume,
label="Final Size Distribution",
color=TAILWIND["rose"]["400"],
linestyle="-",
)
# Set axis scales and labels
ax.set_xscale("log")
ax.set_xlabel("Particle Radius (m)")
ax.set_ylabel("Concentration (#/mâ»Âł)")
ax.legend()
plt.tight_layout()
plt.show()
đ§ Step 8: Track Water Vapor Saturation Over Time¶
To understand how the gas phase responds during the simulation, we track the vapor saturation ratio for water (RH) throughout the 120 seconds.
đ What Is the Saturation Ratio?¶
The saturation ratio is the ratio of actual vapor concentration to the saturation vapor concentration at a given temperature. For water:
$$ \text{Saturation Ratio} = \frac{C_{\text{vapor}}}{C_{\text{sat}}} $$
- 1.0 indicates equilibrium (100% RH)
- < 1.0 indicates subsaturation
- > 1.0 indicates supersaturation
đ What We Learn:¶
In this case, the saturation ratio remains nearly constant, meaning:
- The amount of water vapor in the gas phase is high enough
- The mass transferred into particles (via condensation or evaporation) is too small to significantly disturb the vapor reservoir
This is typical in low mass-loading regimes, such as indoor air with sparse respiratory aerosols.
In [32]:
Copied!
# %% Step 8: Plot water vapor saturation ratio over time
fig, ax = plt.subplots(figsize=(7, 5))
# Plot water vapor saturation ratio (species index 0)
ax.plot(
time_array,
vapor_saturation[:, 0],
color=TAILWIND["sky"]["800"],
label="Vapor Saturation Ratio (Water)",
)
ax.set_ylim(0, 1.1)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Vapor Saturation Ratio")
ax.set_title("Vapor Saturation Ratio for Water Over Time")
ax.legend()
plt.tight_layout()
plt.show()
# %% Step 8: Plot water vapor saturation ratio over time
fig, ax = plt.subplots(figsize=(7, 5))
# Plot water vapor saturation ratio (species index 0)
ax.plot(
time_array,
vapor_saturation[:, 0],
color=TAILWIND["sky"]["800"],
label="Vapor Saturation Ratio (Water)",
)
ax.set_ylim(0, 1.1)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Vapor Saturation Ratio")
ax.set_title("Vapor Saturation Ratio for Water Over Time")
ax.legend()
plt.tight_layout()
plt.show()
đ§Ș Step 9: Track Composition Evolution â Particle Mass Fractions¶
As water evaporates from particles, their composition changes over time. This plot shows the mean mass fraction of each species remaining in the particles during the simulation.
đ What This Plot Tells Us:¶
- Water (sky blue) dominates at the start
- As evaporation progresses, non-volatile species (salts, organics, proteins) become relatively more concentrated
- Supermicron particles (initially 4â40 ÎŒm) retain water much longer, showing water mass persists before stabilizing.
This plot illustrates size-dependent evaporation timescales: larger droplets lose water more slowly due to lower surface-area-to-volume ratios.
In [33]:
Copied!
# %% Step 9: Plot particle mass fraction evolution
# Compute total mass fraction over time
mass_fraction = species_mass / np.sum(species_mass, axis=1, keepdims=True)
# Set color palette (one color per species)
color_list = [
TAILWIND["sky"]["800"], # Water
TAILWIND["rose"]["400"], # NaCl
TAILWIND["red"]["800"], # KHCO3
TAILWIND["red"]["600"], # KH2PO4
TAILWIND["amber"]["500"], # Ca(HCO3)2
TAILWIND["orange"]["400"], # Mg(HCO3)2
TAILWIND["teal"]["900"], # Urea
TAILWIND["teal"]["700"], # Lactate
TAILWIND["green"]["500"], # Glucose
TAILWIND["indigo"]["600"], # Lysozyme
TAILWIND["indigo"]["400"], # IgA
TAILWIND["cyan"]["500"], # Cholesterol
]
# Create stackplot of species mass fraction
fig, ax = plt.subplots(figsize=(8, 5))
ax.stackplot(
time_array,
mass_fraction.T, # Each row = one species
labels=short_name, # Species names
colors=color_list,
alpha=1.0,
)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Mean Mass Fraction in Particles")
ax.set_ylim(0, 1)
ax.grid(True)
ax.legend(
title="Species",
loc="lower left",
framealpha=0.9,
fontsize=10,
)
plt.tight_layout()
plt.show()
# %% Step 9: Plot particle mass fraction evolution
# Compute total mass fraction over time
mass_fraction = species_mass / np.sum(species_mass, axis=1, keepdims=True)
# Set color palette (one color per species)
color_list = [
TAILWIND["sky"]["800"], # Water
TAILWIND["rose"]["400"], # NaCl
TAILWIND["red"]["800"], # KHCO3
TAILWIND["red"]["600"], # KH2PO4
TAILWIND["amber"]["500"], # Ca(HCO3)2
TAILWIND["orange"]["400"], # Mg(HCO3)2
TAILWIND["teal"]["900"], # Urea
TAILWIND["teal"]["700"], # Lactate
TAILWIND["green"]["500"], # Glucose
TAILWIND["indigo"]["600"], # Lysozyme
TAILWIND["indigo"]["400"], # IgA
TAILWIND["cyan"]["500"], # Cholesterol
]
# Create stackplot of species mass fraction
fig, ax = plt.subplots(figsize=(8, 5))
ax.stackplot(
time_array,
mass_fraction.T, # Each row = one species
labels=short_name, # Species names
colors=color_list,
alpha=1.0,
)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Mean Mass Fraction in Particles")
ax.set_ylim(0, 1)
ax.grid(True)
ax.legend(
title="Species",
loc="lower left",
framealpha=0.9,
fontsize=10,
)
plt.tight_layout()
plt.show()
đ Step 10: Size Distribution Evolution Over Time¶
This contour plot shows how the number concentration of aerosol particles changes across time and particle size. It's an excellent way to visualize both evaporation dynamics and distribution narrowing.
đ What to Look For:¶
- The Y-axis (log scale) shows particle radius from $10^{-7}$ to $10^{-4}$ meters
- The X-axis tracks simulation time from 0 to 120 seconds
- The color map shows $\log_{10}$ of number concentration (#/mÂł)
đ§ Interpretation:¶
- Supermicron particles (initially >1 ”m) evaporate rapidly in the first few seconds
- As water is lost, particles shrink and the size distribution stabilizes
- Smaller droplets reach their equilibrium sizes faster due to greater surface-to-volume ratios
This plot captures the transient drying behavior of cough droplets in a controlled environment.
In [34]:
Copied!
# %% Step 10: Contour plot of particle size distribution over time
fig, ax = plt.subplots(figsize=(7, 5))
# Create 2D mesh grid for time and particle radius
X, Y = np.meshgrid(time_array, edges[:-1])
# Convert number concentration to log10 scale (avoid log(0) issues)
log_conc = np.log10(
concentrations,
where=concentrations > 0,
out=np.full_like(concentrations, np.nan), # mask zeros
)
# Contour plot of log concentration
cont = ax.contourf(X, Y, log_conc.T)
# Axis settings
ax.set_yscale("log")
ax.set_ylim(1e-7, 1e-4)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Particle Radius (m)")
# Add colorbar
fig.colorbar(cont, label="Logââ Number Concentration")
plt.tight_layout()
plt.show()
# %% Step 10: Contour plot of particle size distribution over time
fig, ax = plt.subplots(figsize=(7, 5))
# Create 2D mesh grid for time and particle radius
X, Y = np.meshgrid(time_array, edges[:-1])
# Convert number concentration to log10 scale (avoid log(0) issues)
log_conc = np.log10(
concentrations,
where=concentrations > 0,
out=np.full_like(concentrations, np.nan), # mask zeros
)
# Contour plot of log concentration
cont = ax.contourf(X, Y, log_conc.T)
# Axis settings
ax.set_yscale("log")
ax.set_ylim(1e-7, 1e-4)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Particle Radius (m)")
# Add colorbar
fig.colorbar(cont, label="Logââ Number Concentration")
plt.tight_layout()
plt.show()
â Summary: Modeling Cough Droplet Evaporation with Particula¶
In this notebook, we used Particula to simulate the isothermal evaporation of cough droplets in a well-mixed air environment. We walked step-by-step through the setup of both gas-phase and particle-phase systems, created a realistic multicomponent chemical composition, and applied condensation dynamics to understand how droplets evolve over time.
The simulated particles ranged from submicron to supermicron sizes, following a lognormal distribution that reflects real-world respiratory emissions. As expected, supermicron particles took significantly longer for their water mass fraction to decrease, highlighting the importance of size when modeling indoor exposure or transport of pathogens.
We tracked changes in size distribution, vapor saturation ratio, and particle composition over time, and visualized how water loss drives dynamic shifts in mass fraction and droplet behavior.
đ§ Take-Home Points¶
Multicomponent aerosols are easy to build with Particula You can represent realistic droplet mixtures using species-specific mass fractions, thermodynamic properties, and hygroscopicity (Îș values).
Water loss dominates early evaporation dynamics In low humidity environments, water rapidly evaporates from smaller particles, while larger ones persist due to slower mass transfer rates.
Gasâparticle coupling is size- and concentration-dependent In this example, the gas-phase RH stayed nearly constant, showing that low particle loading may have minimal impact on ambient vapor conditions.
Particula supports modular, testable dynamics Processes like condensation, coagulation, and even fixed-RH simulations can be toggled and combined, offering flexibility for exploring a wide range of scenarios.
This notebook provides a foundational example for respiratory aerosol simulations. You can now extend it to test:
- Different relative humidity levels
- Effects of ventilation or air exchange
- Indoor vs. outdoor aerosol transport scenarios
- Inclusion of more complex chemistry or temperature effects