CFD Simulations AC7-01

From KBwiki
Jump to navigation Jump to search

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.


AC7-01 fig10.png
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.

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



Contributed by: *** — ***

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice


© copyright ERCOFTAC 2019