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

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.

Table 5: Characteristics of Meshes 1 – 3. ${\Delta r_{min}}$ is the height of the wall-adjacent cells. ${\lambda }$ the average expansion ratio of the prism layers and ${V_{cell,avg}}$, the average cell volume in the domain.

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}=\left\{{\begin{array}{ll}{\frac {24}{Re_{p}}}\left(1+0.15{Re}_{p}^{0.687}\right)&{\text{if }}Re_{p}\leq 1000\\0.44&{\text{if }}Re_{p}>1000\end{array}}\right.$

${(17)}$

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 ${\Delta t_{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 figure 4 (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 ${\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.

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 a later section, 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, 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 ${1.184\ kg/m^{3}}$ and ${20.13\times 10^{-6}\ kg/ms}$, respectively.

Table 6: Characteristics of Meshes 1 – 2. ${\Delta r_{min}}$ is the height of the wall-adjacent cells, ${\lambda }$ the average expansion ratio of the prism layers and ${V_{cell,avg}}$, the average cell volume in the domain.

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:

${{\frac {d{\vec {x}}_{p}}{dt}}={\vec {u}}_{p}}$	${(21)}$
${m_{p}{\frac {d{\vec {u}}_{p}}{dt}}=\sum _{i}F_{i}}$	${(22)}$

where ${x_{p}}$ are the coordinates of the particle position, ${U_{p}}$ is the particle linear velocity, ${m_{p}}$ is the mass of the particle and ${\sum _{i}F_{i}}$ is the sum of all forces acting on the particles. The time step of the particle tracking calculation is automatically and independently adjusted along the trajectories by considering all relevant time scales, which are also changing throughout the flow field:

• The time required for a particle to cross a control volume ${\Delta t_{CV}}$
• The particle response time ${\tau _{p}}$
• The integral time scale of turbulence ${T_{L}}$

In order to ensure that the particle tracking is solved properly, the Lagrangian time step ATL must be a fraction of the minimum of the time scales detailed above. In the present calculations, the Lagrangian time step is chosen to be five times smaller than the smaller time scale:

${\Delta T_{L}={\frac {min(\Delta T_{CV},\tau _{p},T_{L})}{5}}}$

${(23)}$

The forces acting on the particle in the RAN tests include Drag, Slip-Shear Lift Force, Brownian Motion and Gravitational Force. Drag force dominates the particle motion in most fluid-particle systems and can be expressed by Sommerfeld et al. (2008)

${{\vec {F}}_{D}={\frac {3}{4}}{\frac {\rho _{f}m_{p}}{\rho _{p}D_{p}}}{\frac {C_{D}}{C_{U}}}({\vec {u}}_{f}-{\vec {u}}_{p})|{\vec {u}}_{f}-{\vec {u}}_{p}|}$

${(24)}$

where ${\rho _{f}}$ and ${\rho _{p}}$ are the fluid and particle density, respectively, ${D_{p}}$ is the particle diameter and ${{\vec {u}}_{f}}$ is the instantaneous fluid velocity. The drag coefficient ${C_{D}}$ is given as a function of the particle Reynolds number ${Re_{p}}$:

$Re_{p}={\frac {\rho _{f}D_{p}|{\vec {f}}_{p}-{\vec {u}}_{p}|}{\mu _{f}}}\left\{{\begin{array}{ll}Re_{p}\leq 0.5&C_{d}={\frac {24}{Re_{p}}}\\0.5<Re_{p}\leq 1000&C_{d}={\frac {24}{Re_{p}}}(1+0.15Re_{p}^{0.687})\\Re_{p}>1000&C_{d}\approx 0.44\end{array}}\right.$

${(25)}$

For small particles the drag coefficient may be remarkably reduced due to the slip effect (often also refereed to rarefaction effect). Hence, this reduction of the drag coefficient may be accounted for by a correction function, the so called Cunningham correlation Cu (valid for: ${Re_{p}<0.025,\ 0.1<Kn<1000}$):

${Cu=1+Kn(2.514+0.8e^{-0.55/Kn})}$

${(26)}$

The importance of the slip effect may be estimated based on the ratio of mean free path of the gas molecules to the particle diameter Sommerfeld et al. (2008), which is the Knudsen number Kn (increasing Knudsen number yields a drastic reduction of the drag coeflicient):

${\lambda ={\frac {\mu _{f}}{0.499{\overline {c_{mol}}}\rho _{f}}},\quad {\overline {c_{mol}}}={\sqrt {\frac {8p}{\pi \rho _{f}}}},\quad Kn={\frac {\lambda }{D_{p}}}{\sqrt {2}}={\frac {\mu _{f}}{D_{p}}}{\sqrt {\frac {\pi \rho _{f}}{p}}}}$

${(27)}$

where ${\overline {c_{mol}}}$ is the standard deviation of the mean relative velocity between gas molecules, ${\lambda }$ is the mean free path of the gas molecules, ${\mu _{f}}$ is the fluid absolute viscosity and ${p}$ is the system pressure.

Especially in the airway systems of lungs with wide range of flow Reynolds numbers, strong shear gradients are found, inducing a transverse lift force onto the particles. According to Crowe et al. (2012), the calculation of this slip-shear lift force is based on analytical results of Saffman (1965) and extended for high particle Reynolds number according to Mei (1992):

${{\vec {F}}_{L,S}=1.615D_{p}\mu _{f}Re_{s}^{1/2}{\frac {\left[({\vec {u}}_{f}-{\vec {u}}_{p})\times {\vec {\omega }}\right]}{|{\vec {\omega }}|^{0.5}}}}$

${(28)}$

where ${{\vec {\omega }}={\vec {\nabla }}\times {\vec {u}}_{f}}$ is the fluid rotation, ${Re_{s}=\rho _{f}D_{p}^{2}|{\vec {\omega }}|/\mu _{f}}$ is the particle Reynolds number of the shear flow and the lift coefficient is expressed as ${C_{LS}=f(Re_{p},Re_{s})={\frac {F_{LS}}{F_{LS,Saff}}}}$ being the ratio of the extended lift force to the Saffman force (Mei, 1992) as given below (Saffman, 1965):

$f(Re_{p},Re_{s})=\left\{{\begin{array}{ll}(1-0.3314\beta ^{\,0.5}e^{(-0.1)Re_{p}}+0.331\beta ^{\,0.5}&{\text{for }}Re_{p}\leq 40\\0.0524(\beta Re_{p})^{0.5}&{\text{for }}Re_{p}>40\end{array}}\right.$

${(29)}$

where ${\beta =0.5Re_{s}/Re_{p}}$.

Sub-micrometre particles are also subjected to Brownian motion, caused by random collisions with the gas molecules. The Brownian force can be modelled as a Gaussian white-noise random process. Li & Ahmadi (1992) studied the effects of Brownian diffusion on particle dispersion and deposition in turbulent channel flow. They observed that Brownian motion is the dominant mechanism for dispersion of particles smaller than 0.1${\mu m}$. For particles larger than 0.5${\mu m}$, the effect of Brownian diffusion was negligibly small.

${{\vec {F}}_{B}=\zeta _{i}m_{p}{\sqrt {\frac {216k_{Boltz}\mu _{f}T_{f}}{\pi D_{p}^{5}\rho _{p}^{2}Cu\Delta t_{L}}}},{\text{with }}0\leq \zeta _{i}\leq 1}$

${(30)}$

where ${k_{Boltz}=1.38\times 10^{-23}J/K}$ is the Boltzmann constant, ${T_{f}}$ is the fluid temperature at particle position, ${\zeta _{i}}$ is a random number with equal probability in the range between zero and one and ${\Delta t_{L}}$ is the particle time-step. The gravitational force including buoyancy effects acting on the particle can be expressed as:

${{\vec {F}}_{G}=m_{p}{\vec {g}}\left(1-{\frac {\rho _{f}}{\rho _{p}}}\right)}$

${(31)}$

Turbulent dispersion

In order to model the effect of turbulence on particle transport, the instantaneous fluid velocity has to be used in the fluid forces given above. From the calculation of the fluid field only the mean velocities are available, wherefore the instantaneous fluid fluctuating velocities have to be generated with the help of the calculated turbulence properties. For that purpose, the Langevin equations suggested by Legg & Raupach (1982) were used to obtain the fluid velocity fluctuation experienced by the particle (Sommerfeld et al., 1993):

${u_{i}^{'\,n+1}=R_{p,i}\left(\Delta t_{L},\Delta r\right)u_{i}^{'\,n}+\sigma _{F}{\sqrt {1-R_{p,i}^{2}(\Delta t_{L},\Delta r)}}\zeta _{i}}$

${(32)}$

where the superscripts denote the time step and the subscripts the spatial component. ${\Delta t_{L}}$ is the Lagrangian time step and ${\Delta r}$ is the spatial separation between the virtual fluid element and the particle during the time ${\Delta t_{L}}$. Here, turbulence is considered to be isotropic so that ${\sigma _{F}}$ represents the rms value of the fluid velocity fluctuation and ${\zeta _{i}}$ denote Gaussian distributed random numbers with zero mean and unit variance. The first term on the right hand side of the equation represents the correlated part, the second term the random contribution to the velocity fluctuation, which depend on the degree of correlation and hence the turbulent Stokes number of the particles. The correlation functions ${R_{p,i}(\Delta t_{L},\Delta r)}$ have Lagrangian and Eulerian components:

${R_{p,i}(\Delta t_{L},\Delta r)=R_{L}(\Delta t_{L})R_{E,ij}(\Delta r)}$

${(33)}$

The Lagrangian correlation describes the instantaneous velocity fluctuation along the way of a virtual fluid element and depends on the Lagrangian integral time scale:

${R_{L}(\Delta t_{L})=\exp \left(-{\frac {\Delta t_{L}}{T_{L}}}\right)}$

${(34)}$

On the other hand the Eulerian correlation reflects the deviation of the particle trajectory from the path of the virtual fluid element, the so-called crossing trajectory effect:

${R_{E,ij}(\Delta r)=\left[f(\delta r)-g(\delta r)\right]{\frac {\Delta r_{i}\Delta r_{j}}{\left|\Delta r_{i}\right|^{2}}}+g(\Delta r)\delta _{ij}}$

${(35)}$

where ${f(\Delta r)}$ and ${f(\Delta r)}$ are the longitudinal and transverse two-point correlation functions (Sommerfeld et al., 1993; von Karman & Horwarth, 1938). The required integral time ${T_{L}}$ and the turbulent length scale of turbulence ${L_{E}}$ are estimated with the turbulent kinetic energy ${k}$ and the dissipation rate ${\epsilon }$:

${T_{L}=c_{T}{\frac {\sigma _{F}^{2}}{\epsilon }},\quad L_{E}=c_{L}\sigma _{F}T_{L},\quad \sigma _{F}={\sqrt {2/3k}},\quad c_{T}=0.24,\quad c_{L}=3.0}$

${(36)}$

From numerous studies related to the dispersion of fine particles, it is known that in anisotropic turbulence such fine particle tend to drift out of regions with high turbulence intensity (Sommerfeld et al., 1993). This is why in the framework of RANS, standard dispersion models yield particle accumulation near the wall and hence increased deposition. To solve this problem, the drift correction should compensate such an unphysical spurious drift of very small particles into regions of low turbulence, like in the near wall region. For that, the scaling factor proposed by Matida et al. (2004) is used in the present study, which modifies turbulence properties depending on the particle-wall distance ${y^{+}}$:

	${T_{E}=f_{v}^{2}k/\epsilon }$	${(37)}$
	${\sigma _{F}=f_{v}{\sqrt {2/3k}}}$	${(38)}$
	${f_{v}=1-e^{-0.02y^{+}}}$	${(39)}$

where ${f_{v}}$ is the corrector factor proposed by Matida et al. (2004) and ${y^{+}}$ is the dimensionless distance to the wall.

The particle phase consists of spheres having a density of 914 ${kg/m^{3}}$ (di-ethylhexyl sebacate) and different sizes according to the experiments executed by the Brno University of Technology team (i.e. 0.5—10${\mu m}$). For representing the particle phase, 100,000 particles are injected. The particle injection velocity is the bulk gas velocity at the inlet in the stream-wise direction and zero in the transverse components, plus a random fluctuation drawn from a normal distribution with a rms value of 3% of the gas bulk velocity for all three particle velocity components. Because of the small size of the particles, no angular velocity is considered (no particle rotation). In all calculated cases the implemented turbulent dispersion model and the relevant forces (drag force, slip-shear lift force, gravitational force and Brownian force) are accounted for. The gravitational force is included in the negative z directions of the coordinate system (see Figure 18).

Numerical accuracy

As mentioned before, two different meshes were used to examine the accuracy of the RANS scheme, with emphasis on the near-wall resolution (meshes shown in figure 16). The improvement of the mesh layer resolutions resulted in a better agreement with the experimental data, as well as with the results obtained using Large Edge Simulations (LES), as observed in the deposition fractions as function of particle size (figure 17) for the case of 60 L/min and different particle sizes ranging from 0.5-10${\mu m}$. In figure 18, the deposition pattern of particles with size 4.3${\mu m}$ in the airways system is shown. It is clear that with the improvement of the resolution of the mesh layer near the wall, a remarkable agreement with the experimental data is obtained. This is mainly due the fact that to calculate the particle dispersion, it is necessary to obtain properties from the continuous phase at particle position, such as mean flow velocities and turbulence properties. The implemented models in OpenFOAM 4.1 by the authors take into account the interpolation of such properties for the particle positions; however it is necessary to have this extra refinement in the near wall region for a better solution of these properties. For all further calculations, only the refined mesh is considered.

Figure 17: Deposition Fraction as function of particle size for different numerical methods, different grid size (flow rate 60 L/min).

Figure 18: Colour maps of the deposition pattern in the airways system obtained for a gas flow of 60 L/min and a particle size of 4.3${\mu m}$.

RANS results

In this section, deposition results for the RANS simulations are reported. Figure 19(a) shows deposition fractions in the overall geometry as a function of particle size and for the three flowrates examined. Figure 19(b)—(d) show deposition fractions per segment at the three examined flowrates. Deposition fractions reported in figure 19 can be found in file DF_RANS.xlsx.

Figure 19: RANS results: (a) overall deposition fractions as a function of particle size and for the three flowrates examined: (b)—(d) deposition fractions in the segments.

---

Contributed by: P. Koullapis, F. Lizal, J. Jedelsky, L. Nicolaou, K. Bauer, O. Sgrott, M. Jicha, M. Sommerfeld, S. C. Kassinos — University of Cyprus

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

© copyright ERCOFTAC 2019