CFD Simulations AC7-01: Difference between revisions
Line 438: | Line 438: | ||
The positions and velocities of the rigid and spherical particles are calculated by integrating | The positions and velocities of the rigid and spherical particles are calculated by integrating | ||
the following equations: | the following equations: | ||
{|align="center" border="0" cellpadding="0" | |||
|width=630 align=center|<math>{\frac{d\vec{x}_p}{dt}=\vec{u}_p</math> | |||
|width=20 align=right|<math>{(21)}</math> | |||
|- | |||
|align=center|<math>{\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],}</math> | |||
|align=right|<math>{(22)}</math> | |||
|} | |||
<br/> | <br/> | ||
---- | ---- |
Revision as of 15:22, 8 October 2019
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:
- Computational meshes
- Turbulence modeling
- Different outlet boundary conditions
- 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 , the number of prism layers near the walls, the average expansion ratio of the prism layers (), the total number of computational cells, the average cell volume () and the average and maximum values. The higher 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,
where
and 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 for the cases of 15, 30 L/min and 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 and 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 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, :
where is the particle mass, and and 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:
The drag force acting on the spherical particles is given by,
where is the fluid velocity interpolated at the position of the particle, and is the particle response time, defined as:
with being the particle diameter, the dynamic fluid viscosity and the particle Reynolds number. is the Cunningham correction factor, defined as:
where 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 m. The drag coefficient, , is based on the correlation proposed by Schiller & Naumann (1935):
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),
, |
where is a zero mean variant from a Gaussian probability density function, T=310K is the absolute temperature in the lungs, is the Brownian diffusion coefficient, K is the Boltzmann constant and 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 to , and the particle density was set to , 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 , the advancement of all particles follows for the same time window. Firstly, the particles new position is computed explicitly using the particle velocity calculated at the end of the previous time step ():
The time step 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 , 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 . 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:
The non—linear terms in the above equation (such as ) are linearised using the particles known properties from the previous time step. Spatial interpolation is used to approximate the fluid velocity and the minimum distance vector 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: is obtained at each time step from the solution of flow equations while 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 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 tracheobronchial 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 .
![]() |
![]() |
![]() |
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.
![]() |
Figure 15: Deposition fractions per segment (LES) |
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.
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, the steady—state solver for incompressible and turbulent flow simpleFOAM is used. More details about the solver used, as well the modelling behind the solution, may be found in the OpenFOAM user guide. The gas density and the dynamic Viscosity are set to and , respectively.
A no-slip boundary condition is used at the pipe walls and a zero pressure is assigned to all outlet boundaries. Wall functions are applied in order to solve the turbulence properties in the near wall region, therefore not demanding extra refinements near the wall for a proper solution. A mapped boundary condition is used at the inlet in order to obtain a developed inlet flow without the need of simulating a longer pipe at the entrance of the lung model. This procedure is done setting an averaged value for the desired property (e.g. velocity, turbulent kinetic energy, etc.) at the inlet and an offset distance (0.002m from the inlet in the considered cases) for a reference plane. The mapped BC takes the value of the desired property from this offset plane and uses it as input on the inlet BC for the next iteration step. More information about all the boundary conditions used in the calculations, as well as all the wall functions and properties needed in the calculations are extensively detailed in the OpenFOAM user guide. A summary of the operational conditions and the boundary condition setup is shown in Table 7.
![]() |
Table 7: Configuration of the boundary conditions for the gas calculations considering a gas flow of 60 L/min. |
Solution strategy — Particle transport and deposition
The dispersed phase is treated in the Lagrangian framework, in which the particles are simultaneously and time-dependently tracked through the domain. The particles are treated as point masses and their shape is assumed to be spherical. Only one-way coupling simulations are considered, since the particle size and the volume fraction are very small to have a relevant effect on the continuous phase.
Particle equation of motion
The positions and velocities of the rigid and spherical particles are calculated by integrating the following equations:
Failed to parse (syntax error): {\displaystyle {\frac{d\vec{x}_p}{dt}=\vec{u}_p} | |
Contributed by: *** — ***
© copyright ERCOFTAC 2019