CFD Simulations AC7-01

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

Aerosol deposition in the human upper airways

Application Challenge AC7-01   © copyright ERCOFTAC 2019

CFD Simulations

Overview of CFD Simulations

LES and RANS simulations were carried out in the benchmark geometry. The details of these numerical tests and their predicted deposition are given in the following paragraphs. In summary, the main differences in LES and RANS simulations are:

1. Computational meshes
2. Turbulence modeling
3. Different outlet boundary conditions
4. In RANS simulations a turbulent dispersion model was used to account for the effect of turbulence on particle transport.

Large Eddy Simulations

Computational domain and meshes

The geometry used in the calculations is the same as the one used in the experiments developed by the group in Brno University of Technology (BUT). The computational domain, shown in Figure 10, has one inlet and ten different outlets, for which appropriate boundary conditions must be specified in the simulations.

Figure 10: Computational domain viewed from different angles.

The digital model of the physical geometry was used to generate a proper computational mesh in order to perform the simulations. For the LES simulations, three meshes were generated to allow us to examine the sensitivity of the results on the mesh resolution. The coarser mesh includes 10 million computational cells, the intermediate one 30 million cells and the finer approximately 50 million cells. In these meshes, the near-wall region was resolved with prismatic elements, while the core of the domain was meshed with tetrahedral elements. Cross-sectional Views of these meshes at seven stations are shown in figure 11. A grid convergence analysis was carried out in order to determine the appropriate resolution for the LES simulations. This analysis is presented in section 3.2.4.

Table 5 reports grid characteristics, such as the height of the wall-adjacent cells ${\Delta r_{min}}$, the number of prism layers near the walls, the average expansion ratio of the prism layers (${\lambda }$), the total number of computational cells, the average cell volume (${V}$) and the average and maximum ${y+}$ values. The higher ${y+}$ values (above 1) are found near the glottis constriction and the bifurcation carinas, which are characterised by high wall shear stresses.

Figure 11: Cross—sectional views of the three generated meshes at seven stations (A—G).

Solution strategy and boundary conditions – Airflow

Large Eddy Simulations (LES) are performed using the dynamic version of the Smagorinsky-Lilly subgrid scale model (Lilly, 1992) in order to examine the unsteady flow in the realistic airway geometries. Previous studies have shown that this model performs well in transitional flows in the human airways (Radhakrishnan & Kassinos, 2009; Koullapis et al., 2016). The airflow is described by the filtered set of incompressible Navier-Stokes equations,

${{\frac {\partial {\bar {u}}_{j}}{\partial x_{j}}}=0}$	${(10)}$
${{\frac {\partial {\bar {u}}_{i}}{\partial t}}+{\bar {u}}_{j}{\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}=-{\frac {1}{\rho }}{\frac {\partial {\bar {p}}}{\partial x_{i}}}+{\frac {\partial }{\partial x_{j}}}\left[(\nu +\nu _{sgs}){\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}\right],}$	${(11)}$

where ${{\bar {u}}_{i},{\bar {p}},\rho =1.2kg/m^{3},\nu =1.7\times 10^{-5}m^{2}/s}$ and ${\nu _{sgs}}$ are the velocity component in the i-direction, the pressure, the density and kinematic viscosity of air, and the subgrid-scale (SGS) turbulent eddy viscosity, respectively. The overbar denotes resolved quantities.

The governing equations are discretized using a finite volume method and solved using OpenFOAM, an open-source CFD code. In this framework, unstructured boundary fitted meshes are used with a collocated cell-centred variable arrangement. The finite volume method in OpenFOAM is in general 2nd-order accurate in space, depending on the convection differencing scheme (CDS) used. Whenever possible (usually in the cases with lower Reynolds numbers), the 2nd-order linear CDS is used. The order of accuracy had to be decreased in some cases in order to stabilize the simulation. In these cases, the clippedLinear scheme was used, which provides a good compromise between the accuracy of the (2nd order) linear scheme and the stability of the (1st order) mid-point scheme. The temporal derivative is discretized using backward differencing, which is also second order accurate in time and implicit. To ensure numerical stability, the time steps used are ${5\times 10^{-6}s}$ for the cases of 15, 30 L/min and ${2.5\times 10^{-6}s}$ for the case of 60 L/min. The resulting maximum CFL numbers were 1.37 2.6 and 1.14 (values calculated on the nodes) for the simulations at 15, 30 and 60 L/min, respectively. Although CFL values greater than 1 were recorded, these are limited to very small regions near the distal segments and we expect that the accuracy of the simulation results will not be affected. In particular, the percentage of the number of nodes with CFL greater than 1 over the total number of nodes was ${1.2\times 10^{-5},0.175}$ and ${4.6\times 10^{-4}\%}$ for the computations at 15, 30 and 60 L/min, respectively. Moreover, using smaller time step values than the ones used would require a much larger number of time iterations to reach the desired physical time of the simulations (that is the time when the majority of injected particles in the domain have either deposited or exited). This would result in significant increases of the computational cost.

The non-linearity in the momentum equation is lagged in OpenFOAM (linearization of equation before discretisation). The system of partial differential equations is treated in a segregated way, with each equation being solved separately with explicit coupling between the results. In turbulent flow applications, where the time step is kept small enough to capture the smaller turbulent time scales, the pressure-implicit split-operator algorithm, or PISO—algorithm, is used.

For the lower flowrates of 15 and 30 L/min, the Reynolds number at the inlet of the model is in the laminar regime and thus a parabolic velocity profile is imposed. For the higher flowrate of 60 L/min, the Reynolds number at the inlet, based on the inflow bulk velocity and inlet tube diameter, is ${Re={\frac {U_{in}D_{in}}{\nu }}=3745}$ which lies in the turbulent regime. In order to generate turbulent inflow conditions, a mapped inlet, or recycling, boundary condition was used (Tabor et al., 2004). To apply this boundary condition, the pipe at the inlet was extended by a length equal to ten times its diameter. The pipe section was initially fed with an instantaneous turbulent velocity field generated in a precursor pipe flow LES. During the simulation of the airway geometry, the velocity field from the mid-plane of the pipe domain was mapped to the inlet boundary. Scaling of the velocities was applied to enforce the specified bulk flow rate. In this manner, turbulent flow is sustained in the extended pipe section, and a turbulent velocity profile enters the mouth inlet. The volumetric flowrates at the 10 terminal outlets are prescribed based on the values measured in vitro (Table 3). These outlet conditions result in high asymmetry in the ventilation of the two lungs: the left lung receives 29% of the inhaled air whereas the right lung receives 71%. A no-slip velocity condition is imposed on the airway walls and atmospheric pressure is set at the inlet boundary.

Solution strategy – Particle transport and deposition

Particle equation of motion

A Lagrangian point-particle approach is adopted for the simulation of transport and de— position of particles in the human airways. The motion of each particle is individually computed by solving Newton's equation to determine the particle velocity, ${{\vec {u}}_{p}}$:

${m_{p}{\frac {d{\vec {p}}}{dt}}={\vec {F}}_{D}+{\vec {F}}_{G}+{\vec {F}}_{B}}$

${(12)}$

where ${m_{p}}$ is the particle mass, and ${{\vec {F}}_{D},{\vec {F}}_{G}=m_{p}{\vec {g}}}$ and ${{\vec {F}}_{B}}$ are the drag, gravity and Brownian forces, respectively. Gravity acts in the vertical direction (parallel to the trachea). The position of the particle can be obtained from the kinematic equation:

${{\frac {d{\vec {x}}_{p}}{dt}}={\vec {u}}_{p}}$

${(13)}$

The drag force acting on the spherical particles is given by,

${{\vec {F}}_{D}={\frac {m_{p}}{\tau _{p}/\alpha }}({\vec {u}}-{\vec {u}}_{p})}$

${(14)}$

where ${\vec {u}}$ is the fluid velocity interpolated at the position of the particle, ${\alpha =C_{D}{\frac {Re_{p}}{24}}}$ and ${\tau _{p}}$is the particle response time, defined as:

${\tau _{p}={\frac {\rho _{p}d_{p}^{2}C_{c}}{18\mu _{f}}}}$

${(15)}$

with ${d_{p}}$ being the particle diameter, ${\mu _{f}}$ the dynamic fluid viscosity and ${Re_{p}=d_{p}\left|{\vec {\bar {u}}}-{\vec {u}}_{p}\right|/\nu _{f}}$ the particle Reynolds number. ${C_{c}}$ is the Cunningham correction factor, defined as:

${C_{c}={\frac {2\lambda }{d_{p}}}\left[1.257+0.4\exp(-0.55d_{p}/\lambda )\right]}$

${(16)}$

where ${\lambda =0.070}$m is the mean free path of air. This factor accounts for slip at the particle surface due to non—continuum effects, which appear when the size of the particle becomes comparable to the mean free path of the molecules of the surrounding gas. In our applications this is practically important when ${dp<1\mu }$m. The drag coefficient, ${C_{D}}$, is based on the correlation proposed by Schiller & Naumann (1935):

${C_{d}=}$	${{\frac {24}{Re_{p}}}\left(1+0.15{Re}_{p}^{0.687}\right)\quad {\text{if}}\ Re_{p}\leq 1000}$	${(17)}$
	${0.44\quad {\text{if}}\ Re_{p}>1000}$

The Brownian force is important for submicron particles and causes diffusion due to collisions with the air molecules (Finlay, 2001). The expression for the amplitude of its ith component is based on the correlation proposed by Li & Ahmadi (1992),

${F_{Brownian,i}=\zeta _{i}{\sqrt {{\frac {1}{\tilde {D}}}{\frac {2k_{B}^{2}T^{2}}{\Delta t_{p}}}}}}$,

${(18)}$

where ${\zeta _{i}}$ is a zero mean variant from a Gaussian probability density function, T=310K is the absolute temperature in the lungs, ${{\tilde {D}}=(k_{B}TC_{c})/(3\pi \mu _{f}d_{p})}$ is the Brownian diffusion coefficient, ${k_{B}=1.3806488\times 10^{-23}m^{2}\ kg/s^{2}}$ K is the Boltzmann constant and ${\Delta t}$ is the time step used for integration of the particle equations.

Deposition is assumed once a particle comes into contact with the airway walls. Reflection and re-suspension were not included, since the in vitro experiments used liquid particles which deposit when they hit the surface of the cast. In vivo, the existence of a mucus layer on the inner walls of the airways ensures that particles colliding with the surface deposit. The particle diameters ranged from ${d_{p}=0.5\mu m}$ to ${10\mu m}$, and the particle density was set to ${\rho _{p}=914kg/m^{3}}$, which corresponds to di-ethylhexyl sebacate (DEHS) particles in air at room temperature. At every time step, 10 particles for each size are released from random positions at the mouth inlet. Particles are released over a time period equal to a flow-through in the trachea, and the total number of injected particles is 100,000 for each particle size. The initial velocity of the particles is set to match the air velocity at the inlet. One-way coupling is considered, assuming dilute particle suspensions.

Particle tracking algorithm

Particle tracking is performed using the computationally eflicient and robust tracking algorithm proposed and implemented in OpenFOAM by Macpherson et al. (2009), with some additional improvements that were introduced in the latest versions of OpenFOAM (e.g. tracking of the particles using implicit decomposition of each cell into tetrahedra). A brief description of the numerical procedure follows.

Once the flow field has been advanced for a discrete time step ${|Deltat_{f}}$, the advancement of all particles follows for the same time window. Firstly, the particles new position ${{\vec {x}}_{p}^{n+1}}$ is computed explicitly using the particle velocity calculated at the end of the previous time step (${{\vec {u}}_{p}^{n}}$):

${{\vec {x}}_{p}^{n+1}={\vec {x}}_{p}^{n}+{\vec {u}}_{p}^{n}\Delta t_{p}}$

${(19)}$

The time step ${\Delta t_{p}}$ used in the above equation can differ from the flow time step since the particles are tracked from cell to cell by calculating and identifying face crossings. The advantage of the face crossing approach is the much more eflicient tracking of the particles in complex geometries that comprises of unstructured, arbitrary polyhedral cells, compared to methods that redetermine the hosting grid cell in every iteration. In this manner, for a flow time step ${\Delta t_{f}}$, the motion of a particle can be performed as a series of individual tracking events, each ending when the particle either crosses a face of a cell or arrives at the final destination. Thus, the maximum time step used to track a particle is the one defined for the continuous phase simulation ${(\Delta t_{f})}$. At the particle new position, which is either on a face that has been crossed or the final destination of the particle, eq. 12 is integrated using the implicit Euler scheme to compute the new particle velocity:

${{\vec {u}}_{p}^{n+1}={\frac {{\frac {\Delta t_{p}}{\tau _{p}/\alpha }}{\vec {u}}+{\vec {u}}_{p}^{n}+{\frac {\Delta t_{p}}{m_{p}}}\left({\vec {F}}_{Buoyancy}+{\vec {F}}_{Brownian}\right)}{1+{\frac {\Delta t_{p}}{\tau _{p}/\alpha }}}}}$

${(20)}$

The non—linear terms in the above equation (such as ${\alpha (Re_{p}({\vec {u}}_{p}))}$) are linearised using the particles known properties from the previous time step. Spatial interpolation is used to approximate the fluid velocity ${\vec {u}}$ and the minimum distance vector ${{\vec {r}}_{p,w}}$ at the position of a particle, which are needed in the calculation of drag forces, respectively. Both these vectors are mesh variables, i.e. their values are stored at the cell centres: ${\vec {u}}$ is obtained at each time step from the solution of flow equations while ${{\vec {r}}_{p,w}}$ is computed once in the beginning of the simulation for each cell and fitted to the mesh in order to reduce the computational cost. The spatial interpolation method used decomposes the grid cell that hosts the particle to tetrahedra, using the cell centre-point, face centre-point and two vertices. In this manner, a tetrahedral cell is decomposed into 12 tetrahedra and a prismatic cell into 18 tetrahedra. Then, it searches for the tetrahedron enclosing the position of the particle and when it finds it, linear interpolation is performed based on the 4 vertices of the tetrahedron (Baker, 2003). This strategy leads to an interpolated fluid velocity that is continuous inside each grid cell as well as across the grid cells.

The above procedure, i.e. determination of new position and velocity of the particle is repeated until the particle reaches its final destination at the end of the flow time step.

Numerical accuracy

The sensitivity of the calculations (airflow fields and particle deposition) on the mesh resolution was tested at the lower flowrate of 15 L/min and the results are presented in this section. Figure 12 displays contours of the mean velocity magnitude at cross sections in the mouth-throat/trachea, the main bifurcation and the left/right bronchi for the three different meshes examined. Similar mean velocity contours are observed on the three meshes, suggesting that at the lower flowrate of 15 L/min even the coarser mesh could be used. The main flow features include:

• Recirculation zones in the upper part of the oral cavity, the posterior part of the upper trachea (just downstream the glottis constriction) and in the left main bronchus.
• Flow acceleration at the glottis constriction (formation of laryngeal jet) and development of a separated shear layer at the level of the vocal cords due to the front inclination of trachea.
• In the lower part of the trachea, the high-speed velocity shifts from the anterior to the posterior wall and this leads to the formation of a small region of separation at the anterior wall.

The mean velocity features remain similar as the flowrate increases but larger velocity fluctuations and turbulent kinetic energy levels are observed.

Figure 12: Contours of the mean velocity magnitude at cross sections in the mouth—throat/trachea (1st row), the main bifurcation (2nd row) and the left/right bronchi (3rd/4th rows) for the three different meshes examined at 15 L/min.

Figure 13 shows deposition results on the three meshes. Figure 13(a) and (b) show overall and mouth-throat deposition fractions as a function of particles size, whereas Figure 13(c) and (h) show deposition fractions per segment for six particle sizes. The numbering of the segments is shown in figure4 (deposition in the 10 larger outlet segments is excluded).

Figure 13: Deposition fractions at 15 L/min for the three examined meshes. (a) and (b) show overall and mouth—throat deposition fractions versus particle size. (c)—(h) show deposition fractions per segment. The numbering of the segments is shown in figure 4.

Results on Meshes 2 and 3 show good agreement, whereas deposition on Mesh 1 is overpredicted for the smaller particle sizes. Based on the analysis conducted, Mesh 2 is used for the LES at 15 and 30 L/min and Mesh 3 is used for the simulation at 60 L/min.

LES results

In this section, deposition results for the LES simulations are reported. Figure 14 shows deposition fractions in the overall geometry (a), in the mouth—throat (b) and the tracheo— bronchial regions (c). Deposition fractions are reported as a function of particle size and for the three flowrates examined. The particle sizes examined are: 0.5, 1, 2, 2.5, 4.3, 6, 8, 10 ${\mu m}$.

Figure 14: Deposition fractions as a function of particle size and for the three flowrates examined in (a) the entire airway geometry: (b) the mouth—throat region: and (c) the tracheobronchial tree.

Figure 15 shows deposition fractions per segment for the eight particle sizes mentioned before and at the three examined flowrates. The numbering of the segments is shown in figure4 (deposition in the 10 larger outlet segments is excluded). Deposition reported in figure 15 can be found in file DF_LES.xlsx.

RANS Simulations

The numerical calculations of the airways system are based on the Euler/Lagrange approach and were conducted using the computational open source code OpenFOAM version 4.1. Further details about the modeling of the problem are detailed as follows.

Computational domain and meshes

The geometry used in the RANS calculations is essentially the same geometry as the one used in the LES case (shown in fig. 10).

The meshing process is one of the most important parts in solving the airways system. A mesh with insufficient quality of the computational cells (high mesh non-orthogonality or skewness) could result in numerical diffusion and consequently prediction of the gas phase with significant errors (Jasak, 1996). The geometry of the airways system is extremely complex with several changes of sections, branches, constrictions and expansions. Therefore, it is not possible to generate a hexahedral and structured mesh. The solution found for this problem lies in the use of tetrahedral elements with this configuration being more adaptable to the complexity of the mesh. However this kind of mesh structure can lead to numerical errors mainly within the boundary layers, e.g. near-wall region, making it necessary to create a layer of prismatic elements close to the wall. Two meshes were generated with different near-wall resolutions. Cross-sectional Views of these meshes at five stations are shown in figure 16. In the coarser mesh (Mesh 1), only three layers of prismatic elements were used close to the wall; however the results obtained with such mesh resolution were not satisfactory. This problem was solved by refining the mesh close to the wall. Specifically, the near-wall element was divided three times having the two parts closer to the wall at 25% of the original element thickness and the third one having a thickness of 50%, as can be observed in Figure 16(b). Table 6 reports grid characteristics for meshes 1 and 2. In section 3.3.4, we assess the effect of mesh resolution near the wall on deposition in the benchmark geometry.

Figure 16: (a) Cross-sectional views of the three generated meshes at five stations (A—E). (b) Cross-sectional views at the entrance pipe, showing coarser mesh with three near-wall prism layers and finer mesh with the near-wall element divided by three.

Solution strategy and boundary conditions — Airflow

Fluid-phase calculations are performed for steady—state and are based on the Euler approach by solving the Reynolds-Averaged Navier-Stokes (RANS) equations in connection with the standard k-ω-SST turbulence model. For coupling the velocity and pressure fields, the SIMPLE (semi-Implicit-Method-Of-Pressure-Linked-Equations) algorithm is applied. To solve the conservation and momentum equations,

---

Contributed by: *** — ***

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

© copyright ERCOFTAC 2019