homogeneous nucleation of droplets & expansion chamber examples (Clare's branch)

.gitignore

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

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 .cnt import CNT
+from .constant import Constant
+from .null import Null

PySDM/physics/homogeneous_liquid_nucleation_rate/cnt.py

@@ -0,0 +1,31 @@
+"""
+classical nucleation theory (CNT)
+formulae from Seinfeld and Pandis Eq (11.47) and (11.52)
+"""
+
+import numpy as np
+
+
+class CNT:
+    def __init__(self, _):
+        return
+
+    @staticmethod
+    def j_liq_homo(const, T, S, e_s):
+        m1 = const.Mv / const.N_A  # kg per molecule
+        v1 = m1 / const.rho_w  # m3 per molecule
+        N1 = (S * e_s) / (m1 * const.Rv * T)  # molecules per m3
+        return (
+            ((2 * const.sgm_w) / (np.pi * m1)) ** (1 / 2)
+            * (v1 * N1**2 / S)
+            * np.exp(
+                (-16 * np.pi * v1**2 * const.sgm_w**3)
+                / (3 * const.k_B**3 * T**3 * np.log(S) ** 2)
+            )
+        )
+
+    @staticmethod
+    def r_liq_homo(const, T, S):
+        m1 = const.Mv / const.N_A  # kg per molecule
+        v1 = m1 / const.rho_w  # m3 per molecule
+        return (2 * const.sgm_w * v1) / (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

PySDM/physics/homogeneous_liquid_nucleation_rate/null.py

@@ -0,0 +1,17 @@
+"""
+do-nothing null formulation (needed as other formulations require parameters
+ to be set before instantiation of Formulae)
+"""
+
+
+class Null:  # pylint: disable=unused-argument
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def j_liq_homo(const, T, S, e_s):
+        return 0
+
+    @staticmethod
+    def r_liq_homo(const, T, S):
+        return 0

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"],
+    "title": "Droplet Nucleation In a Rapid Expansion Aerosol Chamber",
+    "label": "Erinin et al. 2024 (arXiv:2501.01467)"
+  },
+  "https://books.google.com/books/about/Atmospheric_Chemistry_and_Physics.html?id=n_RmCgAAQBAJ": {
+    "usages": ["PySDM/physics/homogeneous_liquid_nucleation_rate/cnt.py"],
+    "title": "Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 3rd Edition",
+    "label": "Seinfeld and Pandis, Ed.3"
   }
 }

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,48 @@
+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.Lognormal(
+                    norm_factor=N,
+                    m_mode=158 * si.nm,
+                    s_geom=2,
+                ),
+            },
+        )
+        # mean diameter 316nm, standard deviation 257nm
+        # not sure how to interpret the standard deviation given in the paper
+        # because it looks like it's a lognormal distribution in Fig 2

examples/PySDM_examples/expansion_chamber/expansion_simulation.py

@@ -0,0 +1,161 @@
+from collections import namedtuple
+
+import numpy as np
+from PySDM import Builder
+from PySDM import products as PySDM_products
+from PySDM.backends import CPU
+from PySDM.dynamics import AmbientThermodynamics, Condensation
+from PySDM.environments import Parcel
+from PySDM.initialisation import equilibrate_wet_radii
+from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
+from PySDM.physics import si
+
+
+def run_expansion(
+    formulae,
+    aerosol,
+    n_sd_per_mode,
+    RH0=0.7,
+    T0=296 * si.K,
+    p0=1000 * si.hPa,
+    pf=500 * si.hPa,
+    dz=10 * si.m,
+    t_lift=2 * si.s,
+    t_max=4 * si.s,
+):
+
+    # calculate w, dt, n_steps
+    H = 8500 * si.m  # atmospheric scale height
+    z_lift = -H * np.log(pf / p0)
+    w_lift = z_lift / t_lift
+    w = lambda t: w_lift if t < t_lift else 1e-9
+    dt = dz / w_lift
+    n_steps = int(np.ceil(t_max / dt))
+    # print(f"z_lift={z_lift}")
+    # print(f"w_lift={w_lift}")
+    # print(f"dz={dz}")
+    # print(f"dt={dt}")
+    # print(f"n_steps={n_steps}")
+
+    dry_radius_bin_edges = np.geomspace(50 * si.nm, 2000 * si.nm, 40, endpoint=False)
+    wet_radius_bin_edges = np.geomspace(1 * si.um, 40 * si.um, 40, endpoint=False)
+    products = (
+        PySDM_products.WaterMixingRatio(unit="g/kg", name="liquid_water_mixing_ratio"),
+        PySDM_products.PeakSupersaturation(name="s"),
+        PySDM_products.AmbientRelativeHumidity(name="RH"),
+        PySDM_products.AmbientTemperature(name="T"),
+        PySDM_products.AmbientPressure(name="p"),
+        PySDM_products.AmbientWaterVapourMixingRatio(
+            unit="g/kg", name="water_vapour_mixing_ratio"
+        ),
+        PySDM_products.ParcelDisplacement(name="z"),
+        PySDM_products.Time(name="t"),
+        PySDM_products.ParticleSizeSpectrumPerVolume(
+            name="dry:dN/dR",
+            unit="m^-3 m^-1",
+            radius_bins_edges=dry_radius_bin_edges,
+            dry=True,
+        ),
+        PySDM_products.ParticleSizeSpectrumPerVolume(
+            name="wet:dN/dR",
+            unit="m^-3 m^-1",
+            radius_bins_edges=wet_radius_bin_edges,
+            dry=False,
+        ),
+        PySDM_products.ActivatedEffectiveRadius(
+            name="reff", unit="um", count_activated=True, count_unactivated=False
+        ),
+    )
+
+    const = formulae.constants
+    pv0 = RH0 * formulae.saturation_vapour_pressure.pvs_water(T0)
+
+    env = Parcel(
+        dt=dt,
+        mass_of_dry_air=1e3 * si.kg,
+        p0=p0,
+        initial_water_vapour_mixing_ratio=const.eps * pv0 / (p0 - pv0),
+        w=w,
+        T0=T0,
+    )
+
+    n_sd = n_sd_per_mode * len(aerosol.modes)
+
+    builder = Builder(backend=CPU(formulae), n_sd=n_sd, environment=env)
+    builder.add_dynamic(AmbientThermodynamics())
+    builder.add_dynamic(Condensation())
+    builder.request_attribute("critical supersaturation")
+
+    attributes = {
+        k: np.empty(0) for k in ("dry volume", "kappa times dry volume", "multiplicity")
+    }
+    for i, mode in enumerate(aerosol.modes):
+        kappa, spectrum = mode["kappa"]["Constant"], mode["spectrum"]
+        r_dry, concentration = ConstantMultiplicity(spectrum).sample(n_sd_per_mode)
+        v_dry = builder.formulae.trivia.volume(radius=r_dry)
+        specific_concentration = concentration / builder.formulae.constants.rho_STP
+        attributes["multiplicity"] = np.append(
+            attributes["multiplicity"],
+            specific_concentration * builder.particulator.environment.mass_of_dry_air,
+        )
+        attributes["dry volume"] = np.append(attributes["dry volume"], v_dry)
+        attributes["kappa times dry volume"] = np.append(
+            attributes["kappa times dry volume"], v_dry * kappa
+        )
+
+    r_wet = equilibrate_wet_radii(
+        r_dry=builder.formulae.trivia.radius(volume=attributes["dry volume"]),
+        environment=builder.particulator.environment,
+        kappa_times_dry_volume=attributes["kappa times dry volume"],
+    )
+    attributes["volume"] = builder.formulae.trivia.volume(radius=r_wet)
+
+    particulator = builder.build(attributes, products=products)
+
+    output = {product.name: [] for product in particulator.products.values()}
+    output_attributes = {
+        "multiplicity": tuple([] for _ in range(particulator.n_sd)),
+        "volume": tuple([] for _ in range(particulator.n_sd)),
+        "critical volume": tuple([] for _ in range(particulator.n_sd)),
+        "critical supersaturation": tuple([] for _ in range(particulator.n_sd)),
+    }
+
+    for _ in range(n_steps):
+        particulator.run(steps=1)
+        for product in particulator.products.values():
+            if product.name == "dry:dN/dR" or product.name == "wet:dN/dR":
+                continue
+            value = product.get()
+            if product.name == "t":
+                output[product.name].append(value)
+            else:
+                output[product.name].append(value[0])
+        for key, attr in output_attributes.items():
+            attr_data = particulator.attributes[key].to_ndarray()
+            for drop_id in range(particulator.n_sd):
+                attr[drop_id].append(attr_data[drop_id])
+
+    dry_spectrum = particulator.products["dry:dN/dR"].get()
+    wet_spectrum = particulator.products["wet:dN/dR"].get()
+
+    Output = namedtuple(
+        "Output",
+        [
+            "profile",
+            "attributes",
+            "aerosol",
+            "dry_radius_bin_edges",
+            "dry_spectrum",
+            "wet_radius_bin_edges",
+            "wet_spectrum",
+        ],
+    )
+    return Output(
+        profile=output,
+        attributes=output_attributes,
+        aerosol=aerosol,
+        dry_radius_bin_edges=dry_radius_bin_edges,
+        dry_spectrum=dry_spectrum,
+        wet_radius_bin_edges=wet_radius_bin_edges,
+        wet_spectrum=wet_spectrum,
+    )