homogeneous nucleation of liquid droplets & expansion chamber example

.gitignore

@@ -4,6 +4,9 @@ PySDM_examples/utils/temporary_files/*
 *.pkl
 *.vtu
 *.vts
+*.png
+*.csv
+*.xlsx
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

PySDM/environments/expansion_chamber.py

@@ -0,0 +1,138 @@
+"""
+Zero-dimensional expansion chamber framework
+"""
+
+import numpy as np
+
+from PySDM.environments.impl.moist import Moist
+from PySDM.impl.mesh import Mesh
+from PySDM.initialisation.equilibrate_wet_radii import (
+    default_rtol,
+    equilibrate_wet_radii,
+)
+from PySDM.environments.impl import register_environment
+
+
+@register_environment()
+class ExpansionChamber(Moist):
+    def __init__(
+        self,
+        *,
+        dt,
+        volume: float,
+        p0: float,
+        T0: float,
+        initial_water_vapour_mixing_ratio: float,
+        pf: float,
+        tf: float,
+        variables,
+        mixed_phase=False,
+    ):
+        super().__init__(dt, Mesh.mesh_0d(), variables, mixed_phase=mixed_phase)
+
+        self.volume = volume
+        self.p0 = p0
+        self.T0 = T0
+        self.initial_water_vapour_mixing_ratio = initial_water_vapour_mixing_ratio
+        self.pf = pf
+        self.tf = tf
+
+        self.formulae = None
+        self.delta_liquid_water_mixing_ratio = np.nan
+        self.params = None
+
+    def register(self, builder):
+        self.formulae = builder.particulator.formulae
+        pd0 = self.formulae.trivia.p_d(self.p0, self.initial_water_vapour_mixing_ratio)
+        rhod0 = self.formulae.state_variable_triplet.rhod_of_pd_T(pd0, self.T0)
+        self.mesh.dv = self.volume
+
+        Moist.register(self, builder)
+
+        params = (self.p0, self.T0, self.initial_water_vapour_mixing_ratio, 0)
+        self["p"][:] = params[0]
+        self["T"][:] = params[1]
+        self["water_vapour_mixing_ratio"][:] = params[2]
+        self["t"][:] = params[3]
+
+        self._tmp["water_vapour_mixing_ratio"][:] = params[0]
+        self.sync_chamber_vars()
+        Moist.sync(self)
+        self.notify()
+
+    def init_attributes(
+        self,
+        *,
+        n_in_dv: [float, np.ndarray],
+        kappa: float,
+        r_dry: [float, np.ndarray],
+        rtol=default_rtol,
+        include_dry_volume_in_attribute: bool = True,
+    ):
+        if not isinstance(n_in_dv, np.ndarray):
+            r_dry = np.array([r_dry])
+            n_in_dv = np.array([n_in_dv])
+
+        attributes = {}
+        dry_volume = self.formulae.trivia.volume(radius=r_dry)
+        attributes["kappa times dry volume"] = dry_volume * kappa
+        attributes["multiplicity"] = n_in_dv
+        r_wet = equilibrate_wet_radii(
+            r_dry=r_dry,
+            environment=self,
+            kappa_times_dry_volume=attributes["kappa times dry volume"],
+            rtol=rtol,
+        )
+        attributes["volume"] = self.formulae.trivia.volume(radius=r_wet)
+        if include_dry_volume_in_attribute:
+            attributes["dry volume"] = dry_volume
+        return attributes
+
+    def formulae_drho_dp(const, p, T, water_vapour_mixing_ratio, lv, dql_dp):
+        R_q = const.Rv / (1 / water_vapour_mixing_ratio + 1) + const.Rd / (
+            1 + water_vapour_mixing_ratio
+        )
+        cp_q = const.c_pv / (1 / water_vapour_mixing_ratio + 1) + const.c_pd / (
+            1 + water_vapour_mixing_ratio
+        )
+        return -1 / R_q / T - p / R_q / T**2 * lv / cp_q * dql_dp
+
+    def advance_chamber_vars(self):
+        dt = self.particulator.dt
+        T = self["T"][0]
+        p = self["p"][0]
+        t = self["t"][0]
+
+        dp_dt = (self.pf - self.p0) / (self.tf) * (self.dt / 2)  # mid-point
+        water_vapour_mixing_ratio = (
+            self["water_vapour_mixing_ratio"][0]
+            - self.delta_liquid_water_mixing_ratio / 2
+        )
+
+        drho_dp = self.formulae_drho_dp(
+            p=p,
+            T=T,
+            water_vapour_mixing_ratio=water_vapour_mixing_ratio,
+            lv=self.formulae.latent_heat.lv(T),
+            dql_dp=self.delta_liquid_water_mixing_ratio / dt,
+        )
+        drhod_dp = drho_dp
+
+        self.particulator.backend.explicit_euler(self._tmp["t"], dt, 1)
+        self.particulator.backend.explicit_euler(self._tmp["p"], dt, dp_dt)
+        self.particulator.backend.explicit_euler(
+            self._tmp["rhod"], dt, drhod_dp * dp_dt
+        )
+
+    def sync_chamber_vars(self):
+        self.delta_liquid_water_mixing_ratio = (
+            self._tmp["water_vapour_mixing_ratio"][0]
+            - self["water_vapour_mixing_ratio"][0]
+        )
+        for var in self.variables:
+            self._tmp[var][:] = self[var][:]
+
+    def sync(self):
+        self.sync_chamber_vars()
+        self.advance_chamber_vars()
+        super().sync()

PySDM/formulae.py

@@ -44,6 +44,7 @@ def __init__(  # pylint: disable=too-many-locals
         hydrostatics: str = "Default",
         freezing_temperature_spectrum: str = "Null",
         heterogeneous_ice_nucleation_rate: str = "Null",
+        homogeneous_liquid_nucleation_rate: str = "Null",
         fragmentation_function: str = "AlwaysN",
         isotope_equilibrium_fractionation_factors: str = "Null",
         isotope_meteoric_water_line_excess: str = "Null",
@@ -76,6 +77,7 @@ def __init__(  # pylint: disable=too-many-locals
         self.hydrostatics = hydrostatics
         self.freezing_temperature_spectrum = freezing_temperature_spectrum
         self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate
+        self.homogeneous_liquid_nucleation_rate = homogeneous_liquid_nucleation_rate
         self.fragmentation_function = fragmentation_function
         self.isotope_equilibrium_fractionation_factors = (
             isotope_equilibrium_fractionation_factors

PySDM/physics/__init__.py

@@ -25,6 +25,7 @@
     fragmentation_function,
     freezing_temperature_spectrum,
     heterogeneous_ice_nucleation_rate,
+    homogeneous_liquid_nucleation_rate,
     hydrostatics,
     hygroscopicity,
     impl,

PySDM/physics/constants_defaults.py

@@ -56,6 +56,7 @@
 
 R_str = sci.R * si.joule / si.kelvin / si.mole
 N_A = sci.N_A / si.mole
+k_B = sci.k * si.joule / si.kelvin
 
 D0 = 2.26e-5 * si.metre**2 / si.second
 D_exp = 1.81

PySDM/physics/homogeneous_liquid_nucleation_rate/__init__.py

@@ -0,0 +1,7 @@
+"""
+homogeneous liquid droplet nucleation rate formulations
+"""
+
+from PySDM.impl.null_physics_class import Null
+from .cnt import CNT
+from .constant import Constant

PySDM/physics/homogeneous_liquid_nucleation_rate/cnt.py

@@ -0,0 +1,32 @@
+"""
+classical nucleation theory (CNT)
+formulae from [Seinfeld and Pandis](https://archive.org/details/0237-pdf-atmospheric-chemistry-and-physics-2nd-ed-j.-seinfeld-s.-pandis-wiley-2006-ww)
+Eq (11.47) and (11.52)
+"""  # pylint: disable=line-too-long
+
+import numpy as np
+
+
+class CNT:
+    def __init__(self, _):
+        return
+
+    @staticmethod
+    def j_liq_homo(const, T, S, e_s):
+        mass_per_moleculue = const.Mv / const.N_A
+        volume_per_molecule = mass_per_moleculue / const.rho_w
+        N1 = (S * e_s) / (mass_per_moleculue * const.Rv * T)  # molecules per m3
+        return (
+            ((2 * const.sgm_w) / (np.pi * mass_per_moleculue)) ** (1 / 2)
+            * (volume_per_molecule * N1**2 / S)
+            * np.exp(
+                (-16 * np.pi * volume_per_molecule**2 * const.sgm_w**3)
+                / (3 * const.k_B**3 * T**3 * np.log(S) ** 2)
+            )
+        )
+
+    @staticmethod
+    def r_liq_homo(const, T, S):
+        mass_per_moleculue = const.Mv / const.N_A
+        volume_per_molecule = mass_per_moleculue / const.rho_w
+        return (2 * const.sgm_w * volume_per_molecule) / (const.k_B * T * np.log(S))

PySDM/physics/homogeneous_liquid_nucleation_rate/constant.py

@@ -0,0 +1,19 @@
+"""
+constant rate formulation (for tests)
+"""
+
+import numpy as np
+
+
+class Constant:
+    def __init__(self, const):
+        assert np.isfinite(const.J_LIQ_HOMO)
+        assert np.isfinite(const.R_LIQ_HOMO)
+
+    @staticmethod
+    def j_liq_homo(const, T, S, e_s):  # pylint: disable=unused-argument
+        return const.J_LIQ_HOMO
+
+    @staticmethod
+    def r_liq_homo(const, T, S):  # pylint: disable=unused-argument
+        return const.R_LIQ_HOMO

docs/bibliography.json

@@ -691,5 +691,15 @@
     "usages": ["PySDM/physics/saturation_vapour_pressure/wexler_1976.py"],
     "title": "Vapor Pressure Formulation for Water in Range 0 to f 00 °C. A Revision",
     "label": "Wexler 1976 (J. Res. NBS A Phys. Ch. 80A)"
+  },
+  "https://doi.org/10.48550/arXiv.2501.01467": {
+    "usages": ["examples/PySDM_examples/expansion_chamber/aerosol.py", "examples/PySDM_examples/expansion_chamber/expansion_experiment.ipynb"],
+    "title": "Droplet Nucleation In a Rapid Expansion Aerosol Chamber",
+    "label": "Erinin et al. 2024 (arXiv:2501.01467)"
+  },
+  "https://archive.org/details/0237-pdf-atmospheric-chemistry-and-physics-2nd-ed-j.-seinfeld-s.-pandis-wiley-2006-ww": {
+    "usages": ["PySDM/physics/homogeneous_liquid_nucleation_rate/cnt.py"],
+    "title": "Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 2nd Edition",
+    "label": "Seinfeld & Pandis 2006 (Wiley)"
   }
 }

docs/generate_html.py

@@ -67,7 +67,7 @@ def check_urls(urls_from_json):
                 r"\b(https://doi\.org/10[.][0-9]{4,}(?:[.][0-9]+)*/(?:(?![\"&\'<>^\\])\S)+)\b",
                 r"\b(https://digitallibrary\.un\.org/record/(?:[0-9])+)\b",
                 r"\b(http://mi\.mathnet\.ru/dan(?:[0-9])+)\b",
-                r"\b(https://archive.org/details/(?:[0-9a-z])+)\b",
+                r"\b(https://archive.org/details/(?:[0-9a-z_\.-])+)\b",
             ):
                 urls = re.findall(pattern, text)
                 if urls:

examples/PySDM_examples/expansion_chamber/__init__.py

@@ -0,0 +1,4 @@
+# pylint: disable=invalid-name
+"""
+expansion chamber example imagined as an ascending parcel with updraft velocity matched to dp/dt
+"""

examples/PySDM_examples/expansion_chamber/aerosol.py

@@ -0,0 +1,47 @@
+""" aerosol parameters from [Erinin et al. 2025](https://doi.org/10.48550/arXiv.2501.01467) """
+
+from chempy import Substance
+from pystrict import strict
+
+from PySDM.initialisation import spectra
+from PySDM.initialisation.aerosol_composition import DryAerosolMixture
+from PySDM.physics import si
+
+
+@strict
+class AerosolChamber(DryAerosolMixture):
+    def __init__(
+        self,
+        water_molar_volume: float,
+        N: float = 1e5 / si.cm**3,
+    ):
+        mode1 = {"CaCO3": 1.0}
+
+        super().__init__(
+            compounds=("CaCO3",),
+            molar_masses={
+                "CaCO3": Substance.from_formula("CaCO3").mass * si.g / si.mole
+            },
+            densities={
+                "CaCO3": 2.71 * si.g / si.cm**3,
+            },
+            is_soluble={
+                "CaCO3": True,
+            },
+            ionic_dissociation_phi={
+                "CaCO3": 1,
+            },
+        )
+        self.modes = (
+            {
+                "kappa": self.kappa(
+                    mass_fractions=mode1,
+                    water_molar_volume=water_molar_volume,
+                ),
+                "spectrum": spectra.Gaussian(
+                    norm_factor=N,
+                    loc=316 / 2 * si.nm,
+                    scale=257 / 2 * si.nm,
+                ),
+            },
+        )