diff --git a/app/cara.ipynb b/app/cara.ipynb index 906bf165..931d3197 100644 --- a/app/cara.ipynb +++ b/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 -} +} \ No newline at end of file diff --git a/cara/models.py b/cara/models.py index 82052fd9..31649aa4 100644 --- a/cara/models.py +++ b/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 minutes every 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 minutes every 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) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 72f0c708..206c5f05 100644 --- a/cara/tests/test_known_quantities.py +++ b/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)