Adapted dilution factor formula with exhalation coefficient

caimira/models.py

@@ -1154,28 +1154,42 @@ class ShortRangeModel:
         '''
         The dilution factor for the respective expiratory activity type.
         '''
-        # Average mouth diameter
+        # Mouth opening diameter (m)
         mouth_diameter = 0.02
-        # Convert Breathing rate from m3/h to m3/s
-        BR = np.array(self.activity.exhalation_rate/3600.)
-        # Area of the mouth assuming a perfect circle
-        Am = np.pi*(mouth_diameter**2)/4 
-        # Initial velocity from the division of the Breathing rate with the area
-        u0 = np.array(BR/Am)
 
+        # Breathing rate, from m3/h to m3/s
+        BR = np.array(self.activity.exhalation_rate/3600.)
+
+        # Exhalation coefficient. Ratio between the duration of a breathing cycle and the duration of 
+        # the exhalation. 4 sec breathing cycle assumed.
+        exh_coef = 2
+
+        # Exhalation airflow
+        Q_exh = exh_coef * BR
+
+        # Area of the mouth assuming a perfect circle (m2)
+        Am = np.pi*(mouth_diameter**2)/4 
+
+        # Initial velocity of the exhalation airflow (m/s)
+        u0 = np.array(Q_exh/Am)
+
+        #: Duration of the expiration (s)
         tstar = 2.0
+        
+        #: Streamwise and radial penetration coefficients
         𝛽r1 = 0.18
         𝛽r2 = 0.2
         𝛽x1 = 2.4
 
-        # The expired flow rate during the expiration period, m^3/s
-        Q0 = u0 * np.pi/4*mouth_diameter**2 
         # Parameters in the jet-like stage
+        # Position of virtual origin
         x0 = mouth_diameter/2/𝛽r1
         # Time of virtual origin
-        t0 = (x0/𝛽x1)**2 * (Q0*u0)**(-0.5)
+        t0 = (np.sqrt(np.pi)*(mouth_diameter**3))/(8*(𝛽r1**2)*(𝛽x1**2)*Q_exh)
+        # Aux to test         
+        t0_test = (x0/𝛽x1)**2 * (Am*u0**2)**(-0.5)
         # The transition point, m
-        xstar = np.array(𝛽x1*(Q0*u0)**0.25*(tstar + t0)**0.5 - x0)
+        xstar = np.array(𝛽x1*(Q_exh*u0)**0.25*(tstar + t0)**0.5 - x0)
         # Dilution factor at the transition point xstar
         Sxstar = np.array(2*𝛽r1*(xstar+x0)/mouth_diameter)
 

caimira/tests/models/test_short_range_model.py

@@ -53,11 +53,11 @@ def test_short_range_model_ndarray(concentration_model, short_range_model):
 
 @pytest.mark.parametrize(
     "activity, expected_dilution", [
-        ["Seated", 176.04075727780327],
+        ["Seated", 85.73002264],
-        ["Standing", 157.12965288170005],
+        ["Standing", 76.19303543],
-        ["Light activity", 69.06672998536413],
+        ["Light activity", 32.45103906],
-        ["Moderate activity", 47.165817446310115],
+        ["Moderate activity", 21.79749405],
-        ["Heavy exercise", 23.759992220217875],
+        ["Heavy exercise", 16.372],
     ]
 )
 def test_dilution_factor(activity, expected_dilution):
@@ -67,7 +67,7 @@ def test_dilution_factor(activity, expected_dilution):
                                     distance=0.854).build_model(SAMPLE_SIZE)
     assert isinstance(model.dilution_factor(), np.ndarray)
     np.testing.assert_almost_equal(
-        model.dilution_factor(), expected_dilution, decimal=10
+        model.dilution_factor(), expected_dilution
     )
 
 
@@ -100,9 +100,9 @@ def test_extract_between_bounds(short_range_model, time1, time2,
 @pytest.mark.parametrize(
     "time, expected_short_range_concentration", [
         [8.5, 0.],
-        [10.5, 5.401601371244907],
+        [10.5, 11.266605],
-        [10.6, 5.401601371244907],
+        [10.6, 11.266605],
-        [11.0, 5.401601371244907],
+        [11.0, 11.266605],
         [12.0, 0.],
     ]
 )

caimira/tests/test_full_algorithm.py

@@ -196,6 +196,9 @@ class SimpleShortRangeModel:
     
     #: Breathing rate (m^3/h)
     breathing_rate: _VectorisedFloat = 0.51
+
+    #: Exhalation coefficient
+    exh_coef = 2
     
     #: Tuple with BLO factors
     BLO_factors: typing.Tuple[float, float, float] = (1,0,0)
@@ -207,16 +210,15 @@ class SimpleShortRangeModel:
     diameter_max: float = 100.
     
     #: Mouth opening diameter (m)
-    D: float = 0.02
+    mouth_diameter: float = 0.02
     
     #: Duration of the expiration (s)
     tstar: float = 2.
     
     #: Streamwise and radial penetration coefficients
-    Cr1: float = 0.18
+    𝛽r1: float = 0.18
-    Cx1: float = 2.4
+    𝛽r2: float = 0.2
-    Cr2: float = 0.2
+    𝛽x1: float = 2.4
-    Cx2: float = 2.2
     
     @method_cache
     def dilution_factor(self) -> _VectorisedFloat:
@@ -227,25 +229,24 @@ class SimpleShortRangeModel:
         x = np.array(self.distance)
         dilution = np.empty(x.shape, dtype=np.float64)
         # Expired flow rate during the expiration period, m^3/s
-        Q0 = np.array(self.breathing_rate/3600)
+        Q_exh = self.exh_coef * np.array(self.breathing_rate/3600)
         # The expired flow velocity at the noozle (mouth opening), m/s
-        u0 = np.array(Q0/(np.pi/4. * self.D**2))
+        u0 = np.array(Q_exh/(np.pi/4. * self.mouth_diameter**2))
         # Parameters in the jet-like stage
         # position of virtual origin
-        x01 = self.D/2/self.Cr1
+        x0 = self.mouth_diameter/2/self.𝛽r1
         # Time of virtual origin
-        t01 = (x01/self.Cx1)**2 * (Q0*u0)**(-0.5)
+        t0 = (x0/self.𝛽x1)**2 * (Q_exh*u0)**(-0.5)
         # Transition point (in m)
-        xstar = np.array(self.Cx1*(Q0*u0)**0.25*(self.tstar + t01)**0.5
+        xstar = np.array(self.𝛽x1*(Q_exh*u0)**0.25*(self.tstar + t0)**0.5 - x0)
-                         - x01)
         # Dilution factor at the transition point xstar
-        Sxstar = np.array(2.*self.Cr1*(xstar+x01)/self.D)
+        Sxstar = np.array(2.*self.𝛽r1*(xstar+x0)/self.mouth_diameter)
 
         # Calculate dilution factor at the short-range distance x
-        dilution[x <= xstar] = 2.*self.Cr1*(x[x <= xstar] + x01)/self.D
+        dilution[x <= xstar] = 2.*self.𝛽r1*(x[x <= xstar] + x0)/self.mouth_diameter
-        dilution[x > xstar] = Sxstar[x > xstar]*(1. + self.Cr2*(x[x > xstar]
+        dilution[x > xstar] = Sxstar[x > xstar]*(1. + self.𝛽r2*(x[x > xstar]
                                 - xstar[x > xstar])
-                                /self.Cr1/(xstar[x > xstar] + x01))**3
+                                /self.𝛽r1/(xstar[x > xstar] + x0))**3
 
         return dilution