**2. CFPD models for unconventional lung aerosols transport**

Computational fluid-particle dynamics (CFPD) incorporates computational fluid dynamics (CFD) theories with multiphase flow modeling techniques for particle dynamics in fluid flows, i.e., Euler-Euler and Euler-Lagrange methods [12, 13]. As shown in **Figure 1**, basic steps to proceed simulations of transport and deposition of inhaled aerosol in human respiratory systems using CFPD models include [11–13]: (1) reconstruction of the human respiratory system, (2) mesh generation and independence test, (3) formulation of governing equation system, (4) numerical discretization and solution, (5) model validation and revision, and (6) aerosol transport and deposition visualization and analyses with insights.

**Figure 1.** Basic steps to perform CFPD simulations for transport and deposition of aerosols in human respiratory systems.

The underlying assumptions of conventional CFPD models are that the aerosols are spherical, noninteracting, and constant in particulate diameter [12]. These assumptions induce inaccuracies in local deposition predictions when modeling the dynamics of unconventional aerosol particles, e.g., nonspherical particles and hygroscopic droplets. To accurately capture the transport and deposition of particulate matters, more physical and chemical principles must be considered, such as the rotational motion of nonspherical fibers, particle-particle collisions and coagulations, hygroscopic behaviors of droplets, as well as droplet-vapor phase transition, which strongly impact their trajectories and subsequent deposition patterns in human respiratory systems. Thus, advanced CFPD models have been developed and applied to modeling unconventional particle dynamics [2, 3, 10, 13]. In this section, procedures for CFPD modeling are discussed, followed by introducing governing equation systems of two advanced CFPD models for ellipsoidal fiber dynamics and multicomponent droplet-vapor aerosol transport in human respiratory systems.

#### **2.1. Reconstruction of representative human respiratory systems**

Since geometric variations have strong impacts on the airflow patterns leading to different particle trajectories, the geometric accuracy of human respiratory system is of great importance. Therefore, as the flow domain for CFPD modeling, the human respiratory system geometry must be precise in anatomy and physiology.

#### *2.1.1. Anatomy of human respiratory systems*

i.e., Euler-Euler and Euler-Lagrange methods [12, 13]. As shown in **Figure 1**, basic steps to proceed simulations of transport and deposition of inhaled aerosol in human respiratory systems using CFPD models include [11–13]: (1) reconstruction of the human respiratory system, (2) mesh generation and independence test, (3) formulation of governing equation system, (4) numerical discretization and solution, (5) model validation and revision, and (6)

**Figure 1.** Basic steps to perform CFPD simulations for transport and deposition of aerosols in human respiratory sys-

The underlying assumptions of conventional CFPD models are that the aerosols are spherical, noninteracting, and constant in particulate diameter [12]. These assumptions induce inaccuracies in local deposition predictions when modeling the dynamics of unconventional aerosol particles, e.g., nonspherical particles and hygroscopic droplets. To accurately capture the transport and deposition of particulate matters, more physical and chemical principles must be considered, such as the rotational motion of nonspherical fibers, particle-particle collisions and coagulations, hygroscopic behaviors of droplets, as well as droplet-vapor phase transition, which strongly impact their trajectories and subsequent deposition patterns in human respiratory systems. Thus, advanced CFPD models have been developed and applied to modeling unconventional particle dynamics [2, 3, 10, 13]. In this section, procedures for CFPD modeling are discussed, followed by introducing governing equation systems of two advanced CFPD models for ellipsoidal fiber dynamics and multicomponent droplet-vapor aerosol

Since geometric variations have strong impacts on the airflow patterns leading to different particle trajectories, the geometric accuracy of human respiratory system is of great impor-

tems.

52 Aerosols - Science and Case Studies

transport in human respiratory systems.

**2.1. Reconstruction of representative human respiratory systems**

aerosol transport and deposition visualization and analyses with insights.

Human respiratory system consists of respiratory passages, lung lobes, and respiratory muscles (see **Figure 2**). Respiration includes three separate but related functions: ventilation (breathing), gas exchange (between air and blood in the lungs and between the blood and tissues), and oxygen utilization (by the tissues in the energy-liberating reactions of cell respiration). During the respiration process, airborne particulate matters are also inhaled and may settle in different locations. The internal walls of the human respiratory system are covered by mucus layers. As a filtering mechanism, mucus layers serve to trap and clear the PM carried by the lung airflow.

**Figure 2.** Schematics of human respiratory system (reprinted from [12]).

The human respiratory system ventilation path contains mouth, nose, pharynx, glottis, larynx, trachea, bronchi, bronchioles (including terminal bronchioles (generation 16 (G16)) and respiratory bronchioles (generation 17–19 (G17–G19)), and alveoli [14]. The path can be divided into two functional zones. The conducting zone includes all of the anatomical structures through which air passes before reaching the respiratory zone, i.e., mouth to terminal bronchioles (see **Figure 2**). The respiratory zone is the region where gas exchange occurs, which includes respiratory bronchioles and alveoli.

Human respiratory system can also be divided into upper airway and lower airway. The upper airway is all structures above the glottis; the lower airway is from the vocal cords to the lung. The upper airway include oral cavity, pharynx, larynx, and upper part of the trachea. Specifically, oral cavity can be considered as having the shape of a prolate spheroid that features different radius of curvatures along the sagittal and coronal planes. The shape of the oropharynx is more circular before the glottal region. In the lower airway, the lung consists of 23 generations, encompassing 223 airways plus millions of alveoli of different shapes. There are about 300 million alveoli (0.25–0.5 mm in diameter) in the lung which provide a high surface area, i.e., 60–80 m2 for diffusion of gases. Each alveolus is only one cell-layer thick, so that the total "air-blood barrier" is only two cells across (an alveolar cell and a capillary endothelial cell), which is approximately 2 μm.

Another way is to divide the human respiratory system into three parts [12]: the extrathoracic region (i.e., nose, mouth, and throat), the tracheobronchial part (i.e., trachea and bronchial tree), and the alveolar region (i.e., alveolar ducts and sacs).

### *2.1.2. Idealized human upper airway geometries*

Idealized geometries are built up of simple geometric shapes but possess all the basic anatomical features of real human respiratory systems. At the early age, due to the limitations of computational power and imaging techniques handling subject-specific geometries, the use of idealized representative upper airway geometries is indispensable to investigate airflow and aerosol deposition patterns [15–19]. Idealized tracheobronchial geometries are usually created following airway dimensions measured by either Weibel [20] or Horsfield et al. [21]. They also follow different sets of measurements [22, 23]. Two popular idealized geometries are shown in **Figure 3(a)** and **(b)** [3, 18]. **Figure 3(a)** contains oral cavity, pharynx, larynx, and a symmetric, planar triple-bifurcation lung airway representing G0–G3 after Cheng et al. [22] and Weibel [20]. **Figure 3(b)** presents another idealized geometry with oral airways, nasal airways, and asymmetric triple bifurcation unit (TBU) representing G0 to G3 following Alberta-Finlay model [22] and Horsfield et al. [21]. Newly developed idealized human upper airway geometries emphasize more on matching average aerosol deposition data measured in different airway replicas [24].

Although idealized geometries are still beneficial for extensive fundamental transport mechanism studies and easier model validations with in-vitro measurements [2, 3, 10, 25], subject-specific respiratory system geometries are necessary to be employed for individualized risk assessment or disease treatment. The importance of airway geometry on micron-particle deposition was investigated and confirmed [26].

#### *2.1.3. Subject-specific human upper airway geometries*

Currently, more CAD-like subject-specific human upper airway geometry files are generated via DICOM file conversion using either specialized software or in-house codes. An example is shown in **Figure 4**. The geometry is from a 20-mm opening mouth inlet, to the tracheobronchial airways up to G9, which has been generated based on CT images. The subject for this study is a 47-year-old healthy male volunteer, 174 cm of height and 78 kg of weight. GE's 64 slice CT scanner was used to take images of 500 × 500 mm (i.e., 512 × 512 pixels on the plane) cross-sections with 2.5 mm image-slice thickness. The images were taken from the extracranial skull base to the abdominal region. The geometries in-between slices (slice thickness: 1.25 mm) were created via an interpolation methodology. Compared with idealized geometries (see **Figure 3(a)** and **(b)**), the subject-specific geometry in **Figure 4** maintains more details of local airways and will have more realistic geometric influence on airflow and aerosol transport.

Human respiratory system can also be divided into upper airway and lower airway. The upper airway is all structures above the glottis; the lower airway is from the vocal cords to the lung. The upper airway include oral cavity, pharynx, larynx, and upper part of the trachea. Specifically, oral cavity can be considered as having the shape of a prolate spheroid that features different radius of curvatures along the sagittal and coronal planes. The shape of the oropharynx is more circular before the glottal region. In the lower airway, the lung consists of 23 generations, encompassing 223 airways plus millions of alveoli of different shapes. There are about 300 million alveoli (0.25–0.5 mm in diameter) in the lung which provide a high surface

total "air-blood barrier" is only two cells across (an alveolar cell and a capillary endothelial

Another way is to divide the human respiratory system into three parts [12]: the extrathoracic region (i.e., nose, mouth, and throat), the tracheobronchial part (i.e., trachea and bronchial

Idealized geometries are built up of simple geometric shapes but possess all the basic anatomical features of real human respiratory systems. At the early age, due to the limitations of computational power and imaging techniques handling subject-specific geometries, the use of idealized representative upper airway geometries is indispensable to investigate airflow and aerosol deposition patterns [15–19]. Idealized tracheobronchial geometries are usually created following airway dimensions measured by either Weibel [20] or Horsfield et al. [21]. They also follow different sets of measurements [22, 23]. Two popular idealized geometries are shown in **Figure 3(a)** and **(b)** [3, 18]. **Figure 3(a)** contains oral cavity, pharynx, larynx, and a symmetric, planar triple-bifurcation lung airway representing G0–G3 after Cheng et al. [22] and Weibel [20]. **Figure 3(b)** presents another idealized geometry with oral airways, nasal airways, and asymmetric triple bifurcation unit (TBU) representing G0 to G3 following Alberta-Finlay model [22] and Horsfield et al. [21]. Newly developed idealized human upper airway geometries emphasize more on matching average aerosol deposition data measured in different

Although idealized geometries are still beneficial for extensive fundamental transport mechanism studies and easier model validations with in-vitro measurements [2, 3, 10, 25], subject-specific respiratory system geometries are necessary to be employed for individualized risk assessment or disease treatment. The importance of airway geometry on micron-particle

Currently, more CAD-like subject-specific human upper airway geometry files are generated via DICOM file conversion using either specialized software or in-house codes. An example is shown in **Figure 4**. The geometry is from a 20-mm opening mouth inlet, to the tracheobronchial airways up to G9, which has been generated based on CT images. The subject for this

for diffusion of gases. Each alveolus is only one cell-layer thick, so that the

area, i.e., 60–80 m2

54 Aerosols - Science and Case Studies

airway replicas [24].

cell), which is approximately 2 μm.

*2.1.2. Idealized human upper airway geometries*

deposition was investigated and confirmed [26].

*2.1.3. Subject-specific human upper airway geometries*

tree), and the alveolar region (i.e., alveolar ducts and sacs).

**Figure 3.** Configurations of two idealized human upper airway geometries: (a) symmetric idealized geometry from mouth to G3 with structured hexahedral finite-volume mesh and (b) asymmetric idealized geometry from mouth and nose to G3 (reprinted from [18]).

Creating high-quality boundary-layer mesh is paramount in computational mesh generation for the flow domain. It consists of multilayer hybrid tetrahedral/pentahedral elements near the wall surface to fully contain the viscous sublayers and to resolve any geometric features present there (see **Figures 3** and **4** as examples). Such high local mesh resolution can also help to accurately calculate values of near-wall derivatives, such as the deposition fluxes. Grid independence can be checked via indicators such as flow field solution and vapor deposition fraction (DF).

Realistic 3D imaging of the entire human respiratory system that is 100% accurate is currently unavailable because: (1) the lung consists of 223 airways plus millions of alveoli, while the resolution of CT/MRI is insufficient to capture geometries in small scales, i.e., airways exceeding G9; (2) in-vivo measurements are difficult because the whole respiratory system geometry is time-dependent according to the respiratory movements.

**Figure 4.** Configuration of a subject-specific human upper airway geometry from mouth to G9 with tetrahedral mesh including near-wall prism layers: (a) a subject-specific human upper airway model; (b) mesh details near the glottis; (c) mesh details at an outlet; (d) sagittal contour with regional names.

#### **2.2. Governing equations**

To provide a numerical description of airflow and aerosol behaviors, the governing equation system for CFPD modeling must be formulated. It contains physical and chemical conservation laws as well as supplementary equations (see **Figure 1**). The nonlinear differential governing equations subject to appropriate initial conditions and boundary conditions (BCs) will be discretized into algebraic equations and solved using different numerical techniques, e.g., finite difference or finite volume method. To promote the understanding of the physical and chemical characteristics of airflow and aerosol motions in lung, this section introduces the procedure to formulate equation systems for CFPD models for both conventional and unconventional aerosol transport and deposition.

#### *2.2.1. Airflow field equations*

The airflow dynamics in the respiratory tract is intrinsically unsteady, driven by the pressure differences under the action of cyclic breathing process. The incompressible Navier-Stokes (N-S) equation is employed to characterize airflow in the human respiratory tract, accompanied by continuity equation, energy equation, and constitutive equations [27], which can be written in tensor form as follows:

#### Continuity equation

$$\frac{\partial u\_j}{\partial \mathbf{x}\_j} = \mathbf{0} \tag{1}$$

in which *uj* represents the fluid velocity and *ρ* is the air density.

Navier-Stokes (N-S) equation

$$
\rho \frac{\partial u\_i}{\partial t} + \rho \frac{\partial (u \mu\_j)}{\partial \mathbf{x}\_j} = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial \tau\_y}{\partial \mathbf{x}\_j} + B\_i \tag{2}
$$

where *p* is the pressure, *Bi* is the body force including gravity, electromagnetic force, etc. The viscous stress tensor *τij* in Eq. (2) is given by:

$$\mathbf{r}\_{\circ} = \mu \left( \frac{\partial \mathbf{u}\_{i}}{\partial \mathbf{x}\_{j}} + \frac{\partial \mathbf{u}\_{j}}{\partial \mathbf{x}\_{i}} \right) \tag{3}$$

In Eq. (3), *μ* is the air viscosity.

Energy equation

**Figure 4.** Configuration of a subject-specific human upper airway geometry from mouth to G9 with tetrahedral mesh including near-wall prism layers: (a) a subject-specific human upper airway model; (b) mesh details near the glottis; (c)

To provide a numerical description of airflow and aerosol behaviors, the governing equation system for CFPD modeling must be formulated. It contains physical and chemical conservation laws as well as supplementary equations (see **Figure 1**). The nonlinear differential governing equations subject to appropriate initial conditions and boundary conditions (BCs) will be discretized into algebraic equations and solved using different numerical techniques, e.g., finite difference or finite volume method. To promote the understanding of the physical and chemical characteristics of airflow and aerosol motions in lung, this section introduces the procedure to formulate equation systems for CFPD models for both conventional and uncon-

The airflow dynamics in the respiratory tract is intrinsically unsteady, driven by the pressure differences under the action of cyclic breathing process. The incompressible Navier-Stokes (N-S) equation is employed to characterize airflow in the human respiratory tract, accompanied by continuity equation, energy equation, and constitutive equations [27], which can be written

> 0 *<sup>j</sup> j u x*

¶ <sup>=</sup> ¶ (1)

mesh details at an outlet; (d) sagittal contour with regional names.

ventional aerosol transport and deposition.

**2.2. Governing equations**

56 Aerosols - Science and Case Studies

*2.2.1. Airflow field equations*

in tensor form as follows:

Continuity equation

$$
\rho \frac{\partial \left(c\_{\rho} T\right)}{\partial t} + \rho \frac{\partial (c\_{\rho} \mu\_{j} T)}{\partial \mathbf{x}\_{j}} = \frac{\partial}{\partial \mathbf{x}\_{j}} \left[k \frac{\partial T}{\partial \mathbf{x}\_{j}}\right] + \Phi + S\_{\mathrm{T}} \tag{4}
$$

where *k* is thermal conductivity and *cp* is specific heat capacity. Furthermore, *ST* is the thermal sink or source term due to radiation, chemical reaction, phase change, etc. *Φ* is the dissipation function which can be given as:

$$
\Phi = \pi\_{\psi} \frac{\partial u\_i}{\partial \mathbf{x}\_j} \tag{5}
$$

The equation of state as stated by the ideal gas law should complement the above equations for model closure, i.e., *pVm* = *RT*, where *Vm* is the molar volume and *R* is the ideal gas constant.

At typical breathing rates, airflow through the oral airway region and first few generations is incipient turbulent. It becomes laminar at G6-8 and remains so thereafter. In order to capture the airflow structures in the laminar-to-turbulent flow regimes in human upper airways under common inhalation flow rates, the low-Reynolds-number (LRN) k-ω model and shear stress transport (SST) transition model are selected and adapted [27–32], based on their overall performance in predicting the onset of "laminar-to-turbulent" transition, their computational efficiency and reasonable accuracy as compared with large eddy simulations (LESs) [30]. For these Reynolds averaged Navier-Stokes (RANS) models, only averaged values of velocity, pressure, and other turbulence variables are calculated. The instantaneous fluctuation components can be recovered, being approximated by the eddy-interaction model (EIM) [33, 34].

#### *2.2.2. Particle dynamics equations (Lagrangian phase)*

For any given inlet concentration of effectively spherical particles (i.e., conventional aerosol particles), a Lagrangian frame of reference for the trajectory-computations can be employed. In light of the large particle-to-air density ratio and negligible thermophoretic forces, the generalized particle trajectory equation for particle *i* can be described as [35]:

$$\frac{d}{dt}(m\_p \vec{u}\_i^{\rho}) = \vec{F}\_i^D + \vec{F}\_i^L + \vec{F}\_i^{BM} + \vec{F}\_i^G + \vec{F}\_i^{\text{inensation}} \tag{6}$$

Here, and *mp* are the velocity and mass of the particle, respectively; represents the drag force which can be given by:

$$\vec{F}\_i^D = \frac{1}{2} \rho C\_D A\_p (\vec{u}\_i - \vec{u}\_i^p) \mid \vec{u}\_i - \vec{u}\_i^p \mid \tag{7}$$

in which *Ap* is the projected particle area and *CD* is the particle drag coefficient [3, 4, 10]. In Eq. (6), is the gravity [35], is the Brownian motion-induced force [3, 10, 36], is the Saffman lift force [37], and interaction is the force due to particle-particle interactions. The interactions may consist of van der Waals interaction and electrostatic interaction. The effect of the lubrication force, or near-wall drag modifications, may also need to be considered. All forces in Eq. (6) are point forces acting at the center of the particle, which employ specific empirical correlations corresponding to different particle aerodynamic properties.

#### *2.2.3. Advection-diffusion equation for vapor (Eulerian phase)*

For vapors and nanoparticles less than 50 nm in diameter, their inertial effects can be neglected. Thus, advection-diffusion equation can be used to describe vapor mass transfer which can be written as:

$$\frac{\partial Y}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} \left( \mu\_j Y \right) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \tilde{D}\_a + \tilde{D}\_{\text{subbalance}} \right) \frac{\partial Y}{\partial \mathbf{x}\_j} \right] + S\_Y \tag{8}$$

where *Y* is the vapor mass fraction, is the vapor molecular diffusivity in air, is the turbulence diffusivity which in laminar flow is zero, and *SY* is the source term.

#### *2.2.4. Supplementary equations for modeling unconventional aerosol particles*

components can be recovered, being approximated by the eddy-interaction model (EIM) [33,

For any given inlet concentration of effectively spherical particles (i.e., conventional aerosol particles), a Lagrangian frame of reference for the trajectory-computations can be employed. In light of the large particle-to-air density ratio and negligible thermophoretic forces, the

> interaction ( ) == + + + + rrr rr <sup>r</sup> *p D L BM G pi i i i i i <sup>d</sup> mu F F F F F*

> > <sup>1</sup> ( - )| - | <sup>2</sup> *<sup>D</sup> p p F CA u u u u i Dp i i i i* =

in which *Ap* is the projected particle area and *CD* is the particle drag coefficient [3, 4, 10].

actions. The interactions may consist of van der Waals interaction and electrostatic interaction. The effect of the lubrication force, or near-wall drag modifications, may also need to be considered. All forces in Eq. (6) are point forces acting at the center of the particle, which employ specific empirical correlations corresponding to different particle aerody-

For vapors and nanoparticles less than 50 nm in diameter, their inertial effects can be neglected. Thus, advection-diffusion equation can be used to describe vapor mass transfer which can be

> ( *<sup>j</sup>* ) ( *a turbulence* ) *<sup>Y</sup> j j j Y Y uY D D <sup>S</sup> tx x x* ¶¶ ¶ é ù ¶ + =+ + ê ú ¶¶ ¶ ¶ ë û

where *Y* is the vapor mass fraction, is the vapor molecular diffusivity in air, is

the turbulence diffusivity which in laminar flow is zero, and *SY* is the source term.

*dt* (6)

<sup>r</sup> rr rr (7)

is the Brownian motion-induced force [3, 10, 36],

interaction is the force due to particle-particle inter-

% % (8)

represents the drag

generalized particle trajectory equation for particle *i* can be described as [35]:

and *mp* are the velocity and mass of the particle, respectively;

r

*2.2.2. Particle dynamics equations (Lagrangian phase)*

is the gravity [35],

*2.2.3. Advection-diffusion equation for vapor (Eulerian phase)*

is the Saffman lift force [37], and

34].

58 Aerosols - Science and Case Studies

Here,

In Eq. (6),

namic properties.

written as:

force which can be given by:

The governing equations for conventional lung aerosol dynamics discussed above are for spherical particles with constant diameters. They are not sufficiently accurate to predict the transport and deposition of unconventional aerosol particles, i.e., nonspherical fiber and hygroscopic multicomponent droplets interacting with ambient vapors. Therefore, supplementary equations need to be introduced to encompass more physical and chemical mechanisms based on first principles [2–4, 10, 38].

#### *2.2.4.1. Translational and rotational motion equations for ellipsoidal fiber*

Based on the point-mass and point-force assumptions adapted by the conventional CFPD models, particles are in isotropic shape and their rotational motions are neglected. However, for nonspherical fibers with anisotropic shapes, the orientation determined by their rotational motions greatly influences the magnitudes and directions of the forces they experience in translational motion (see Eq. (6)). These forces will differ from those exerted on spherical particles with the same aerodynamic diameters. In order to capture this effect and provide accurate prediction of ellipsoidal fiber trajectory (see **Figure 5(a)** and **(b)**) based on realistic physical and chemical mechanisms, the conventional CFPD model was modified as follows [10, 13, 27]:


**Figure 5.** Sketches for ellipsoidal fiber modeling: (a) sketch of geometric parameters: semi-minor axis *ap*, semi-major axis *bp*, and aspect ratio *β*; (b) sketch of coordinate systems: space-fixed frame *xyz*, body-fixed frame *x*′*y*′*z*′, and co-moving frame *x″y″z″*.

The force and torque correlations for ellipsoidal fibers are documented in [10, 13]. Euler's rotation equations for ellipsoidal fiber orientation calculation in the body-fixed *x*′*y*′*z*′-frame (see **Figure 5(b)**) are given by:

$$I\_{\mathbf{x}^\*} \frac{d\alpha\_{\mathbf{x}^\*}}{dt} - \alpha\_{\mathbf{y}^\*} \alpha\_{\mathbf{z}^\*} \left(I\_{\mathbf{y}^\*} - I\_{\mathbf{z}^\*}\right) = T\_{\mathbf{x}^\*} \tag{9a}$$

$$I\_{\mathbf{y}^\*} \frac{d\alpha\_{\mathbf{y}^\*}}{dt} - \alpha\_{\mathbf{z}^\*} \alpha\_{\mathbf{x}^\*} (I\_{\mathbf{z}^\*} - I\_{\mathbf{x}^\*}) = T\_{\mathbf{y}^\*} \tag{9b}$$

$$I\_{z^\cdot} \frac{d\,\alpha\_{z^\cdot}}{dt} - \alpha\_{z^\cdot} \alpha\_{y^\cdot} \left(I\_{x^\cdot} - I\_{y^\cdot}\right) = T\_{z^\cdot} \tag{9c}$$

Here, (*Ix* ′, *Iy* ′, *Iz* ′) are the particle moments of inertia about the principal axes *x*′, *y*′, *z*′; (*ωx* ′, *ωy* ′, *ωz* ′) are the particle angular velocities with respect to the principal axes *x*′, *y*′, *z*′; and (*Tx* ′, *Ty* ′, *Tz* ′) are the torques acting on the particle with respect to the principal axes *x*′, *y*′, *z*′. For ellipsoidal particles, (*Ix* ′, *Iy* ′, *Iz* ′) can be written as:

$$I\_{\mathbf{x}^\*} = I\_{\mathbf{y}^\*} = \frac{\left(1 + \boldsymbol{\beta}^2\right) \cdot \boldsymbol{a}\_p^2}{\mathfrak{S}} \boldsymbol{m}\_p \tag{10a,b}$$

$$I\_{z.} = \frac{2a\_p^2}{\mathfrak{S}} m\_p \tag{10c}$$

where *ap*, *bp*, and *β* are semi-minor axis, semi-major axis, and aspect ratio of the ellipsoidal fiber, respectively (see **Figure 5(a)**). Details of this advanced CFPD model, i.e., Euler-Lagrange method enhanced by Euler's rotation equation (EL-ER), are discussed by Feng and Kleinstreuer [10]. With different correlations of forces and torques acting on particles with specific anisotropic shapes, EL-ER method can be extended to simulate the transport and deposition of other nonspherical aerosol particles.

#### *2.2.4.2. Governing equations for transport and phase change of droplet-vapor mixture*

For liquid aerosols (i.e., droplet-vapor mixture), the phase change between droplet and vapor (see **Figure 6**) significantly affects their dispersion and deposition because the induced inertia changes during their transport. Therefore, the multicomponent mixture plus discrete-droplet (MCM-DD) model has been developed to simultaneously simulate multicomponent dropletvapor interactions with evaporation and condensation effects during their transport in human respiratory systems [2, 3].

Computational Fluid-Particle Dynamics Modeling for Unconventional Inhaled Aerosols in Human Respiratory Systems http://dx.doi.org/10.5772/65361 61

**Figure 6.** Schematic of droplet-vapor interactions and phase change dynamics.

The force and torque correlations for ellipsoidal fibers are documented in [10, 13]. Euler's rotation equations for ellipsoidal fiber orientation calculation in the body-fixed *x*′*y*′*z*′-frame

( ) '

(9a)

(9b)

(9c)

= = (10a,b)

*<sup>a</sup> I m* <sup>=</sup> (10c)

 - -= w w

'' ' ' '

( ) '

 - -= w w

'' ' ' '

( ) '

 - -= w w

'' ' ' '

Here, (*Ix* ′, *Iy* ′, *Iz* ′) are the particle moments of inertia about the principal axes *x*′, *y*′, *z*′; (*ωx* ′, *ωy* ′, *ωz* ′) are the particle angular velocities with respect to the principal axes *x*′, *y*′, *z*′; and (*Tx* ′, *Ty* ′, *Tz* ′) are the torques acting on the particle with respect to the principal axes *x*′, *y*′, *z*′. For

( ) 2 2

*p*

5

2

where *ap*, *bp*, and *β* are semi-minor axis, semi-major axis, and aspect ratio of the ellipsoidal fiber, respectively (see **Figure 5(a)**). Details of this advanced CFPD model, i.e., Euler-Lagrange method enhanced by Euler's rotation equation (EL-ER), are discussed by Feng and Kleinstreuer [10]. With different correlations of forces and torques acting on particles with specific anisotropic shapes, EL-ER method can be extended to simulate the transport and deposition

For liquid aerosols (i.e., droplet-vapor mixture), the phase change between droplet and vapor (see **Figure 6**) significantly affects their dispersion and deposition because the induced inertia changes during their transport. Therefore, the multicomponent mixture plus discrete-droplet (MCM-DD) model has been developed to simultaneously simulate multicomponent dropletvapor interactions with evaporation and condensation effects during their transport in human

2 5 *p z p*

1

'

*2.2.4.2. Governing equations for transport and phase change of droplet-vapor mixture*

*x y p <sup>a</sup> II m* + × b

*z zy x y z <sup>d</sup> I II T*

*y zx z x y*

*I II T*

*x yz y z x <sup>d</sup> I II T*

'

'

'

ellipsoidal particles, (*Ix* ′, *Iy* ′, *Iz* ′) can be written as:

of other nonspherical aerosol particles.

respiratory systems [2, 3].

*d*

*x*

*y*

*z*

' '

*dt* w

*dt* w

*dt* w

(see **Figure 5(b)**) are given by:

60 Aerosols - Science and Case Studies

#### *2.2.4.3. Governing equations for continuous phase (air-vapor mixture)*

Air-vapor mixture is described as a single continuous phase. Therefore, the conservation laws for mass, momentum, and energy are used to describe the air-vapor mixture transport phenomena and the advection-diffusion equation for vapor species transport (see Eqs. (1)–(5) and Eq. (8)). Droplet-vapor interaction will be realized by: (a) introducing source terms into the energy equation and vapor species transport equation (see Eqs. (4) and (8)) and (b) employing the local vapor mass fraction in the droplet mass conservation equation (see Eq. (15)). Specifically, the energy equation of air-vapor mixture can be stated as:

$$\mathbf{S}\_{v \sim d, s}^{(T)} = \int\_{t\_{i, \text{max}}}^{t\_{i, \text{max}}} \left( \sum\_{i=1}^{N\_{d, \text{all}}} \left( \overline{j\_s} A\_d \right)\_i \right) dt\_d \wedge \left( V\_{\text{call}} \Delta t\_f \right) \tag{11}$$

where subscript *a−v* means "air-vapor" and *N* represents the number of volatile components in droplets. The energy source term is due to the latent heat of evaporation or condensation which is released or absorbed by the droplets per local mesh cell:

$$\mathbf{S}\_{v-d}^{(E)} = \sum\_{s=1}^{N} \mathbf{S}\_{v-d,s}^{(E)} = \left\{ \sum\_{i=1}^{N\_{d,sd}} \left[ \left( \sum\_{s=1}^{N} L\_s \overline{j\_s} \right) A\_d \right] \right\} / V\_{coll} \tag{12}$$

where *Nd*,*cell* is the total droplet number in a specified mesh cell, subscript *v−d* means "between vapor and droplet", is the average evaporation/condensation mass flux normal to the droplet surface, and *Ls* is the latent heat of liquid-vapor phase transition of the *s*th species.

Furthermore, the advection-diffusion equation of the *s*th vapor species can be given as:

$$\frac{\partial(\rho\_{a-v}Y\_{v,s})}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_{j}}(\rho\_{a-v}\boldsymbol{\mu}\_{j}Y\_{v,s}) = \frac{\partial}{\partial \mathbf{x}\_{j}}\left|\,\rho\_{a-v}\tilde{D}\_{a-v,s}\frac{\partial Y\_{v,s}}{\partial \mathbf{x}\_{j}}\right| + S\_{v=d,s}^{(Y)}\tag{13}$$

Here, the local vaporized/condensed vapor mass flow rate of the aerosol components is added to its advection-diffusion equation as a source term , (kg m−3 s−1) which is given by:

$$\frac{\partial \left( \rho c\_{\rho} T \right)\_{a-\upsilon}}{\partial t} + \frac{\partial \left( \rho c\_{\rho} u\_{j} T \right)\_{a-\upsilon}}{\partial \mathbf{x}\_{j}} = \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ k\_{a-\upsilon} \frac{\partial T\_{a-\upsilon}}{\partial \mathbf{x}\_{j}} \right] + \Phi + \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \sum\_{s=1}^{N} h\_{s} \mathbf{p}\_{a-\upsilon} \tilde{D}\_{a-\upsilon} \frac{\partial Y\_{s}}{\partial \mathbf{x}\_{j}} \right] + S^{(E)}\_{\upsilon=4} \tag{14}$$

in which *Ad* is the droplet surface area. Additionally, *dtd* represents the droplet phase time differential and Δ*tf* is the flow time-step length. It is worth mentioning that Eq. (14) describes the major physics of liquid-vapor mass change.

#### *2.2.4.4. Governing equations for multicomponent droplets*

The Lagrangian approach is chosen to track multicomponent droplets, neglecting rotational motion. The governing equations for discrete droplets are the translational equation (see Eq. (6)) as well as the mass and energy conservations of droplets and the supplementary equations. The mass conservation of droplets can be expressed as:

$$\frac{dm\_d}{dt} = -\sum\_{s=1}^{m} \int\_{surf} j\_s dA \approx -\sum\_{s=1}^{m} (\overline{j\_s} \cdot A\_d) \tag{15}$$

where the average evaporation/condensation mass flux is given by [39]:

$$\overline{j\_s} = \rho\_{a-v} C\_m \left( Sh \cdot \tilde{D}\_{a-v,s} / d\_d \right) \ln \left[ \left( 1 - Y\_{s,\text{coll}} \right) / \left( 1 - Y\_{s,\text{surf}} \right) \right] \tag{16}$$

in which *Ys*,*vsurf* and *Ys*,*cell* are the mass fractions of the *s*th vapor phase at the droplet surface and the cell center that the droplet is currently in. Eq. (16) assumes that the droplet radius is much smaller than the distance between the droplet mass center and the mesh cell center. Specifically, *Ys*,*cell* is determined by Eq. (13), which realizes the droplet-vapor interaction modeling.

The energy conservation of droplets is given by:

$$m\_d \mathcal{C}\_{p,d} \left( dT\_d / dt \right) = C\_T \left( \frac{k\_{a-v} \cdot Nu}{d\_d} \right) \left( T\_{a-v, \text{all}} - T\_d \right) A\_d - \left( \sum\_{s=1}^N L\_s \overline{j\_s} \right) A\_d \tag{17}$$

where *Nu* is the Nusselt number.

,

surface, and *Ls* is the latent heat of liquid-vapor phase transition of the *s*th species.

Furthermore, the advection-diffusion equation of the *s*th vapor species can be given as:

( ) *a v vs v s <sup>Y</sup>*

*Y Y*

¶ ¶ ¶ ¶ é ù += + ê ú ¶¶ ¶ ¶ ê ú ë û

*tx x x*

r


( ) , , ( )

*j jj*

Here, the local vaporized/condensed vapor mass flow rate of the aerosol components is added

to its advection-diffusion equation as a source term , (kg m−3 s−1) which is given by:

( ) ( ) ( )

*j jj j s cT cuT <sup>T</sup> k hS*

¶r ¶r ¶¶ ¶ éù é ù <sup>+</sup> <sup>=</sup> êú ê ú +F+ r + ¶ ¶ ¶¶ ¶ û ë ¶ ë û

*<sup>N</sup> <sup>p</sup> p j a v a v a v <sup>E</sup>*

in which *Ad* is the droplet surface area. Additionally, *dtd* represents the droplet phase time

The Lagrangian approach is chosen to track multicomponent droplets, neglecting rotational motion. The governing equations for discrete droplets are the translational equation (see Eq. (6)) as well as the mass and energy conservations of droplets and the supplementary equations.

1 1


*j C Sh D d Y Y s a v m a vs d*

*<sup>s</sup> s d surf s s dm j dA j A dt* = =

*m m*

*v d v ds s s d cell s is i S S* - - *Lj A V* = ==

where *Nd*,*cell* is the total droplet number in a specified mesh cell, subscript *v−d* means "between

*N N Nd cell*

/

is the average evaporation/condensation mass flux normal to the droplet

, ,,

1


*a v s av v*

*<sup>Y</sup> <sup>D</sup> <sup>x</sup> xx*

is the flow time-step length. It is worth mentioning that Eq. (14) describes

( )

=- »- × å å ò (15)

is given by [39]:

ë û % (16)

% (13)

*s a v j*

¶ å % (14)


*d*

*a v j vs a v a vs v ds*

*u Y D S*


 r

ì ü ï ï é ù æ ö = = í ý ê ú ç ÷ ï ï î þ ë û è ø å åå (12)

() ()

*E E*

vapor and droplet",

62 Aerosols - Science and Case Studies

differential and Δ*tf*

r

*t xx*

the major physics of liquid-vapor mass change.

*2.2.4.4. Governing equations for multicomponent droplets*

The mass conservation of droplets can be expressed as:

*d*

where the average evaporation/condensation mass flux

r


, 1 11

#### *2.2.5. Lung deposited dose calculation for pharmacokinetic studies*

Different mechanisms can induce aerosol deposition in pulmonary route, i.e., diffusion, sedimentation, inertial impaction, interception, etc. Local deposition of particles in human airways can be quantified in terms of the deposition fraction (DF) or deposition efficiency (DE) in a specific region (e.g., oral airway, first, second, and third bifurcations). The DF is calculated as:

$$DF\_i = \frac{\text{mass of aerosols deposited at region i}}{\text{total mass of aerosols in haled}} \tag{18}$$

The DE can be calculated as:

$$DE\_i = \frac{\text{mass of aerosols deposited at region } i}{\text{mass of aerosols entering region } i} \tag{19}$$

In addition to the above two traditional deposition parameters DE and DF, a deposition enhancement factor (DEF) can quantify local particle deposition patterns, which is defined as:

$$DEF = \left( DF\_i / A\_i \right) / \left( \sum\_{i=1}^{n} DF\_i / \sum\_{i=1}^{n} A\_i \right) \tag{20}$$

where *Ai* is the area of a local wall cell *i*, *n* is the number of wall cells in one specific airway region, and *DFi* is the local deposition fraction in wall cell *i*. Clearly, high DEF-values suggest the presence of "hot spots" of local particle depositions as compared to the regional average. See Refs. [3, 10] on how to calculate depositions for nonspherical particles, droplets, and vapors.
