Merge branch 'ventilation' into 'master'

Ventilation

See merge request cara/cara!2

app/cara.ipynb

@@ -42,7 +42,8 @@
     "    \"\"\"\n",
     "    model = cara.models.Model(\n",
     "        room=cara.models.Room(volume=volume),\n",
-    "        ventilation=cara.models.Ventilation(),\n",
+    "        ventilation=cara.models.PeriodicWindow(period=120, duration=120, inside_temp=293, outside_temp=283,\n",
+    "                                               window_height=1.6, opening_length=0.6, cd_b=0.6),\n",
     "        infected=cara.models.InfectedPerson(\n",
     "            virus=cara.models.Virus.types['SARS_CoV_2'],\n",
     "            present_times=((0, 4), (5, 8)),\n",
@@ -145,7 +146,7 @@
    ],
    "source": [
     "ventilation_widgets = {\n",
-    "    'Natural': widgets.Label('Currently hard-coded to 514.74 / volume'),\n",
+    "    'Natural': widgets.Label('Currently hard-coded to window-example from mathematica notebook'),\n",
     "    'other': widgets.Label('Not yet implemented.')\n",
     "}\n",
     "for name, widget in ventilation_widgets.items():\n",
@@ -288,4 +289,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 4
-}
+}

cara/models.py

@@ -1,6 +1,7 @@
 import functools
 import numpy as np
 import typing
+from abc import abstractmethod
 
 
 from dataclasses import dataclass
@@ -8,15 +9,67 @@ from dataclasses import dataclass
 
 @dataclass(frozen=True)
 class Room:
+    # The total volume of the room
     volume: int
 
 
 @dataclass(frozen=True)
 class Ventilation:
-    QairNat: float = 514.74
+    """
+    An abstract class for ventilation schemes
+    """
 
-    def air_change_per_hour(self, room: Room):
-        return self.QairNat / room.volume
+    @abstractmethod
+    def air_exchange(self, room: Room, time: float) -> float:
+        # Returns the rate at which air is being exchanged in the given room per cubic meter at a given time
+        pass
+
+
+@dataclass(frozen=True)
+class PeriodicWindow(Ventilation):
+
+    # The window is opened for <duration> minutes every <period> minutes
+    period: int   #: How often the window is opened (minutes)
+    duration: int   #: How long the window remains opened for (minutes)
+
+    inside_temp: float   #: The temperature inside the room (Kelvin)
+    outside_temp: float   #: The temperature outside of the window (Kelvin)
+
+    window_height: float   #: The height of the window
+
+    opening_length: float   #: The length of the opening-gap when the window is open
+
+    cd_b: float = 0.6   #: Discharge coefficient: what portion effective area is used to exchange air (0 <= cd_b <= 1)
+
+    def air_exchange(self, room: Room, time: float) -> float:
+        # Returns the rate at which air is being exchanged in the given room per cubic meter at a given time
+
+        # If the window is closed, no air is being exchanged
+        if time % self.period < (self.period - self.duration):
+            return 0
+
+        root = np.sqrt(9.81 * self.window_height * (abs(self.inside_temp - self.outside_temp)) / self.outside_temp)
+
+        return (3600 / (3 * room.volume)) * self.cd_b * self.window_height * self.opening_length * root
+
+
+@dataclass(frozen=True)
+class PeriodicHEPA(Ventilation):
+
+    # The HEPA is switched on for <duration> minutes every <period> minutes
+    period: int   #: How often the HEPA is switched on (minutes)
+    duration: int   #: How long the HEPA remains switched on for (minutes)
+
+    q_air_mech: float   #: The rate at which the HEPA exchanges air (when switched on)
+
+    def air_exchange(self, room: Room, time: float) -> float:
+        # Returns the rate at which air is being exchanged in the given room per cubic meter at a given time
+
+        # If the HEPA is off, no air is being exchanged
+        if time % self.period < (self.period - self.duration):
+            return 0
+
+        return self.q_air_mech / room.volume
 
 
 @dataclass(frozen=True)
@@ -195,8 +248,7 @@ class Model:
     def virus(self):
         return self.infected.virus
 
-    @property
-    def infectious_virus_removal_rate(self):
+    def infectious_virus_removal_rate(self, time: float) -> float:
         # Particle deposition on the floor
         vg = 1 * 10 ** -4
         # Height of the emission source to the floor - i.e. mouth/nose (m)
@@ -204,12 +256,12 @@ class Model:
         # Deposition rate (h^-1)
         k = (vg * 3600) / h
 
-        return k + self.virus.decay_constant + self.ventilation.air_change_per_hour(self.room)
+        return k + self.virus.decay_constant + self.ventilation.air_exchange(self.room, time)
 
     @functools.lru_cache()
     def concentration(self, time: float) -> float:
         t = time
-        IVRR = self.infectious_virus_removal_rate
+        IVRR = self.infectious_virus_removal_rate(time)
         V = self.room.volume
         Ni = self.infected_occupants
         ER = self.infected.emission_rate(time)

cara/tests/test_known_quantities.py

@@ -1,5 +1,6 @@
 import numpy.testing as npt
 import pytest
+from numpy import linspace
 
 import cara.models as models
 
@@ -26,7 +27,8 @@ def test_no_mask_emission_rate(baseline_model):
 def baseline_model():
     model = models.Model(
         room=models.Room(volume=75),
-        ventilation=models.Ventilation(),
+        ventilation=models.PeriodicWindow(period=120, duration=120, inside_temp=293, outside_temp=283, cd_b=0.6,
+                                          window_height=1.6, opening_length=0.6),
         infected=models.InfectedPerson(
             virus=models.Virus.types['SARS_CoV_2'],
             present_times=((0, 4), (5, 8)),
@@ -41,12 +43,42 @@ def baseline_model():
     return model
 
 
+@pytest.fixture
+def baseline_periodic_window():
+    return models.PeriodicWindow(period=120, duration=15, inside_temp=293, outside_temp=283, cd_b=0.6,
+                                 window_height=1.6, opening_length=0.6)
+
+
+@pytest.fixture
+def baseline_room():
+    return models.Room(volume=75)
+
+
+@pytest.fixture
+def baseline_periodic_hepa():
+    return models.PeriodicHEPA(period=120, duration=15, q_air_mech=514.74)
+
+
 def test_r0(baseline_model):
     saturated = 2.909312e-01
     ts = [0, 4, 5, 7, 10]
     concentrations = [baseline_model.concentration(t) for t in ts]
     npt.assert_allclose(
         concentrations,
-        [0.000000e+00, saturated, 1.274225e-04, saturated, 5.580870e-08],
+        [0.000000e+00, 2.909211e-01, 1.273836e-04, 2.909210e-01, 5.577662e-08],
         rtol=1e-5
     )
+
+
+def test_periodic_window(baseline_periodic_window, baseline_room):
+    ts = linspace(0, 9 * 60, 9)
+    aes = [baseline_periodic_window.air_exchange(baseline_room, t) for t in ts]
+    answers = [0, 0, 0, 0, 0, 0, 0, 514.76 / 75, 0]
+    npt.assert_allclose(aes, answers, rtol=1e-5)
+
+
+def test_periodic_hepa(baseline_periodic_hepa, baseline_room):
+    ts = linspace(0, 9 * 60, 9)
+    aes = [baseline_periodic_hepa.air_exchange(baseline_room, t) for t in ts]
+    answers = [0, 0, 0, 0, 0, 0, 0, 514.74 / 75, 0]
+    npt.assert_allclose(aes, answers, rtol=1e-5)