Merge branch 'feature/dilution_update' into 'master'

Dilution factor update

See merge request cara/caimira!408
This commit is contained in:
Andre Henriques 2022-11-09 15:34:13 +01:00
commit dedf5f2bf8
4 changed files with 58 additions and 42 deletions

View file

@ -142,10 +142,14 @@ In addition, for each individual interaction, the expiration type may be differe
To calculate the short-range component, we first need to calculate what is the **dilution factor**, that depends on the distance :math:`x` as a random variable, from a log normal distribution in :meth:`caimira.monte_carlo.data.short_range_distances`.
This factor is calculated in a two-stage expiratory jet model, with its transition point defined as follows:
:math:`\mathrm{xstar}=𝛽_{\mathrm{x1}} (Q_{0} \cdot u_{0})^\frac{1}{4} \cdot (\mathrm{tstar} + t_{0})^\frac{1}{2} - x_{0}`,
:math:`\mathrm{xstar}=𝛽_{\mathrm{x1}} (Q_{\mathrm{exh}} \cdot u_{0})^\frac{1}{4} \cdot (\mathrm{tstar} + t_{0})^\frac{1}{2} - x_{0}`,
where the :math:`Q_{0}` is the expired flow rate during the expiration period, in :math:`m^{3} s^{-1}`, :math:`u_{0}` is the expired jet speed (in :math:`m s^{-1}`) given by :math:`u_{0}=\frac{Q_{0}}{A_{m}}`, :math:`A_{m}` being the area of the mouth assuming a perfect circle (average `mouth_diameter` of `0.02m`).
The time of the transition point :math:`\mathrm{tstar}` is defined as `2s` and corresponds to the end of the exhalation period, i.e. when the jet is interrupted. The distance of the virtual origin of the puff-like stage is defined by :math:`x_{0}=\frac{\textrm{mouth_diameter}}{2𝛽_{\mathrm{r1}}}` (in m), and the corresponding time is given by :math:`t_{0} = \frac{\sqrt{\pi}D^3}{8𝛽_{\mathrm{r1}}^2𝛽_{\mathrm{x1}}^2Q_{0}}` (in s).
where :math:`Q_{\mathrm{exh}}= φ \mathrm{BR}` is the expired flow rate during the expiration period, in :math:`m^{3} s^{-1}`, `φ` is the exhalation coefficient
(dimensionless) and represents the ratio between the total period of a breathing cycle and the duration of the exhalation alone.
Assuming the duration of the inhalation part is equal to the exhalation and one starts immediately after the other, `φ` will always be equal to `2` no matter what is the breating cycle time. :math:`\mathrm{BR}` is the given exhalation rate.
:math:`u_{0}` is the expired jet speed (in :math:`m s^{-1}`) given by :math:`u_{0}=\frac{Q_{\mathrm{exh}}}{A_{m}}`, :math:`A_{m}` being the area of the mouth assuming a perfect circle (average `mouth_diameter` of `0.02m`).
The time of the transition point :math:`\mathrm{tstar}` is defined as `2s` and corresponds to the end of the exhalation period, i.e. when the jet is interrupted. The distance of the virtual origin of the puff-like stage is defined by
:math:`x_{0}=\frac{\textrm{mouth_diameter}}{2𝛽_{\mathrm{r1}}}` (in m), and the corresponding time is given by :math:`t_{0} = \frac{\sqrt{\pi} \cdot \textrm{mouth_diameter}^3}{8𝛽_{\mathrm{r1}}^2𝛽_{\mathrm{x1}}^2Q_{exh}}` (in s).
Having the distance for the transition point, we can calculate the dilution factor at the transition point, defined as follows:
:math:`\mathrm{Sxstar}=2𝛽_{\mathrm{r1}}\frac{(xstar + x_{0})}{\textrm{mouth_diameter}}`.

View file

@ -1154,33 +1154,44 @@ class ShortRangeModel:
'''
The dilution factor for the respective expiratory activity type.
'''
# Average mouth diameter
# Average 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.
φ = 2
# Exhalation airflow, as per Jia et al. (2022)
Q_exh = φ * 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 period(s), assuming a 4s breath-cycle
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)
# 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)
distances = np.array(self.distance)
factors = np.empty(distances.shape, dtype=np.float64)
factors[distances < xstar] = 2*𝛽r1*(distances[distances < xstar]
+ x0)/mouth_diameter

View file

@ -53,11 +53,11 @@ def test_short_range_model_ndarray(concentration_model, short_range_model):
@pytest.mark.parametrize(
"activity, expected_dilution", [
["Seated", 176.04075727780327],
["Standing", 157.12965288170005],
["Light activity", 69.06672998536413],
["Moderate activity", 47.165817446310115],
["Heavy exercise", 23.759992220217875],
["Seated", 85.73002264],
["Standing", 76.19303543],
["Light activity", 32.45103906],
["Moderate activity", 21.79749405],
["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.6, 5.401601371244907],
[11.0, 5.401601371244907],
[10.5, 11.266605],
[10.6, 11.266605],
[11.0, 11.266605],
[12.0, 0.],
]
)

View file

@ -196,6 +196,9 @@ class SimpleShortRangeModel:
#: Breathing rate (m^3/h)
breathing_rate: _VectorisedFloat = 0.51
#: Exhalation coefficient
φ = 2
#: Tuple with BLO factors
BLO_factors: typing.Tuple[float, float, float] = (1,0,0)
@ -206,17 +209,16 @@ class SimpleShortRangeModel:
#: Maximum diameter for integration (short-range only) (microns)
diameter_max: float = 100.
#: Mouth opening diameter (m)
D: float = 0.02
#: Average mouth opening diameter (m)
mouth_diameter: float = 0.02
#: Duration of the expiration (s)
#: Duration of the expiration period(s), assuming a 4s breath-cycle
tstar: float = 2.
#: Streamwise and radial penetration coefficients
Cr1: float = 0.18
Cx1: float = 2.4
Cr2: float = 0.2
Cx2: float = 2.2
𝛽r1: float = 0.18
𝛽r2: float = 0.2
𝛽x1: float = 2.4
@method_cache
def dilution_factor(self) -> _VectorisedFloat:
@ -226,26 +228,25 @@ 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)
# Exhalation airflow, as per Jia et al. (2022), m^3/s
Q_exh = self.φ * 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
- x01)
xstar = np.array(self.𝛽x1*(Q_exh*u0)**0.25*(self.tstar + t0)**0.5 - x0)
# 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] = Sxstar[x > xstar]*(1. + self.Cr2*(x[x > xstar]
dilution[x <= xstar] = 2.*self.𝛽r1*(x[x <= xstar] + x0)/self.mouth_diameter
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