Composition of different EC droplets (mass fraction) used for modeling.

## Abstract

The awareness is growing of health hazards and pharmaceutical benefits of micro-/nano-aerosol particles which are mostly nonspherical and hygroscopic, and categorized as “unconventional” vs. solid spheres. Accurate and realistic numerical models will significantly contribute to answering public health questions. In this chapter, fundamentals and future trends of computational fluid-particle dynamics (CFPD) models for lung aerosol dynamics are discussed, emphasizing the underlying physics to simulate unconventional inhaled aerosols such as fibers, droplets, and vapors. Standard simulation procedures are presented, including reconstruction of the human respiratory system, CFPD model formulation, finite-volume mesh generation, etc. Case studies for fiber and droplet transport and deposition in lung are also provided. Furthermore, challenges and future directions are discussed to develop next-generation models. The ultimate goal is to establish a roadmap to link different numerical models, and to build the framework of a new multiscale numerical model, which will extend exposure and lung deposition predictions to health endpoints, e.g., tissue and delivered doses, by calculating absorption and translocation into alveolar regions and systemic regions using discrete element method (DEM), lattice Boltzmann method (LBM), and/or physiologically based pharmacokinetic (PBPK) models. It will enable simulations of extremely complex airflow-vapor-particle-structure dynamics in the entire human respiratory system at detailed levels.

### Keywords

- computational fluid-particle dynamics (CFPD) model
- lung aerosol dynamics
- unconventional inhaled aerosol
- fiber
- droplet
- vapor
- whole lung model
- lattice Boltzmann method (LBM)
- compartmental modeling
- physiologically based pharmacokinetic (PBPK) modeling

## 1. Introduction

Inhalable aerosols consist of multiple nano-to-micro-scale solid or liquid particulate matters (PMs) with dissolved or embedded compounds, as well as associated vapors. Instead of being spheres with constant diameters, most of these airborne PMs are nonspherical in shape, with everchanging sizes due to evaporation/condensation. Therefore, they are defined as “unconventional aerosol particles”. Examples include fiber-like carbon nanotubes (CNTs) [1], hygroscopic cigarette particles [2–4], and pulmonary drug powders [5].

Both active and passive human exposures to ultrafine airborne aerosols are growing. Contributing factors include the increase in nanomaterial generation during product/device manufacturing, handling and use [6], consumptions of conventional and new tobacco products [2–4], as well as nanoparticle pulmonary drug development [7, 8]. Of great concern are the toxic health risks and therapeutic benefits induced by the inhalation. This is mainly because the small sizes facilitate ultrafine PM transport into deeper lung airways, translocation to the bloodstream crossing body membrane barriers, thus affecting organs and tissues in both pulmonary routes and systemic regions. For example, the transmission of respiratory infectious disease begins from the deposition of pathogens carried by airborne PM into the respiratory tract [9]. Since infection initiation is usually region specific, local dosages (i.e., lung deposited, tissue, and delivered doses) in pulmonary routes and systemic regions are essential for precise risk assessments. Another example is the worker exposure to airborne nanomaterials in the workplace [6]. Furthermore, pulmonary targeted drug delivery is favorable for both lung and systemic disease treatments, owing to the strong capability of the lung to absorb pharmaceuticals. Using ultrafine unconventional particles as drug carriers has also been promoted in medical device design, encapsulation, and drug delivery [7, 8, 10], because of their extraordinary capabilities to penetrate into deeper lung airways by virtue of the nonspherical shape or the hygroscopic growth characteristics. In this case, accurate local dosage information is also essential for the evaluation of drug delivery efficacy to targeted sites.

In light of the above discussion, the primary goal of relevant researches is to provide high-resolution data and enhance the fundamental understanding of the physical and chemical mechanisms for the dynamics of unconventional inhaled aerosols in human respiratory systems and systemic regions. Data of interests include: transport, deposition, dissolution, absorption, distribution, and clearance of unconventional aerosols in lung airways and systemic regions under various exposure conditions. These data will also help to answer the key question: “*What type of inhaled aerosols deposits where at what surface concentrations in subject-specific human respiratory systems under distinct breathing conditions, and what delivered dose into the site of interest in systemic regions*?”

Experiments or clinical tests are not able to provide such high-resolution data on human subjects, because of ethical reasons and lack of reliable measurements. Hence, developing and applying an accurate and realistic computer simulation model is desired. Different types of numerical models have been developed specifically to determine the fate of inhaled aerosols in human bodies, as has been reviewed recently by Phalen and Raabe [11]. To simulate particle trajectories, interactions, and lung deposition of unconventional inhaled aerosols, computational fluid-particle dynamics (CFPD) models based on Euler-Lagrange and Euler-Euler methods are the best options compared to other numerical models such as empirical and semiempirical ones [11]. Specifically, accurate prediction of the fate of such aerosols requires the model to be able to capture the effects of particle shape and size change on their trajectories before deposition. Based on more underlying physics and less extrapolations and simplifications, CFPD models are able to describe the fundamental behaviors of fluids and particles with high-resolution local data by solving conservation laws with constitutive equations [3, 4, 10]. Hence, these validated noninvasive, cost-effective, and accurate CFPD models will complement in vitro and in vivo studies on unconventional lung aerosol dynamics by simulating fluid-particle flows that are difficult to be reproduced experimentally. They will be able to deliver new perspectives of scientific knowledge and healthcare data for medical device improvements for effective and targeted drug deliveries, novel lung therapeutics, noninvasive disease diagnostic methodologies, and exposure health risk evaluations. Accordingly, CFPD models are expected to contribute to the following aspects:

Understand fundamental, physical, and chemical mechanisms that influence transport phenomena of different inhaled aerosol.

Virtualize corrective surgery leading to improved breathing patterns in abnormal airways.

Enhance evaluating capabilities of toxicologists on adverse health effects of inhaled particulate matters.

Design and engineer novel carriers for pulmonary drug delivery.

Establish guidelines for limiting workplace exposure to new nanoparticles and micron fibers, considering both short-term and long-term exposures.

This chapter is intended to serve as an overview of advanced CFPD models for unconventional inhaled aerosols and a demonstration of how to build and apply CFPD models on lung aerosol dynamics simulations. It also shows a roadmap to the next-generation multiscale lung dosimetry models. Specifically, fundamental mathematics of advanced CFPD modeling techniques is introduced in Section 2, focusing on predicting unconventional aerosol transport and deposition in human upper airways, i.e., nonspherical fiber, hygroscopic droplet, vapor mixture, etc. Basic steps of CFPD modeling are also discussed in this section. Section 3 presents applications of different advanced CFPD models with parametric sensitivity analyses on particle transport and deposition characteristics in airways. Challenges and future directions are discussed in Section 4 on developing next-generation multiscale numerical models, which will be able to simulate unconventional aerosol transport, deposition, and translocation from the complete pulmonary routes into systemic regions. Structured literature reviews on the above topics are also integrated in the following sections. Other reviews can be found in [12, 13].

## 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.

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

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.

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 2^{23} 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 m^{2} 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.

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 2^{23} 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.

### 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

in which *u*_{j} represents the fluid velocity and *ρ* is the air density.

Navier-Stokes (N-S) equation

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

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

Energy equation

where *k* is thermal conductivity and *c*_{p} is specific heat capacity. Furthermore, *S*_{T} is the thermal sink or source term due to radiation, chemical reaction, phase change, etc. *Φ* is the dissipation function which can be given as:

The equation of state as stated by the ideal gas law should complement the above equations for model closure, i.e., *pV*_{m} = *RT*, where *V*_{m} 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]:

Here, *m*_{p} are the velocity and mass of the particle, respectively;

in which *A*_{p} is the projected particle area and *C*_{D} is the particle drag coefficient [3, 4, 10]. In Eq. (6),

#### 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:

where *Y* is the vapor mass fraction, *S*_{Y} is the source term.

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

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]:

The translational equation (see Eq. (6)) was revised with force correlations specifically for ellipsoidal fibers;

Euler’s rotation equations were introduced with the Euler quaternions to calculate fiber rotational motion and orientation at each time step;

The hydrodynamic torques acting on fibers are employed with correlations pertinent to ellipsoids.

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:

Here, (*I*_{x ′}, *I*_{y ′}, *I*_{z ′}) 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 (*T*_{x ′}, *T*_{y ′}, *T*_{z ′}) are the torques acting on the particle with respect to the principal axes *x*′, *y*′, *z*′. For ellipsoidal particles, (*I*_{x ′}, *I*_{y ′}, *I*_{z ′}) can be written as:

where *a*_{p}, *b*_{p}, 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 droplet-vapor interactions with evaporation and condensation effects during their transport in human respiratory systems [2, 3].

#### 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:

where subscript *a−v* means “air-vapor” and *N* represents the number of volatile components in droplets. The energy source term

where *N*_{d,cell} is the total droplet number in a specified mesh cell, subscript *v−d* means “between vapor and droplet”, *L*_{s} 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:

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

in which *A*_{d} is the droplet surface area. Additionally, *dt*_{d} represents the droplet phase time differential and Δ*t*_{f} 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:

where the average evaporation/condensation mass flux

in which *Y*_{s,vsurf} and *Y*_{s,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, *Y*_{s,cell} is determined by Eq. (13), which realizes the droplet-vapor interaction modeling.

The energy conservation of droplets is given by:

where *Nu* is the Nusselt number.

#### 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:

The DE can be calculated as:

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:

where *A*_{i} is the area of a local wall cell *i*, *n* is the number of wall cells in one specific airway region, and *DF*_{i} 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.

## 3. Case studies

Employing the advanced CFPD models for unconventional aerosol particle dynamics discussed in Section 2.4, two case studies are briefly shown in this section: (a) the transport and deposition of ellipsoidal fiber as a new pulmonary drug carrier in a subject-specific human upper airway model and (b) the transport, hygroscopic behavior, and deposition of multicomponent electronic cigarette (EC) droplets in an idealized triple bifurcation unit (TBU).

### 3.1. Transport and deposition of ellipsoidal fibers in a human upper airway

Due to the toxicity of airborne ellipsoidal fibers to human and hence the lack of experimental data using human volunteers because of the severe ethical constraints, it is essential to numerically investigate the transport and deposition of fibers in human respiratory systems. On the other hand, multifunctional micron fibers are also being used as pulmonary drug carriers for disease treatment. Therefore, the numerical investigations will provide data for estimations of dosimetry, safety, and the efficacy of drugs in the lungs which are also critical factors in the development of inhaled medicines.

Considering the significant impact the anisotropic shape has on the accuracy of trajectory prediction, the one-way coupled EL-ER method discussed in Section 2.2.4.1 was employed [10] to calculate the translational and rotational motion of submicron ellipsoidal fibers of different aspect ratios (1 < *β* < 20). Brownian motion was also considered in both Poiseuille flow (see Figure 7(a)–(d)) and complex flows inside a subject-specific human upper airway (see Figure 8).

To gain basic insight for ellipsoidal particle dynamics, Poiseuille flow with an averaged Reynolds number of Re_{ave} = 173 is considered where ellipsoidal fibers with randomly initialized incidence angles were released at different inlet-plane positions, their trajectories computed and visualized (see Figure 7(a)–(d)). Figure 7(a) and (b) show the rotational/translational motion of the ellipsoidal fibers (*β* = 14 and *a*_{p} = 0.5 μm) released at two locations, i.e., off-center and center. Released at an off-center position (see Figure 7(a)), the coupling of rotation and translation of the ellipsoidal particle can be clearly observed with periodical sudden turns of 180° in a “clockwise” manner. Released at the center (see Figure 7(b)), the flow velocity gradients were zero, i.e., no torques were affecting the particle. As a result, the particle kept its initial orientation till it reached a location where the velocity gradients and resulting torques were sufficiently large to rotate the particle. The frequency of the sudden rotation is influenced by the stability of the fibers at different locations of the flow field. It must be noted that using conventional CFPD models for spherical particles, the rotational motion is neglected and not able to be visualized.

Parametric sensitivity analyses were also performed for the aspect ratio effect on the translational and rotational motion for ellipsoidal fibers with the same equivalent volume diameter (see Figure 7(c) and (d)). Specifically, Figure 7(c) shows different *y*-direction velocities of fibers with the same volume but different aspect ratios. With the increase in aspect ratio *β*, the *y*-direction sedimentation velocity decreases and the rotation frequency of the fibers decreases. Indeed, with different aspect ratios, the Stokes resistance correction factors in the three principal directions of the ellipsoidal particle are all changing, which leads to the variation of the sedimentation velocity in *y*-direction. Figure 7(d) provides different trajectories of particles with the same volume and different aspect ratios. Particles with larger *β*-values have a stronger tendency to follow the mainstream and hence travel farther than a particle with a smaller *β*. Therefore, it can be conjectured that fibers are more likely to migrate deeper into the lung airways when compared to spherical particles of the same volume. In other words, when quasi-aligned to the flow, a fiber experiences a larger drag force in the gravitational direction (i.e., *y*-direction), so that particles with larger aspect ratios may not deposit in parallel flow.

To further investigate the transport characteristics of ellipsoidal fibers with different aspect ratios and the same volume, parametric studies on fiber transport in a subject-specific human upper airway with 20-mm inlet are performed. Specifically, a comparison of total deposition efficiencies of spheres and ellipsoidal fibers for two inhalation flow rates is shown in Figure 8. The total deposition data indicate that slender toxic ellipsoidal fibers (i.e., those with higher aspect ratios) are potentially more harmful than thicker ones due to their ability to penetrate into deeper lung regions when somewhat aligned with the mainstream flow. In contrast, to treat deeper lung diseases using pulmonary drugs, slender ellipsoidal fibers as drug carriers may enhance the drug delivery efficacy compared to conventional spherical ones.

### 3.2. Size change dynamics of multicomponent droplets in a G3-G6 TBU

As tobacco-smoking is addictive with potentially detrimental effects and hence related to high health-care cost, alternatives such as electronic cigarettes (ECs) have been developed. EC aerosols consist of multicomponent droplets and vapors, e.g., water, nicotine, glycerol, and PG. Although ECs can look like and ultimately should “taste” like conventional cigarettes, their design and function as well as product quality are very different. Debates are still going on for the use of ECs, with reported positive and negative aspects of EC consumptions [4]. Therefore, investigations on the differences between conventional cigarette smoke and EC aerosols are of great interest for EC product design as well as government regulations.

In this study, as the first step, the focus is on hygroscopic growth of nanosize multicomponent droplets and droplet-vapor interactions during transport with subsequent deposition in a triple bifurcating unit (TBU) representing G3–G6 lung airway (see Figure 9). Using the MCM-DD model derived in Section 2.2.4.2, electronic cigarette (EC) droplets with different compositions were compared with conventional smoke particles (CSPs) as well as solid particles on their different behaviors of evaporation/condensation dynamics during their transport and deposition in the TBU, with body temperature *T* = 310.15 K and relative humidity RH = 99.5% [2].

Figure 10(a) shows the comparisons of hygroscopic behavior between three EC droplets with different compositions and two conventional cigarette smoke particles (CSPs) with different compositions (see Table 1 for compositions). At RH = 99.5%, both EC-droplets and CSPs grow due to condensation effect. However, CSPs #1 and #2 grow less than EC-droplets #B1 to #B3. This suggests that with initial droplet diameter and ambient environment (humidity) being the same, EC-droplets #B1 to #B3 will be always larger in size than the liquid CSPs #1 and #2. Further, these EC-droplets with larger diameters with have higher inertia, and thus being less affected by the local airflow structure than CSPs.

Water | Glycerol | PG | Nicotine | |
---|---|---|---|---|

EC Droplet #B1 | 0.25 | 0.398 | 0.332 | 0.02 |

EC Droplet #B2 | 0.15 | 0.498 | 0.332 | 0.02 |

EC Droplet #B3 | 0.05 | 0.598 | 0.332 | 0.02 |

Figure 10(b) shows the total deposition of EC-droplets #B2, CSPs and solid particles in the G3-G6 TBU, assuming the same initial particle diameters (i.e., 264 nm ≤ *d*_{d,ini} ≤ 3200 nm). It can be observed that with the same initial diameter, the total deposition fraction for the three particles is EC-droplets #B2 < CSPs #1 < solid particles. Clearly, the hygroscopic growth of droplets and CSPs limited their Brownian diffusion, thereby reducing TDE in G3–G6. Since this effect is more prominent for EC-droplets than for CSPs (see Figure 10), the latter therefore has higher total deposition efficiency values. In other words, CSPs #1 deposit more in the G3–G6 TBU than #B2 EC-droplets due to lower hygroscopic growth and stronger Brownian diffusion effect. It can thus be conjectured that EC-droplets may deposit more than CSPs #1 up to the first bifurcation, due to more pronounced inertial impaction. In contrast, those EC-droplets tend to deposit more after the first bifurcation in G3–G6 TBUs and deeper lung regions where Brownian diffusion dominates deposition. Results also indicate the important effect of hygroscopic behaviors of droplets on the accuracy of deposited-dose predictions.

## 4. Challenges and future directions

Challenges faced by the current CFPD models and numerical studies for unconventional lung aerosol dynamics are mostly induced by the point-mass and point-force assumptions, the resolution of medical imaging techniques, the computational hardware capabilities, as well as the gap between CFPD prediction capabilities and the needs of toxicologists and health-care providers. Specifically,

The point-mass and point-force assumptions raise the difficulty to apply CFPD modeling techniques on aerosol transport and deposition simulation in the alveolar region, where the particle size will not be much smaller than a single mesh element and the assumptions are not valid.

The current resolution of medical imaging techniques is not sufficiently high to identify small lung airways beyond G8. Therefore, most CFPD models are still for human upper airways, and it is impossible to reconstruct a complete subject-specific human respiratory system using CT/MRI imaging data. This restricts CFPD models from simulating aerosol transport in a full inspiratory-expiratory breathing cycle, i.e., the quantitative prediction of exhalation clearance of aerosols.

Due to the computational resources limitation, most CFPD models still use simplified mechanisms to reach good compromises between accuracy and computational efficiency, e.g., one-way coupling for dilute particle suspensions. Therefore, CFPD models are not suitable for dense particle suspensions in which the particle-particle interactions are not negligible.

Although CFPD models are able to provide high-resolution lung deposition patterns for unconventional aerosols, the deposited dose cannot be used directly for health risk assessment or drug delivery efficacy evaluation. Improvements need to be made for the current CFPD models to calculate the tissue and delivered doses.

This section shows the roadmap to bridge the current simulation techniques to emerging toxicological/pharmaceutical paradigms to overcome the limitations of current numerical simulations. The long-term goal is to build an accurate and cost-effective digital respiratory simulation framework of subject-specific exposure to inhalable aerosols under realistic breathing cycles. A proposed framework and kinetic process of this multiscale model, i.e., a hybrid CFPD-LBM-PBPK model (see Sections 4.1–4.3), is shown in Figure 11. Details of future researches for extensions and improvements of the current CFPD model and the development of the next-generation model are discussed below.

### 4.1. Reconstructing and modeling the entire conducting zone of human lung

In order to provide full-scale regional lung deposited doses in the oral/nasal cavity, conducting and respiratory zone, it is necessary to find methodologies to represent or reconstruct the complete human respiratory system geometry, including the smaller airways, which are not visible via CT/MRI imaging. Accordingly, two methods can be further developed: (1) coupling methods with single-path or multiple-path lung airway trees and (2) virtual lung model with idealized airways for higher generations covering the entire lung conducting zone (G0–G17). This is crucial to create a full-scale lung aerosol dynamics simulation featuring full breathing cycles.

#### 4.1.1. Coupling methods with single-path or multiple-path lung airway trees

Coupling methods combine coupled boundary conditions with truncated lung airway trees to represent the whole conducting zone. Realization of coupling methods can be performed in the following steps:

Divide the whole lung airway trees into segments [40] or truncate the airway tree with single path or multiple paths left [41–43];

Simulate each segment in series or in parallel;

The boundary conditions applied to each truncated airway ends will be obtained from the previous simulation of the next larger section (inhalation) or the next smaller section (exhalation).

Kleinstreuer and Zhang [40] used the TBU of Weibel type A and constructed the symmetric and in-plane tracheobronchial airways from G0 to G15, and used coupling method to simulate airflow and particle deposition in the lung airways with steady-state inhalation flow rate *Q*_{in} = 30 L/min. Walters and Luke [42] created an 8-generation airway tree with random connection angles following Weibel’s lung morphology. They applied coupled pressure boundary conditions for inhalation cases, and validated the coupling method for inhalation cases. They extended the application of this coupling method on particle transport in lung airways [43]. Tian et al. [41] developed the stochastic individual path (SIP) and multiple stochastic individual path (MSIP) approaches to represent the whole-lung geometry. Compared to the coupling method applied on detached airway geometries [41], this model is able to simulate transient breathing scenarios. Recently, similar models have also been developed [44, 45].

Since modeling the whole lung is computationally expensive, using the coupling method is potentially advantageous for fast-solution purposes. However, accuracy of the coupling methods still needs to be tested compared with using virtual lung model containing the entire conducting zone without any truncations.

#### 4.1.2. Virtual lung model for the entire conducting zone

To encompass the entire conducting zone, the main objective of virtual lung model reconstruction is to develop a stochastic algorithm that generates 3D asymmetric human lower lung airways, which can then be connected to the subject-specific upper lung airways [46–48]. To establish such a virtual lung model, the upper airway (i.e., mouth and nasal airways, throat, and trachea) will be exactly reconstructed using image processing software that converts CT/MRI images. Conducting airways, from the trachea down to the level of the terminal bronchioles, will be generated via a stochastic algorithm [46–48]. Specifically, by generating geometric variables (i.e., airway diameters, branching angles, and rotation angles of successive branching planes) stochastically with morphometric restrictions [46], following fractal model (i.e., Hess-Murray law), or grow from 1-D centerline airway tree [48], a complete tracheobronchial airway tree starting from the trachea can be constructed. The accuracy of the geometry can be judged in two ways: (1) the airways follow the centerlines of a 1D lung airway tree and (2) it conforms to surface data of the airways segmented from medical imaging. The combined geometry can be exported as a stereolithography (STL) file that can be used for mesh generation as needed for computer modeling and simulations.

Compared with coupling methods, constructing a virtual lung model will enable CFPD simulations with more realistic boundary conditions. It will be computationally achievable to simulate the multiphase inhaled aerosols in the entire lung conducting zone, with foreseeable increase in computational power.

### 4.2. Lattice Boltzmann method (LBM) for aerosol transport in alveoli

A portion of the inhaled aerosol will enter the respiratory zone (G17–G23) of human lung, which is mainly the alveolar region. In this case, the difficulty in using CFPD model is that the point-mass and point-force assumptions do not hold and particle-particle interactions are not negligible. Lattice Boltzmann method (LBM) can be employed for the two-way fluid-particle interaction and many-particle interaction in alveoli, due to its advantage in resolving particles that are comparative or larger than the mesh cells in size.

#### 4.2.1. Background and literature review

Lattice Boltzmann method (LBM) solves the Lattice Boltzmann equation (LBE) to obtain the probability density distribution, which is then used to obtain the pressure and velocity field through integration. Details of the LBE, discrete speeds, discretization, and boundary conditions are given in [49] and avoided here for brevity. LBM consists of two main steps: (1) propagation: particle distribution moves along the lattice bonds to neighboring lattice nodes and (2) collision: particles on the same node interact and adjust velocities to conserve momentum and mass. Communication between nodes in the lattice is required only in the propagation step requiring message passing or update of the shared memory location.

LBM has been applied to the airflow and particle transport simulations in the human airways due to the advantages mentioned above compared with CFPD models. Many studies focused on the flow structures in nasal cavities. For example, Finck et al. [50] first employed LBM for steady simulations of laminar nasal flow. The model predictions agree well with those of finite-volume based model using structured grids. However, the LBM has high flexibility in fast grid generation and easy boundary-condition implementation for such complex geometries compared to conventional Navier-Stokes solvers. Hörschler et al. [51] compared steady and unsteady flows in real nasal cavity, and concluded that the transient effects are only important at Re < 1500. Henn et al. [52] coupled Lagrangian particle tracking method to LBM to simulate aerosol dynamics in a human nasal cavity. They concluded that transient simulations are necessary to study particle dynamics in nasal cavities. Other studies modeled airflow in the nasal cavities with diagnosed pathologies [53].

There are also existing researches focusing on air flow simulation in human upper airway using LBM. Ball et al. [54] showed that LBM yields higher fidelity than the RANS models for laminar-to-turbulent transition flow regime predictions inside an idealized human extrathoracic airway, as compared to experimental results. Schroeder et al. [55] presented LB simulation results of pulsatile flow at two breathing frequencies in a human respiratory system down to G6. They showed that inspiratory flow is characterized by flow separation and secondary flow, while expiratory flow is characterized by more homogeneous structure but higher vorticity. This is in contrast with the study by Eite et al. [56], which showed that less vorticity is generated at exhalation than at inhalation in a nasal cavity. Krause [53] investigated respirations in the upper part of the human lungs using LBM coupled with a model for the rest of lung, thereby constructing a two-scale whole lung model. This latter model assumed that the lower part of the bronchial tree can be described by a symmetrical dyadic tree of pipes, where Poiseuille type flow holds. Imai et al. [57] developed a LBM with GPU acceleration using an adaptive meshing method. The model was used to study the effect of breath holding on the deposition of microparticles in pulmonary airways during inhalation. The subject-specific airway model contains thirteen generations of airways. The study found that breath holding can effectively increase particle deposition; yet the effect varies for different particle sizes. Wang and Elghobashi [58] analyzed flow structures in human upper airways including the nasal cavity, the pharynx, the larynx, and the upper part of the trachea for both healthy and diseased subjects. The laminar-transitional-turbulent flow in the airway was captured in inspiration, and laminar flow was found during expiration. The authors proposed that the maximum positive value of the local time-averaged second derivative of pressure can be used to locate the obstructions due to disease in the upper airway. It should be noted that most studies claimed that LBM provides better accuracy for air flow simulation in respiratory systems. However, Chen and Gutmark [59] did a comparison study of airflow in an idealized human extrathoracic airway. They claimed that LES is more accurate in predicting root-mean-square flow field compared with Reynolds stress model and LBM.

In contrast to many existing LBM simulations on airflow patterns in human respiratory systems, there are only a few studies using LBM for lung aerosol dynamics. Henn et al. [52] used LBM coupled with Lagrangian tracking method to study microparticle transport in human nasal cavities. Trunk et al. [60] developed an Euler-Euler LBM method to investigate the particle dynamics in an idealized bifurcation geometry representing the split from trachea to main bronchi. Particle capturing BC and numerical stabilization were discussed for lung aerosol dynamics simulations. Li and Kleinstreuer [61] analyzed transient airflow in two-dimensional alveoli and bifurcating alveolar ducts. The results highlighted the importance of alveolar geometries and shapes, as well as the effect of the moving alveolus walls. Implications of these effects on particle deposition in the alveolar sac were discussed. However, simulation of particle suspension has not been performed.

#### 4.2.2. LBM vs. CFPD models

As compared with the CFPD models, LBM has advantages in the following aspects [54]: the convection operator is linear in LBM; the conservation laws are exactly preserved in LBM; the calculations are arithmetic operations in LBM; and the boundary conditions (BCs) like no-slip or moving boundaries are easy to implement. In addition, the entropic lattice Boltzmann method is attractive in modeling transitional and fully turbulent flows for high-Re simulations [62]. It should be noted that when the particle size is much larger than the lattice size, LBM can readily account for the two-way fluid-particle interaction and many body interaction by accurately calculating the hydrodynamic force and torque exerted on the particle by the fluids [63–65]. Also, Brownian motion of particles arises spontaneously from stochastic fluctuations in the fluid stress tensor, avoiding application of a random force on particles. Therefore, LBM is suitable for the simulation of aerosol dynamics, especially in the alveoli sacs. Furthermore, LBM can also be used to simulate the mucociliary clearance process when coupled with immersed boundary method or finite difference method [66, 67].

In conclusion, LBM is a rapidly developing method in lung aerosol dynamics simulation. Many researches have been done concerning airflow dynamics. On the other hand, LBM simulation including particle transport and deposition in alveoli is still largely unexplored. Much more efforts are demanded to couple particle transport with airflow simulation using LBM, with both sub-grid-size particles, and macro-size particles.

### 4.3. Physiologically based pharmacokinetic (PBPK) models for aerosol translocation into systemic regions

#### 4.3.1. Background and literature review

Carried by ultrafine aerosols, toxicants and therapeutic components can easily translocate from lung deposition sites into the systemic circulation, and accumulate in extrapulmonary organs. Therefore, after obtaining deposited dose of environmental pollutants or pulmonary drugs in the human respiratory systems, the judgment of potential adverse effects or curing effectiveness depends on agents’ toxicokinetic/pharmacokinetic information with more insights into their in-vivo behavior, and then toxicology, i.e., absorption, distribution, metabolism, and elimination. In order to provide comprehensive regional data for risk assessment, it is necessary to link a whole-body PBPK model to CFPD models, which extends exposure and deposition modeling to health endpoints. The integration of preclinical data generated by CFPD-PBPK models is also a key to optimize the drug delivery processes for clinical practice.

Physiologically based pharmacokinetic (PBPK) models are based on differential equations that are designed to simulate the absorption, distribution, metabolism, and excretion (ADME) in compartments which are connected biological regions, usually in whole animals or humans [11, 66]. Realistic descriptions of physiology and relevant pathways of metabolism can be incorporated into PBPK models. The history and details of PBPK models have been well discussed in [67, 68].

Focusing on the development of hybrid CFPD-PBPK models, the respiratory system serves as one compartment, but has significant influence on the whole body as it is the main barrier between inhaled aerosols and systemic circulation (see Figure 11). As direct inputs to the PBPK model, regional airway-tissue doses must be accurately calculated from the lung deposited doses acquired from the CFPD simulations. Therefore, the aim of incorporating CFPD models with PBPK models is to develop a new boundary condition for site-specific airway-tissue dosimetry for PBPK modeling that includes essential mechanisms, such as clearance, tissue reaction, tissue metabolism, diffusion, and blood perfusion. Assuming a linear first-order constitutive relationship, the material conservation law for compartment *i* of the airway tissue can be expressed as:

where *c*_{i} is the amount of pathogen per unit volume accumulating in compartment *i*, *k*_{ij} is the transfer rate from compartment *j* to *i*, and *S*_{ij} is a general source term representing exterior inputs, reactions, clearance, etc. The transfer rate, *k*_{ij}, can be determined from existing experimental data, or in vivo deposition/retention and mass transfer data obtained from laboratory animals. Other biological material properties and site-specific boundary conditions for the PBPK model will be obtained from the open literature.

Efforts have been made to link CFPD and PBPK/PBTK models for simulations of inhaled drugs or toxicants, from deposition to systemic regions of interests [69–75]. For example, Rigg et al. [69] established a numerical model for particle dissolution, absorption, and clearance (DAC) in the nasal cavity. Their DAC model can estimate the tissue dose for particles after deposition, which are the more accurate inputs than deposited dose for PBPK/PBTK calculation. Using subject-specific respiratory system geometries of different species, Corley et al. [70] compared the parameter difference between humans and animals and adapted the nasal CFD-PBPK model [71] to perform interspecies variability analyses and sensitivity analyses for vapor exposure. Anatomic, physiological, and biochemical parameter values were also listed [70]. Parameters such as airway geometry, inhalation flow rates, vapor concentration, partition coefficient between air and tissue, tissue thickness, and metabolism rate were identified as the significant factors that influence local vapor uptakes. Recently, they extended the CFD-PBPK model to transient simulation for toxic vapors [72]. Furthermore, PBPK models with different compartmental divisions have been established to estimate pulmonary drug delivery [73, 74]. Parameter values were also provided in the context.

#### 4.3.2. Comparisons among CFPD-PBPK, CFPD, and PBPK models

Compared with PBPK models, CFPD-PBPK models have the advantage in considering the subject variabilities induced by interindividual morphologies of human respiratory systems, which can be reflected by the different local deposition patterns provided by CFPD calculations. Additionally, CFPD-PBPK models extend exposure and deposition modeling of CFPD to health endpoints, i.e., delivered dosage estimations at site of interests in the whole human body.

Although research efforts have been made to develop hybrid CFPD-PBPK models as discussed in Section 4.3.1, none of the existing models display the realism needed for accurate dosimetry calculations, because they ignored the influence of essential physical and chemical mechanisms such as hygroscopic characteristics on the dispersion and deposition of droplets in the respiratory system, and used incomplete human lung airway geometries.

Therefore, future directions to develop more accurate and computationally efficient CFPD-PBPK models are to make more realistic compartment divisions, mechanism selections, and parameter value determinations to reflect the underlying aerosol dynamics in systemic regions [76]. For example, the compartmental division of cellular-level PBPK model should be based on specific anatomy and mass transfer characteristics of different lung regions. It is also important to collaborate closely with experimentalists and clinical doctors to find reliable data sets reflecting responses and pharmacokinetics for evaluating the performance of the CFPD-PBPK model [77]. Moreover, the determination of human physiological and biochemical parameters requires strong collaborations. These parameters include breathing rate, skin surface area, cardiac output, tissue blood flow rates, tissue volumes, partition coefficients, and biochemical rate constants [68]. Subject-variability must be considered too, i.e., physiological and biochemical parameter differences due to the differences in interindividual response to the biological effects of inhaled aerosols. It will not only contribute to the evaluation of the confidence in PBPK models on the lack of precise knowledge of the parameter values, but also contribute to the uncertainty analyses for the extrapolations between animals and humans. The error source must be identified and reduced to enhance the accuracy of the numerical results [75].

### 4.4. Others

Other numerical methods that can be integrated or employed to improve or replace current CFPD models are direct numerical simulations (DNSs) for turbulence [78], complete numerical simulations (CNSs) for multiphase flow [79], discrete element methods (DEMs) [80], etc. They will be able to extend the CFPD model capabilities by encompassing more underlying physics (e.g., particle-particle interactions) and less simplifications (e.g., point-mass and point-force assumptions).

#### 4.4.1. Complete numerical simulations (CNSs) for multiphase flow

Complete numerical simulation (CNS) methods are numerical techniques where the Navier-Stokes equations are applied to finite-size particles instead of introducing point-mass and point-force assumptions. Particle motion is governed by the Newton’s second law, based on the hydrodynamic forces calculated from the fluid equations. For CNS-based multiphase flow, all surface and collision forces due to fluid-particle and particle-particle interactions should be accurately calculated to obtain the velocity, pressure, and stress fields surrounding each particle [35]. Specifically, mesh scales for CNS methods are much lower than particle sizes; hence, surface forces such as drag force and lift force, as well as interaction force between particles can be directly integrated from the shear stress and pressure distribution along the particle surfaces. This avoids empirical correlations, such as drag and lift coefficients.

Compared with CFPD models which are based on Euler-Euler or Euler-Lagrange methods, CNS is a mesoscale model that can accurately describe motions of particles with arbitrary shapes. As such, CNS of multiphase flows can promote fundamental understanding of the underlying physics, e.g., the nonlinear and geometrically complicated phenomena of particle-particle and particle-wall interactions.

With the enhancement of the computational power in the future, CFPD models will be replaced by DNS and CNS methods to simulate unconventional lung aerosol dynamics. However, the computational cost of DNS and CNS methods is still too high for engineering application. Applying these methods to a large number of particles is not realistic according to the current computational resources [79].

#### 4.4.2. Discrete element method (DEM)

Discrete element method (DEM) can provide transient forces and torques acting on individual particles. DEM was first proposed by Cundall and Strack [81] based on molecular dynamics. The most attractive feature of DEM method is its relatively efficient algorithms for the contact detection and contact force calculation between arbitrary shaped particles [82]. Soft-sphere models are the most common type of models used in DEM, which simplify particle-particle interaction force with spring-dashpot-slider systems [27]. Determinations of the stiffness coefficients and damping coefficients are essential in DEM interaction models. Contact forces are calculated in two steps in DEM methods [83]:

Contact detection: identify which particles are contacting;

Contact resolution: the actual contact forces are calculated. This step involves the accurate calculation of the contact geometry and kinematics, typically characterized by the overlap depth/separation and sometimes also relative tangential motion.

Additionally, for nonspherical particles, representation of the particle shape is also essential in DEM simulations. Two major treatments of the particle shape representation are widely used for the contact detection and contact force calculation between nonspherical particles [84]: superquadric method and multisphere method. The algorithm to integrate DEM with Euler-Euler and Euler-Lagrange models (see Figure 12 for a standard flow chart) for dense particulate flow is straightforward [38]. However, there are two key challenges:

Spring-dashpot-slider systems are not designed for micro-/nano-particles. Coefficients, i.e., equivalent stiffness and damping coefficient, need to be determined via experiments specifically for ultrafine particles. Or else, the validity of using soft-sphere models needs to be tested.

DEM is not computationally efficient on handling aerosol transport and deposition in human respiratory systems, in which the difference in characteristic dimensions between particles and flow domains is too large. Therefore, efficient algorithms need to be developed.

## 5. Summary

In this chapter, CFPD models have been systematically overviewed, with emphases on simulating the transport and deposition of unconventional ultrafine aerosol particles in human upper airways. Compared to other numerical models for inhaled aerosol dose [1], CFPD models provide more detailed characterizations of airflow and aerosol deposition patterns with a reasonable degree of accuracy and computational cost.

However, in order to provide comprehensive tissue and delivered dosage data for risk assessment or drug delivery efficacy evaluation, it is necessary to improve current CFPD models and extend lung deposition calculation by integrating with other simulation techniques to build a multiscale dosimetry model (see Figure 11 for the framework). This multidisciplinary effort requires expertise in medical imaging, airway geometric reconstruction, computational techniques, pulmonary physiology and medicine, fluid mechanics, etc. From a physics viewpoint, the future numerical models for lung aerosol dynamics will provide comprehensive analysis of the fate of inhaled particulate matter under realistic conditions, namely, actual particle characteristics, inhalation flow rates, and human respiratory systems. The anticipated results will greatly increase understanding of distinct particle transport/deposition phenomena under realistic inhalation conditions as well as quantify the knowledge base for dosimetry and health-effect studies in the case of toxic particles, and optimal drug-aerosol targeting in the case of therapeutic particles. Relevant numerical researches will lead to understanding and quantification of dosimetry effects of inhalation hazards, more accurate evaluation of toxic particle fate in the human lung, and improved estimates of drug delivery points for treating lung and other diseases.