(Created page with "<!-- metadata commented in wiki content ==== ''', , Rainald Löhner <math>\cdot </math>Harbir Antil <math>\cdot </math>Sergio Idelsohn <math>\cdot </math> Eugenio Oñ...")
 
Line 1: Line 1:
 
<!-- metadata commented in wiki content
 
<!-- metadata commented in wiki content
====
+
==DETAILED SIMULATION OF VIRAL PROPAGATION IN THE BUILT ENVIRONMENT==
  
''', , Rainald Löhner <math>\cdot </math>Harbir Antil <math>\cdot </math>Sergio Idelsohn <math>\cdot </math>        Eugenio Oñate'''
+
'''Rainald Löhner Harbir Antil Sergio Idelsohn         Eugenio Oñate'''
 
-->
 
-->
  
Line 9: Line 9:
 
A summary is given of the mechanical characteristics of virus contaminants and the transmission via droplets and aerosols. The ordinary and partial differential equations describing the physics of these processes with high fidelity are presented, as well as appropriate numerical schemes to solve them. Several examples taken from recent evaluations of the built environment are shown, as well as the optimal placement of sensors.
 
A summary is given of the mechanical characteristics of virus contaminants and the transmission via droplets and aerosols. The ordinary and partial differential equations describing the physics of these processes with high fidelity are presented, as well as appropriate numerical schemes to solve them. Several examples taken from recent evaluations of the built environment are shown, as well as the optimal placement of sensors.
  
'''keywords''' Covid-19 <math>\cdot </math>Particle Methods <math>\cdot </math>Finite Elements <math>\cdot </math>          Computational Fluid Dynamics
+
'''keywords''' Covid-19 Particle Methods Finite Elements           Computational Fluid Dynamics
  
==2 The Covid-19 Crisis==
+
==1 The Covid-19 Crisis==
  
 
Starting in Wuhan, China, in the fall of 2019, the Covid-19 pandemic has claimed and will continue to claim millions of infected patients and hundreds of thousands of deaths. The lockdowns that followed its outbreak have led to mass unemployment, stalled economic development and loss of productivity that will take years to recover. Some changes in habits and lifestyles may be permanent: in the future, working from home or in a `socially distanced manner' may be the prevalent modus operandi for large segments of society.  
 
Starting in Wuhan, China, in the fall of 2019, the Covid-19 pandemic has claimed and will continue to claim millions of infected patients and hundreds of thousands of deaths. The lockdowns that followed its outbreak have led to mass unemployment, stalled economic development and loss of productivity that will take years to recover. Some changes in habits and lifestyles may be permanent: in the future, working from home or in a `socially distanced manner' may be the prevalent modus operandi for large segments of society.  
Line 17: Line 17:
 
The present paper gives a short description of computational techniques that can elucidate the flow and propagation of viruses or other contaminants in built environments in order to mitigate or avoid infections.
 
The present paper gives a short description of computational techniques that can elucidate the flow and propagation of viruses or other contaminants in built environments in order to mitigate or avoid infections.
  
==3 Virus Infection==
+
==2 Virus Infection==
  
 
Before addressing the requirements for the numerical simulation of virus propagation a brief description of virus propagation and lifetime are given. Covid-19 is one of many corona-viruses. The virus is usually present in the air or some surface, and makes its way into the body either via inhalation (nose, mouth), ingestion (mouth) or attachment (eyes, hands, clothes). In many cases the victim inadvertedly touches an infected surface or viruses are deposited on its hands, and then the hands touch either the nose, the eyes or the mouth, thus allowing the virus to enter the body.  
 
Before addressing the requirements for the numerical simulation of virus propagation a brief description of virus propagation and lifetime are given. Covid-19 is one of many corona-viruses. The virus is usually present in the air or some surface, and makes its way into the body either via inhalation (nose, mouth), ingestion (mouth) or attachment (eyes, hands, clothes). In many cases the victim inadvertedly touches an infected surface or viruses are deposited on its hands, and then the hands touch either the nose, the eyes or the mouth, thus allowing the virus to enter the body.  
  
An open question of great importance for all that will follow is how many viruses it takes to overwhelm the body's natural defense mechanism and trigger an infection. This number, which is sometimes called the ''viral load'' or the ''infectious dose'' will depend on numerous factors, among them the state of immune defenses of the individual, the timing of viral entry (all at once, piece by piece), and the amount of hair and mucous in the nasal vessels. In principle, a single organism in a favourable environment may replicate sufficiently to cause disease <span id='citeF-98'></span>[[#cite-98|[98]]]. Data from research performed on biological warfare agents <span id='citeF-37'></span>[[#cite-37|[37]]] suggests that both bacteria and viruses can produce disease with as few as 1-100 organisms (e.g. brucellosis 10-100, Q fever 1-10, tularaemia 10-50, smallpox 10-100, viral haemorrhagic fevers 1-10 organisms, tuberculosis 1). Compare these numbers and consider that as many as 3,000 organisms can be produced by talking for 5&nbsp;minutes or a single cough, with sneezing producing many more <span id='citeF-80'></span><span id='citeF-55'></span><span id='citeF-102'></span><span id='citeF-87'></span><span id='citeF-108'></span>[[#cite-80|[80,55,102,87,108]]]. Figure&nbsp;1, reproduced from <span id='citeF-102'></span>[[#cite-102|[102]]], shows a typical number and size distribution.
+
An open question of great importance for all that will follow is how many viruses it takes to overwhelm the body's natural defense mechanism and trigger an infection. This number, which is sometimes called the ''viral load'' or the ''infectious dose'' will depend on numerous factors, among them the state of immune defenses of the individual, the timing of viral entry (all at once, piece by piece), and the amount of hair and mucous in the nasal vessels. In principle, a single organism in a favourable environment may replicate sufficiently to cause disease <span id='citeF-98'></span>[[#cite-98|[98]]]. Data from research performed on biological warfare agents <span id='citeF-37'></span>[[#cite-37|[37]]] suggests that both bacteria and viruses can produce disease with as few as 1-100 organisms (e.g. brucellosis 10-100, Q fever 1-10, tularaemia 10-50, smallpox 10-100, viral haemorrhagic fevers 1-10 organisms, tuberculosis 1). Compare these numbers and consider that as many as 3,000 organisms can be produced by talking for 5&nbsp;minutes or a single cough, with sneezing producing many more <span id='citeF-80'></span><span id='citeF-55'></span><span id='citeF-102'></span><span id='citeF-87'></span><span id='citeF-108'></span>[[#cite-80|[80,55,102,87,108]]]. Figure [[#img-1|1]], reproduced from <span id='citeF-102'></span>[[#cite-102|[102]]], shows a typical number and size distribution.
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">  figures/general/droplet_distribution.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 1&nbsp;&nbsp; Counts of Particles of Various Diameters in Air Expelled by (90) Coughs <span id='citeF-80'></span>[[#cite-80|[80]]] </div>
+
<div id='img-1'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_514971305-droplet_distribution.png|600px|Counts of Particles of Various Diameters in Air Expelled by (90) Coughs <span id='citeF-80'></span>[[#cite-80|[80]]]]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="1" | '''Figure 1:''' Counts of Particles of Various Diameters in Air Expelled by (90) Coughs <span id='citeF-80'></span>[[#cite-80|[80]]]
 +
|}
  
==4 Virus Lifetime Outside the Body==
+
==3 Virus Lifetime Outside the Body==
  
 
Current evidence points to lifetimes outside the body that can range from 1-2 hours in air to several days on particular surfaces <span id='citeF-51'></span><span id='citeF-106'></span>[[#cite-51|[51,106]]]. There has also been some documentation of lifetime variation depending on humidity.
 
Current evidence points to lifetimes outside the body that can range from 1-2 hours in air to several days on particular surfaces <span id='citeF-51'></span><span id='citeF-106'></span>[[#cite-51|[51,106]]]. There has also been some documentation of lifetime variation depending on humidity.
  
==5 Virus Transmission==
+
==4 Virus Transmission==
  
===5.1 Human Sneezing and Coughing===
+
===4.1 Human Sneezing and Coughing===
  
 
In the sequel, we consider human sneezing and coughing as the main conduits of virus transmission. Clearly, breathing and talking will lead to the exhalation of air, and, consequently the exhalation of viruses for infected victims <span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-3|[3,4]]]. However, it stands to reason that the size and amount of particles released -&nbsp;and hence the amount of viruses in them&nbsp;- is much higher and much more concentrated when sneezing or coughing <span id='citeF-36'></span><span id='citeF-102'></span><span id='citeF-50'></span><span id='citeF-56'></span><span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-36|[36,102,50,56,3,4]]].  
 
In the sequel, we consider human sneezing and coughing as the main conduits of virus transmission. Clearly, breathing and talking will lead to the exhalation of air, and, consequently the exhalation of viruses for infected victims <span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-3|[3,4]]]. However, it stands to reason that the size and amount of particles released -&nbsp;and hence the amount of viruses in them&nbsp;- is much higher and much more concentrated when sneezing or coughing <span id='citeF-36'></span><span id='citeF-102'></span><span id='citeF-50'></span><span id='citeF-56'></span><span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-36|[36,102,50,56,3,4]]].  
  
The velocity of air at a person's mouth during sneezing and coughing has been a source of heated debate, particularly in the media. The experimental evidence points to exit velocities of the order of 2-14&nbsp;m/sec <span id='citeF-40'></span><span id='citeF-41'></span><span id='citeF-100'></span><span id='citeF-101'></span>[[#cite-40|[40,41,100,101]]]. A typical amount and size of particles can be seen in Figure&nbsp;1.
+
The velocity of air at a person's mouth during sneezing and coughing has been a source of heated debate, particularly in the media. The experimental evidence points to exit velocities of the order of 2-14&nbsp;m/sec <span id='citeF-40'></span><span id='citeF-41'></span><span id='citeF-100'></span><span id='citeF-101'></span>[[#cite-40|[40,41,100,101]]]. A typical amount and size of particles can be seen in Figure [[#img-1|1]].
  
===5.2 Sink Velocities===
+
===4.2 Sink Velocities===
  
 
If, for the sake of argument, we consider Stoke's law for the drag of spherical particles, valid below Reynolds numbers of <math display="inline">Re=1</math>, the terminal sink velocity (also known as the settling velocity) of particles will be given by:
 
If, for the sake of argument, we consider Stoke's law for the drag of spherical particles, valid below Reynolds numbers of <math display="inline">Re=1</math>, the terminal sink velocity (also known as the settling velocity) of particles will be given by:
Line 131: Line 137:
 
|}
 
|}
  
===5.3 Evaporation===
+
===4.3 Evaporation===
  
 
Depending on the relative humidity and the temperature of the ambient air, the smaller particles can evaporate in milliseconds. However, as the mucous and saliva evaporate, they build a gel-like structure that surrounds the virus, allowing it to survive. This implies that extremely small particles with possible viruses will remain infectious for extended periods of times - up to an hour according to some studies.  
 
Depending on the relative humidity and the temperature of the ambient air, the smaller particles can evaporate in milliseconds. However, as the mucous and saliva evaporate, they build a gel-like structure that surrounds the virus, allowing it to survive. This implies that extremely small particles with possible viruses will remain infectious for extended periods of times - up to an hour according to some studies.  
Line 137: Line 143:
 
An important question is then whether a particle/droplet will first reach the ground or evaporate. Figure&nbsp;2, taken from <span id='citeF-109'></span>[[#cite-109|[109]]], shows that below 120&nbsp;<math display="inline">\mu m</math> the particles evaporate before falling 2&nbsp;m (i.e. reaching the ground).
 
An important question is then whether a particle/droplet will first reach the ground or evaporate. Figure&nbsp;2, taken from <span id='citeF-109'></span>[[#cite-109|[109]]], shows that below 120&nbsp;<math display="inline">\mu m</math> the particles evaporate before falling 2&nbsp;m (i.e. reaching the ground).
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/general/evap_fall_droplets.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 2&nbsp;&nbsp; Evaporation Time and Falling Time of Droplets of Varying Diameter</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 2&nbsp;&nbsp; Evaporation Time and Falling Time of Droplets of Varying Diameter</div>
  
 
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">(<math>T_{p0}=33^oC, T_{\infty }=18^oC, RH=0%</math>) From <span id='citeF-109'></span>[[#cite-109|[109]]] </div>
 
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">(<math>T_{p0}=33^oC, T_{\infty }=18^oC, RH=0%</math>) From <span id='citeF-109'></span>[[#cite-109|[109]]] </div>
  
===5.4 Viral Load===
+
===4.4 Viral Load===
  
 
A central question that requires an answer is then: how many viruses are in these small particles&nbsp;? An approximate answer may be obtained from the experiments that are being carried out on animals to trace and monitor infections. For ferrets <span id='citeF-53'></span>[[#cite-53|[53]]] <math display="inline">O(10^5-10^6)</math> have been used to infect via intranasal swabs, while for mice <span id='citeF-104'></span>[[#cite-104|[104]]] <math display="inline">O(10^4)</math> seem to suffice. Viral titers can vary a lot, but one may assume on the order of <math display="inline">O(10^6)</math> viruses/ml for a nasopharyngeal swab <span id='citeF-53'></span><span id='citeF-104'></span>[[#cite-53|[53,104]]]. Table&nbsp;2 lists the number of viruses per droplet and the number of droplets needed to contain just 1&nbsp;virus. Note that while for a droplet with a diameter of 1&nbsp;mm one can expect <math display="inline">O(500)</math> viruses, only every 2,000th particle of diameter 10&nbsp;<math display="inline">\mu m</math> does contain a single virus.
 
A central question that requires an answer is then: how many viruses are in these small particles&nbsp;? An approximate answer may be obtained from the experiments that are being carried out on animals to trace and monitor infections. For ferrets <span id='citeF-53'></span>[[#cite-53|[53]]] <math display="inline">O(10^5-10^6)</math> have been used to infect via intranasal swabs, while for mice <span id='citeF-104'></span>[[#cite-104|[104]]] <math display="inline">O(10^4)</math> seem to suffice. Viral titers can vary a lot, but one may assume on the order of <math display="inline">O(10^6)</math> viruses/ml for a nasopharyngeal swab <span id='citeF-53'></span><span id='citeF-104'></span>[[#cite-53|[53,104]]]. Table&nbsp;2 lists the number of viruses per droplet and the number of droplets needed to contain just 1&nbsp;virus. Note that while for a droplet with a diameter of 1&nbsp;mm one can expect <math display="inline">O(500)</math> viruses, only every 2,000th particle of diameter 10&nbsp;<math display="inline">\mu m</math> does contain a single virus.
Line 178: Line 184:
 
Similar numbers are seen in field studies as well. The size of viruses varies from 0.02-0.3 <math display="inline">\mu m</math>, while the size of bacteria varies from 0.5-10 <math display="inline">\mu m</math>. The influenza virus RNA detected by quantitative polymerase chain reaction in human exhaled breath suggests that it may be contained in fine particles generated during tidal breathing and not only coughs <span id='citeF-36'></span><span id='citeF-55'></span><span id='citeF-56'></span><span id='citeF-87'></span>[[#cite-36|[36,55,56,87]]]. Influenza RNA and Mycobacterium tuberculosis have been reported in particles that range in size from 0.5-4.0 <math display="inline">\mu m</math> (<span id='citeF-36'></span><span id='citeF-55'></span><span id='citeF-56'></span><span id='citeF-87'></span>[[#cite-36|[36,55,56,87]]] and references cited therein).
 
Similar numbers are seen in field studies as well. The size of viruses varies from 0.02-0.3 <math display="inline">\mu m</math>, while the size of bacteria varies from 0.5-10 <math display="inline">\mu m</math>. The influenza virus RNA detected by quantitative polymerase chain reaction in human exhaled breath suggests that it may be contained in fine particles generated during tidal breathing and not only coughs <span id='citeF-36'></span><span id='citeF-55'></span><span id='citeF-56'></span><span id='citeF-87'></span>[[#cite-36|[36,55,56,87]]]. Influenza RNA and Mycobacterium tuberculosis have been reported in particles that range in size from 0.5-4.0 <math display="inline">\mu m</math> (<span id='citeF-36'></span><span id='citeF-55'></span><span id='citeF-56'></span><span id='citeF-87'></span>[[#cite-36|[36,55,56,87]]] and references cited therein).
  
==6 Physical Modeling of Aerosol Propagation==
+
==5 Physical Modeling of Aerosol Propagation==
  
 
When solving the two-phase equations, the air, as a continuum, is best represented by a set of partial differential equations (the Navier-Stokes equations) that are numerically solved on a mesh. Thus, the gas characteristics are calculated at the mesh points within the flowfield. However, as the droplets/particles may be relatively sparse in the flowfield, they can be modeled by either:
 
When solving the two-phase equations, the air, as a continuum, is best represented by a set of partial differential equations (the Navier-Stokes equations) that are numerically solved on a mesh. Thus, the gas characteristics are calculated at the mesh points within the flowfield. However, as the droplets/particles may be relatively sparse in the flowfield, they can be modeled by either:
Line 189: Line 195:
 
While a Lagrangian approach offers many potential advantages, this method also creates problems that need to be addressed. For instance, large numbers of particles may cause a Lagrangian analysis to be memory intensive. This problem is circumvented by treating parcels of particles, i.e. doing the detailed analysis for one particle and then applying the effect to many. In addition, continuous mapping and remapping of particles to their respective elements may increase computational requirements, particularly for unstructured grids.
 
While a Lagrangian approach offers many potential advantages, this method also creates problems that need to be addressed. For instance, large numbers of particles may cause a Lagrangian analysis to be memory intensive. This problem is circumvented by treating parcels of particles, i.e. doing the detailed analysis for one particle and then applying the effect to many. In addition, continuous mapping and remapping of particles to their respective elements may increase computational requirements, particularly for unstructured grids.
  
===6.1 Equations Describing the Motion of the Air===
+
===5.1 Equations Describing the Motion of the Air===
  
 
As seen from the experimental evidence, the velocities of air encountered during coughing and sneezing never exceed a Mach-number of <math display="inline">Ma=0.1</math>. Therefore, the air may be assumed as a Newtonian, incompressible liquid, where buoyancy effects are modeled via the Boussinesq approximation. The equations describing the conservation of momentum, mass and energy for incompressible, Newtonian flows may be written as
 
As seen from the experimental evidence, the velocities of air encountered during coughing and sneezing never exceed a Mach-number of <math display="inline">Ma=0.1</math>. Therefore, the air may be assumed as a Newtonian, incompressible liquid, where buoyancy effects are modeled via the Boussinesq approximation. The equations describing the conservation of momentum, mass and energy for incompressible, Newtonian flows may be written as
Line 198: Line 204:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> \rho {\boldsymbol v}_{,t} + \rho {\boldsymbol v}\cdot \nabla {\boldsymbol v}+ \nabla p =    \nabla \cdot \mu \nabla {\boldsymbol v}+ \rho {\boldsymbol g}+ \beta \rho {\boldsymbol g}(T - T_0)      + {\boldsymbol s}_v  ~~,                                                              </math>
+
| style="text-align: center;" | <math> \rho \vvec _{,t} + \rho \vvec \cdot \Grd \vvec + \Grd p =    \nabla \cdot \mu \nabla \vvec + \rho \gvec + \beta \rho \gvec (T - T_0)      + \svec _v  ~~,                                                              </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.1)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.1)
Line 208: Line 214:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>                            \nabla \cdot {\boldsymbol v}= 0      ~~.      </math>
+
| style="text-align: center;" | <math>                            \Div \vvec = 0      ~~.      </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.2)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.2)
Line 218: Line 224:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> \rho c_p T_{,t} + \rho c_p {\boldsymbol v}\cdot \nabla T =        \nabla \cdot k \nabla T + s_e ~~,                                                              </math>
+
| style="text-align: center;" | <math> \rho c_p T_{,t} + \rho c_p \vvec \cdot \Grd T =        \nabla \cdot k \nabla T + s_e ~~,                                                              </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.3)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.3)
 
|}
 
|}
  
Here <math display="inline">\rho , {\boldsymbol v}, p, \mu , {\boldsymbol g}, \beta , T, T_0, c_p, k</math> denote the density, velocity vector, pressure, viscosity, gravity vector, coefficient of thermal expansion, temperature, reference temperature, specific heat coefficient and conductivity respectively, and <math display="inline">{\boldsymbol s}_v, s_e</math> momentum and energy source terms (e.g. due to particles or external forces/heat sources). For turbulent flows both the viscosity and the conductivity are obtained either from additional equations or directly via a large eddy simulation (LES) assumption through monotonicity induced LES (MILES) <span id='citeF-20'></span><span id='citeF-38'></span><span id='citeF-39'></span><span id='citeF-47'></span>[[#cite-20|[20,38,39,47]]].
+
Here <math display="inline">\rho , \vvec , p, \mu , \gvec , \beta , T, T_0, c_p, k</math> denote the density, velocity vector, pressure, viscosity, gravity vector, coefficient of thermal expansion, temperature, reference temperature, specific heat coefficient and conductivity respectively, and <math display="inline">\svec _v, s_e</math> momentum and energy source terms (e.g. due to particles or external forces/heat sources). For turbulent flows both the viscosity and the conductivity are obtained either from additional equations or directly via a large eddy simulation (LES) assumption through monotonicity induced LES (MILES) <span id='citeF-20'></span><span id='citeF-38'></span><span id='citeF-39'></span><span id='citeF-47'></span>[[#cite-20|[20,38,39,47]]].
  
===6.2 Equations Describing the Motion of Particles/Droplets===
+
===5.2 Equations Describing the Motion of Particles/Droplets===
  
 
In order to describe the interaction of particles/droplets with the flow, the mass, forces and energy/work exchanged between the flowfield and the particles must be defined. As before, we denote for '''fluid (air)''' by <math display="inline">\rho , p, T, k, v_i, \mu </math> and <math display="inline">c_p</math> the density, pressure, temperature, conductivity, velocity in direction <math display="inline">x_i</math>, viscosity, and the specific heat at constant pressure. For the '''particles''', we denote by <math display="inline">\rho _p, T_p, v_{pi}, d, c_{pp}</math> and <math display="inline">Q</math> the density, temperature, velocity in direction <math display="inline">x_i</math>, equivalent diameter, and heat transferred per unit volume. In what follows, we will refer to droplet and particles, collectively as particles.  
 
In order to describe the interaction of particles/droplets with the flow, the mass, forces and energy/work exchanged between the flowfield and the particles must be defined. As before, we denote for '''fluid (air)''' by <math display="inline">\rho , p, T, k, v_i, \mu </math> and <math display="inline">c_p</math> the density, pressure, temperature, conductivity, velocity in direction <math display="inline">x_i</math>, viscosity, and the specific heat at constant pressure. For the '''particles''', we denote by <math display="inline">\rho _p, T_p, v_{pi}, d, c_{pp}</math> and <math display="inline">Q</math> the density, temperature, velocity in direction <math display="inline">x_i</math>, equivalent diameter, and heat transferred per unit volume. In what follows, we will refer to droplet and particles, collectively as particles.  
  
Making the classical assumptions that the particles may be represented by an equivalent sphere of diameter <math display="inline">d</math>, the drag forces <math display="inline">{\boldsymbol D}</math> acting on the particles will be due to the difference of fluid and particle velocity:
+
Making the classical assumptions that the particles may be represented by an equivalent sphere of diameter <math display="inline">d</math>, the drag forces <math display="inline">\Dvec </math> acting on the particles will be due to the difference of fluid and particle velocity:
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 236: Line 242:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {\boldsymbol D}= {{\pi d^2} \over 4} \cdot c_D \cdot          { 1 \over 2} \rho | {\boldsymbol v}- {\boldsymbol v}_p | ( {\boldsymbol v}- {\boldsymbol v}_p )                                                  ~~.  </math>
+
| style="text-align: center;" | <math> \Dvec = {{\pi d^2} \over 4} \cdot c_D \cdot          { 1 \over 2} \rho | \vvec - \vvec _p | ( \vvec - \vvec _p )                                                  ~~.  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
Line 248: Line 254:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> Re = {{\rho | {\boldsymbol v}- {\boldsymbol v}_p | d } \over { \mu }}  </math>
+
| style="text-align: center;" | <math> Re = {{\rho | \vvec - \vvec _p | d } \over { \mu }}  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
Line 320: Line 326:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> \rho _p {{\pi d^3} \over 6 } \cdot {{ d{\boldsymbol v}_p} \over {dt}} = {\boldsymbol D}    + \rho _p {{\pi d^3} \over 6 } {\boldsymbol g}~~.                                                        </math>
+
| style="text-align: center;" | <math> \rho _p {{\pi d^3} \over 6 } \cdot {{ d\vvec _p} \over {dt}} = \Dvec    + \rho _p {{\pi d^3} \over 6 } \gvec ~~.                                                        </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
Line 332: Line 338:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {{ d{\boldsymbol v}_p} \over {dt}} = {{3 \rho } \over {4 \rho _p d}} \cdot c_d                                | {\boldsymbol v}- {\boldsymbol v}_p | ( {\boldsymbol v}- {\boldsymbol v}_p )                    + {\boldsymbol g}                    = \alpha _v | {\boldsymbol v}- {\boldsymbol v}_p | ( {\boldsymbol v}- {\boldsymbol v}_p )                    + {\boldsymbol g} ~~.                                                          </math>
+
| style="text-align: center;" | <math> {{ d\vvec _p} \over {dt}} = {{3 \rho } \over {4 \rho _p d}} \cdot c_d                                | \vvec - \vvec _p | ( \vvec - \vvec _p )                    + \gvec                      = \alpha _v | \vvec - \vvec _p | ( \vvec - \vvec _p )                    + \gvec  ~~.                                                          </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
Line 344: Line 350:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {{ d{\boldsymbol x}_p} \over {dt}} = {\boldsymbol v}_p ~~.    </math>
+
| style="text-align: center;" | <math> {{ d\xvec _p} \over {dt}} = \vvec _p ~~.    </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
Line 380: Line 386:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {{d{\boldsymbol u}_p} \over {dt}} = {\boldsymbol r}({\boldsymbol u}_p, {\boldsymbol x}, {\boldsymbol u}_f) ~~,                                                          </math>
+
| style="text-align: center;" | <math> {{d\uvec _p} \over {dt}} = \rvec (\uvec _p, \xvec , \uvec _f) ~~,                                                          </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
 
|}
 
|}
  
where <math display="inline">{\boldsymbol u}_p, {\boldsymbol x}, {\boldsymbol u}_f</math> denote the particle unknowns, the position of the particle and the fluid unknowns at the position of the particle.
+
where <math display="inline">\uvec _p, \xvec , \uvec _f</math> denote the particle unknowns, the position of the particle and the fluid unknowns at the position of the particle.
  
===6.3 Numerical Integration of the Motion of the Air===
+
===5.3 Numerical Integration of the Motion of the Air===
  
 
The last six decades have seen a large number of schemes that may be used to solve numerically the incompressible Navier-Stokes equations given by Eqns.(6.1-6.3). In the present case, the following design criteria were implemented:
 
The last six decades have seen a large number of schemes that may be used to solve numerically the incompressible Navier-Stokes equations given by Eqns.(6.1-6.3). In the present case, the following design criteria were implemented:
Line 401: Line 407:
 
The resulting discretization in time is given by the following projection scheme <span id='citeF-68'></span><span id='citeF-71'></span>[[#cite-68|[68,71]]]:
 
The resulting discretization in time is given by the following projection scheme <span id='citeF-68'></span><span id='citeF-71'></span>[[#cite-68|[68,71]]]:
  
* Advective-Diffusive Prediction: <math display="inline">{\boldsymbol v}^n, p^n \rightarrow {\boldsymbol v}^{*}</math>
+
* Advective-Diffusive Prediction: <math display="inline">\vvec ^n, p^n \rightarrow \vvec ^{*}</math>
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 408: Line 414:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> \svec' = - \nabla p^n + \rho {\boldsymbol g}          + \beta \rho {\boldsymbol g}(T^n - T_0) + {\boldsymbol s}_v ~~, </math>
+
| style="text-align: center;" | <math> \svec' = - \Grd p^n + \rho \gvec            + \beta \rho \gvec (T^n - T_0) + \svec _v ~~, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
Line 418: Line 424:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {\boldsymbol v}^i = {\boldsymbol v}^n + \alpha ^i \gamma {\Delta t}\left(  - {\boldsymbol v}^{i-1} \cdot \nabla {\boldsymbol v}^{i-1}    \nabla \cdot \mu \nabla {\boldsymbol v}^{i-1} + \svec' \right) ~~; ~~i=1,k-1~~; </math>
+
| style="text-align: center;" | <math> \vvec ^i = \vvec ^n + \alpha ^i \gamma \dt \left(  - \vvec ^{i-1} \cdot \Grd \vvec ^{i-1}    \nabla \cdot \mu \nabla \vvec ^{i-1} + \svec' \right) ~~; ~~i=1,k-1~~; </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (21a)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (21a)
Line 428: Line 434:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>  \left[{ 1 \over {\Delta t}} - \theta \nabla \cdot \mu \nabla \right]    \left({\boldsymbol v}^{k} - {\boldsymbol v}^n \right)  + {\boldsymbol v}^{k-1} \cdot \nabla {\boldsymbol v}^{k-1} =    \nabla \cdot \mu \nabla {\boldsymbol v}^{k-1} + \svec' ~~.  </math>
+
| style="text-align: center;" | <math>  \left[{ 1 \over \dt } - \theta \nabla \cdot \mu \nabla \right]    \left(\vvec ^{k} - \vvec ^n \right)  + \vvec ^{k-1} \cdot \Grd \vvec ^{k-1} =    \nabla \cdot \mu \nabla \vvec ^{k-1} + \svec' ~~.  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (21b)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (21b)
Line 440: Line 446:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>  \nabla \cdot {\boldsymbol v}^{n+1} = 0                          ~~; </math>
+
| style="text-align: center;" | <math>  \Div \vvec ^{n+1} = 0                          ~~; </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
Line 450: Line 456:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>  {{ {\boldsymbol v}^{n+1} - {\boldsymbol v}^{*} }\over {\Delta t}} + \nabla ( p^{n+1} - p^n )    = 0                                          ~~; </math>
+
| style="text-align: center;" | <math>  {{ \vvec ^{n+1} - \vvec ^{*} }\over \dt } + \Grd ( p^{n+1} - p^n )    = 0                                          ~~; </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
Line 462: Line 468:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>  \nabla ^2 ( p^{n+1} - p^n ) = {{\nabla \cdot {\boldsymbol v}^{*} }\over {\Delta t}} ~~; </math>
+
| style="text-align: center;" | <math>  \nabla ^2 ( p^{n+1} - p^n ) = {{\Div \vvec ^{*} }\over \dt } ~~; </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
 
|}
 
|}
  
* Velocity Correction: <math display="inline">{\boldsymbol v}^{*} \rightarrow {\boldsymbol v}^{n+1}</math>
+
* Velocity Correction: <math display="inline">\vvec ^{*} \rightarrow \vvec ^{n+1}</math>
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 474: Line 480:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math>  {\boldsymbol v}^{n+1} = {\boldsymbol v}^{*} - {\Delta t}\nabla ( p^{n+1} - p^n ) ~~. </math>
+
| style="text-align: center;" | <math>  \vvec ^{n+1} = \vvec ^{*} - \dt \Grd ( p^{n+1} - p^n ) ~~. </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
 
|}
 
|}
  
<math display="inline">\theta </math> denotes the implicitness-factor for the viscous terms (<math display="inline">\theta=1</math>: 1st order, fully implicit, <math display="inline">\theta=0.5</math>: 2nd order, Crank-Nicholson). <math display="inline">\alpha ^i</math> are the standard low-storage Runge-Kutta coefficients <math display="inline">\alpha ^i=1/(k+1-i)</math>. The <math display="inline">k-1</math> stages of Eqn.(21a) may be seen as a predictor (or replacement) of <math display="inline">{\boldsymbol v}^n</math> by <math display="inline">{\boldsymbol v}^{k-1}</math>. The original right-hand side has not been modified, so that at steady-state <math display="inline">{\boldsymbol v}^n={\boldsymbol v}^{k-1}</math>, preserving the requirement that the steady-state be independent of the timestep <math display="inline">{\Delta t}</math>. The factor <math display="inline">\gamma </math> denotes the local ratio of the stability limit for explicit timestepping for the viscous terms versus the timestep chosen. Given that the advective and viscous timestep limits are proportional to:
+
<math display="inline">\theta </math> denotes the implicitness-factor for the viscous terms (<math display="inline">\theta=1</math>: 1st order, fully implicit, <math display="inline">\theta=0.5</math>: 2nd order, Crank-Nicholson). <math display="inline">\alpha ^i</math> are the standard low-storage Runge-Kutta coefficients <math display="inline">\alpha ^i=1/(k+1-i)</math>. The <math display="inline">k-1</math> stages of Eqn.(21a) may be seen as a predictor (or replacement) of <math display="inline">\vvec ^n</math> by <math display="inline">\vvec ^{k-1}</math>. The original right-hand side has not been modified, so that at steady-state <math display="inline">\vvec ^n=\vvec ^{k-1}</math>, preserving the requirement that the steady-state be independent of the timestep <math display="inline">\dt </math>. The factor <math display="inline">\gamma </math> denotes the local ratio of the stability limit for explicit timestepping for the viscous terms versus the timestep chosen. Given that the advective and viscous timestep limits are proportional to:
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 486: Line 492:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {\Delta t}_a \approx {h \over {|{\boldsymbol v}|}} ~~;~~    {\Delta t}_v \approx {{\rho h^2} \over \mu } ~~, </math>
+
| style="text-align: center;" | <math> \dt _a \approx {h \over {|\vvec |}} ~~;~~    \dt _v \approx {{\rho h^2} \over \mu } ~~, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
Line 498: Line 504:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> \gamma = {{{\Delta t}_v} \over {{\Delta t}_a}}    \approx {{\rho |{\boldsymbol v}| h }\over \mu } \approx Re_h  ~~, </math>
+
| style="text-align: center;" | <math> \gamma = {{\dt _v} \over {\dt _a}}    \approx {{\rho |\vvec | h }\over \mu } \approx Re_h  ~~, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
Line 517: Line 523:
 
In regions away from boundary layers, this factor is <math display="inline">O(1)</math>, implying that a high-order Runge-Kutta scheme is recovered. Conversely, for regions where <math display="inline">Re_h=O(0)</math>, the scheme reverts back to the usual 1-stage Crank-Nicholson scheme. Besides higher accuracy, an important benefit of explicit multistage advection schemes is the larger timestep one can employ. The increase in allowable timestep is roughly proportional to the number of stages used (and has been exploited extensively for compressible flow simulations <span id='citeF-49'></span>[[#cite-49|[49]]]). Given that for an incompressible solver of the projection type given by Eqns.(20-25) most of the CPU time is spent solving the pressure-Poisson system Eqn.(24), the speedup achieved is also roughly proportional to the number of stages used.  
 
In regions away from boundary layers, this factor is <math display="inline">O(1)</math>, implying that a high-order Runge-Kutta scheme is recovered. Conversely, for regions where <math display="inline">Re_h=O(0)</math>, the scheme reverts back to the usual 1-stage Crank-Nicholson scheme. Besides higher accuracy, an important benefit of explicit multistage advection schemes is the larger timestep one can employ. The increase in allowable timestep is roughly proportional to the number of stages used (and has been exploited extensively for compressible flow simulations <span id='citeF-49'></span>[[#cite-49|[49]]]). Given that for an incompressible solver of the projection type given by Eqns.(20-25) most of the CPU time is spent solving the pressure-Poisson system Eqn.(24), the speedup achieved is also roughly proportional to the number of stages used.  
  
At steady state, <math display="inline">{\boldsymbol v}^{*}={\boldsymbol v}^n={\boldsymbol v}^{n+1}</math> and the residuals of the pressure correction vanish, implying that the result does not depend on the timestep <math display="inline">{\Delta t}</math>.  
+
At steady state, <math display="inline">\vvec ^{*}=\vvec ^n=\vvec ^{n+1}</math> and the residuals of the pressure correction vanish, implying that the result does not depend on the timestep <math display="inline">\dt </math>.  
  
 
The spatial discretization of these equations is carried out via linear finite elements. The resulting matrix system is re-written as an edge-based solver, allowing the use of consistent numerical fluxes to stabilize the advection and divergence operators <span id='citeF-71'></span>[[#cite-71|[71]]].  
 
The spatial discretization of these equations is carried out via linear finite elements. The resulting matrix system is re-written as an edge-based solver, allowing the use of consistent numerical fluxes to stabilize the advection and divergence operators <span id='citeF-71'></span>[[#cite-71|[71]]].  
Line 523: Line 529:
 
The energy (temperature) equation (Eqn.(6.3)) is integrated in a manner similar to the advective-diffusive prediction (Eqn(21)), i.e. with an explicit, high order Runge-Kutta scheme for the advective parts and an implicit, 2nd order Crank-Nicholson scheme for the conductivity.
 
The energy (temperature) equation (Eqn.(6.3)) is integrated in a manner similar to the advective-diffusive prediction (Eqn(21)), i.e. with an explicit, high order Runge-Kutta scheme for the advective parts and an implicit, 2nd order Crank-Nicholson scheme for the conductivity.
  
===6.4 Numerical Integration of the Motion of Particles/Droplets===
+
===5.4 Numerical Integration of the Motion of Particles/Droplets===
  
 
The equations describing the position, velocity and temperature of a particle (Eqns. 15-19) may be formulated as a system of nonlinear Ordinary Differential Equations of the form:
 
The equations describing the position, velocity and temperature of a particle (Eqns. 15-19) may be formulated as a system of nonlinear Ordinary Differential Equations of the form:
Line 532: Line 538:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {{d{\boldsymbol u}_p} \over {dt}} = {\boldsymbol r}({\boldsymbol u}_p, {\boldsymbol x}, {\boldsymbol u}_f) ~~.                                                          </math>
+
| style="text-align: center;" | <math> {{d\uvec _p} \over {dt}} = \rvec (\uvec _p, \xvec , \uvec _f) ~~.                                                          </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
Line 544: Line 550:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {\boldsymbol u}^{n+i}_p = {\boldsymbol u}^n_p + \alpha ^i \Delta t \cdot    {\boldsymbol r}({\boldsymbol u}^{n+i-1}_p, {\boldsymbol x}^{n+i-1}, {\boldsymbol u}^{n+i-1}_f) ~~, ~~ i=1,k  ~~. </math>
+
| style="text-align: center;" | <math> \uvec ^{n+i}_p = \uvec ^n_p + \alpha ^i \Delta t \cdot    \rvec (\uvec ^{n+i-1}_p, \xvec ^{n+i-1}, \uvec ^{n+i-1}_f) ~~, ~~ i=1,k  ~~. </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
Line 570: Line 576:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {\boldsymbol M}\Delta {\boldsymbol u}= {\boldsymbol r}~~,    </math>
+
| style="text-align: center;" | <math> \Mmat \duvec = \rvec ~~,    </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
 
|}
 
|}
  
where <math display="inline">{\boldsymbol M}, \Delta {\boldsymbol u}</math> and <math display="inline">{\boldsymbol r}</math> denote, respectively, the consistent mass matrix, increment of the unknowns vector and right-hand side vector. Given the `host element' of each particle, i.e. the fluid mesh element that contains the particle, the forces and heat transferred to <math display="inline">{\boldsymbol r}</math> are added as follows:
+
where <math display="inline">\Mmat , \duvec </math> and <math display="inline">\rvec </math> denote, respectively, the consistent mass matrix, increment of the unknowns vector and right-hand side vector. Given the `host element' of each particle, i.e. the fluid mesh element that contains the particle, the forces and heat transferred to <math display="inline">\rvec </math> are added as follows:
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 582: Line 588:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {\boldsymbol r}^i_D = \sum _{el~surr~i} N^i({\boldsymbol x}_p) {\boldsymbol D}_p ~~.  </math>
+
| style="text-align: center;" | <math> \rvec ^i_D = \sum _{el~surr~i} N^i(\xvec _p) \Dvec _p ~~.  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
 
|}
 
|}
  
Here <math display="inline">N^i({\boldsymbol x}_p)</math> denotes the shape-function values of the host element for the point coordinates <math display="inline">{\boldsymbol x}_p</math>, and the sum extends over all elements that surround node <math display="inline">i</math>. As the sum of all shape-function values is unity at every point:
+
Here <math display="inline">N^i(\xvec _p)</math> denotes the shape-function values of the host element for the point coordinates <math display="inline">\xvec _p</math>, and the sum extends over all elements that surround node <math display="inline">i</math>. As the sum of all shape-function values is unity at every point:
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 594: Line 600:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> \sum N^i({\boldsymbol x}) = 1 ~~\forall {\boldsymbol x}~~,    </math>
+
| style="text-align: center;" | <math> \sum N^i(\xvec ) = 1 ~~\forall \xvec ~~,    </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
Line 608: Line 614:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
| style="text-align: center;" | <math> {\boldsymbol f}_p =  \rho _p {{\pi d^3}\over 6}              {{\left({\boldsymbol v}^{n+1}_p - {\boldsymbol v}^n_p \right)} \over {\Delta t}}                                            ~~,    </math>
+
| style="text-align: center;" | <math> \fvec _p =  \rho _p {{\pi d^3}\over 6}              {{\left(\vvec ^{n+1}_p - \vvec ^n_p \right)} \over {\Delta t}}                                            ~~,    </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
Line 625: Line 631:
 
These quantities are multiplied by the number of particles in a packet in order to obtain the final values transmitted to the fluid. Before going on, we summarize the basic steps required in order to update the particles one timestep:
 
These quantities are multiplied by the number of particles in a packet in order to obtain the final values transmitted to the fluid. Before going on, we summarize the basic steps required in order to update the particles one timestep:
  
* Initialize Fluid Source-Terms: <math display="inline">{\boldsymbol r}=0</math>
+
* Initialize Fluid Source-Terms: <math display="inline">\rvec=0</math>
 
* '''DO''' : For Each Particle:
 
* '''DO''' : For Each Particle:
 
* - DO: For Each Runge-Kutta Stage:
 
* - DO: For Each Runge-Kutta Stage:
* ~&nbsp;&nbsp;- Find Host Element of Particle: IELEM, <math display="inline">N^i({\boldsymbol x})</math>
+
* ~&nbsp;&nbsp;- Find Host Element of Particle: IELEM, <math display="inline">N^i(\xvec )</math>
 
* ~&nbsp;&nbsp;- Obtain Fluid Variables Required
 
* ~&nbsp;&nbsp;- Obtain Fluid Variables Required
 
* ~&nbsp;&nbsp;- Update Particle: Velocities, Position, Temperature, ...
 
* ~&nbsp;&nbsp;- Update Particle: Velocities, Position, Temperature, ...
Line 635: Line 641:
 
* '''ENDDO'''
 
* '''ENDDO'''
  
====6.4.1 Particle Parcels====
+
====5.4.1 Particle Parcels====
  
 
For a large number of very small particles, it becomes impossible to carry every individual particle in a simulation. The solution is to:
 
For a large number of very small particles, it becomes impossible to carry every individual particle in a simulation. The solution is to:
Line 645: Line 651:
 
Beyond a reasonable number of particles per element (typically <math display="inline">> 8</math>), this procedure produces accurate results without any deterioration in physical fidelity.
 
Beyond a reasonable number of particles per element (typically <math display="inline">> 8</math>), this procedure produces accurate results without any deterioration in physical fidelity.
  
====6.4.2 Other Numerics====
+
====5.4.2 Other Numerics====
  
 
In order to achieve a robust particle integrator, a number of additional precautions and algorithms need to be implemented. The most important of these are:
 
In order to achieve a robust particle integrator, a number of additional precautions and algorithms need to be implemented. The most important of these are:
Line 654: Line 660:
 
* Particle Tracking: A common feature of all particle-grid applications is that the particles do not move far between timesteps. This makes physical sense: if a particle jumped ten gridpoints during one timestep, it would have no chance to exchange information with the points along the way, leading to serious errors. Therefore, the assumption that the new host elements of the particles are in the vicinity of the current ones is a valid one. For this reason, the most efficient way to search for the new host elements is via the vectorized neighbour-to-neighbour algorithm described in <span id='citeF-59'></span><span id='citeF-71'></span>[[#cite-59|[59,71]]].
 
* Particle Tracking: A common feature of all particle-grid applications is that the particles do not move far between timesteps. This makes physical sense: if a particle jumped ten gridpoints during one timestep, it would have no chance to exchange information with the points along the way, leading to serious errors. Therefore, the assumption that the new host elements of the particles are in the vicinity of the current ones is a valid one. For this reason, the most efficient way to search for the new host elements is via the vectorized neighbour-to-neighbour algorithm described in <span id='citeF-59'></span><span id='citeF-71'></span>[[#cite-59|[59,71]]].
  
==7 Examples==
+
==6 Examples==
  
 
The techniques described above were implemented in FEFLO, a general-purpose computational fluid dynamics (CFD) code based on the following general principles:
 
The techniques described above were implemented in FEFLO, a general-purpose computational fluid dynamics (CFD) code based on the following general principles:
Line 736: Line 742:
 
The examples given show clearly the dangers of droplet- and aerosol- based infections in the built environment.
 
The examples given show clearly the dangers of droplet- and aerosol- based infections in the built environment.
  
===7.1 Sneezing in Transportation Security Agency (TSA) Queues===
+
===6.1 Sneezing in Transportation Security Agency (TSA) Queues===
  
 
One of the obvious vectors for viral contamination and spread are security and passport examination queues in airports. Air flow is moderate, passengers are in very close proximity, and in some airports queues wind back and forth in narrow lanes. Figure&nbsp;3a,b show the arrangement of pedestrians, as well as the discretization chosen. Note the smaller elements close to the bodies and in the region of interest between the two pedestrians in the middle row. This particular mesh had 12.74Mels. The distribution of particles and the absolute value of the velocity in the centerplane over time can be discerned from Figures&nbsp;3c-e. One can see that the large (red) particles follow a ballistic path and have some influence on the flow (e.g. at time <math display="inline">t=0.20</math>). This `ballistic phase' ends at about <math display="inline">t=1~sec</math>. The (green) particles of size <math display="inline">d=0.1~mm</math> are quickly stopped by the air, and then sink slowly towards the floor in close proximity to the individual sneezing. The even smaller (cyan, blue) particles rise with the cloud of warmer air exhaled by the sneezing individual, and disperse much further at later times.
 
One of the obvious vectors for viral contamination and spread are security and passport examination queues in airports. Air flow is moderate, passengers are in very close proximity, and in some airports queues wind back and forth in narrow lanes. Figure&nbsp;3a,b show the arrangement of pedestrians, as well as the discretization chosen. Note the smaller elements close to the bodies and in the region of interest between the two pedestrians in the middle row. This particular mesh had 12.74Mels. The distribution of particles and the absolute value of the velocity in the centerplane over time can be discerned from Figures&nbsp;3c-e. One can see that the large (red) particles follow a ballistic path and have some influence on the flow (e.g. at time <math display="inline">t=0.20</math>). This `ballistic phase' ends at about <math display="inline">t=1~sec</math>. The (green) particles of size <math display="inline">d=0.1~mm</math> are quickly stopped by the air, and then sink slowly towards the floor in close proximity to the individual sneezing. The even smaller (cyan, blue) particles rise with the cloud of warmer air exhaled by the sneezing individual, and disperse much further at later times.
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/tsa_queue/wire_far.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3a&nbsp;&nbsp;TSA Queue: Arrangement of Pedestrians and Surface Mesh</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3a&nbsp;&nbsp;TSA Queue: Arrangement of Pedestrians and Surface Mesh</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/tsa_queue/wire_close.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3b&nbsp;&nbsp;TSA Queue: Surface Mesh and Cut Plane</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3b&nbsp;&nbsp;TSA Queue: Surface Mesh and Cut Plane</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/tsa_queue/que2_0.02.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3c&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.02~sec</math></div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3c&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.02~sec</math></div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/tsa_queue/que2_0.10.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3d&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.10~sec</math></div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3d&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.10~sec</math></div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/tsa_queue/que2_0.20.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3e&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.20~sec</math></div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3e&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.20~sec</math></div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/tsa_queue/que2_0.40.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3f&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.40~sec</math></div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3f&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.40~sec</math></div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/tsa_queue/que2_0.80.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3g&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.80~sec</math></div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3g&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=0.80~sec</math></div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/tsa_queue/que2_1.60.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3h&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=1.60~sec</math></div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 3h&nbsp;&nbsp;TSA Queue: Particle Distribution at <math>t=1.60~sec</math></div>
  
===7.2 Sneezing in a Generic Hospital Room===
+
===6.2 Sneezing in a Generic Hospital Room===
  
 
This case considers a typical hospital room. Of interest here was the dispersion of particles in the '''first minute''' after coughing, in particular the reach into neighbouring halls and the amount of `negative pressure' needed to keep all contaminants in the room. Figure&nbsp;4a shows the arrangement of the room, with patient and caregiver clearly visible. This particular mesh had 2.25Mels. The distribution of particles over time can be discerned from Figures&nbsp;4c-l. As before, one can see that the large (red) particles follow a ballistic path. This `ballistic phase' ends at about <math display="inline">t=1~sec</math>. The (green) particles of size <math display="inline">d=0.1~mm</math> are quickly stopped by the air, and then sink slowly towards the patient. The even smaller (cyan, blue) particles rise with the cloud of warmer air exhaled by the sneezing individual, and disperse much further at later times, covering almost the entire room. The velocity distribution in the room may be inferred from Figure&nbsp;4m.
 
This case considers a typical hospital room. Of interest here was the dispersion of particles in the '''first minute''' after coughing, in particular the reach into neighbouring halls and the amount of `negative pressure' needed to keep all contaminants in the room. Figure&nbsp;4a shows the arrangement of the room, with patient and caregiver clearly visible. This particular mesh had 2.25Mels. The distribution of particles over time can be discerned from Figures&nbsp;4c-l. As before, one can see that the large (red) particles follow a ballistic path. This `ballistic phase' ends at about <math display="inline">t=1~sec</math>. The (green) particles of size <math display="inline">d=0.1~mm</math> are quickly stopped by the air, and then sink slowly towards the patient. The even smaller (cyan, blue) particles rise with the cloud of warmer air exhaled by the sneezing individual, and disperse much further at later times, covering almost the entire room. The velocity distribution in the room may be inferred from Figure&nbsp;4m.
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/wire.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4a&nbsp;&nbsp;Hospital Room: Surface Mesh</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4a&nbsp;&nbsp;Hospital Room: Surface Mesh</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0000.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4b&nbsp;&nbsp;Hospital Room: Particle Distribution at t=0.0 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4b&nbsp;&nbsp;Hospital Room: Particle Distribution at t=0.0 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0001.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4c&nbsp;&nbsp;Hospital Room: Particle Distribution at t=0.2 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4c&nbsp;&nbsp;Hospital Room: Particle Distribution at t=0.2 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0002.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4d&nbsp;&nbsp;Hospital Room: Particle Distribution at t=0.4 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4d&nbsp;&nbsp;Hospital Room: Particle Distribution at t=0.4 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0003.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4e&nbsp;&nbsp;Hospital Room: Particle Distribution at t=0.6 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4e&nbsp;&nbsp;Hospital Room: Particle Distribution at t=0.6 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0005.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4f&nbsp;&nbsp;Hospital Room: Particle Distribution at t=1.0 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4f&nbsp;&nbsp;Hospital Room: Particle Distribution at t=1.0 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0010.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4g&nbsp;&nbsp;Hospital Room: Particle Distribution at t=2.0 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4g&nbsp;&nbsp;Hospital Room: Particle Distribution at t=2.0 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0020.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4h&nbsp;&nbsp;Hospital Room: Particle Distribution at t=4.0 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4h&nbsp;&nbsp;Hospital Room: Particle Distribution at t=4.0 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0032.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4i&nbsp;&nbsp;Hospital Room: Particle Distribution at t=10.0 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4i&nbsp;&nbsp;Hospital Room: Particle Distribution at t=10.0 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0052.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4j&nbsp;&nbsp;Hospital Room: Particle Distribution at t=20.0 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4j&nbsp;&nbsp;Hospital Room: Particle Distribution at t=20.0 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0072.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4k&nbsp;&nbsp;Hospital Room: Particle Distribution at t=30.0 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4k&nbsp;&nbsp;Hospital Room: Particle Distribution at t=30.0 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/just_parts.0112.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4l&nbsp;&nbsp;Hospital Room: Particle Distribution at t=50.0 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4l&nbsp;&nbsp;Hospital Room: Particle Distribution at t=50.0 sec</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/hospital_room/velo_and_parts.0132.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4m&nbsp;&nbsp;Hospital Room: Surface Velocities at t=60.0 sec</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 4m&nbsp;&nbsp;Hospital Room: Surface Velocities at t=60.0 sec</div>
  
==8 Reopening After the Crisis==
+
==7 Reopening After the Crisis==
  
 
A lingering question facing all levels of society is how and when to reopen facilities where people congregate in close proximity. One key technology that would allow opening is testing and sensing. We consider sensing in the sequel. Several vendors have announced measuring devices for Covid-19 in the next half year. Given that these sensors are expensive, and that a hospital or university many need hundreds of these, the question becomes how best to deploy them. In other words: given an arbitrary number of contamination or infection scenarios, which is the minimum number of sensors needed to detect them, and where should they be placed&nbsp;? A partial answer to this non-trivial question was given in <span id='citeF-67'></span><span id='citeF-110'></span>[[#cite-67|[67,110]]]. If we assume a given number of sensors, every contaminant/infection scenario (location and amount of release, flow conditions, etc.) will lead to a sensor input. The data recorded from all the possible release scenarios at all possible sensor locations allows the identification of the best or optimal sensor locations. Clearly, if only one sensor is to be placed, it should be at the location that recorded the highest number of releases. This argument can be used recursively by removing from further consideration all releases already recorded by sensors previously placed. The procedure is repeated recursively until no undetected release cases are left, or the available sensors have been exhausted.   
 
A lingering question facing all levels of society is how and when to reopen facilities where people congregate in close proximity. One key technology that would allow opening is testing and sensing. We consider sensing in the sequel. Several vendors have announced measuring devices for Covid-19 in the next half year. Given that these sensors are expensive, and that a hospital or university many need hundreds of these, the question becomes how best to deploy them. In other words: given an arbitrary number of contamination or infection scenarios, which is the minimum number of sensors needed to detect them, and where should they be placed&nbsp;? A partial answer to this non-trivial question was given in <span id='citeF-67'></span><span id='citeF-110'></span>[[#cite-67|[67,110]]]. If we assume a given number of sensors, every contaminant/infection scenario (location and amount of release, flow conditions, etc.) will lead to a sensor input. The data recorded from all the possible release scenarios at all possible sensor locations allows the identification of the best or optimal sensor locations. Clearly, if only one sensor is to be placed, it should be at the location that recorded the highest number of releases. This argument can be used recursively by removing from further consideration all releases already recorded by sensors previously placed. The procedure is repeated recursively until no undetected release cases are left, or the available sensors have been exhausted.   
Line 792: Line 798:
 
See <span id='citeF-21'></span><span id='citeF-45'></span>[[#cite-21|[21,45]]] for an in-depth analysis of robust sensor placement under uncertainty.
 
See <span id='citeF-21'></span><span id='citeF-45'></span>[[#cite-21|[21,45]]] for an in-depth analysis of robust sensor placement under uncertainty.
  
===8.1 Hospital Room===
+
===7.1 Hospital Room===
  
 
This case considers the same hospital room as shown before. The boundary conditions determining the flow are assumed as steady, with air entering the room through vents&nbsp;1-3 and exiting the room through the bathroom exhaust or the door. Figures&nbsp;5a-c show the outlay of the room, average velocities and the `age of air' after 5&nbsp;minutes. Note the high values for the age of air in the corners and the back of the room. This particular mesh had 2.2Mels. Four contaminant release scenarios were considered: cases&nbsp;1-3 assumed contaminant coming in through each of the vents (separately) during the first minute, while case&nbsp;4 assumed virus production from the patient for a period of 10&nbsp;seconds. The case was run for 5&nbsp;minutes of real time, and the contaminant concentration was measured on all walls/ceilings. The maximum concentrations measured have been summarized in Figure&nbsp;5d. Note the different areas covered depending on the release scenario. It was assumed that sensors should only be allowed above a certain height, and should be located on a wall or the ceiling. Table&nbsp;4 summarizes the points that measured data above a set threshold. As one can see, none of the possible sensor locations is able to measure/detect all 4&nbsp;cases, and many possible sensor locations do not detect even a single case. There are many possible pairs of sensors that can detect all 4&nbsp;cases. The pair selected is the one that achieves the highest relative measurement values, and is shown in Figure&nbsp;5e. Note that this makes good sense: one sensor close the HVAC exits, and one close to the patient.
 
This case considers the same hospital room as shown before. The boundary conditions determining the flow are assumed as steady, with air entering the room through vents&nbsp;1-3 and exiting the room through the bathroom exhaust or the door. Figures&nbsp;5a-c show the outlay of the room, average velocities and the `age of air' after 5&nbsp;minutes. Note the high values for the age of air in the corners and the back of the room. This particular mesh had 2.2Mels. Four contaminant release scenarios were considered: cases&nbsp;1-3 assumed contaminant coming in through each of the vents (separately) during the first minute, while case&nbsp;4 assumed virus production from the patient for a period of 10&nbsp;seconds. The case was run for 5&nbsp;minutes of real time, and the contaminant concentration was measured on all walls/ceilings. The maximum concentrations measured have been summarized in Figure&nbsp;5d. Note the different areas covered depending on the release scenario. It was assumed that sensors should only be allowed above a certain height, and should be located on a wall or the ceiling. Table&nbsp;4 summarizes the points that measured data above a set threshold. As one can see, none of the possible sensor locations is able to measure/detect all 4&nbsp;cases, and many possible sensor locations do not detect even a single case. There are many possible pairs of sensors that can detect all 4&nbsp;cases. The pair selected is the one that achieves the highest relative measurement values, and is shown in Figure&nbsp;5e. Note that this makes good sense: one sensor close the HVAC exits, and one close to the patient.
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/sensor_placement/room_outlay.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5a&nbsp;&nbsp;Hospital Room: Outlay of Room and Boundary Conditions</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5a&nbsp;&nbsp;Hospital Room: Outlay of Room and Boundary Conditions</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/sensor_placement/absvelo.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5b&nbsp;&nbsp;Hospital Room: Average Velocities (5&nbsp;mins)</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5b&nbsp;&nbsp;Hospital Room: Average Velocities (5&nbsp;mins)</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/sensor_placement/ageofair.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5c&nbsp;&nbsp;Hospital Room: Age of Air (5&nbsp;mins)</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5c&nbsp;&nbsp;Hospital Room: Age of Air (5&nbsp;mins)</div>
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/sensor_placement/maxvalues5mins.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5d&nbsp;&nbsp;Hospital Room: Maximum Contaminant Concentration Over 5 Minutes</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5d&nbsp;&nbsp;Hospital Room: Maximum Contaminant Concentration Over 5 Minutes</div>
  
  
Line 828: Line 834:
 
|}
 
|}
  
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> figures/sensor_placement/optimal_place.eps </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5e&nbsp;&nbsp;Hospital Room: Optimal Sensor Locations</div>
+
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">   </div> <div class="center" style="width: auto; margin-left: auto; margin-right: auto;">Figure 5e&nbsp;&nbsp;Hospital Room: Optimal Sensor Locations</div>
  
==9 Conclusions and Outlook==
+
==8 Conclusions and Outlook==
  
 
The present paper has summarized some of the mechanical characteristics of virus contaminants and the transmission via droplets and aerosols. The ordinary and partial differential equations describing the physics of these processes with high fidelity were given, as well as appropriate numerical schemes to solve them. Several examples taken from recent evaluations of the built environment were given, as well as the optimal placement of sensors.
 
The present paper has summarized some of the mechanical characteristics of virus contaminants and the transmission via droplets and aerosols. The ordinary and partial differential equations describing the physics of these processes with high fidelity were given, as well as appropriate numerical schemes to solve them. Several examples taken from recent evaluations of the built environment were given, as well as the optimal placement of sensors.

Revision as of 12:24, 30 June 2020


Abstract

A summary is given of the mechanical characteristics of virus contaminants and the transmission via droplets and aerosols. The ordinary and partial differential equations describing the physics of these processes with high fidelity are presented, as well as appropriate numerical schemes to solve them. Several examples taken from recent evaluations of the built environment are shown, as well as the optimal placement of sensors.

keywords Covid-19 Particle Methods Finite Elements Computational Fluid Dynamics

1 The Covid-19 Crisis

Starting in Wuhan, China, in the fall of 2019, the Covid-19 pandemic has claimed and will continue to claim millions of infected patients and hundreds of thousands of deaths. The lockdowns that followed its outbreak have led to mass unemployment, stalled economic development and loss of productivity that will take years to recover. Some changes in habits and lifestyles may be permanent: in the future, working from home or in a `socially distanced manner' may be the prevalent modus operandi for large segments of society.

The present paper gives a short description of computational techniques that can elucidate the flow and propagation of viruses or other contaminants in built environments in order to mitigate or avoid infections.

2 Virus Infection

Before addressing the requirements for the numerical simulation of virus propagation a brief description of virus propagation and lifetime are given. Covid-19 is one of many corona-viruses. The virus is usually present in the air or some surface, and makes its way into the body either via inhalation (nose, mouth), ingestion (mouth) or attachment (eyes, hands, clothes). In many cases the victim inadvertedly touches an infected surface or viruses are deposited on its hands, and then the hands touch either the nose, the eyes or the mouth, thus allowing the virus to enter the body.

An open question of great importance for all that will follow is how many viruses it takes to overwhelm the body's natural defense mechanism and trigger an infection. This number, which is sometimes called the viral load or the infectious dose will depend on numerous factors, among them the state of immune defenses of the individual, the timing of viral entry (all at once, piece by piece), and the amount of hair and mucous in the nasal vessels. In principle, a single organism in a favourable environment may replicate sufficiently to cause disease [98]. Data from research performed on biological warfare agents [37] suggests that both bacteria and viruses can produce disease with as few as 1-100 organisms (e.g. brucellosis 10-100, Q fever 1-10, tularaemia 10-50, smallpox 10-100, viral haemorrhagic fevers 1-10 organisms, tuberculosis 1). Compare these numbers and consider that as many as 3,000 organisms can be produced by talking for 5 minutes or a single cough, with sneezing producing many more [80,55,102,87,108]. Figure 1, reproduced from [102], shows a typical number and size distribution.

Counts of Particles of Various Diameters in Air Expelled by (90) Coughs [80]
Figure 1: Counts of Particles of Various Diameters in Air Expelled by (90) Coughs [80]

3 Virus Lifetime Outside the Body

Current evidence points to lifetimes outside the body that can range from 1-2 hours in air to several days on particular surfaces [51,106]. There has also been some documentation of lifetime variation depending on humidity.

4 Virus Transmission

4.1 Human Sneezing and Coughing

In the sequel, we consider human sneezing and coughing as the main conduits of virus transmission. Clearly, breathing and talking will lead to the exhalation of air, and, consequently the exhalation of viruses for infected victims [3,4]. However, it stands to reason that the size and amount of particles released - and hence the amount of viruses in them - is much higher and much more concentrated when sneezing or coughing [36,102,50,56,3,4].

The velocity of air at a person's mouth during sneezing and coughing has been a source of heated debate, particularly in the media. The experimental evidence points to exit velocities of the order of 2-14 m/sec [40,41,100,101]. A typical amount and size of particles can be seen in Figure 1.

4.2 Sink Velocities

If, for the sake of argument, we consider Stoke's law for the drag of spherical particles, valid below Reynolds numbers of , the terminal sink velocity (also known as the settling velocity) of particles will be given by:

(1)

where denote the density of the particles (essentially water in the present case), density of the gas (air), gravity, dynamic viscosity of the gas and diameter of the particle respectively. The equivalent Reynolds' number is:

(2)

With this yields a limiting diameter for of

(3)

i.e. approximately 1/10th of a millimeter in diameter - a particle size that would still be perceived by the human eye. The corresponding sink velocity is given by:

(4)

with in , i.e. for

(5)

Note the quadratic dependency of the sink velocity with diameter. Table 1 lists the sink velocities for water droplets of different diameters in air. One can see that below diameters of the sink velocity is very low, implying that these particles remain in and move with the air for considerable time (and possibly distances).


Table. 1 Sink Velocities and Reynolds Number For Water Particles in Air
Diameter [mm] sink velocity [m/sec] Re
1.00E-00 3.01E+01 1.99E+03
1.00E-01 3.01E-01 1.99E+00
1.00E-02 3.01E-03 1.99E-03
1.00E-03 3.01E-05 1.99E-06
1.00E-04 3.01E-07 1.99E-09

4.3 Evaporation

Depending on the relative humidity and the temperature of the ambient air, the smaller particles can evaporate in milliseconds. However, as the mucous and saliva evaporate, they build a gel-like structure that surrounds the virus, allowing it to survive. This implies that extremely small particles with possible viruses will remain infectious for extended periods of times - up to an hour according to some studies.

An important question is then whether a particle/droplet will first reach the ground or evaporate. Figure 2, taken from [109], shows that below 120  the particles evaporate before falling 2 m (i.e. reaching the ground).

Figure 2   Evaporation Time and Falling Time of Droplets of Varying Diameter
() From [109]

4.4 Viral Load

A central question that requires an answer is then: how many viruses are in these small particles ? An approximate answer may be obtained from the experiments that are being carried out on animals to trace and monitor infections. For ferrets [53] have been used to infect via intranasal swabs, while for mice [104] seem to suffice. Viral titers can vary a lot, but one may assume on the order of viruses/ml for a nasopharyngeal swab [53,104]. Table 2 lists the number of viruses per droplet and the number of droplets needed to contain just 1 virus. Note that while for a droplet with a diameter of 1 mm one can expect viruses, only every 2,000th particle of diameter 10  does contain a single virus.


Table. 2 Estimated Number of Viruses For Different Particle Diameters
Droplet Diameter [mm] Volume [mm] Viruses/droplet Droplets Needed for 1 Virus
1.00E+00 5.24E-01 5.24E+02 1.00E+00
1.00E-01 5.24E-04 5.24E-01 1.91E+00
1.00E-02 5.24E-07 5.24E-04 1.91E+03
1.00E-03 5.24E-10 5.24E-07 1.91E+06

Similar numbers are seen in field studies as well. The size of viruses varies from 0.02-0.3 , while the size of bacteria varies from 0.5-10 . The influenza virus RNA detected by quantitative polymerase chain reaction in human exhaled breath suggests that it may be contained in fine particles generated during tidal breathing and not only coughs [36,55,56,87]. Influenza RNA and Mycobacterium tuberculosis have been reported in particles that range in size from 0.5-4.0 ([36,55,56,87] and references cited therein).

5 Physical Modeling of Aerosol Propagation

When solving the two-phase equations, the air, as a continuum, is best represented by a set of partial differential equations (the Navier-Stokes equations) that are numerically solved on a mesh. Thus, the gas characteristics are calculated at the mesh points within the flowfield. However, as the droplets/particles may be relatively sparse in the flowfield, they can be modeled by either:

  • A continuum description, i.e. in the same manner as the fluid flow, or
  • A particle (or Lagrangian) description, where individual particles (or groups of particles) are monitored and tracked in the flow.

Although the continuum (so-called multi-fluid) method has been used for some applications, the inherent assumptions of the continuum approach lead to several disadvantages which may be countered with a Lagrangian treatment for dilute flows. The continuum assumption cannot robustly account for local differences in particle characteristics, particularly if the particles are polydispersed. In addition, the only boundary conditions that can be considered in a straightforward manner are slipping and sticking, whereas reflection boundary conditions, such as specular and diffuse reflection, may be additionally considered with a Lagrangian approach. Numerical diffusion of the particle density is eliminated by employing Lagrangian particles due to their pointwise spatial accuracy.

While a Lagrangian approach offers many potential advantages, this method also creates problems that need to be addressed. For instance, large numbers of particles may cause a Lagrangian analysis to be memory intensive. This problem is circumvented by treating parcels of particles, i.e. doing the detailed analysis for one particle and then applying the effect to many. In addition, continuous mapping and remapping of particles to their respective elements may increase computational requirements, particularly for unstructured grids.

5.1 Equations Describing the Motion of the Air

As seen from the experimental evidence, the velocities of air encountered during coughing and sneezing never exceed a Mach-number of . Therefore, the air may be assumed as a Newtonian, incompressible liquid, where buoyancy effects are modeled via the Boussinesq approximation. The equations describing the conservation of momentum, mass and energy for incompressible, Newtonian flows may be written as

Failed to parse (unknown function "\vvec"): \rho \vvec _{,t} + \rho \vvec \cdot \Grd \vvec + \Grd p = \nabla \cdot \mu \nabla \vvec + \rho \gvec + \beta \rho \gvec (T - T_0) + \svec _v ~~,
(6.1)
Failed to parse (unknown function "\Div"): \Div \vvec = 0 ~~.
(6.2)
Failed to parse (unknown function "\vvec"): \rho c_p T_{,t} + \rho c_p \vvec \cdot \Grd T = \nabla \cdot k \nabla T + s_e ~~,
(6.3)

Here Failed to parse (unknown function "\vvec"): {\textstyle \rho , \vvec , p, \mu , \gvec , \beta , T, T_0, c_p, k}

denote the density, velocity vector, pressure, viscosity, gravity vector, coefficient of thermal expansion, temperature, reference temperature, specific heat coefficient and conductivity respectively, and Failed to parse (unknown function "\svec"): {\textstyle \svec _v, s_e}
momentum and energy source terms (e.g. due to particles or external forces/heat sources). For turbulent flows both the viscosity and the conductivity are obtained either from additional equations or directly via a large eddy simulation (LES) assumption through monotonicity induced LES (MILES) [20,38,39,47].

5.2 Equations Describing the Motion of Particles/Droplets

In order to describe the interaction of particles/droplets with the flow, the mass, forces and energy/work exchanged between the flowfield and the particles must be defined. As before, we denote for fluid (air) by and the density, pressure, temperature, conductivity, velocity in direction , viscosity, and the specific heat at constant pressure. For the particles, we denote by and the density, temperature, velocity in direction , equivalent diameter, and heat transferred per unit volume. In what follows, we will refer to droplet and particles, collectively as particles.

Making the classical assumptions that the particles may be represented by an equivalent sphere of diameter , the drag forces Failed to parse (unknown function "\Dvec"): {\textstyle \Dvec }

acting on the particles will be due to the difference of fluid and particle velocity:
Failed to parse (unknown function "\Dvec"): \Dvec = {{\pi d^2} \over 4} \cdot c_D \cdot { 1 \over 2} \rho | \vvec - \vvec _p | ( \vvec - \vvec _p ) ~~.
(7)

The drag coefficient is obtained empirically from the Reynolds-number :

Failed to parse (unknown function "\vvec"): Re = {{\rho | \vvec - \vvec _p | d } \over { \mu }}
(8)

as (see, e.g. [95]):

(9)

The lower bound of is required to obtain the proper limit for the Euler equations, when . The heat transferred between the particles and the fluid is given by

(10)

where is the film coefficient and the radiation coefficient. For the class of problems considered here, the particle temperature and kinetic energy are such that the radiation coefficient may be ignored. The film coefficient is obtained from the Nusselt-Number :

(11)

where is the Prandtl-number of the gas

(12)

as

(13)

Having established the forces and heat flux, the particle motion and temperature are obtained from Newton's law and the first law of thermodynamics. For the particle velocities, we have:

Failed to parse (unknown function "\vvec"): \rho _p {{\pi d^3} \over 6 } \cdot {{ d\vvec _p} \over {dt}} = \Dvec + \rho _p {{\pi d^3} \over 6 } \gvec ~~.
(14)

This implies that:

Failed to parse (unknown function "\vvec"): {{ d\vvec _p} \over {dt}} = {{3 \rho } \over {4 \rho _p d}} \cdot c_d | \vvec - \vvec _p | ( \vvec - \vvec _p ) + \gvec = \alpha _v | \vvec - \vvec _p | ( \vvec - \vvec _p ) + \gvec ~~.
(15)

where . The particle positions are obtained from:

Failed to parse (unknown function "\xvec"): {{ d\xvec _p} \over {dt}} = \vvec _p ~~.
(16)

The temperature change in a particle is given by:

(17)

which may be expressed as:

(18)

with . Equations (15, 16, 18) may be formulated as a system of Ordinary Differential Equations (ODEs) of the form:

Failed to parse (unknown function "\uvec"): {{d\uvec _p} \over {dt}} = \rvec (\uvec _p, \xvec , \uvec _f) ~~,
(19)

where Failed to parse (unknown function "\uvec"): {\textstyle \uvec _p, \xvec , \uvec _f}

denote the particle unknowns, the position of the particle and the fluid unknowns at the position of the particle.

5.3 Numerical Integration of the Motion of the Air

The last six decades have seen a large number of schemes that may be used to solve numerically the incompressible Navier-Stokes equations given by Eqns.(6.1-6.3). In the present case, the following design criteria were implemented:

  • Spatial discretization using unstructured grids (in order to allow for arbitrary geometries and adaptive refinement);
  • Spatial approximation of unknowns with simple linear finite elements (in order to have a simple input/output and code structure);
  • Edge-based data structures (for reduced access to memory and indirect addressing);
  • Temporal approximation using implicit integration of viscous terms and pressure (the interesting scales are the ones associated with advection);
  • Temporal approximation using explicit, high-order integration of advective terms;
  • Low-storage, iterative solvers for the resulting systems of equations (in order to solve large 3-D problems); and
  • Steady results that are independent from the timestep chosen (in order to have confidence in convergence studies).

The resulting discretization in time is given by the following projection scheme [68,71]:

  • Advective-Diffusive Prediction: Failed to parse (unknown function "\vvec"): {\textstyle \vvec ^n, p^n \rightarrow \vvec ^{*}}


Failed to parse (unknown function "\svec"): \svec' = - \Grd p^n + \rho \gvec + \beta \rho \gvec (T^n - T_0) + \svec _v ~~,
(20)
Failed to parse (unknown function "\vvec"): \vvec ^i = \vvec ^n + \alpha ^i \gamma \dt \left( - \vvec ^{i-1} \cdot \Grd \vvec ^{i-1} \nabla \cdot \mu \nabla \vvec ^{i-1} + \svec' \right) ~~; ~~i=1,k-1~~;
(21a)
Failed to parse (unknown function "\dt"): \left[{ 1 \over \dt } - \theta \nabla \cdot \mu \nabla \right] \left(\vvec ^{k} - \vvec ^n \right) + \vvec ^{k-1} \cdot \Grd \vvec ^{k-1} = \nabla \cdot \mu \nabla \vvec ^{k-1} + \svec' ~~.
(21b)
  • Pressure Correction:
Failed to parse (unknown function "\Div"): \Div \vvec ^{n+1} = 0 ~~;
(22)
Failed to parse (unknown function "\vvec"): {{ \vvec ^{n+1} - \vvec ^{*} }\over \dt } + \Grd ( p^{n+1} - p^n ) = 0 ~~;
(23)
  • which results in
Failed to parse (unknown function "\Div"): \nabla ^2 ( p^{n+1} - p^n ) = {{\Div \vvec ^{*} }\over \dt } ~~;
(24)
  • Velocity Correction: Failed to parse (unknown function "\vvec"): {\textstyle \vvec ^{*} \rightarrow \vvec ^{n+1}}


Failed to parse (unknown function "\vvec"): \vvec ^{n+1} = \vvec ^{*} - \dt \Grd ( p^{n+1} - p^n ) ~~.
(25)

denotes the implicitness-factor for the viscous terms (: 1st order, fully implicit, : 2nd order, Crank-Nicholson). are the standard low-storage Runge-Kutta coefficients . The stages of Eqn.(21a) may be seen as a predictor (or replacement) of Failed to parse (unknown function "\vvec"): {\textstyle \vvec ^n}

by Failed to parse (unknown function "\vvec"): {\textstyle \vvec ^{k-1}}

. The original right-hand side has not been modified, so that at steady-state Failed to parse (unknown function "\vvec"): {\textstyle \vvec ^n=\vvec ^{k-1}} , preserving the requirement that the steady-state be independent of the timestep Failed to parse (unknown function "\dt"): {\textstyle \dt } . The factor denotes the local ratio of the stability limit for explicit timestepping for the viscous terms versus the timestep chosen. Given that the advective and viscous timestep limits are proportional to:

Failed to parse (unknown function "\dt"): \dt _a \approx {h \over {|\vvec |}} ~~;~~ \dt _v \approx {{\rho h^2} \over \mu } ~~,
(26)

we immediately obtain

Failed to parse (unknown function "\dt"): \gamma = {{\dt _v} \over {\dt _a}} \approx {{\rho |\vvec | h }\over \mu } \approx Re_h ~~,
(27)

or, in its final form:

(28)

In regions away from boundary layers, this factor is , implying that a high-order Runge-Kutta scheme is recovered. Conversely, for regions where , the scheme reverts back to the usual 1-stage Crank-Nicholson scheme. Besides higher accuracy, an important benefit of explicit multistage advection schemes is the larger timestep one can employ. The increase in allowable timestep is roughly proportional to the number of stages used (and has been exploited extensively for compressible flow simulations [49]). Given that for an incompressible solver of the projection type given by Eqns.(20-25) most of the CPU time is spent solving the pressure-Poisson system Eqn.(24), the speedup achieved is also roughly proportional to the number of stages used.

At steady state, Failed to parse (unknown function "\vvec"): {\textstyle \vvec ^{*}=\vvec ^n=\vvec ^{n+1}}

and the residuals of the pressure correction vanish, implying that the result does not depend on the timestep Failed to parse (unknown function "\dt"): {\textstyle \dt }

.

The spatial discretization of these equations is carried out via linear finite elements. The resulting matrix system is re-written as an edge-based solver, allowing the use of consistent numerical fluxes to stabilize the advection and divergence operators [71].

The energy (temperature) equation (Eqn.(6.3)) is integrated in a manner similar to the advective-diffusive prediction (Eqn(21)), i.e. with an explicit, high order Runge-Kutta scheme for the advective parts and an implicit, 2nd order Crank-Nicholson scheme for the conductivity.

5.4 Numerical Integration of the Motion of Particles/Droplets

The equations describing the position, velocity and temperature of a particle (Eqns. 15-19) may be formulated as a system of nonlinear Ordinary Differential Equations of the form:

Failed to parse (unknown function "\uvec"): {{d\uvec _p} \over {dt}} = \rvec (\uvec _p, \xvec , \uvec _f) ~~.
(29)

They can be integrated numerically in a variety of ways. Due to its speed, low memory requirements and simplicity, we have chosen the following k-step low-storage Runge-Kutta procedure to integrate them:

Failed to parse (unknown function "\uvec"): \uvec ^{n+i}_p = \uvec ^n_p + \alpha ^i \Delta t \cdot \rvec (\uvec ^{n+i-1}_p, \xvec ^{n+i-1}, \uvec ^{n+i-1}_f) ~~, ~~ i=1,k ~~.
(30)

For linear ODEs the choice

(31)

leads to a scheme that is -th order accurate in time. Note that in each step the location of the particle with respect to the fluid mesh needs to be updated in order to obtain the proper values for the fluid unknowns. The default number of stages used is . This would seem unnecessarily high, given that the flow solver is of second-order accuracy, and that the particles are integrated separately from the flow solver before the next (flow) timestep, i.e. in a staggered manner. However, it was found that the 4-stage particle integration preserves very well the motion in vortical structures and leads to less `wall sliding' close to the boundaries of the domain [78]. The stability/ accuracy of the particle integrator should not be a problem as the particle motion will always be slower than the maximum wave speed of the fluid (fluid velocity).

The transfer of forces and heat flux between the fluid and the particles must be accomplished in a conservative way, i.e. whatever is added to the fluid must be subtracted from the particles and vice-versa. The finite element discretization of the fluid equations will lead to a system of ODE's of the form:

Failed to parse (unknown function "\Mmat"): \Mmat \duvec = \rvec ~~,
(32)

where Failed to parse (unknown function "\Mmat"): {\textstyle \Mmat , \duvec }

and Failed to parse (unknown function "\rvec"): {\textstyle \rvec }
denote, respectively, the consistent mass matrix, increment of the unknowns vector and right-hand side vector. Given the `host element' of each particle, i.e. the fluid mesh element that contains the particle, the forces and heat transferred to Failed to parse (unknown function "\rvec"): {\textstyle \rvec }
are added as follows:
Failed to parse (unknown function "\rvec"): \rvec ^i_D = \sum _{el~surr~i} N^i(\xvec _p) \Dvec _p ~~.
(33)

Here Failed to parse (unknown function "\xvec"): {\textstyle N^i(\xvec _p)}

denotes the shape-function values of the host element for the point coordinates Failed to parse (unknown function "\xvec"): {\textstyle \xvec _p}

, and the sum extends over all elements that surround node . As the sum of all shape-function values is unity at every point:

Failed to parse (unknown function "\xvec"): \sum N^i(\xvec ) = 1 ~~\forall \xvec ~~,
(34)

this procedure is strictly conservative.

From Eqns.(15-18) and their equivalent numerical integration via Eqn.(30), the change in momentum and energy for one particle is given by:

Failed to parse (unknown function "\fvec"): \fvec _p = \rho _p {{\pi d^3}\over 6} {{\left(\vvec ^{n+1}_p - \vvec ^n_p \right)} \over {\Delta t}} ~~,
(35)
(36)

These quantities are multiplied by the number of particles in a packet in order to obtain the final values transmitted to the fluid. Before going on, we summarize the basic steps required in order to update the particles one timestep:

  • Initialize Fluid Source-Terms: Failed to parse (unknown function "\rvec"): {\textstyle \rvec=0}
  • DO : For Each Particle:
  • - DO: For Each Runge-Kutta Stage:
  • ~  - Find Host Element of Particle: IELEM, Failed to parse (unknown function "\xvec"): {\textstyle N^i(\xvec )}
  • ~  - Obtain Fluid Variables Required
  • ~  - Update Particle: Velocities, Position, Temperature, ...
  • - ENDDO
  • - Transfer Loads to Element Nodes
  • ENDDO

5.4.1 Particle Parcels

For a large number of very small particles, it becomes impossible to carry every individual particle in a simulation. The solution is to:

  • Agglomerate the particles into so-called packets of particles;
  • Integrate the governing equations for one individual particle; and
  • Transfer back to the fluid times the effect of one particle.

Beyond a reasonable number of particles per element (typically ), this procedure produces accurate results without any deterioration in physical fidelity.

5.4.2 Other Numerics

In order to achieve a robust particle integrator, a number of additional precautions and algorithms need to be implemented. The most important of these are:

  • Agglomeration/Subdivision of Particle Parcels: As the fluid mesh may be adaptively refined and coarsened in time, or the particle traverses elements of different sizes, it may be important to adapt the parcel concentrations as well. This is necessary to ensure that there is sufficient parcel representation in each element and yet, that there are not too many parcels as to constitute an inefficient use of CPU and memory.
  • Limiting During Particle Updates: As the particles are integrated independently from the flow solver, it is not difficult to envision situations where for the extreme cases of very light or very heavy particles physically meaningless or unstable results may be obtained. In order to prevent this, the changes in particle velocities and temperatures are limited in order not to exceed the differences in velocities and temperature between the particles and the fluid [78].
  • Particle Contact/Merging: In some situations, particles may collide or merge in a certain region of space.
  • Particle Tracking: A common feature of all particle-grid applications is that the particles do not move far between timesteps. This makes physical sense: if a particle jumped ten gridpoints during one timestep, it would have no chance to exchange information with the points along the way, leading to serious errors. Therefore, the assumption that the new host elements of the particles are in the vicinity of the current ones is a valid one. For this reason, the most efficient way to search for the new host elements is via the vectorized neighbour-to-neighbour algorithm described in [59,71].

6 Examples

The techniques described above were implemented in FEFLO, a general-purpose computational fluid dynamics (CFD) code based on the following general principles:

  • Use of unstructured grids (automatic grid generation and mesh refinement);
  • Finite element discretization of space;
  • Separate flow modules for compressible and incompressible flows;
  • Edge-based data structures for speed;
  • Optimal data structures for different architectures;
  • Bottom-up coding from the subroutine level to assure an open-ended, expandable architecture.

The code has had a long history of relevant applications involving compressible flow simulations in the areas of transonic flow [81,82,83,84,63,85], store separation [7,10,12,14,15], blast-structure interaction [6,8,11,13,16], [66,72,94,105,97], incompressible flows [91,93,65,68,5,76], free-surface hydrodynamics [62,69,70], dispersion [22,23,67,24], patient-based haemodynamics [26,63,27,1,73] and aeroacoustics [57]. The code has been ported to vector [64], shared memory [61,96], distributed memory [91,60,92,74] and GPU-based [30,31,32,33,75] machines.

The cases shown all simulate sneezing/coughing in different environments. The ambient temperature was assumed to be C. In order to simulate a sneeze/cough, the velocity and temperature in a spherical region of radius () near the patient's mouth was reset at the beginning of each timestep according to the following triangular function:

Failed to parse (unknown function "\cases"): f(t)= \cases { {t \over t_{mid}} ~~if: 0 \le t \le t_{mid}
Failed to parse (syntax error): 0 ~~if: 2 t_{mid} \le t } ~~.
(37)
(38)

The droplets were initialized with 4 different sizes and different velocities, and released every seconds during seconds. This resulted in a final number of particle packets of . The temperature was set to and the velocity to . Table 3 summarizes the diameters and resulting mass distribution.


Table. 3 Initial Conditions for Particles
Droplet Diameter [mm] Mass [gr] Nr. of Packets Nr. of Particles
1.00E+00 5.50E+00 1.05E+03 1.05E+04
1.00E-01 0.11E+00 2.10E+03 2.10E+05
1.00E-02 0.58E-02 1.12E+04 1.12E+07
1.00E-03 0.58E-05 1.12E+04 1.12E+08

In the cases shown different temporal scales appear:

  • The fast, ballistic drop of the larger () particles, occurring in the range of ;
  • The slower drop of particles of diameter , occurring in the range of ; and
  • The transport of the even smaller particles through the air, occurring in the range of .

We have attempted to show these phases in the results, and for this reason the results are not displayed at equal time intervals. Unless otherwise noted, the particles have been colored according to the logarithm of the diameter, with red colors representing the largest and blue the smallest particles.

The examples given show clearly the dangers of droplet- and aerosol- based infections in the built environment.

6.1 Sneezing in Transportation Security Agency (TSA) Queues

One of the obvious vectors for viral contamination and spread are security and passport examination queues in airports. Air flow is moderate, passengers are in very close proximity, and in some airports queues wind back and forth in narrow lanes. Figure 3a,b show the arrangement of pedestrians, as well as the discretization chosen. Note the smaller elements close to the bodies and in the region of interest between the two pedestrians in the middle row. This particular mesh had 12.74Mels. The distribution of particles and the absolute value of the velocity in the centerplane over time can be discerned from Figures 3c-e. One can see that the large (red) particles follow a ballistic path and have some influence on the flow (e.g. at time ). This `ballistic phase' ends at about . The (green) particles of size are quickly stopped by the air, and then sink slowly towards the floor in close proximity to the individual sneezing. The even smaller (cyan, blue) particles rise with the cloud of warmer air exhaled by the sneezing individual, and disperse much further at later times.

Figure 3a  TSA Queue: Arrangement of Pedestrians and Surface Mesh
Figure 3b  TSA Queue: Surface Mesh and Cut Plane
Figure 3c  TSA Queue: Particle Distribution at
Figure 3d  TSA Queue: Particle Distribution at
Figure 3e  TSA Queue: Particle Distribution at
Figure 3f  TSA Queue: Particle Distribution at
Figure 3g  TSA Queue: Particle Distribution at
Figure 3h  TSA Queue: Particle Distribution at

6.2 Sneezing in a Generic Hospital Room

This case considers a typical hospital room. Of interest here was the dispersion of particles in the first minute after coughing, in particular the reach into neighbouring halls and the amount of `negative pressure' needed to keep all contaminants in the room. Figure 4a shows the arrangement of the room, with patient and caregiver clearly visible. This particular mesh had 2.25Mels. The distribution of particles over time can be discerned from Figures 4c-l. As before, one can see that the large (red) particles follow a ballistic path. This `ballistic phase' ends at about . The (green) particles of size are quickly stopped by the air, and then sink slowly towards the patient. The even smaller (cyan, blue) particles rise with the cloud of warmer air exhaled by the sneezing individual, and disperse much further at later times, covering almost the entire room. The velocity distribution in the room may be inferred from Figure 4m.

Figure 4a  Hospital Room: Surface Mesh
Figure 4b  Hospital Room: Particle Distribution at t=0.0 sec
Figure 4c  Hospital Room: Particle Distribution at t=0.2 sec
Figure 4d  Hospital Room: Particle Distribution at t=0.4 sec
Figure 4e  Hospital Room: Particle Distribution at t=0.6 sec
Figure 4f  Hospital Room: Particle Distribution at t=1.0 sec
Figure 4g  Hospital Room: Particle Distribution at t=2.0 sec
Figure 4h  Hospital Room: Particle Distribution at t=4.0 sec
Figure 4i  Hospital Room: Particle Distribution at t=10.0 sec
Figure 4j  Hospital Room: Particle Distribution at t=20.0 sec
Figure 4k  Hospital Room: Particle Distribution at t=30.0 sec
Figure 4l  Hospital Room: Particle Distribution at t=50.0 sec
Figure 4m  Hospital Room: Surface Velocities at t=60.0 sec

7 Reopening After the Crisis

A lingering question facing all levels of society is how and when to reopen facilities where people congregate in close proximity. One key technology that would allow opening is testing and sensing. We consider sensing in the sequel. Several vendors have announced measuring devices for Covid-19 in the next half year. Given that these sensors are expensive, and that a hospital or university many need hundreds of these, the question becomes how best to deploy them. In other words: given an arbitrary number of contamination or infection scenarios, which is the minimum number of sensors needed to detect them, and where should they be placed ? A partial answer to this non-trivial question was given in [67,110]. If we assume a given number of sensors, every contaminant/infection scenario (location and amount of release, flow conditions, etc.) will lead to a sensor input. The data recorded from all the possible release scenarios at all possible sensor locations allows the identification of the best or optimal sensor locations. Clearly, if only one sensor is to be placed, it should be at the location that recorded the highest number of releases. This argument can be used recursively by removing from further consideration all releases already recorded by sensors previously placed. The procedure is repeated recursively until no undetected release cases are left, or the available sensors have been exhausted.

See [21,45] for an in-depth analysis of robust sensor placement under uncertainty.

7.1 Hospital Room

This case considers the same hospital room as shown before. The boundary conditions determining the flow are assumed as steady, with air entering the room through vents 1-3 and exiting the room through the bathroom exhaust or the door. Figures 5a-c show the outlay of the room, average velocities and the `age of air' after 5 minutes. Note the high values for the age of air in the corners and the back of the room. This particular mesh had 2.2Mels. Four contaminant release scenarios were considered: cases 1-3 assumed contaminant coming in through each of the vents (separately) during the first minute, while case 4 assumed virus production from the patient for a period of 10 seconds. The case was run for 5 minutes of real time, and the contaminant concentration was measured on all walls/ceilings. The maximum concentrations measured have been summarized in Figure 5d. Note the different areas covered depending on the release scenario. It was assumed that sensors should only be allowed above a certain height, and should be located on a wall or the ceiling. Table 4 summarizes the points that measured data above a set threshold. As one can see, none of the possible sensor locations is able to measure/detect all 4 cases, and many possible sensor locations do not detect even a single case. There are many possible pairs of sensors that can detect all 4 cases. The pair selected is the one that achieves the highest relative measurement values, and is shown in Figure 5e. Note that this makes good sense: one sensor close the HVAC exits, and one close to the patient.

Figure 5a  Hospital Room: Outlay of Room and Boundary Conditions
Figure 5b  Hospital Room: Average Velocities (5 mins)
Figure 5c  Hospital Room: Age of Air (5 mins)
Figure 5d  Hospital Room: Maximum Contaminant Concentration Over 5 Minutes


Table. 4 Data Measurement Summary
Cases Measured Number
0 4308
1 3377
2 1010
3 0
4 0
Figure 5e  Hospital Room: Optimal Sensor Locations

8 Conclusions and Outlook

The present paper has summarized some of the mechanical characteristics of virus contaminants and the transmission via droplets and aerosols. The ordinary and partial differential equations describing the physics of these processes with high fidelity were given, as well as appropriate numerical schemes to solve them. Several examples taken from recent evaluations of the built environment were given, as well as the optimal placement of sensors.

Current efforts are directed at increasing the realism of the physical processes modeled (e.g. by adding the effect of moving pedestrians [79]), streamlining the simulation toolbox and workflow, and fielding these tools so that the post-pandemic opening can occur as smoothly as possible.

Acknowledgements

H. Antil is partially supported by NSF grants DMS-1818772, DMS-1913004, the Air Force Office of Scientific Research (AFOSR) under Award NO: FA9550-19-1-0036, and Department of Navy, Naval PostGraduate School under Award NO: N00244-20-1-0005

BIBLIOGRAPHY

[1] S. Appanaboyina, F. Mut, R. Löhner, C. M. Putman and J. R. Cebral - Computational Fluid Dynamics of Stented Intracranial Aneurysms Using Adaptive Embedded Unstructured Grids; Int. J. Num. Meth. Fluids 57, 5, 475-493 (2008).

[2] S. Asadi, A.S. Wexler, C.D. Cappa, S. Barreda, N.M. Bouvier and W. Ristenpart - Aerosol Emission and Superemission During Human Speech Increase with Voice Loudness; Nature Scientific Reports 9, (1):2348 (2019). www.nature.com/scientificreports/ https://doi.org/10.1038/s41598-019-38808-z

[3] S. Asadi, A.S. Wexler, C.D. Cappa, S. Barreda, N.M. Bouvier and W. Ristenpart - Effect of Voicing and Articulation Manner on Aerosol Particle Emission During Human Speech PLoS ONE 15(1):e0227699 (2020). https://doi.org/10.1371/journal.pone.0227699

[4] S. Asadi, N.M. Bouvier, A.S. Wexler and W. Ristenpart - The Coronavirus Pandemic and Aerosols: Does COVID-19 Transmit via Expiratory Particles ? Aerosol Science and Technology (2020). https://doi.org/10.1080/02786826.2020.1749229

[5] R. Aubry, F. Mut, R. Löhner and J. R. Cebral - Deflated Preconditioned Conjugate Gradient Solvers for the Pressure-Poisson Equation; J. Comp. Phys. 227, 24, 10196-10208 (2008).

[6] J.D. Baum and R. Löhner - Numerical Simulation of Shock Interaction with a Modern Main Battlefield Tank; AIAA-91-1666 (1991).

[7] J.D. Baum and R. Löhner - Numerical Simulation of Pilot/Seat Ejection from an F-16; AIAA-93-0783 (1993).

[8] J.D. Baum. H. Luo and R. Löhner - Numerical Simulation of a Blast Inside a Boeing 747; AIAA-93-3091 (1993).

[9] J.D. Baum, H. Luo and R. Löhner - Numerical Simulation of a Blast Withing a Multi-Room Shelter; pp. 451-463 in Proc. MABS-13 Conf. The Hague, Netherlands, September (1993).

[10] J.D. Baum, H. Luo and R. Löhner - A New ALE Adaptive Unstructured Methodology for the Simulation of Moving Bodies; AIAA-94-0414 (1994).

[11] J.D. Baum, H. Luo and R. Löhner - Numerical Simulation of Blast in the World Trade Center; AIAA-95-0085 (1995).

[12] J.D. Baum, H. Luo and R. Löhner - Validation of a New ALE, Adaptive Unstructured Moving Body Methodology for Multi-Store Ejection Simulations; AIAA-95-1792 (1995).

[13] J.D. Baum, H. Luo, R. Löhner, C. Yang, D. Pelessone and C. Charman - Coupled Fluid/Structure Modeling of Shock Interaction with a Truck; AIAA-96-0795 (1996).

[14] J.D. Baum, H. Luo, R. Löhner, E. Goldberg and A. Feldhun - Application of Unstructured Adaptive Moving Body Methodology to the Simulation of Fuel Tank Separation from an F-16 c/d Fighter; AIAA-97-0166 (1997).

[15] J.D. Baum, R. Löhner, T.J. Marquette and H. Luo - Numerical Simulation of Aircraft Canopy Trajectory; AIAA-97-1885 (1997).

[16] J.D. Baum, H. Luo, E. Mestreau, R. Löhner, D. Pelessone and C. Charman - Coupled CFD/CSD Methodology for Modeling Weapon Detonation and Fragmentation; AIAA-99-0794 (1999).

[17] J.D. Baum, E. Mestreau, H. Luo, R. Löhner, D. Pelessone, M.E. Giltrud and J.K. Gran - Modeling of Near-Field Blast Wave Evolution; AIAA-06-0191 (2006).

[18] K. Balakrishnan and S. Menon - On the Role of Ambient Reactive Particles in the Mixing and Afterburn Behind Explosive Blast Waves; Combust. Sci. and Tech. 182, 186–214 (2010).

[19] K. Benkiewicz and K. Hayashi - Two-Dimensional Numerical Simulations of Multi-Headed Detonations in Oxygen-Aluminum Mixtures Using an Adaptive Mesh Refinement; Shock Waves , 12, 5, 385-402 (2003).

[20] J.P. Boris, F.F. Grinstein, E.S. Oran, and R.J. Kolbe - New Insights Into Large Eddy Simulation; Fluid Dynamics Research 10, 199-228 (1992).

[21] J.A. Burns and C.N. Rautenberg - The Infinite-Dimensional Optimal Filtering Problem with Mobile and Stationary Sensor Networks; Numerical Functional Analysis and Optimization 36, 2, 181-224 (2015). https://doi.org/10.1080/01630563.2014.970647

[22] F. Camelli and R. Löhner - Assessing Maximum Possible Damage for Contaminant Release Events; Engineering Computations 21, 7, 748-760 (2004).

[23] F. Camelli, R. Löhner, W.C. Sandberg and R. Ramamurti - VLES Study of Ship Stack Gas Dynamics; AIAA-04-0072 (2004).

[24] F. Camelli and R. Löhner - VLES Study of Flow and Dispersion Patterns in Heterogeneous Urban Areas; AIAA-06-1419 (2006).

[25] F. Camelli, J. Lien, D. Dayong, D. W. Wong, M. Rice, R. Löhner and C. Yang - Generating Seamless Surfaces for Transport and Dispersion Modeling in GIS; submitted to GeoInformatica 16, 2, 207-327 (2012).

[26] J.R. Cebral and R. Löhner - From Medical Images to Anatomically Accurate Finite Element Grids; Int. J. Num. Meth. Eng. 51, 985-1008 (2001).

[27] J.R. Cebral and R. Löhner - Efficient Simulation of Blood Flow Past Complex Endovascular Devices Using an Adaptive Embedding Technique; IEEE Transactions on Medical Imaging 24, 4, 468-476 (2005).

[28] C. Chao, M.P. Wan, L. Morawska, G. Johnson, R. Graham, Z. Ristovski M. Hargreaves, K. Mengersen, L. Kerrie C. Steve, Y. Li, X. Xie and S. Katoshevski - Characterization of Expiration Air Jets and Droplet Size Distributions Immediately at the Mouth Opening; J. of Aerosol Science 40, 2, 122-133 (2009).

[29] R. Clift, J.R. Grace and M.E. Weber - Bubbles, Drops and Particles; Academic Press, New York (1978).

[30] A. Corrigan, F. Camelli and R. Löhner - Porting Of An Edge-Based CFD Solver to GPUs; AIAA-10-0523 (2010).

[31] A. Corrigan, F. Camelli, R. Löhner and F. Mut - Porting of FEFLO to GPUs; Proc. ECCOMAS CFD 2010 Conf. Lisbon, Portugal, June 14-17 (2010).

[32] A. Corrigan, F.F. Camelli, R. Löhner and J. Wallin - Running Unstructured Grid Based CFD Solvers on Modern Graphics Hardware; Int. J. Num. Meth. Fluids 66, 221-229 (2011).

[33] A. Corrigan and R. Löhner - Porting of FEFLO to Multi-GPU Clusters; AIAA-11-0948 (2011).

[34] N.G. Deen, M. v.Sint Annaland and J.A.M. Kuipers - Direct Numerical Simulation of Particle Mixing in Dispersed Gas-Liquid-Solid Flows Using a Combined Volume of Fluid and Discrete Particle Approach; Proc. Fifth Int. Conf. on CFD in the Process Industries, CSIRO, Melbourne, Australia, 13-15 December (2006).

[35] L. Dietz, P.F. Horve, D.A. Coil, M. Fretz, J.A. Eisen and L. van den Wymelenberg - 2019 Novel Coronavirus (COVID-19) Pandemic: Built Environment Considerations to Reduce Transmission; mSystems 5:e00245-20 (2020). doi:10.1128/mSystems.00245-20.

[36] P. Fabian, J.J. McDevitt, W.H. Dehaan, R.O.P. Fung, B.J. Cowling, K.H. Chan, et al. - Influenza Virus in Human Exhaled Breath: An Observational Study; PLoS ONE 3:e2691 (2008).

[37] D.R. Franz, P.B. Jahrling, A.M. Friedlander, et al. - Clinical Recognition and Management of Patients Exposed to Biological Warfare Agents; JAMA 278:399e411 (1997).

[38] C. Fureby and F. Grinstein - Monotonically Integrated Large Eddy Simulation of Free Shear Flows; AIAA J. 37, 5, 544-556 (1999).

[39] F.F. Grinstein and C. Fureby - Recent Progress on MILES for High-Reynolds-Number Flows; J. Fluids Engineering 124, 848-861 (2002).

[40] J.K. Gupta, C-H. Lin and Q. Chen - Flow Dynamics and Characterization of a Cough; Indoor Air 19, 517–525 (2009).

[41] J.K. Gupta, C-H. Lin and Q. Chen - Characterizing Exhaled Airflow from Breathing and Talking; Indoor Air 20, 31-39 (2010).

[42] J.K. Gupta, C-H. Lin and Q. Chen - Inhalation of Expiratory Droplets in Aircraft Cabins; Indoor Air 21, 341-350 (2011). doi:10.1111/j.1600-0668.2011.00709.x

[43] J.K. Gupta, C-H. Lin and Q. Chen - Transport of Expiratory Droplets in an Aircraft Cabin; Indoor Air 21, 3-11 (2011).

[44] S.K. Halloran, A.S. Wexler and W.D. Ristenpart - A Comprehensive Breath Plume Model for Disease Transmission via Expiratory Aerosols; PLoS ONE 7(5):e37088 (2012). https://doi.org/10.1371/journal.pone.0037088

[45] M. Hintermüller, C.N. Rautenberg, M. Mohammadi and M. Kanitsar - Optimal Sensor Placement: A Robust Approach; SIAM J. Control Optim. 55(6), 3609-3639 (2017). https://doi.org/10.1137/16M1088867

[46] J. Hoberock and N. Bell - Thrust: Parallel Template Library, Version 1.3 (2010).

[47] S.R. Idelsohn, N. Nigro, A. Larreteguy, J.M. Gimenez and P. Ryshakov - A Pseudo-DNS Method for the Simulation of Incompressible Fluid Flows with Instabilities at Different Scales; Int. J. Comp. Particle Mechanics (2019). https://doi.org/10.1007/s40571-019-00264-x

[48] M. Ip, J.W. Tang, D.S.C. Hui, A.L.N. Wong, M.T.V. Chan, G.M. Joynt, A.T.P. So, S.D. Hall, P.K.S. Chan and J.J.Y. Sung - Airflow and Droplet Spreading Around Oxygen Masks: A Simulation Model for Infection Control Research; AJIC 35, 10, 684-689 (2007).

[49] A. Jameson, W. Schmidt and E. Turkel - Numerical Solution of the Euler Equations by Finite Volume Methods using Runge-Kutta Time-Stepping Schemes; AIAA-81-1259 (1981).

[50] G.R. Johnson, L. Morawska, Z.D. Ristovski, M. Hargreaves, K. Mengersen, C.Y.H. Chao, M.P. Wan, Y. Li , X. Xie, D. Katoshevski, S. Corbette - Modality of Human Expired Aerosol Size Distributions J. of Aerosol Science 42, 839-851 (2011).

[51] G. Kampf, D. Todt, S. Pfaender, E. Steinmann - Persistence of Coronaviruses on Inanimate Surfaces and Their Inactivation With Biocidal Agents; J. of Hospital Infection 104, 3, 246-251, March 01 (2020). https://doi.org/10.1016/j.jhin.2020.01.022

[52] C.K. Kim, J.G. Moon, J.S. Hwang, M.C. Lai and K.S. Im - Afterburning of TNT Explosive Products in Air With Aluminum Particles; AIAA-2008-1029 (2008).

[53] Y.-I. Kim et al. - Infection and Rapid Transmission of SARS-CoV-2 in Ferrets; Cell Host and Microbe 27, 1-6 (2020). https://doi.org/10.1016/j.chom.2020.03.023

[54] Y. Li et al. - Role of Ventilation in Airborne Transmission of Infactious Agents in the Built Environment - A Multidisciplinary Systematic Review; Indoor Air 17, 2-18 (2007).

[55] W.G. Lindsley, F.M. Blachere, R.E. Thewlis, A. Vishnu, K.A. Davis, G. Cao, et al. - Measurements of Airborne Influenza Virus in Aerosol Particles from Human Coughs; PLoS ONE 5:e15100 (2010).

[56] W.G. Lindsley, T.A. Pearce, J.B. Hudnall, K.A. Davis, S.M. Davis, M.A. Fisher, et al. - Quantity and Size Distribution of Cough-Generated Aerosol Particles Produced by Influenza Patients During and After Illness; J. Occup. Environ. Hyg. 9, 443-9. (2012).

[57] J. Liu, K. Kailasanath, R. Ramamurti, D. Munday, E. Gutmark and R. Löhner - Large-Eddy Simulations of a Supersonic Jet and Its Near-Field Acoustic Properties; AIAA J. 47, 8, 1849-1864 (2009).

[58] R. Löhner, K. Morgan, J. Peraire and M. Vahdati - Finite Element Flux-Corrected Transport (FEM-FCT) for the Euler and Navier-Stokes Equations; Int. J. Num. Meth. Fluids 7, 1093-1109 (1987).

[59] R. Löhner and J. Ambrosiano - A Vectorized Particle Tracer for Unstructured Grids; J. Comp. Phys. 91, 1, 22-31 (1990).

[60] R. Löhner and R. Ramamurti - A Load Balancing Algorithm for Unstructured Grids; Comp. Fluid Dyn. 5, 39-58 (1995).

[61] R. Löhner - Renumbering Strategies for Unstructured-Grid Solvers Operating on Shared-Memory, Cache-Based Parallel Machines; Comp. Meth. Appl. Mech. Eng. 163, 95-109 (1998).

[62] R. Löhner, C. Yang and E. Oñate - Viscous Free Surface Hydrodynamics Using Unstructured Grids; Proc. 22nd Symp. Naval Hydrodynamics, Washington, D.C., August (1998).

[63] R. Löhner, Chi Yang, J. Cebral, O. Soto, F. Camelli, J.D. Baum, H. Luo, E. Mestreau, D. Sharov, R. Ramamurti, W. Sandberg and Ch. Oh - Advances in FEFLO; AIAA-01-0592 (2001).

[64] R. Löhner and M. Galle - Minimization of Indirect Addressing for Edge-Based Field Solvers; Comm. Num. Meth. Eng. 18, 335-343 (2002).

[65] R. Löhner - Multistage Explicit Advective Prediction for Projection-Type Incompressible Flow Solvers; J. Comp. Phys. 195, 143-152 (2004).

[66] R. Löhner, J.D. Baum and D. Rice - Comparison of Coarse and Fine Mesh 3-D Euler Predictions for Blast Loads on Generic Building Configurations; Proc. MABS-18 Conf., Bad Reichenhall, Germany, September (2004).

[67] R. Löhner and F. Camelli - Optimal Placement of Sensors for Contaminant Detection Based on Detailed 3-D CFD Simulations; Engineering Computations 22, 3, 260-273 (2005).

[68] R. Löhner, Chi Yang, J.R. Cebral, F. Camelli, O. Soto and J. Waltz - Improving the Speed and Accuracy of Projection-Type Incompressible Flow Solvers; Comp. Meth. Appl. Mech. Eng. 195, 23-24, 3087-3109 (2006).

[69] R. Löhner, Chi Yang and E. Oñate - On the Simulation of Flows with Violent Free Surface Motion; Comp. Meth. Appl. Mech. Eng. 195, 5597-5620 (2006).

[70] R. Löhner, Chi Yang and E. Oñate - Simulation of Flows With Violent Free Surface Motion and Moving Objects Using Unstructured Grids; Int. J. Num. Meth. Fluids 53, 1315-1338 (2007).

[71] R. Löhner - Applied CFD Techniques, Second Edition; J. Wiley & Sons (2008).

[72] R. Löhner, H. Luo, J.D. Baum and D. Rice - Improvements in Speed for Explicit, Transient Compressible Flow Solvers; Int. J. Num. Meth. Fluids 56, 12, 2229-2244 (2008).

[73] R. Löhner, J.R. Cebral, F.F. Camelli, S. Appanaboyina, J.D. Baum, E.L. Mestreau and O. Soto - Adaptive Embedded and Immersed Unstructured Grid Techniques; Comp. Meth. Appl. Mech. Eng. 197, 2173-2197 (2008).

[74] R. Löhner, F. Mut and F.F. Camelli - Timings OF FEFLO on the SGI-ICE Machines; AIAA-11-1064 (2011).

[75] R. Löhner and A. Corrigan - Semi-Automatic Porting if a General Fortran CFD Code to GPUs: The Difficult Modules; AIAA-11-3219 (2011).

[76] R. Löhner, F. Mut, J.R. Cebral, R. Aubry and G. Houzeaux; Deflated Preconditioned Conjugate Gradient Solvers for the Pressure-Poisson Equation: Extensions and Improvements; Int. J. Num. Meth. Eng. 87, 1-5, 2-14 (2011).

[77] R. Löhner - F2GPU - A General Fortran to GPU Translator; Proc. NVIDIA GTC Conf. , San Jose, CA, May (2012).

[78] R. Löhner, F. Camelli, J.D. Baum, F. Togashi and O. Soto - On Mesh-Particle Techniques; Comp. Part. Mech. 1, 199-209 (2014).

[79] R. Löhner and F. Camelli - Tightly Coupled Computational Fluid and Crowd Dynamics; pp. 505-509 in Proc. Pedestrian and Evacuation Dynamics 2016 (PED 2016), (W. Song, J. Ma and L. Fu eds.), University of Science and Technology Press, Hefei, China, Oct 17-21 (2016).

[80] R.G. Loudon and R.M. Roberts - Droplet Expulsion from the Respiratory Tract; Am. Rev. Respir. Dis. 95, 3, 435–442 (1967).

[81] H. Luo, J.D. Baum and R. Löhner - Edge-Based Finite Element Scheme for the Euler Equations; AIAA J. 32, 6, 1183-1190 (1994).

[82] H. Luo, J.D. Baum, R. Löhner and J. Cabello - Implicit Finite Element Schemes and Boundary Conditions for Compressible Flows on Unstructured Grids; AIAA-94-0816 (1994).

[83] H. Luo, J.D. Baum and R. Löhner - An Accurate, Fast, Matrix-Free Implicit Method for Computing Unsteady Flows on Unstructured Grids; AIAA-99-0937 (1999).

[84] H. Luo, D. Sharov, J.D. Baum and R. Löhner - A Class of Matrix-free Implicit Methods for Compressible Flows on Unstructured Grids; First International Conference on Computational Fluid Dynamics, Kyoto, Japan, July 10-14 (2000).

[85] H. Luo, J.D. Baum and R. Löhner - A Fast, Matrix-Free Implicit Method for Computing Low Mach Number Flows on Unstructured Grids; Int. J. CFD 14, 133-157 (2001).

[86] D. Merrill and A. Grimshaw - Revisiting Sorting for GPGPU Stream Architectures; UVA CS Rep. CS2010-03, Charlottesville, VA (2010).

[87] D.K. Milton, M.P. Fabian, B.J. Cowling, M.L. Grantham, J.J. McDevitt - Influenza Virus Aerosols in Human Exhaled Breath: Particle Size, Culturability, and Effect of Surgical Masks; PLoS Pathog. 9:e1003205 (2013).

[88] NVIDIA Corporation. NVIDIA CUDA 3.2 Programming Guide (2010).

[89] J. Peraire, M. Vahdati, K. Morgan and O.C. Zienkiewicz - Adaptive Remeshing for Compressible Flow Computations; J. Comp. Phys. 72, 449-466 (1987).

[90] P. Peterson - F2PY: Tool for Connecting Fortran and Python Programs; Int. J. Computational Science and Engineering 4, 296-305 (2009).

[91] R. Ramamurti and R. Löhner - Simulation of Flow Past Complex Geometries Using a Parallel Implicit Incompressible Flow Solver; pp. 1049,1050 in Proc. 11th AIAA CFD Conf. , Orlando, FL, July (1993).

[92] R. Ramamurti and R. Löhner - A Parallel Implicit Incompressible Flow Solver Using Unstructured Meshes; Computers and Fluids 5, 119-132 (1996).

[93] R. Ramamurti, W.C. Sandberg and R. Löhner - Computation of Unsteady Flow Past Deforming Geometries; Int. J. Comp. Fluid Dyn. , 83-99 (1999).

[94] D.L. Rice, J.D. Baum, F. Togashi, R. Löhner and A. Amini - First-Principles Blast Diffraction Simulations on a Notebook: Accuracy, Resolution and Turn-Around Issues; Proc. MABS-20 Conf., Oslo, Norway, September (2008).

[95] H. Schlichting - Boundary Layer Theory; McGraw-Hill (1979).

[96] D. Sharov, H. Luo, J.D. Baum and R. Löhner - Implementation of Untructured Grid GMRES+LU-SGS Method on Shared-Memory, Cache-Based Parallel Computers; AIAA-00-0927 (2000).

[97] A. Stück, F. Camelli and R. Löhner - Adjoint-Based Design of Shock Mitigation Devices; Int. J. Num. Meth. Fluids 64, 443-472 (2010).

[98] J.W. Tang, Y. Li, I. Eames, P.K.S. Chan and G.L. Ridgway Factors Involved in the Aerosol Transmission of Infection and Control of Ventilation in Healthcare Premises; J. of Hospital Infection 64, 100-114 (2006).

[99] J.W. Tang, C.J. Noakes, P.V. Nielsen, I. Eames, A. Nicolle, Y. Li and G.S. Settles - Observing and Quantifying Airflows in the Infection Control of Aerosol- and Airborne-Transmitted Diseases: An Overview of Approaches; J. of Hospital Infection 77 213-222 (2011).

[100] J.W. Tang, A.D. Nicolle, J. Pantelic, G.C. Koh, L. Wang, M. Amin, C.A. Klettner, D.K.W. Cheong, C. Sekhar and K.W. Tham - Airflow Dynamics of Coughing in Healthy Human Volunteers by Shadowgraph Imaging: An Aid to Aerosol Infection Control; PLoS ONE 7, 4: e34818 (2012). doi:10.1371/journal.pone.0034818

[101] J.W. Tang, A.D. Nicolle, C.A. Klettner, J. Pantelic, L. Wang, A. Bin Suhaimi, A.Y.L. Tan, G.W.X. Ong, R. Su, C. Sekhar, D.K.W. Cheong and K.W. Tham - Airflow Dynamics of Human Jets: Sneezing and Breathing - Potential Sources of Infectious Aerosols; PLoS ONE 8, 4: e59970 (2013). doi:10.1371/journal.pone.0059970

[102] P.F.M. Teunis, N. Brienen, M.E.E. Kretzschmar - High Infectivity and Pathogenicity of Influenza A Virus Via Aerosol and Droplet Transmission; Epidemics 2, 215–222 (2010).

[103] R. Tilch, A. Tabbal, M. Zhu, F. Decker and R. Löhner - Combination of Body-Fitted and Embedded Grids for External Vehicle Aerodynamics; Engineering Computations 25, 1, 28-41 (2008).

[104] K. K.-W. To et al. - Temporal Profiles of Viral Load in Posterior Oropharyngeal Saliva Samples and Serum Antibody Responses During Infection by SARS-CoV-2: An Observational Cohort Study; Lancet Infect. Dis. (online) (2020). https://doi.org/10.1016/S1473-3099(20)30196-1

[105] F. Togashi, J.D. Baum, E. Mestreau, R. Löhner, and D. Sunshine; Numerical Modeling of Long-Duration Blast Wave Evolution in Confined Facilities; AIAA-09-1531 (2009).

[106] N. van Doremalen, T. Bushmaker, D.H. Morris, M.G. Holbrook, A. Gamble, B.N. Williamson, A. Tamin, J.L. Harcourt, N.J. Thornburg, S.I. Gerber, J.O. Lloyd-Smith, E. de Wit, V.J. Munster - Aerosol and Surface Stability of SARS-CoV-2 as Compared with SARS-CoV-1; The New England Journal of Medicine 382, 16, April 16 (2020). DOI: 10.1056/NEJMc2004973

[107] A.W. Vreman, B.J. Geurts, N.G. Deen and J.A.M. Kuipers - Large-Eddy Simulation of a Particle Laden Turbulent Channel Flow; Proc. Direct and Large-Eddy Simulation V, ERCOFTAC Series Volume 9, pp 271-278 (2004).

[108] J. Wei, Y. Li - Airborne Spread of Infectious Agents in the Indoor Environment; American J. of Infection Control 44, S102-S108 (2016). http://dx.doi.org/10.1016/j.ajic.2016.06.003

[109] X. Xie, Y. Li, A.T.Y. Chwang, P.L. Ho, W.H. Seto - How Far Droplets Can Move in Indoor Environments - Revisiting the wells Evaporation-Falling Curve; Indoor Air 17, 211-225 (2007). doi:10.1111/j.1600-0668.2006.00469.x

[110] T. Zhang, Q. Chen and C.-H. Lin - Optimal Sensor Placement for Airborne Contaminant Detection in an Aircraft Cabin; HVAC&R Research 13, 5, 683-696 (2007).

[111] Y. Zhang, G. Feng, Z. Kang, Y. Bi and Y. Cai - Numerical Simulation of Coughed Droplets in Conference Room; 10th International Symposium on Heating, Ventilation and Air Conditioning, ISHVAC2017, October, 19-22 Jinan, China (2017), Procedia Engineering 205, 302–308 (2017).

Back to Top

Document information

Published on 01/01/2020

DOI: 10.1007/s00466-020-01881-7
Licence: CC BY-NC-SA license

Document Score

0

Views 33
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?