Isotopic relaxation timescales comparison

PySDM/physics/isotope_relaxation_timescale/bolin_1958.py

@@ -10,8 +10,60 @@ 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 tau(const, radius, r_dr_dt):
+    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_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
+
+    @staticmethod
+    # pylint: disable=too-many-arguments
+    def tau_without_assumptions(
+        const, radius, alpha, D, vent_coeff, rho_env, M, temperature, pvs_water, pvs
+    ):
+        return (
+            alpha
+            * radius**2
+            / 3
+            * const.rho_w
+            / D
+            / vent_coeff
+            / const.Mv
+            * pvs
+            / pvs_water
+            / (rho_env / M - pvs / const.R_STR / temperature)
+        )

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

examples/PySDM_examples/Bolin_1958/common.py

@@ -0,0 +1,46 @@
+from functools import partial
+from PySDM import Formulae
+
+
+class IsotopeTimescaleCommon:
+    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)
+
+    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_fun(self, K):
+        return partial(
+            self.formulae.drop_growth.r_dr_dt,
+            T=self.temperature,
+            pvs=self.formulae.saturation_vapour_pressure.pvs_water(self.temperature),
+            D=self.D,
+            K=K,
+            ventilation_factor=self.vent_coeff_fun(),
+        )

examples/PySDM_examples/Bolin_1958/table_1.ipynb

@@ -21,8 +21,8 @@
    "id": "5e4e58050cc01f64",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2024-12-05T21:30:20.214060Z",
-     "start_time": "2024-12-05T21:30:20.209298Z"
+     "end_time": "2025-02-07T14:51:42.815887Z",
+     "start_time": "2025-02-07T14:51:42.812266Z"
     }
    },
    "source": [
@@ -40,113 +40,84 @@
    "id": "95f360dc62f373f9",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2024-12-05T21:30:21.622672Z",
-     "start_time": "2024-12-05T21:30:20.218383Z"
+     "end_time": "2025-02-07T14:51:44.523180Z",
+     "start_time": "2025-02-07T14:51:42.819689Z"
     }
    },
    "source": [
-    "import numpy as np\n",
     "import pandas\n",
+    "import numpy as np\n",
     "from PySDM.physics import in_unit, si\n",
-    "from PySDM import Formulae"
+    "from PySDM import Formulae\n",
+    "from PySDM_examples.Bolin_1958.common import IsotopeTimescaleCommon"
    ],
    "outputs": [],
    "execution_count": 2
   },
   {
-   "cell_type": "code",
-   "id": "e29a92d6680c4e51",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2024-12-05T21:30:21.689116Z",
-     "start_time": "2024-12-05T21:30:21.687254Z"
+     "end_time": "2025-02-07T14:51:44.630655Z",
+     "start_time": "2025-02-07T14:51:44.613229Z"
     }
    },
-   "source": "radii = np.asarray([0.005, 0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.20]) * si.cm",
+   "cell_type": "code",
+   "source": [
+    "any_non_zero_value = 44.0\n",
+    "radii = np.asarray([0.005, 0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.20]) * si.cm\n",
+    "formulae = Formulae(\n",
+    "    terminal_velocity=\"RogersYau\",\n",
+    "    drop_growth=\"Mason1951\",\n",
+    "    diffusion_thermics=\"Neglect\",\n",
+    "    saturation_vapour_pressure=\"AugustRocheMagnus\",\n",
+    "    ventilation=\"Froessling1938\",\n",
+    "    particle_shape_and_density=\"LiquidSpheres\",\n",
+    "    air_dynamic_viscosity=\"ZografosEtAl1987\",\n",
+    "    constants={\"BOLIN_ISOTOPE_TIMESCALE_COEFF_C1\": 1.63},\n",
+    "    isotope_relaxation_timescale=\"Bolin1958\",\n",
+    ")\n",
+    "temperature = formulae.constants.T0 + 10 * si.K\n",
+    "iso_fun = IsotopeTimescaleCommon(settings=settings, temperature=temperature, radii=radii)\n",
+    "r_dr_dt = iso_fun.r_dr_dt_fun(K=any_non_zero_value)\n",
+    "adjustment_time = formulae.isotope_relaxation_timescale.tau_of_rdrdt(radius = radii, r_dr_dt = r_dr_dt(RH=0, RH_eq=1, lv=0))"
+   ],
+   "id": "4a8c2bd612c892f2",
    "outputs": [],
    "execution_count": 3
   },
   {
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2024-12-05T21:30:24.420354Z",
-     "start_time": "2024-12-05T21:30:21.698737Z"
+     "end_time": "2025-02-07T14:51:44.647096Z",
+     "start_time": "2025-02-07T14:51:44.636800Z"
     }
    },
    "cell_type": "code",
    "source": [
-    "formulae = Formulae(\n",
-    "    terminal_velocity='RogersYau',\n",
-    "    drop_growth='Mason1951',\n",
-    "    diffusion_thermics='Neglect',\n",
-    "    saturation_vapour_pressure='AugustRocheMagnus',\n",
-    "    ventilation='Froessling1938',\n",
-    "    particle_shape_and_density='LiquidSpheres',\n",
-    "    air_dynamic_viscosity='ZografosEtAl1987',\n",
-    "    constants={'BOLIN_ISOTOPE_TIMESCALE_COEFF_C1': 1.63},\n",
-    "    isotope_relaxation_timescale='Bolin1958',\n",
-    ")\n",
-    "const = formulae.constants\n",
-    "any_non_zero_value = 44.\n",
-    "temperature = const.T0 + 10 * si.K\n",
-    "pressure = const.p_STP\n",
-    "pvs = formulae.saturation_vapour_pressure.pvs_water(temperature)\n",
-    "v_term = formulae.terminal_velocity.v_term(radii)\n",
-    "eta_air=formulae.air_dynamic_viscosity.eta_air(temperature)\n",
-    "D=formulae.diffusion_thermics.D(T=temperature, p=pressure)\n",
-    "\n",
-    "air_density = pressure/const.Rd/temperature\n",
-    "assert abs(air_density - 1)/air_density <.3\n",
-    "Re = formulae.particle_shape_and_density.reynolds_number(\n",
-    "     radius=radii,\n",
-    "     velocity_wrt_air=v_term,\n",
-    "     dynamic_viscosity=eta_air,\n",
-    "     density=air_density,\n",
-    ")\n",
-    "Sc = formulae.trivia.air_schmidt_number(\n",
-    "    dynamic_viscosity=eta_air, \n",
-    "    diffusivity=D, \n",
-    "    density=air_density,\n",
-    ")\n",
-    "F = formulae.ventilation.ventilation_coefficient(sqrt_re_times_cbrt_sc=Re**(1/2) * Sc**(1/3))\n",
-    "\n",
-    "r_dr_dt = formulae.drop_growth.r_dr_dt(\n",
-    "    RH_eq=1,\n",
-    "    T=temperature,\n",
-    "    RH=0,\n",
-    "    lv=0,\n",
-    "    pvs=pvs,\n",
-    "    D=D,\n",
-    "    K=any_non_zero_value,\n",
-    "    ventilation_factor=F\n",
-    ")\n",
-    "adjustment_time = formulae.isotope_relaxation_timescale.tau(radius = radii, r_dr_dt = r_dr_dt)\n",
-    "\n",
-    "\n",
     "pandas.options.display.float_format = '{:>,.2g}'.format\n",
     "data = pandas.DataFrame({\n",
     "    'radius [cm]': in_unit(radii, si.cm),\n",
     "    'adjustment time [s]': adjustment_time,\n",
-    "    'terminal velocity [m/s]': v_term,\n",
-    "    'distance [m]': v_term * adjustment_time,\n",
+    "    'terminal velocity [m/s]': iso_fun.v_term,\n",
+    "    'distance [m]': iso_fun.v_term * adjustment_time,\n",
     "})\n",
     "\n",
     "data # pylint: disable=pointless-statement"
    ],
-   "id": "314f42c310883cfa",
+   "id": "340812de267c4290",
    "outputs": [
     {
      "data": {
       "text/plain": [
        "   radius [cm]  adjustment time [s]  terminal velocity [m/s]  distance [m]\n",
-       "0        0.005                  1.7                      0.4          0.69\n",
-       "1         0.01                  5.4                      0.8           4.3\n",
-       "2        0.025                   20                        2            40\n",
-       "3         0.05                   48                        4       1.9e+02\n",
-       "4        0.075                   81                      5.5       4.4e+02\n",
-       "5          0.1              1.2e+02                      6.4       7.6e+02\n",
-       "6         0.15                2e+02                      7.8       1.6e+03\n",
-       "7          0.2                3e+02                        9       2.7e+03"
+       "0        0.005                  4.6                      0.4           1.8\n",
+       "1         0.01                   14                      0.8            11\n",
+       "2        0.025                   54                        2       1.1e+02\n",
+       "3         0.05              1.3e+02                        4       5.1e+02\n",
+       "4        0.075              2.1e+02                      5.5       1.2e+03\n",
+       "5          0.1              3.2e+02                      6.4         2e+03\n",
+       "6         0.15              5.4e+02                      7.8       4.2e+03\n",
+       "7          0.2              7.9e+02                        9       7.1e+03"
       ],
       "text/html": [
        "<div>\n",
@@ -177,58 +148,58 @@
        "    <tr>\n",
        "      <th>0</th>\n",
        "      <td>0.005</td>\n",
-       "      <td>1.7</td>\n",
+       "      <td>4.6</td>\n",
        "      <td>0.4</td>\n",
-       "      <td>0.69</td>\n",
+       "      <td>1.8</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>1</th>\n",
        "      <td>0.01</td>\n",
-       "      <td>5.4</td>\n",
+       "      <td>14</td>\n",
        "      <td>0.8</td>\n",
-       "      <td>4.3</td>\n",
+       "      <td>11</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2</th>\n",
        "      <td>0.025</td>\n",
-       "      <td>20</td>\n",
+       "      <td>54</td>\n",
        "      <td>2</td>\n",
-       "      <td>40</td>\n",
+       "      <td>1.1e+02</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>3</th>\n",
        "      <td>0.05</td>\n",
-       "      <td>48</td>\n",
+       "      <td>1.3e+02</td>\n",
        "      <td>4</td>\n",
-       "      <td>1.9e+02</td>\n",
+       "      <td>5.1e+02</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>4</th>\n",
        "      <td>0.075</td>\n",
-       "      <td>81</td>\n",
+       "      <td>2.1e+02</td>\n",
        "      <td>5.5</td>\n",
-       "      <td>4.4e+02</td>\n",
+       "      <td>1.2e+03</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>5</th>\n",
        "      <td>0.1</td>\n",
-       "      <td>1.2e+02</td>\n",
+       "      <td>3.2e+02</td>\n",
        "      <td>6.4</td>\n",
-       "      <td>7.6e+02</td>\n",
+       "      <td>2e+03</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>6</th>\n",
        "      <td>0.15</td>\n",
-       "      <td>2e+02</td>\n",
+       "      <td>5.4e+02</td>\n",
        "      <td>7.8</td>\n",
-       "      <td>1.6e+03</td>\n",
+       "      <td>4.2e+03</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>7</th>\n",
        "      <td>0.2</td>\n",
-       "      <td>3e+02</td>\n",
+       "      <td>7.9e+02</td>\n",
        "      <td>9</td>\n",
-       "      <td>2.7e+03</td>\n",
+       "      <td>7.1e+03</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
@@ -241,6 +212,19 @@
     }
    ],
    "execution_count": 4
+  },
+  {
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2025-02-07T14:51:44.668885Z",
+     "start_time": "2025-02-07T14:51:44.666532Z"
+    }
+   },
+   "cell_type": "code",
+   "source": "",
+   "id": "eeec7e163ba8c4ea",
+   "outputs": [],
+   "execution_count": null
   }
  ],
  "metadata": {