Isotopic relaxation timescales comparison

PySDM/physics/isotope_diffusivity_ratios/grahams_law.py

@@ -12,3 +12,9 @@ def ratio_2H(const, temperature):  # pylint: disable=unused-argument
         return (
             (2 * const.M_1H + const.M_16O) / (const.M_2H + const.M_1H + const.M_16O)
         ) ** const.ONE_HALF
+
+    @staticmethod
+    def ratio_3H(const, temperature):  # pylint: disable=unused-argument
+        return (
+            (2 * const.M_1H + const.M_16O) / (const.M_3H + const.M_1H + const.M_16O)
+        ) ** const.ONE_HALF

PySDM/physics/isotope_relaxation_timescale/bolin_1958.py

@@ -10,8 +10,42 @@ class Bolin1958:  # pylint: disable=too-few-public-methods
     def __init__(self, const):
         assert np.isfinite(const.BOLIN_ISOTOPE_TIMESCALE_COEFF_C1)
 
+    @staticmethod
+    # pylint: disable=too-many-arguments unused-argument
+    def tau_of_rdrdt(const, radius, r_dr_dt, alpha=0):
+        """timescale for evaporation of a falling drop with tritium"""
+        return -(radius**2) / 3 / r_dr_dt * const.BOLIN_ISOTOPE_TIMESCALE_COEFF_C1
+
+    @staticmethod
+    # pylint: disable=too-many-arguments
+    def c1_coeff(
+        const,
+        vent_coeff_iso,
+        vent_coeff,
+        D_iso,
+        D,
+        alpha,
+        rho_env_iso,
+        rho_env,
+        M_iso,
+        pvs_iso,
+        pvs_water,
+        temperature,
+    ):
+        return (
+            vent_coeff_iso
+            * D_iso
+            / vent_coeff
+            / D
+            / alpha
+            / pvs_iso
+            * pvs_water
+            * (rho_env_iso / M_iso - pvs_iso / const.R_str / temperature)
+            / (rho_env / const.Mv - pvs_water / const.R_str / temperature)
+        )
+
     @staticmethod
     # pylint: disable=too-many-arguments
-    def tau(const, radius, r_dr_dt):
+    def tau_of_rdrdt_c1(radius, r_dr_dt, c1_coeff):
         """timescale for evaporation of a falling drop with tritium"""
-        return (-3 / radius**2 * r_dr_dt * const.BOLIN_ISOTOPE_TIMESCALE_COEFF_C1) ** -1
+        return -(radius**2) / 3 / r_dr_dt / c1_coeff

PySDM/physics/isotope_relaxation_timescale/jouzel_et_al_1975.py

@@ -0,0 +1,19 @@
+"""eq. 7 in [Jouzel et al 1975](https://doi.org/10.1029/JC080i036p05015)"""
+
+
+class JouzelEtAl1975:
+    @staticmethod
+    # pylint: disable=too-many-arguments
+    def tau(const, e_s, D_iso, M_iso, vent_coeff_iso, radius, alpha, temperature):
+        """relaxation time for stationary droplet??"""
+        return (
+            radius**2
+            * alpha
+            * const.rho_w
+            * temperature
+            / 3
+            / vent_coeff_iso
+            / D_iso
+            / e_s
+            / M_iso
+        )

PySDM/physics/isotope_relaxation_timescale/miyake_et_al_1968.py

@@ -11,3 +11,8 @@ def tau(const, e_s, D, M, vent_coeff, radius, alpha, temperature):
         return (radius**2 * alpha * const.rho_w * const.R_str * temperature) / (
             3 * e_s * D * M * vent_coeff
         )
+
+    @staticmethod
+    # pylint: disable=too-many-arguments unused-argument
+    def tau_of_rdrdt(const, radius, r_dr_dt, alpha):
+        return -(radius**2) / 3 / r_dr_dt * alpha

PySDM/physics/isotope_relaxation_timescale/stewart_1975.py

@@ -0,0 +1,19 @@
+"""eq. 7 in [Stewart 1975](https://doi.org/10.1029/JC080i009p01133)"""
+
+
+class Stewart1975:
+    @staticmethod
+    # pylint: disable=too-many-arguments
+    def tau(const, e_s, Dn, M_iso, vent_coeff, radius, alpha, temperature):
+        """relaxation time for stationary droplet; Dn denotes D^n"""
+        return (
+            radius**2
+            * alpha
+            * const.rho_w
+            * temperature
+            / 3
+            / vent_coeff
+            / Dn
+            / e_s
+            / M_iso
+        )

PySDM/physics/isotope_relaxation_timescale/zaba_and_arabas_2025.py

@@ -0,0 +1,28 @@
+"""isotope relaxation timescale"""
+
+
+class ZabaAndArabas2025:
+    @staticmethod
+    def tau(
+        const,
+        radius,
+        alpha,
+        D_iso,
+        vent_coeff_iso,
+        rho_env,
+        M_iso,
+        temperature,
+        Rv,
+        pvs_iso,
+    ):
+        return (
+            alpha
+            * radius**2
+            / 3
+            * const.rho_w
+            / D_iso
+            / vent_coeff_iso
+            / const.Mv
+            * Rv
+            / (rho_env / M_iso - pvs_iso / const.R_STR / temperature)
+        )

examples/PySDM_examples/Bolin_1958/common.py

@@ -0,0 +1,71 @@
+from PySDM import Formulae
+
+
+class IsotopeTimescale:
+    def __init__(self, *, settings, temperature, radii):
+        self.radii = radii
+        self.formulae = Formulae(**settings)
+        self.temperature = temperature
+        self.pressure = self.formulae.constants.p_STP
+        self.v_term = self.formulae.terminal_velocity.v_term(radii)
+        self.D = self.formulae.diffusion_thermics.D(T=self.temperature, p=self.pressure)
+        self.D_iso = self.formulae.isotope_diffusivity_ratios.ratio_3H(self.temperature)
+        self.K = 44.0  # any non-zero value
+        self.pvs_water = self.formulae.saturation_vapour_pressure.pvs_water(
+            self.temperature
+        )
+        self.alpha = self.formulae.isotope_equilibrium_fractionation_factors.alpha_i_3H(
+            self.temperature
+        )  # check i/l
+        self.M_iso = self.formulae.constants.M_3H
+
+    def vent_coeff_fun(self):
+        eta_air = self.formulae.air_dynamic_viscosity.eta_air(self.temperature)
+        air_density = self.pressure / self.formulae.constants.Rd / self.temperature
+
+        assert abs(air_density - 1) / air_density < 0.3
+
+        Re = self.formulae.particle_shape_and_density.reynolds_number(
+            radius=self.radii,
+            velocity_wrt_air=self.v_term,
+            dynamic_viscosity=eta_air,
+            density=air_density,
+        )
+        Sc = self.formulae.trivia.air_schmidt_number(
+            dynamic_viscosity=eta_air,
+            diffusivity=self.D,
+            density=air_density,
+        )
+
+        return self.formulae.ventilation.ventilation_coefficient(
+            sqrt_re_times_cbrt_sc=self.formulae.trivia.sqrt_re_times_cbrt_sc(
+                Re=Re, Sc=Sc
+            )
+        )
+
+    def r_dr_dt(self, RH, RH_eq, lv):
+        return self.formulae.drop_growth.r_dr_dt(
+            T=self.temperature,
+            pvs=self.pvs_water,
+            D=self.D,
+            K=self.K,
+            ventilation_factor=self.vent_coeff_fun(),
+            RH=RH,
+            RH_eq=RH_eq,
+            lv=lv,
+        )
+
+    def c1_coeff(self, *, vent_coeff_iso, rho_env, pvs_iso):
+        return self.formulae.isotope_relaxation_timescale.c1_coeff(
+            vent_coeff_iso=vent_coeff_iso,
+            vent_coeff=self.vent_coeff_fun(),
+            D_iso=self.D_iso,
+            D=self.D,
+            alpha=self.alpha,
+            rho_env_iso=self.formulae.constants.VSMOW_R_3H,
+            rho_env=rho_env,
+            M_iso=self.M_iso,
+            pvs_iso=pvs_iso,  # any number
+            pvs_water=self.pvs_water,
+            temperature=self.temperature,
+        )