# Fluid Dynamics and Numerical Simulation of Exhaled Droplets Containing Infectious Viruses

Authors: Jiawen Han
Journal: Security and Communication Networks
 (2022)

---

## Abstract

Understanding the movement and transmission patterns of airborne particles is very important to understand the diseases they carry. One of the most common sources of viral infections is sneezing and coughing. This topic explores the theoretical properties of airborne microbes and their dispersal. It aims at developing a model that can predict the movement and transmission of these particles. The model is formulated using the Runge–Kutta (RK) algorithm, which is a 4th-order standard for solving differential equations. It is used to study the effects of various factors such as wind speed and jet velocity on the movement of droplets. The model is compared to the well-known Maxey–Riley equation. It then explains the various factors that influence the dispersal of airborne microbes. The evaporation of airborne microbes has a significant effect on the movement of smaller particles.


---

## Body

## 1. Introduction

Most of the time, small particles originating from the mouths of humans stay in the air for a long time due to their sizes. They can also stay suspended in the air for a long time through the flow of air. When people cough or sneeze, the larger particles evaporate and become airborne microorganisms. These microorganisms then spread around in the air due to the wind movement.Due to their small size, microbial particles released from the population do not settle quickly and remain suspended in the air with the airflow, usually remaining in the air for a long time through dust and airflow. The largest droplets from sneezing or coughing, due to their own evaporation, may become droplet nuclei or microorganisms containing pathogens, which are spread into the air with the wind, resulting in droplet transmission.Studies show that the coronaviruses can be transported through three routes. One of these is contact transmission, which involves the direct contact of a pathogen with an object. The second type of transmission is droplet transmission, which involves the transfer of airborne microorganisms from one part of the body to another. This occurs when droplets from a sick patient are dispersed over a short distance. The spread of disease by airborne particles is known as airborne transmission. These particles are composed of various microorganisms [1]. The particles are designed to be rigid and have zero fluid velocity on their surface. They can reach varying horizontal distances depending on their mode of propagation. The combined effect between gravity, inertial forces, drag, and environmental forces determines the fate of saliva droplets.In this article, the RK algorithm will be applied to predict the farthest distances that particles of diverse sizes can reach in air and to consider the effect of evaporation of droplets in air on their movement. The particle will accelerate as a result of the variable force's effect. The uniform and accelerated motion of particles will be discussed separately in this section.The article is organized as follows:Section2 discusses the particle will accelerate as a result of the varied force acting on it. The uniform and accelerated motion of particles will be discussed separately in this section. The Stokes law (SL) is to determine the velocity. They evaluate the acceleration of forces under gravity, and drag forces have a significant impact on the particle as it settles, and we will look at how these two forces affect the particle's motion in this section. Section 3 explains the comprehensive force analysis of droplets motion. Section 4 defines the movement phenomena of microbial aerosols and diffusion of particles. Section 5 analyzes the prediction of microbial aerosol movement and evaluate the different result and discussion. Section 6 completes the article.


## 2. Uniform Motion of Particles and Acceleration

The constant motion of a particle is often the result of the combined action of a constant external force like gravity and the resistance of the gas to the particle’s motion. In most cases, the exhaled gas particles reach a relatively stable state in a split second. Under the action of the variable force, the particle will undergo accelerated motion. In this section, we will discuss the uniform and accelerated motion of particles disjointedly.


### 2.1. Stokes’ Law (SL)

In multiphase fluid mechanics, the problem of tiny solid particles, droplets, or bubbles moving in a viscous fluid is often encountered. In situations where the flow rate is not high and the pressure change is small, the compressibility of the gas can be ignored because of the small change in density caused by its movement. We treat air as an incompressible fluid. In spherical coordinates (r,θ,φ), the continuity equation and Navier–Stokes equations of motion for an incompressible fluid are(1)∂∂rurr2sinθ+∂∂θuθrsinθ+∂∂ϕruϕ=0,(2)ρDurDt−uθ2+uϕ2r=ρfr−∂p∂r+μ∇2ur−2urr2−2r2sinθ∂uθsinθ∂θ−2r2sinθ∂uϕ∂ϕ,(3)ρDuθDt+uruθr−uϕ2cotθr=ρfθ−1r∂p∂θ+μ∇2uθ+2r2∂ur∂θ−uθr2sinθ−2r2sin2θ∂uϕ∂ϕ,(4)ρDuϕDt+uϕurr−uϕuθcotθr=ρfϕ−1rsinθ∂p∂ϕ+μ∇2uϕ+2r2sinθ∂ur∂θ+2cosθr2sin2θ∂uθ∂ϕ−uϕr2sin2θ,where(5)DDt=∂∂t+ur∂∂r+uθr∂∂θ+uϕrsinθ∂∂ϕ,(6)∇2=1r2∂∂rr2∂∂r+1r2sinθ∂∂θsinθ∂∂θ+1r2sin2θ∂2∂ϕ2=∂2∂r2+2r∂∂r+1r2∂2∂θ2+cotθr2∂∂θ+1r2sin2θ∂2∂ϕ2.And ur, uθ,uϕ are the velocities in the r, θ, ϕ directions, µ is the fluid viscosity, p is the pressure, ρ is the fluid density, and fr,fθ,fϕ are the body force components.As can be seen from the definition of the Reynolds number (Re=ρvdp/μ, where v, ρ, and µ are the velocity, density, and viscosity coefficients of the fluid, respectively, and dp is the diameter of particles), the Reynolds number for such problems is very small because of the tiny size of the aerosol particles and their low velocity. The inertial force of the fluid is negligible compared to the viscous force when analyzed from a force point of view. In 1851, George Gabriel Stokes derived an expression for the frictional force exerted by a spherical body with a small Reynolds number in a viscous fluid, which is now known as SL [2]. It is assumed that the particle is a rigid sphere and that the fluid velocity on the surface of the particle is zero.Since the flow field has symmetry, we have∂/∂ϕ=0, uϕ=0. Then, the continuity (1) becomes(7)∂ur∂r+1r∂uθ∂θ+2urr+uθcotθr=0.Combining (2) and (3), the equations for r and θ are(8)∂p∂r=μ∂2ur∂r2+1r2∂2ur∂θ2+2r∂ur∂r+cotθr2∂ur∂θ−2r2∂uθ∂θ−2urr2−2cotθr2uθ,(9)1r∂p∂θ=μ∂2uθ∂r2+1r2∂2uθ∂θ2+2r∂uθ∂r+cotθr2∂uθ∂θ+2r2∂ur∂θ−uθr2sin2θ.And the boundary conditions on the sphere and at infinity are(10)ur|r=a=0,uθ|r=a=0,(11)ur|∞=V∞cosθ,uθ|∞=−V∞sinθ,p|∞=p∞,where the radius of the sphere is a and the velocity of the fluid flowing through the sphere is V∞.The three equations (7)–(9) are solved for the three unknowns ur, uθ, and p, respectively, under boundary conditions.Eqs. (10) and (11) are obtained using the separation of variables method. Here, we omit the specific solving steps and give the following results:(12)urr,θ=V∞cos1−32ar+12a3r3,(13)uθr,θ=−V∞sin1−34ar−14a3r3,(14)pr,θ=−32μV∞ar2cosθ+p∞.Since the flow field is symmetrical aboutφ, the force acting on the sphere by the flow field has only one normal stress and one shear stress component, and the specific form is as follows:(15)prr=−p+2μ∂ur∂r,prθ=μ1r∂ur∂θ+∂uθ∂r−uθr.On the spherer = a, from (10) we have(16)∂ur∂θ=0,∂uθ∂θ=0.Solving the equations obtains the following:(17)prr=−p=32μV∞acosθ−p∞,prθ=−3μV∞2asinθ.Due to thex-axis flow being symmetric, all forces acting on the sphere perpendicular to it are zero. As a result, the drag FD is a straight line perpendicular to the x-axis. Integrating prr and prθ over the surface of the sphere S, its magnitude can be determined:(18)FD=∫∫Sprrcosθ−prθsinθdS=∫02πdϕ∫0πprrcosθ−prθsinθa2sinθdθ,FD=∫∫Sprrcosθ−prθsinθdS=∫02πdϕ∫0πprrcosθ−prθsinθa2sinθdθ,=3πμV∞a∫0πsinθdθ=6πμV∞a.The following equation is the SL. Specifically, the total drag of the air to a spherical particle (moving with velocityup relative to the air with viscosity µa) with diameter dp is(19)FD=3πμaupdp.


### 2.2. Settling Velocity (SV)

An important application of SL is to determine the velocity at which aerosol particles settle by gravity in stationary air. When a particle of aerosol is released into the air, the drag pulling on it is exactly equal to gravity FG and in the opposite direction. In this case, the particles will quickly reach a final SV under the following conditions:(20)FD=FG=mg,(21)3πμaUpdp=ρp−ρa16πdp3g,where g is the gravitational acceleration, Up is the SV, and ρp and ρa are the densities of the particle and air, respectively. When considering the effect of buoyancy, the density of air needs to be taken into account.However, the density of air can be neglected when it is relatively small compared to the density of the particle. For example, an exhaled aerosol settling in air has a density ratioρp/ρa = 1000/1.184 ≈ 845(we approximate the density of the exhaled particles as the density of water and at a standard atmospheric pressure, and the density of air is 1.184 kg/m3 at 25 degrees Celsius [3].), so we omit ρa in the following analysis. For low Reynolds numbers, we obtain(22)Up=ρpdp2g18μa.In the subsequent analysis, we find that the particle will reach its final SV in a very short time. It is therefore reasonable to describe the motion of the particle byUp.


### 2.3. Cunningham Correction Factor

When the size of the particles is small enough to approach the mean free path of the gas molecules, unlike what SL assumes, there is slippage on the surface of the particle. In 1910, Ebenezer Cunningham derived a correction factor for SL and experimentally determined the parameters in air by Davies [4].This factor is called the Cunningham correction factor and is always greater than 1. The empirical equation for the Cunningham correction factor, as measured experimentally by slip, is the following:(23)Cc=1+2λdp1.257+0.4⋅e−0.55dp/λ,where λ is the mean free path for air and dp is the particle diameter.In this way, when considering small aerosol particles, Stokes drag (20) can be written as(24)FD=3πμaupdpCc,so that the final SV is(25)up=ρpdp2gCc18μa.It can be seen that the particle SV is one factorCc faster than the value calculated by the uncorrected SL. The equation for the mean free path(metre) of air molecules at standard atmospheric pressure is [5](26)λ=π8⋅μa0.4987⋅1ρaP=6.536×10−8,where P = 1 × 10−5 Pa and µa = 17.9 × 10−6 Pa s. Substituting (27) into (24), we can obtain that the correction factor for a particle with a diameter of 1 micron is 1.1643, which means that it settles 16% faster than the velocity calculated using the uncorrected SL (2.23) equation. The correction factor increases rapidly with decreasing particle size when the particle diameter is less than 1 micron (as shown in Figure 1), at which point the Cunningham correction factor is used. Also, the slip of gas molecules on the surface of particles is large when their diameter is between 0.1 and 0.01 µm.Figure 1
The variation in the Cunningham correction factor with particle sizes from 0.01 to 0.1µm, 0.1 to 1 µm, and 1 to 10 µm.


### 2.4. Relaxation Time

A particle that settles naturally in air, regardless of its initial velocity, will eventually reach a constant SV without the influence of other forces. The relaxation time of the particle is the amount of time that passes before this constant SV is achieved. The particle's mobility must be understood before the expression can be given. In SL, expressed in (25), the drag is proportional to the velocity. We express the mobility of tiny particles in terms of the velocity B caused by a unit force:(27)B=upFD=Cc3πμadp.By definition, the final velocity of an aerosol particle can then be stated simply as the product of the force and the mobility, for example, when the external force is gravity:(28)up=FG⋅B=mg⋅Cc3πμadp=16πdp3ρpg⋅Cc3πμadp.We now define the relaxation time as the product of the mass of the particle and the mobility, and denote it asτ:(29)τ=mB=16πdp3ρpCc3πμadp=ρpdp2Cc18μa.As seen in the above equation, the relaxation time depends on the mass and mobility of the particle and is not affected by the magnitude of the external forces acting on the particle. It is affected by the viscosity of air and Cunningham correction factor.For gravity, substituting (34) into (33) yields that(30)up=τ⋅g.In general, for any external forceF acting invariably on a particle, the final velocity can be expressed as(31)up=τ⋅Fm,where m is the mass of the particle.


### 2.5. Acceleration of Particles

In Subsection2.2, we made the forces on the pellet balanced such that the pellet moved at a uniform speed. In this subsection, the velocity of the particle V (t) is considered as a function of time, and eventually, the particle will accelerate to equilibrium. Gravity and drag forces have a comparatively large effect on the particle as it settles, and in this section, we focus on the effect of these two forces on the motion of the particle. According to Newton’s second law,(32)∑F=dmVtdt=mdVtdt.To simplify the calculations, this section does not consider the evaporation of droplets, that is, the mass of the particle remains constant, and, at the same time, the slip of air molecules on the surface of the particle is neglected. Then, we have(33)FG−FD=mg−3πμaVtdp=mdVtdt.The boundary conditions of the problem are(34)Vt=0whent=0,(35)Vt⟶upast≫τ.Multiplying both sides of the (38) by the mechanical mobility of the particle B gives(36)Bmg−3πμaVtdp=B⋅mdVtdt.Sinceτ = mB and omitting the effect of the correction factor, B=1/3πμadp, (39) becomes(37)τg−Vt=τdVtdt.Meanwhile, substituting SVup = τg into the above equation gives(38)∫0tdtτ=∫0VtdVtup−Vt.The above equation is integrated to give that(39)tτ=−lnup−Vt+lnup,e−t/τ=up−Vtup.The final solution is that(40)Vt=up1−e−t/τ.SupposeV0 is the initial velocity for the particle, and then, more generally, (45) can be written as(41)Vt=up−up−V0e−t/τ.Observing the above equation, the instantaneous velocityV (t) of the particle is never reached, but can be infinitely close to the SV up. If the particle is released from rest, at t = τ, the particle’s velocity at this point differs from the SV by 1/eup. The velocity of the particle will be infinitely close to its final SV when t = 4τ from Figure 2 [6].Figure 2
Velocity versus time for an accelerating particle.


### 2.6. Distance Travelled by the Particles

The displacement of a particle subjected to gravity and drag can be found from the instantaneous velocity obtained in the previous subsection. Assume that the distance of the particle isx(t), and then,(42)dxtdt=Vt.Substituting (41) into the above equation and integrating yields(43)dx=up−up−V0e−t/τ⋅dt,(44)∫0xtdx=∫0tupdt−∫0tup−V0e−t/τdt.This results in the following:(45)xt=upt−τup−V01−e−t/τ.When a particle settles to the ground, it reaches the farthest distance it has experienced, that is, the distance at which the particle’s velocity drops to zero. At this point,up = 0 and(46)xmax=τV0,since t≫τ. We also refer to xmax as the stopping distance. From the definition of the relaxation time, rewriting equation (51) gives(47)xmax=mBV0.For larger diameter particles subject only to gravity and drag, they will eventually reach their farthest distances. Particles of both big and tiny sizes will swiftly slow down to near-stop when the starting velocity is high, resulting in a short stopping distance. A particle with a diameter of 100 microns travelling at a speed of 10 metres per second, for example, will stop in only 12.7 centimetres. What is clear is that if the particles are in an environment with an air wind velocity, they are subject to more ground forces affecting their trajectory and eventually the motion of the particles before settling to the ground can even be dominated by the motion of the ambient gas.


## 2.1. Stokes’ Law (SL)

In multiphase fluid mechanics, the problem of tiny solid particles, droplets, or bubbles moving in a viscous fluid is often encountered. In situations where the flow rate is not high and the pressure change is small, the compressibility of the gas can be ignored because of the small change in density caused by its movement. We treat air as an incompressible fluid. In spherical coordinates (r,θ,φ), the continuity equation and Navier–Stokes equations of motion for an incompressible fluid are(1)∂∂rurr2sinθ+∂∂θuθrsinθ+∂∂ϕruϕ=0,(2)ρDurDt−uθ2+uϕ2r=ρfr−∂p∂r+μ∇2ur−2urr2−2r2sinθ∂uθsinθ∂θ−2r2sinθ∂uϕ∂ϕ,(3)ρDuθDt+uruθr−uϕ2cotθr=ρfθ−1r∂p∂θ+μ∇2uθ+2r2∂ur∂θ−uθr2sinθ−2r2sin2θ∂uϕ∂ϕ,(4)ρDuϕDt+uϕurr−uϕuθcotθr=ρfϕ−1rsinθ∂p∂ϕ+μ∇2uϕ+2r2sinθ∂ur∂θ+2cosθr2sin2θ∂uθ∂ϕ−uϕr2sin2θ,where(5)DDt=∂∂t+ur∂∂r+uθr∂∂θ+uϕrsinθ∂∂ϕ,(6)∇2=1r2∂∂rr2∂∂r+1r2sinθ∂∂θsinθ∂∂θ+1r2sin2θ∂2∂ϕ2=∂2∂r2+2r∂∂r+1r2∂2∂θ2+cotθr2∂∂θ+1r2sin2θ∂2∂ϕ2.And ur, uθ,uϕ are the velocities in the r, θ, ϕ directions, µ is the fluid viscosity, p is the pressure, ρ is the fluid density, and fr,fθ,fϕ are the body force components.As can be seen from the definition of the Reynolds number (Re=ρvdp/μ, where v, ρ, and µ are the velocity, density, and viscosity coefficients of the fluid, respectively, and dp is the diameter of particles), the Reynolds number for such problems is very small because of the tiny size of the aerosol particles and their low velocity. The inertial force of the fluid is negligible compared to the viscous force when analyzed from a force point of view. In 1851, George Gabriel Stokes derived an expression for the frictional force exerted by a spherical body with a small Reynolds number in a viscous fluid, which is now known as SL [2]. It is assumed that the particle is a rigid sphere and that the fluid velocity on the surface of the particle is zero.Since the flow field has symmetry, we have∂/∂ϕ=0, uϕ=0. Then, the continuity (1) becomes(7)∂ur∂r+1r∂uθ∂θ+2urr+uθcotθr=0.Combining (2) and (3), the equations for r and θ are(8)∂p∂r=μ∂2ur∂r2+1r2∂2ur∂θ2+2r∂ur∂r+cotθr2∂ur∂θ−2r2∂uθ∂θ−2urr2−2cotθr2uθ,(9)1r∂p∂θ=μ∂2uθ∂r2+1r2∂2uθ∂θ2+2r∂uθ∂r+cotθr2∂uθ∂θ+2r2∂ur∂θ−uθr2sin2θ.And the boundary conditions on the sphere and at infinity are(10)ur|r=a=0,uθ|r=a=0,(11)ur|∞=V∞cosθ,uθ|∞=−V∞sinθ,p|∞=p∞,where the radius of the sphere is a and the velocity of the fluid flowing through the sphere is V∞.The three equations (7)–(9) are solved for the three unknowns ur, uθ, and p, respectively, under boundary conditions.Eqs. (10) and (11) are obtained using the separation of variables method. Here, we omit the specific solving steps and give the following results:(12)urr,θ=V∞cos1−32ar+12a3r3,(13)uθr,θ=−V∞sin1−34ar−14a3r3,(14)pr,θ=−32μV∞ar2cosθ+p∞.Since the flow field is symmetrical aboutφ, the force acting on the sphere by the flow field has only one normal stress and one shear stress component, and the specific form is as follows:(15)prr=−p+2μ∂ur∂r,prθ=μ1r∂ur∂θ+∂uθ∂r−uθr.On the spherer = a, from (10) we have(16)∂ur∂θ=0,∂uθ∂θ=0.Solving the equations obtains the following:(17)prr=−p=32μV∞acosθ−p∞,prθ=−3μV∞2asinθ.Due to thex-axis flow being symmetric, all forces acting on the sphere perpendicular to it are zero. As a result, the drag FD is a straight line perpendicular to the x-axis. Integrating prr and prθ over the surface of the sphere S, its magnitude can be determined:(18)FD=∫∫Sprrcosθ−prθsinθdS=∫02πdϕ∫0πprrcosθ−prθsinθa2sinθdθ,FD=∫∫Sprrcosθ−prθsinθdS=∫02πdϕ∫0πprrcosθ−prθsinθa2sinθdθ,=3πμV∞a∫0πsinθdθ=6πμV∞a.The following equation is the SL. Specifically, the total drag of the air to a spherical particle (moving with velocityup relative to the air with viscosity µa) with diameter dp is(19)FD=3πμaupdp.


## 2.2. Settling Velocity (SV)

An important application of SL is to determine the velocity at which aerosol particles settle by gravity in stationary air. When a particle of aerosol is released into the air, the drag pulling on it is exactly equal to gravity FG and in the opposite direction. In this case, the particles will quickly reach a final SV under the following conditions:(20)FD=FG=mg,(21)3πμaUpdp=ρp−ρa16πdp3g,where g is the gravitational acceleration, Up is the SV, and ρp and ρa are the densities of the particle and air, respectively. When considering the effect of buoyancy, the density of air needs to be taken into account.However, the density of air can be neglected when it is relatively small compared to the density of the particle. For example, an exhaled aerosol settling in air has a density ratioρp/ρa = 1000/1.184 ≈ 845(we approximate the density of the exhaled particles as the density of water and at a standard atmospheric pressure, and the density of air is 1.184 kg/m3 at 25 degrees Celsius [3].), so we omit ρa in the following analysis. For low Reynolds numbers, we obtain(22)Up=ρpdp2g18μa.In the subsequent analysis, we find that the particle will reach its final SV in a very short time. It is therefore reasonable to describe the motion of the particle byUp.


## 2.3. Cunningham Correction Factor

When the size of the particles is small enough to approach the mean free path of the gas molecules, unlike what SL assumes, there is slippage on the surface of the particle. In 1910, Ebenezer Cunningham derived a correction factor for SL and experimentally determined the parameters in air by Davies [4].This factor is called the Cunningham correction factor and is always greater than 1. The empirical equation for the Cunningham correction factor, as measured experimentally by slip, is the following:(23)Cc=1+2λdp1.257+0.4⋅e−0.55dp/λ,where λ is the mean free path for air and dp is the particle diameter.In this way, when considering small aerosol particles, Stokes drag (20) can be written as(24)FD=3πμaupdpCc,so that the final SV is(25)up=ρpdp2gCc18μa.It can be seen that the particle SV is one factorCc faster than the value calculated by the uncorrected SL. The equation for the mean free path(metre) of air molecules at standard atmospheric pressure is [5](26)λ=π8⋅μa0.4987⋅1ρaP=6.536×10−8,where P = 1 × 10−5 Pa and µa = 17.9 × 10−6 Pa s. Substituting (27) into (24), we can obtain that the correction factor for a particle with a diameter of 1 micron is 1.1643, which means that it settles 16% faster than the velocity calculated using the uncorrected SL (2.23) equation. The correction factor increases rapidly with decreasing particle size when the particle diameter is less than 1 micron (as shown in Figure 1), at which point the Cunningham correction factor is used. Also, the slip of gas molecules on the surface of particles is large when their diameter is between 0.1 and 0.01 µm.Figure 1
The variation in the Cunningham correction factor with particle sizes from 0.01 to 0.1µm, 0.1 to 1 µm, and 1 to 10 µm.


## 2.4. Relaxation Time

A particle that settles naturally in air, regardless of its initial velocity, will eventually reach a constant SV without the influence of other forces. The relaxation time of the particle is the amount of time that passes before this constant SV is achieved. The particle's mobility must be understood before the expression can be given. In SL, expressed in (25), the drag is proportional to the velocity. We express the mobility of tiny particles in terms of the velocity B caused by a unit force:(27)B=upFD=Cc3πμadp.By definition, the final velocity of an aerosol particle can then be stated simply as the product of the force and the mobility, for example, when the external force is gravity:(28)up=FG⋅B=mg⋅Cc3πμadp=16πdp3ρpg⋅Cc3πμadp.We now define the relaxation time as the product of the mass of the particle and the mobility, and denote it asτ:(29)τ=mB=16πdp3ρpCc3πμadp=ρpdp2Cc18μa.As seen in the above equation, the relaxation time depends on the mass and mobility of the particle and is not affected by the magnitude of the external forces acting on the particle. It is affected by the viscosity of air and Cunningham correction factor.For gravity, substituting (34) into (33) yields that(30)up=τ⋅g.In general, for any external forceF acting invariably on a particle, the final velocity can be expressed as(31)up=τ⋅Fm,where m is the mass of the particle.


## 2.5. Acceleration of Particles

In Subsection2.2, we made the forces on the pellet balanced such that the pellet moved at a uniform speed. In this subsection, the velocity of the particle V (t) is considered as a function of time, and eventually, the particle will accelerate to equilibrium. Gravity and drag forces have a comparatively large effect on the particle as it settles, and in this section, we focus on the effect of these two forces on the motion of the particle. According to Newton’s second law,(32)∑F=dmVtdt=mdVtdt.To simplify the calculations, this section does not consider the evaporation of droplets, that is, the mass of the particle remains constant, and, at the same time, the slip of air molecules on the surface of the particle is neglected. Then, we have(33)FG−FD=mg−3πμaVtdp=mdVtdt.The boundary conditions of the problem are(34)Vt=0whent=0,(35)Vt⟶upast≫τ.Multiplying both sides of the (38) by the mechanical mobility of the particle B gives(36)Bmg−3πμaVtdp=B⋅mdVtdt.Sinceτ = mB and omitting the effect of the correction factor, B=1/3πμadp, (39) becomes(37)τg−Vt=τdVtdt.Meanwhile, substituting SVup = τg into the above equation gives(38)∫0tdtτ=∫0VtdVtup−Vt.The above equation is integrated to give that(39)tτ=−lnup−Vt+lnup,e−t/τ=up−Vtup.The final solution is that(40)Vt=up1−e−t/τ.SupposeV0 is the initial velocity for the particle, and then, more generally, (45) can be written as(41)Vt=up−up−V0e−t/τ.Observing the above equation, the instantaneous velocityV (t) of the particle is never reached, but can be infinitely close to the SV up. If the particle is released from rest, at t = τ, the particle’s velocity at this point differs from the SV by 1/eup. The velocity of the particle will be infinitely close to its final SV when t = 4τ from Figure 2 [6].Figure 2
Velocity versus time for an accelerating particle.


## 2.6. Distance Travelled by the Particles

The displacement of a particle subjected to gravity and drag can be found from the instantaneous velocity obtained in the previous subsection. Assume that the distance of the particle isx(t), and then,(42)dxtdt=Vt.Substituting (41) into the above equation and integrating yields(43)dx=up−up−V0e−t/τ⋅dt,(44)∫0xtdx=∫0tupdt−∫0tup−V0e−t/τdt.This results in the following:(45)xt=upt−τup−V01−e−t/τ.When a particle settles to the ground, it reaches the farthest distance it has experienced, that is, the distance at which the particle’s velocity drops to zero. At this point,up = 0 and(46)xmax=τV0,since t≫τ. We also refer to xmax as the stopping distance. From the definition of the relaxation time, rewriting equation (51) gives(47)xmax=mBV0.For larger diameter particles subject only to gravity and drag, they will eventually reach their farthest distances. Particles of both big and tiny sizes will swiftly slow down to near-stop when the starting velocity is high, resulting in a short stopping distance. A particle with a diameter of 100 microns travelling at a speed of 10 metres per second, for example, will stop in only 12.7 centimetres. What is clear is that if the particles are in an environment with an air wind velocity, they are subject to more ground forces affecting their trajectory and eventually the motion of the particles before settling to the ground can even be dominated by the motion of the ambient gas.


## 3. Comprehensive Force Analysis of Droplets Motion

The forces on these particles are studied in order to determine their movement patterns. However, the wide variety of microorganisms present in the air can make the analysis challenging. To simplify the process, the spherical particles were analyzed, and we do not consider the rotation of droplets.The droplet is exhaled from the mouth by a combination of the following eight forces: air resistanceFd (drag) [7], thermophoretic force Fth [8], gravity Fg and buoyancy Fb, added drag Fm [9], Basset force FB10, Saffman lift force FS11, and Brownian motion force FBi12.The motion of the droplet in air is controlled by a number of dominant forces, and thus, we can distinguish exactly which forces are dominant by calculating the order of magnitude of each force. Calculation is based on the following assumptions:(1)
At 36.5 degrees Celsius, the diameter of the virus-containing droplets ejected by coughing is 10 microns [13], and its density and air density are 993 kg/m31.184 kg/m3, respectively.(2)
Whent = 20°C, the thermodynamic temperature is T = 273.15 + t = 293.15 K.(3)
The initial coughing droplet velocity is assumed to be 28 m/s, and the ambient wind speed is 0.4 m/s to simplify the calculations. The aerodynamic viscosity is 17.9 × 10−6 N · s/m2, and the droplet dynamic viscosity is 7.2 × 10−4 N · s/m2.The following table illustrates the order of magnitudes of each force.The amount of the various forces exerted on microbial aerosols is not the same, as can be observed from the preceding results. Because the thermophoretic force has such a small impact on the particles in comparison to the gravity and air resistance forces, it will be omitted from the numerical computations.


## 4. Movement Phenomena of Microbial Aerosols

Ordinary aerosols' motion was previously examined and analyzed in the previous section. Based on their trajectory study, it can be determined that they are subjected to a wide range of forces, making their motion extremely complicated. Physical events such as evaporation, sedimentation, and diffusion accompany the movement of microbiological particles due to their inherent features. This section examines and studies the many physical mechanisms of microbial aerosol propagation in the air, which in turn determine the trajectory of their motion.


### 4.1. Evaporation Nucleation

The human body is one of the sources of microbial aerosol emission. Transmission of microbial aerosols takes place mainly through talking, cough, and sneezing. When microbial aerosols leave the body, they are often in the form of droplets. The human droplets containing microorganisms evaporate continuously during their movement and eventually change into droplet nuclei, also known as bacterial/viral particles. As a result, analyzing droplet evaporation is important in order to understand the movement of microorganisms. Since the 1940s, researchers have examined the distribution of droplet sizes and the total number of droplets. Jennison [14] measured the size of droplets produced by the human nose and mouth when sneezing by means of high-speed photography. The results showed that 40% to 80% of the droplets were under 100 µm, 20% to 40% were under 50 µm, and the speed of the droplets could reach 46 m/s. It was noted that the evaporation time of these droplets was very short, and therefore, most of the droplets would evaporate into droplet nuclei in the air.Duguid [15] also carried out a measurement experiment on human droplets by placing a glass slide in front of the mouth, and after the droplets from human speech, coughing and sneezing had impacted on the slide, and the slide was then placed under a microscope for observation and counting. In his research, he discovered that the particle size of droplets produced by talking, coughing, and sneezing ranged from 1 to 2000 µm, with 95% of the droplets falling between 2 and 100 µm in size. Due to the limitations of the experimental apparatus at the time, only droplets larger than 1 µm could be observed by the microscope, and after studying various breathing behaviors, Duguid found a similar distribution of droplets produced by talking, coughing (with mouth open/closed), and sneezing, and Figure 3 shows the particle size and number of droplets produced by the four different behaviors.Figure 3
Droplet sizes and quantity produced by different respiratory behaviors.Following a series of investigations, Duguid came to the conclusion that the geometric particle size distribution of droplet nuclei spanned from 0.25 to 42µm, with roughly 97% falling between 0.5 and 12 µm and the majority falling between 1 and 2 µm.The use of laser particle counters has made it easier to measure the number of particles produced by coughing, sneezing, talking, etc. Papineni’s laser particle counter tests [16] found that the majority of particles exhaled from the mouth were in the 2 µm range, with very few particles larger than 8 µm. This conclusion may seem inconsistent with the above, but it shows that large droplets change into droplet nuclei in a very short period of time.Irving proposed an equation [17] to calculate problem relevant to evaporation, which is given as follows:(48)dmdt=s∫Ddρ.The mass loss rate is represented asdm/dt, while the shape factor s, the diffusion coefficient D, and the evaporated substance's vapour density ρ are given.According to the gas law, we can write(49)ρ=MRp∞T∞−pdTd.Here, the molecular weightM, the gas constant R = 8.3145 J/(K · mol), and ambient and droplet surface temperatures T∞ and Td are taken into account, with pd being the pressure of vapour at the surface of the particle, and p∞ being the partial pressure of the vapour much far away from the droplet.The shape factor for a spherical particle is(50)s=2πdp⋅bb−dp,where b is the diameter of the outside of the film of gas. If we assume that b is very large compared to particle diameter dp, we simply obtain(51)s=2πdp.Substituting this together with (49) in (48), we find(52)dmdt=12ρpπdp2⋅ddpdt=2πDMdpRp∞T∞−pdTd,(53)ddpdt=4DMRdpρpp∞T∞−pdTd.Calculation of the time for a droplet of diameterd1 evaporating to diameter d2(d1> d2) using the equation above is given as(54)∫d1d2dpddp=∫0t4DMRρpp∞T∞−pdTddt,(55)t=Rρp8DM⋅d22−d12p∞T∞−pdTd−1.A 100µm exhaled droplet, for example, is expelled from the mouth at a temperature ranging from 37°C to 20°C. To get the mass diffusivity of 0.282 cm/s2 and the saturated vapour pressures pd = 6.2795 kPa and p∞ = 2.3388 kPa, we assume the droplets have the same molar mass as water (M = 18 × 10−3kg/mol). When d2 = 0, the droplet has a lifetime of(56)t=8.3145×9938×0.282×10−4×18×10−3⋅0−100×10−622.3388×10320+273.15−6.2795×10337+273.15−1,=1.65723s.It takes 1.66 seconds for a 100µm droplet to evaporate fully in the air, based on the aforementioned calculations. Using these equations, we can see how air molecules' mean free path is smaller than the particle diameter. When the particle diameter is less than 0.1 µm, the evaporation rate of the droplet is influenced by the motion of the gas molecules. Davies [18] proposes an integrated expression for the mass transport rate, where the particle size reduction rate becomes(57)ddpdt=4DMRdpρpp∞T∞−pdTd2λ+dpdp+5.33σ/dp+3.42λ,where λ is the mean free path of the air molecules and σ is the Boltzmann constant. The solution is implicitly given by(58)t=Rρp5.33σ−2.84λ24DMp∞/T∞−pd/Td⋅log2λ+dp+1.42λdp+12dp2.As in the case of the example above, the specific values of the parameters are substituted into equation (57) and the equation is solved using the fourth-order standard RK method (discuss later) to obtain the evaporation times for different diameter droplets, as shown below.Observation of the images shows that the rate of evaporation and shrinkage of the particles is slow at the beginning and accelerates as the particle size decreases. In Figure4, a particle with a diameter of 100 µm takes 1.67 s to evaporate completely in 20° air, four times the time required for a particle with a diameter of 50 µm. As seen in Figure 5, the emission of smaller particles into air at different temperatures has little effect on their evaporation nucleation time. Both graphs illustrate that the time to complete evaporation is very short for small exhaled particles.Figure 4
The variation in diameter with time equation (53).Figure 5
The variation in diameter with time for 15, 50, and 100µm particles, respectively, according to 0.6 µm particle emission at different ambient temperatures (10°, 20°, and 30°) (equation (57)).


### 4.2. Diffusion of Particles

The microscopic pulsating motion of the microbial particles themselves due to Brownian forces is also known as Brownian motion. It is the result of the interaction between microbial aerosol particles and air molecules. Brownian motion can be described by Fick’s laws of diffusion [19]. The relationship is as follows:(59)JB=−DB⋅∇C,where JB is the diffusion flux and ∇C is the concentration gradient. A negative sign indicates a change in concentration from a high to a low concentration. DB is the Brownian diffusion coefficient of a droplet, which characterizes how violently the particle diffuses. It can be studied analytically according to the Stokes–Einstein relation [20]:(60)DB=σBT=σCcT3πμadp,where B is the mobility, or the ratio of the terminal velocity of the particle to the applied force, and T is the thermodynamic temperature. (58) shows that Brownian motion can only occur when there is a concentration gradient of particles in space. It is the main mechanism for the movement of ultrafine particles in space. However, for particles larger than 0.1 µm, Brownian motion is not the main mechanism of propagation. Brownian motion causes aerosol particles to move in an undirected random motion, but in the case of an individual particle, the probability that it will return to its original position after a certain amount of time is extremely small. Over a time, interval t larger than the relaxation time τ, each particle will move a net distance X. William’s book [6] replaces each displacement with the root mean square of the displacement, which can be estimated by the following equation:(61)X=2DBt.In general,t > 10τ. For simplicity, take t = 1 s. Based on the above analysis, the net displacement X of a 0.1 µm diameter droplet in air is given as 2.3239 × 10−5 m, that is, 23.24 µm. The net distance of the particle by Brownian diffusion alone is actually very small, but the overall tendency is to diffuse from higher to lower concentrations.


## 4.1. Evaporation Nucleation

The human body is one of the sources of microbial aerosol emission. Transmission of microbial aerosols takes place mainly through talking, cough, and sneezing. When microbial aerosols leave the body, they are often in the form of droplets. The human droplets containing microorganisms evaporate continuously during their movement and eventually change into droplet nuclei, also known as bacterial/viral particles. As a result, analyzing droplet evaporation is important in order to understand the movement of microorganisms. Since the 1940s, researchers have examined the distribution of droplet sizes and the total number of droplets. Jennison [14] measured the size of droplets produced by the human nose and mouth when sneezing by means of high-speed photography. The results showed that 40% to 80% of the droplets were under 100 µm, 20% to 40% were under 50 µm, and the speed of the droplets could reach 46 m/s. It was noted that the evaporation time of these droplets was very short, and therefore, most of the droplets would evaporate into droplet nuclei in the air.Duguid [15] also carried out a measurement experiment on human droplets by placing a glass slide in front of the mouth, and after the droplets from human speech, coughing and sneezing had impacted on the slide, and the slide was then placed under a microscope for observation and counting. In his research, he discovered that the particle size of droplets produced by talking, coughing, and sneezing ranged from 1 to 2000 µm, with 95% of the droplets falling between 2 and 100 µm in size. Due to the limitations of the experimental apparatus at the time, only droplets larger than 1 µm could be observed by the microscope, and after studying various breathing behaviors, Duguid found a similar distribution of droplets produced by talking, coughing (with mouth open/closed), and sneezing, and Figure 3 shows the particle size and number of droplets produced by the four different behaviors.Figure 3
Droplet sizes and quantity produced by different respiratory behaviors.Following a series of investigations, Duguid came to the conclusion that the geometric particle size distribution of droplet nuclei spanned from 0.25 to 42µm, with roughly 97% falling between 0.5 and 12 µm and the majority falling between 1 and 2 µm.The use of laser particle counters has made it easier to measure the number of particles produced by coughing, sneezing, talking, etc. Papineni’s laser particle counter tests [16] found that the majority of particles exhaled from the mouth were in the 2 µm range, with very few particles larger than 8 µm. This conclusion may seem inconsistent with the above, but it shows that large droplets change into droplet nuclei in a very short period of time.Irving proposed an equation [17] to calculate problem relevant to evaporation, which is given as follows:(48)dmdt=s∫Ddρ.The mass loss rate is represented asdm/dt, while the shape factor s, the diffusion coefficient D, and the evaporated substance's vapour density ρ are given.According to the gas law, we can write(49)ρ=MRp∞T∞−pdTd.Here, the molecular weightM, the gas constant R = 8.3145 J/(K · mol), and ambient and droplet surface temperatures T∞ and Td are taken into account, with pd being the pressure of vapour at the surface of the particle, and p∞ being the partial pressure of the vapour much far away from the droplet.The shape factor for a spherical particle is(50)s=2πdp⋅bb−dp,where b is the diameter of the outside of the film of gas. If we assume that b is very large compared to particle diameter dp, we simply obtain(51)s=2πdp.Substituting this together with (49) in (48), we find(52)dmdt=12ρpπdp2⋅ddpdt=2πDMdpRp∞T∞−pdTd,(53)ddpdt=4DMRdpρpp∞T∞−pdTd.Calculation of the time for a droplet of diameterd1 evaporating to diameter d2(d1> d2) using the equation above is given as(54)∫d1d2dpddp=∫0t4DMRρpp∞T∞−pdTddt,(55)t=Rρp8DM⋅d22−d12p∞T∞−pdTd−1.A 100µm exhaled droplet, for example, is expelled from the mouth at a temperature ranging from 37°C to 20°C. To get the mass diffusivity of 0.282 cm/s2 and the saturated vapour pressures pd = 6.2795 kPa and p∞ = 2.3388 kPa, we assume the droplets have the same molar mass as water (M = 18 × 10−3kg/mol). When d2 = 0, the droplet has a lifetime of(56)t=8.3145×9938×0.282×10−4×18×10−3⋅0−100×10−622.3388×10320+273.15−6.2795×10337+273.15−1,=1.65723s.It takes 1.66 seconds for a 100µm droplet to evaporate fully in the air, based on the aforementioned calculations. Using these equations, we can see how air molecules' mean free path is smaller than the particle diameter. When the particle diameter is less than 0.1 µm, the evaporation rate of the droplet is influenced by the motion of the gas molecules. Davies [18] proposes an integrated expression for the mass transport rate, where the particle size reduction rate becomes(57)ddpdt=4DMRdpρpp∞T∞−pdTd2λ+dpdp+5.33σ/dp+3.42λ,where λ is the mean free path of the air molecules and σ is the Boltzmann constant. The solution is implicitly given by(58)t=Rρp5.33σ−2.84λ24DMp∞/T∞−pd/Td⋅log2λ+dp+1.42λdp+12dp2.As in the case of the example above, the specific values of the parameters are substituted into equation (57) and the equation is solved using the fourth-order standard RK method (discuss later) to obtain the evaporation times for different diameter droplets, as shown below.Observation of the images shows that the rate of evaporation and shrinkage of the particles is slow at the beginning and accelerates as the particle size decreases. In Figure4, a particle with a diameter of 100 µm takes 1.67 s to evaporate completely in 20° air, four times the time required for a particle with a diameter of 50 µm. As seen in Figure 5, the emission of smaller particles into air at different temperatures has little effect on their evaporation nucleation time. Both graphs illustrate that the time to complete evaporation is very short for small exhaled particles.Figure 4
The variation in diameter with time equation (53).Figure 5
The variation in diameter with time for 15, 50, and 100µm particles, respectively, according to 0.6 µm particle emission at different ambient temperatures (10°, 20°, and 30°) (equation (57)).


## 4.2. Diffusion of Particles

The microscopic pulsating motion of the microbial particles themselves due to Brownian forces is also known as Brownian motion. It is the result of the interaction between microbial aerosol particles and air molecules. Brownian motion can be described by Fick’s laws of diffusion [19]. The relationship is as follows:(59)JB=−DB⋅∇C,where JB is the diffusion flux and ∇C is the concentration gradient. A negative sign indicates a change in concentration from a high to a low concentration. DB is the Brownian diffusion coefficient of a droplet, which characterizes how violently the particle diffuses. It can be studied analytically according to the Stokes–Einstein relation [20]:(60)DB=σBT=σCcT3πμadp,where B is the mobility, or the ratio of the terminal velocity of the particle to the applied force, and T is the thermodynamic temperature. (58) shows that Brownian motion can only occur when there is a concentration gradient of particles in space. It is the main mechanism for the movement of ultrafine particles in space. However, for particles larger than 0.1 µm, Brownian motion is not the main mechanism of propagation. Brownian motion causes aerosol particles to move in an undirected random motion, but in the case of an individual particle, the probability that it will return to its original position after a certain amount of time is extremely small. Over a time, interval t larger than the relaxation time τ, each particle will move a net distance X. William’s book [6] replaces each displacement with the root mean square of the displacement, which can be estimated by the following equation:(61)X=2DBt.In general,t > 10τ. For simplicity, take t = 1 s. Based on the above analysis, the net displacement X of a 0.1 µm diameter droplet in air is given as 2.3239 × 10−5 m, that is, 23.24 µm. The net distance of the particle by Brownian diffusion alone is actually very small, but the overall tendency is to diffuse from higher to lower concentrations.


## 5. Prediction of Microbial Aerosol Movement

Because of the complexity of CFD simulations, a straightforward technique of motion calculation is presented in this section in order to make the analysis more straightforward. Based on the results of the previous analysis, we would disregard the problem of gas-liquid coupling as well as the collisional adhesion phenomenon between the particles in the current scenario. For an individual droplet, the equation of motion could be written as follows using Newton's second law and above thorough examination of forces:(62)mpdVp⟶dt=∑F=Fd+Fg+Fb+Fm+FB+FS+FBi,where V⟶p is the initial velocity vector of the particle, V = (u,v,w) is the ambient wind speed, and V⟶p=up,vp,wp) is the particle velocity. Then, the component form of the equation is(63)dupdt=Fdx+Fmx+FBi+FBx,(64)dvpdt=Fdy+Fmy+FBi+FBy,(65)dwpdt=Fg+Fb+Fdz+Fmz+FBi+FBz+FS.The velocity components up, vp, and wp are the three unknowns in these three equations. Integration of a droplet particle's velocity function yields its displacement. The Reynolds number Rep=dpu−upρa/μa is used to determine the drag coefficient Cd. Clearly, the differential equation for the motion of a droplet in air is difficult to solve analytically because of its complexity. For the computation of differential equations, numerical methods are typically employed.


### 5.1. RK Method

An ordinary differential equation's initial value problem can be summarized as follows:(66)y′=fx,y,x0=a≤x≤b,yx0=y0.Over a series of discrete nodes, the so-called numerical solution method seeks an approximationyi(i = 1,2, ..., n) to the value y(xi) of the solution of the following equation:(67)x0=a<x1<x2<…<xn−1<xn=b.The distancehi = xi − xi−1 between two adjacent nodes is known as the step from node xi−1 to node xi. The steps hi (i = 1,2, …, n) can either be equal or not equal at all.The method is “stepwise,” which means the solution process proceeds step by step in the order in which the nodes are arranged. To describe this type of algorithm, the recursive formula foryi+1 can be calculated given the known information yi,yi−1,yi−2,...,y0. The recursive formulas are usually divided into two categories: those that use only the value of the previous point yi in the calculation of yi+1, called the single-step method and those that use the value of the k points yi,yi−1,yi−2,...,yi−k+1 before yi+1, called the k-step method.The first step in the numerical solution is to try to eliminate its derivative term, a process known as discretization. As the difference is an approximation to the differentiation, the basic way to actually discretize is to replace the derivative with the difference quotient. Examining the difference quotientyxi+1−yxi/h, according to the mean value theorem, there exists a point ξ ∈ (xi,xi+1) such that(68)yxi+1−yxih=y′ξ,thus using the differential equation y0 = f(x,y) to obtain the following relation:(69)yxi+1=yxi+hfξ,yξ,where k∗≜fξ,yξ is called the average slope on the interval [xi,xi+1], so that whenever an algorithm is provided for the average slope, a format for calculating it is derived accordingly from equation (70). Numerical integration of ordinary differential equations can be simplified using the Euler method [21], which is the simplest RK method. The famous Euler formula,(70)yi+1=yi+hfxi,yi,i=0,1,2,…,n,simply takes the slope value k1 = f(xi,yi) at a point xi as the average slope k∗. The Taylor expansion is usually used to discuss local truncation errors in numerical methods. For the Euler method, the local truncation error is the error that is made in a single step, and it is denoted by O(h2).The modified Euler’s scheme is(71)k1=fxi,yi,k2=fxi+h,yi+hk1,yi+1=yi+h2k1+k2.It takes the arithmetic average of the slope values k1 and k2 at two pointsxi and xi+1 as an arithmetic average for the slope k, while the slope value k2 at xi+1 is forecast using known information yi by the Eulerian approach. The local truncation error of (72) is O(h3).Euler’s method and the modified Euler’s method reveal that if one manages to take the slope values of a few more points, saym points, in the interval [xi, xi+1] and then, weight them as an average slope k∗, it is possible to construct a computational format with higher accuracy such that the local truncation error is O(hm+1). This is the fundamental concept behind the RK method.For the numerical calculations, the fourth-order RK format with fourth-order accuracy is used as the primary numerical calculation method. In order to solve the original differential equations, the ode45 function in MATLAB was called. Rather than directly calling the ode45 function in this article, we write the code for each step of the fourth-order RK method (fixed step). Prior to solving, the initial iteration values, such as the droplet displacement and velocity at rest and the wind speed, need to be determined. It is also necessary to determine the step size of the iterative process(h) as well as the time span of the droplets in the process. The following are the specific equations for the fourth-order RK method:(72)ki1=fitj,y1j,y2j,…,ynj,ki2=fitj+h2,y1j+hk112,…,ynj+hkn12,ki3=fitj+h2,y1j+hk122,…,ynj+hkn22,ki4=fitj+h,y1j+hk13,…,ynj+hkn3,yi,j+1=yi,j+h6ki1+2ki2+2ki3+ki4,i=1,2,...,n,where i denotes the equation entry in the set of differential equations and j denotes the time step.


### 5.2. Results and Discussion

To simulate the dynamic mechanisms of exhaled droplets, we use initial particle (V⟶p) and environmental (V) velocities, as well as initial size distributions (dp), which were derived from reports in the experimental literature. We looked at the effect of different droplet diameters, jet angles, and initial velocities on the trajectory of a sneezing droplet in three directions [14]. The comparison between the Maxley–Riley equation and our model is the focus of this article.Modified Maxey–Riley equations [22] model the motion of a rigid sphere in an uncompressible flow:(73)ρp⋅dVp⟶dt=ρa⋅DV⟶Dt+ρp−ρag,−18μaρadp2Vp⟶−V⟶−124dp2⋅∇2V⟶−12ρadVp⟶dt−DDtV⟶+140dp2⋅∇2V⟶,−9ρadpμaπ∫0t1t−sdVp⟶sds−ddsV⟶+124dp2⋅∇2V⟶ds.Here, the material derivative of velocity is denoted as(74)Dy⟶dt=∂y⟶∂t+V⟶∂y⟶∂x⟶.The individual force listed in separate terms on the right-hand side of (73) has the following physical meaning: force on particles from undisturbed flow, gravity and buoyancy force, Stokes drag, added mass term, and the Basset–Boussinesq memory term. Similarly, equation (73) was modelled using the fourth-order RK method. Substituting the parameters mentioned earlier in this article into the material derivative and Laplace operator gives(75)V⟶=0.400.44z100.2,∇2V⟶=∂2V⟶∂2x+∂2V⟶∂2y+∂2V⟶∂2z=00−0.444194z−1.8,(76)DV⟶Dt=∂V⟶∂t+∇V⟶V⟶=000000000.0555242z−0.60.400.44z100.2=000.0154147z−0.6,(77)DDtV⟶+140dp2⋅∇2V⟶=DDt00−0.00111049dp2z−1.8+0.277621z0.2=000.00055493dp2z−2.6+0.0154147z−0.6.The terms includingdp2∇2V⟶ are referred to as the Faux´en corrections [23]. To simplify the operation, we omit consideration of the terms of the particles when their diameters are very small, which is obviously satisfied for saliva droplets.Given the ambient velocity and the initial velocity of the particle, we find that the velocity and displacement of the particle in they-direction under the Maxey–Riley equation are both 0. For a better comparison, the same particle with a diameter of 380 µm and an initial velocity of 30 m/s is taken. The displacement of the droplet in the x and z-directions is shown below.Figure6 illustrates that the particle reaches a maximum distance of 2.3 m horizontally under the Maxey–Riley equation, which is 0.8 m further than the maximum distance obtained in our model. According to Table 1, drag force (Fd), gravity (Fg), and Saffman lift (FS) are significant in the motion of particles in equations (63)–(65), while gravity, Stokes drag, and Basset force have the greatest impact on particles in the Maxey–Riley equation. The two equations describing the particle trajectory lead to a difference of nearly 1m in the final(maximum) displacement of the particle. Even for moderate sneezing, the result given by the Maxey–Riley equation has exceeded the safe social distance(2m).Figure 6
Maximum displacement of particle motion under the Maxey–Riley equation.Table 1
An overview of the relative magnitudes of the forces acting on moving droplets.
ForcesFdFthFgFbFmFBFSFBiOrders of magnitude(N)10−810−1710−1210−1510−17<10−12<10−1110−15Two equations ((62) and (73)) give different results due to the different perspectives considered for the force analysis of the particle and the different expressions for the same force. The calculation of the drag force (Stokes force) is found to be the primary cause of the equations' differences. The drag force is proportional to the square of the relative velocity of the particle and the surrounding air, such that the relative velocity is the main factor affecting the force on particles. For Maxey–Riley equation, if the Stokes drag is multiplied by |ua − up|, the horizontal distance for particle settling at 380 µm is close to 1.5 m. The following figure shows the displacement and time difference of the two equations in the horizontal direction.As can be seen from Figure7, the 380 µm particle settles to the ground within 1 second under both equations of the simulation. At approximately 0.62 s, the horizontal direction will arrive at 1.45 m according to the estimate of the Maxey–Riley equation, while the particles under the action of model 1 are still floating in the air at this time and will eventually arrive at 1.5 m at 0.75 s. It can also be seen that the particle at the same initial velocity has greater inertia under the action of model 1 and therefore travels a relatively greater distance just after being ejected. At about 0.52 s, the horizontal displacement of the particle under the Maxey–Riley equation exceeds the particles of model 1, but it is finally surpassed. This physical model 1 is validated by comparison with the Maxey–Riley model, which yields relatively small and acceptable errors. As to which model is more accurate, more experimental data and numerical results will be needed. However, with this comparison, we can conclude that the present model can describe the motion of a virus-containing droplet.Figure 7
Variation of the horizontal maximum displacement with time for different models. Model 1 indicates equation (63), and model 2 represents the Maxey–Riley equation (5.12).
(a)(b)


## 5.1. RK Method

An ordinary differential equation's initial value problem can be summarized as follows:(66)y′=fx,y,x0=a≤x≤b,yx0=y0.Over a series of discrete nodes, the so-called numerical solution method seeks an approximationyi(i = 1,2, ..., n) to the value y(xi) of the solution of the following equation:(67)x0=a<x1<x2<…<xn−1<xn=b.The distancehi = xi − xi−1 between two adjacent nodes is known as the step from node xi−1 to node xi. The steps hi (i = 1,2, …, n) can either be equal or not equal at all.The method is “stepwise,” which means the solution process proceeds step by step in the order in which the nodes are arranged. To describe this type of algorithm, the recursive formula foryi+1 can be calculated given the known information yi,yi−1,yi−2,...,y0. The recursive formulas are usually divided into two categories: those that use only the value of the previous point yi in the calculation of yi+1, called the single-step method and those that use the value of the k points yi,yi−1,yi−2,...,yi−k+1 before yi+1, called the k-step method.The first step in the numerical solution is to try to eliminate its derivative term, a process known as discretization. As the difference is an approximation to the differentiation, the basic way to actually discretize is to replace the derivative with the difference quotient. Examining the difference quotientyxi+1−yxi/h, according to the mean value theorem, there exists a point ξ ∈ (xi,xi+1) such that(68)yxi+1−yxih=y′ξ,thus using the differential equation y0 = f(x,y) to obtain the following relation:(69)yxi+1=yxi+hfξ,yξ,where k∗≜fξ,yξ is called the average slope on the interval [xi,xi+1], so that whenever an algorithm is provided for the average slope, a format for calculating it is derived accordingly from equation (70). Numerical integration of ordinary differential equations can be simplified using the Euler method [21], which is the simplest RK method. The famous Euler formula,(70)yi+1=yi+hfxi,yi,i=0,1,2,…,n,simply takes the slope value k1 = f(xi,yi) at a point xi as the average slope k∗. The Taylor expansion is usually used to discuss local truncation errors in numerical methods. For the Euler method, the local truncation error is the error that is made in a single step, and it is denoted by O(h2).The modified Euler’s scheme is(71)k1=fxi,yi,k2=fxi+h,yi+hk1,yi+1=yi+h2k1+k2.It takes the arithmetic average of the slope values k1 and k2 at two pointsxi and xi+1 as an arithmetic average for the slope k, while the slope value k2 at xi+1 is forecast using known information yi by the Eulerian approach. The local truncation error of (72) is O(h3).Euler’s method and the modified Euler’s method reveal that if one manages to take the slope values of a few more points, saym points, in the interval [xi, xi+1] and then, weight them as an average slope k∗, it is possible to construct a computational format with higher accuracy such that the local truncation error is O(hm+1). This is the fundamental concept behind the RK method.For the numerical calculations, the fourth-order RK format with fourth-order accuracy is used as the primary numerical calculation method. In order to solve the original differential equations, the ode45 function in MATLAB was called. Rather than directly calling the ode45 function in this article, we write the code for each step of the fourth-order RK method (fixed step). Prior to solving, the initial iteration values, such as the droplet displacement and velocity at rest and the wind speed, need to be determined. It is also necessary to determine the step size of the iterative process(h) as well as the time span of the droplets in the process. The following are the specific equations for the fourth-order RK method:(72)ki1=fitj,y1j,y2j,…,ynj,ki2=fitj+h2,y1j+hk112,…,ynj+hkn12,ki3=fitj+h2,y1j+hk122,…,ynj+hkn22,ki4=fitj+h,y1j+hk13,…,ynj+hkn3,yi,j+1=yi,j+h6ki1+2ki2+2ki3+ki4,i=1,2,...,n,where i denotes the equation entry in the set of differential equations and j denotes the time step.


## 5.2. Results and Discussion

To simulate the dynamic mechanisms of exhaled droplets, we use initial particle (V⟶p) and environmental (V) velocities, as well as initial size distributions (dp), which were derived from reports in the experimental literature. We looked at the effect of different droplet diameters, jet angles, and initial velocities on the trajectory of a sneezing droplet in three directions [14]. The comparison between the Maxley–Riley equation and our model is the focus of this article.Modified Maxey–Riley equations [22] model the motion of a rigid sphere in an uncompressible flow:(73)ρp⋅dVp⟶dt=ρa⋅DV⟶Dt+ρp−ρag,−18μaρadp2Vp⟶−V⟶−124dp2⋅∇2V⟶−12ρadVp⟶dt−DDtV⟶+140dp2⋅∇2V⟶,−9ρadpμaπ∫0t1t−sdVp⟶sds−ddsV⟶+124dp2⋅∇2V⟶ds.Here, the material derivative of velocity is denoted as(74)Dy⟶dt=∂y⟶∂t+V⟶∂y⟶∂x⟶.The individual force listed in separate terms on the right-hand side of (73) has the following physical meaning: force on particles from undisturbed flow, gravity and buoyancy force, Stokes drag, added mass term, and the Basset–Boussinesq memory term. Similarly, equation (73) was modelled using the fourth-order RK method. Substituting the parameters mentioned earlier in this article into the material derivative and Laplace operator gives(75)V⟶=0.400.44z100.2,∇2V⟶=∂2V⟶∂2x+∂2V⟶∂2y+∂2V⟶∂2z=00−0.444194z−1.8,(76)DV⟶Dt=∂V⟶∂t+∇V⟶V⟶=000000000.0555242z−0.60.400.44z100.2=000.0154147z−0.6,(77)DDtV⟶+140dp2⋅∇2V⟶=DDt00−0.00111049dp2z−1.8+0.277621z0.2=000.00055493dp2z−2.6+0.0154147z−0.6.The terms includingdp2∇2V⟶ are referred to as the Faux´en corrections [23]. To simplify the operation, we omit consideration of the terms of the particles when their diameters are very small, which is obviously satisfied for saliva droplets.Given the ambient velocity and the initial velocity of the particle, we find that the velocity and displacement of the particle in they-direction under the Maxey–Riley equation are both 0. For a better comparison, the same particle with a diameter of 380 µm and an initial velocity of 30 m/s is taken. The displacement of the droplet in the x and z-directions is shown below.Figure6 illustrates that the particle reaches a maximum distance of 2.3 m horizontally under the Maxey–Riley equation, which is 0.8 m further than the maximum distance obtained in our model. According to Table 1, drag force (Fd), gravity (Fg), and Saffman lift (FS) are significant in the motion of particles in equations (63)–(65), while gravity, Stokes drag, and Basset force have the greatest impact on particles in the Maxey–Riley equation. The two equations describing the particle trajectory lead to a difference of nearly 1m in the final(maximum) displacement of the particle. Even for moderate sneezing, the result given by the Maxey–Riley equation has exceeded the safe social distance(2m).Figure 6
Maximum displacement of particle motion under the Maxey–Riley equation.Table 1
An overview of the relative magnitudes of the forces acting on moving droplets.
ForcesFdFthFgFbFmFBFSFBiOrders of magnitude(N)10−810−1710−1210−1510−17<10−12<10−1110−15Two equations ((62) and (73)) give different results due to the different perspectives considered for the force analysis of the particle and the different expressions for the same force. The calculation of the drag force (Stokes force) is found to be the primary cause of the equations' differences. The drag force is proportional to the square of the relative velocity of the particle and the surrounding air, such that the relative velocity is the main factor affecting the force on particles. For Maxey–Riley equation, if the Stokes drag is multiplied by |ua − up|, the horizontal distance for particle settling at 380 µm is close to 1.5 m. The following figure shows the displacement and time difference of the two equations in the horizontal direction.As can be seen from Figure7, the 380 µm particle settles to the ground within 1 second under both equations of the simulation. At approximately 0.62 s, the horizontal direction will arrive at 1.45 m according to the estimate of the Maxey–Riley equation, while the particles under the action of model 1 are still floating in the air at this time and will eventually arrive at 1.5 m at 0.75 s. It can also be seen that the particle at the same initial velocity has greater inertia under the action of model 1 and therefore travels a relatively greater distance just after being ejected. At about 0.52 s, the horizontal displacement of the particle under the Maxey–Riley equation exceeds the particles of model 1, but it is finally surpassed. This physical model 1 is validated by comparison with the Maxey–Riley model, which yields relatively small and acceptable errors. As to which model is more accurate, more experimental data and numerical results will be needed. However, with this comparison, we can conclude that the present model can describe the motion of a virus-containing droplet.Figure 7
Variation of the horizontal maximum displacement with time for different models. Model 1 indicates equation (63), and model 2 represents the Maxey–Riley equation (5.12).
(a)(b)


## 6. Conclusion

Coughing and sneezing produce airborne droplets, and our primary goal is to investigate the transport mechanisms and fluid dynamics of these droplets in an airborne environment. A simple force analysis of aerosol particles in air, based on SL, is first carried out to determine their SV and stopping distance in air, setting the stage for the complex motion of the particles in the fluid that follows. Next, this article’s focus is on virus-carrying individual droplets in the air, so the forces that be applied to aerosol particles after they have been dispersed in the air are analyzed and their orders of magnitudes are calculated. To help determine how long it takes for microbial aerosols to evaporate into nuclei, the time it takes for particles of different sizes to evaporate is calculated. In order to calculate the horizontal displacement that can be reached by the particles at various starting velocities, diameters, and jet angles, we used the Runge–Kutta method to solve for the contaminated region using the trial data. Then, the discrepancies between our model and the well-known Maxey–Riley equation are then discussed. The two equations' results are combined to calculate safe settling distances and provide theoretical justification for infectious illness control. Finally, the evaporation effect is incorporated into the particle movement model, resulting in a more realistic outcome.



---

*Source: 6708401-2022-05-17.xml*