CFD Simulations AC7-02

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

Airflow in the human upper airways

Application Challenge AC7-02   © copyright ERCOFTAC 2020

Overview of CFD Simulations

Three LES (LES1-3) and one RANS simulation were carried out in the benchmark geometry at an air flowrate of 60 L/min. This flowrate results in a Reynolds number of 4921 in the model’s trachea. This is slightly higher than the Reynolds number in the experiments, due to the reasons discussed in section 2.4. The details of the numerical tests are given in the following paragraphs. In summary, the main differences between the simulations are:
1. Computational meshes: LES1&2 employ the same comp. meshes, whereas LES3 and RANS use different meshes.
2. Turbulence modeling: LES1 uses the dynamic version of the Smagorinsky-Lilly subgrid scale model, whereas LES2&3 employ the Wall-adapting eddy viscosity (WALE) SGS model. In RANS simulations, the k-ω SST, standard k-ε and RNG k-ε turbulence model have been tested.
3. Discretisation method: LES1&2 and RANS use the Finite Volume approach whereas LES3 employs the Finite Element Method.

Large Eddy Simulations — Case LES1

Computational domain and meshes

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

Figure 8: 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 LES1, two meshes were generated to allow us to examine the sensitivity of the results to the mesh resolution. The coarser mesh includes 10 million computational 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 Fig. 9. A grid convergence analysis was carried out in order to determine the appropriate resolution for the simulations. This analysis is presented in section 3.2.3. Table 3 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 9: Cross-sectional views of the three generated meshes employed in LES1&2 at seven stations (A-G).

Table 3: 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.

Solution strategy and boundary conditions — Airflow

LES were performed using the dynamic version of the Smagorinsky-Lilly subgrid scale model with localized filtering (Lilly, 1992; Zang et al., 1993) 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 {\overline {u}}_{j}}{\partial x_{j}}}=0\qquad (2)$

${\frac {\partial {\overline {u}}_{i}}{\partial t}}+{\overline {u}}_{j}{\frac {\partial {\overline {u}}_{i}}{\partial x_{j}}}=-{\frac {1}{\rho }}{\frac {\partial {\overline {p}}}{\partial x_{i}}}+{\frac {\partial }{\partial x_{j}}}{\bigg [}(\nu +\nu _{sgs}){\frac {\partial {\overline {u}}_{i}}{\partial x_{j}}}{\bigg ]}\qquad (3),$

where ${\overline {u}}_{i}$, p, $\rho =1.2kg/m^{3}$, $\nu =1.59\times 10^{-5}m^{2}/s$ and $\nu _{sgs}$ are the filtered velocity component in the i-direction, the filtered 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 (OpenFOAM Foundation, 2013b,a). 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 (OpenFOAM Foundation, 2013b). The temporal derivative is discretized using backward differencing, which is is also second order accurate in time and implicit. To ensure numerical stability the time step used is $2.5\times 10^{-6}s$ at inhalation flowrate of 60 L/min. 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.

At the inhalation flowrate of 60 L/min, the Reynolds number at the inlet, based on the inflow bulk velocity and inlet tube diameter, is $Re_{in}={\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, as shown in fig. 10. 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 extended pipe’s inlet boundary (fig. 10). 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. To match the experimental conditions, uniform pressure is prescribed in the simulations at the 10 terminal outlets. A no-slip velocity condition is imposed on the airway walls.

Figure 10: Explanatory figure for the mapped inlet, or recycling, boundary condition applied at the inlet of the model.

Numerical accuracy

The sensitivity of the calculations to the mesh resolution was tested and the results are presented in this section. Fig. 11 displays contours and profiles of the mean velocity magnitude at cross sections in the mouth-throat/trachea, the main bifurcation and the left/right bronchi for the two different meshes examined. Good agreement is observed for the mean velocities between the coarse and fine meshes. Slight deviations are found at station D1-D2, located just downstream of the glottis. Airflow recirculation is evident at this station that is not well captured on the coarse mesh.

Fig. 12 displays contours and profiles of the turbulent kinetic energy (TKE) at cross sections in the mouth-throat/trachea, the main bifurcation and the left/right bronchi for the two different meshes examined. TKE is calculated as:

$TKE={\frac {1}{2}}<u_{i}'u_{i}'>$

Higher levels of TKE are recorded on the finer mesh. These are mainly generated in the shear layers. As observed from the 2D profiles, TKE peaks are underpredicted on the coarse mesh. In conclusion, although the coarse mesh yields results that are in reasonable agreement with the finer mesh predictions, in the following comparisons between CFD and PIV measurements, the finer LES predictions are considered.

Figure 11: LES1: Comparison of mean velocity contours and profiles at selected stations between coarse (Mesh 1) and fine (Mesh 2) meshes.

Figure 12: LES1: Comparison of turbulent kinetic energy contours and profiles at selected stations between coarse (Mesh 1) and fine (Mesh 2) meshes.

Large Eddy Simulations — Case LES2

Computational domain and meshes

The computational domain and the meshes employed in these simulations are the ones described in Section 3.2.1.

Solution strategy and boundary conditions — Airflow

In LES2, an eddy-viscosity-type model following the Boussinesq hypothesis (Sagaut, 2006) is used to close the system of filtered Navier-Stokes equations, (2&3). The Wall-adapting eddy viscosity model (WALE) SGS model (Nicoud & Ducros, 1999) has been employed. This model is based on the square of the velocity gradient tensor. The SGS viscosity obtained with this model takes into account the strain and the rotation rate of the smallest resolved turbulent fluctuations. Some features of this model are its capability of switching off in two-dimensional flows and in laminar flows. It also has a cubic behaviour near walls with respect to the normal direction of the wall.

The CFD simulations presented in this section have been carried out using the in-house CFD software TermoFluids (Lehmkuhl et al., 2007), based on the finite volume method (FVM). This CFD code is parallel, highly scalable and designed to work in both structured and unstructured meshes. The convective term in the momentum equation (3) is discretized using a second order Symmetry-Preserving (SP) scheme (Verstappen & Veldman, 2003). This discretization scheme constructs an anti-symmetric discrete convective operator. This operator does not introduce artificial dissipation in the momentum equation, and therefore the kinetic-energy is preserved. The diffusive operator is discretized by means of a second-order Central Differencing Scheme (CDS). Both convective and diffusive operators are integrated explicitly by means of a one-leg second-order time-integration scheme (Trias & Lehmkuhl, 2011b). This scheme includes a free-parameter allowing to adapt its stability region. The free-parameter is dynamically selected maximizing the stability region in function of the eigenvalues of the discrete operators. Moreover, this time-integration strategy also selects the optimal time-step at each iteration. The pressure-velocity coupling is solved by means of the Fractional Step projection method (Kim & Moin, 1985). The idea behind this technique is to split the momentum within two steps: a first explicit step where an intermediate velocity is obtained, followed by a second step where the pressure is solved implicitly and the intermediate velocity is corrected obtaining the physical divergence-free velocity field. The Poisson equation is solved by means of the iterative Conjugate-Gradient (CG) method with Jacobi diagonal scaling.

The presented case has an inlet flowrate Q=60 L/min, which falls within the turbulent regime. Therefore, in order to generate the turbulent velocity profile for the inlet section, an auxiliary simulation of a pipe with the same diameter as the inlet and periodic boundary condition in the streamwise direction have been employed. The velocity field generated at the mid-plane of the pipe is saved and later imposed at the inlet section of the simulation during running time. At the respiratory airways walls, the boundary condition is defined as no-slip. Regarding the outlets, two different conditions were employed. For the cases presented in section 3.3.3.1, a constant flowrate was fixed en each one of the outlets. On the other hand, for the cases solved in section 3.3.3.2 and in Section 4, the pressure is fixed at the ten outlets with a value equal to atmospheric pressure.

Numerical accuracy

Two of the main numerical aspects that will determine the accuracy of a CFD simulation employing a LES approach are the mesh and the turbulence model employed for the subgrid scales. Therefore, in the present section these two parameters have been analyzed in detail.

Mesh convergence analysis

The effects of the mesh on the simulation results are studied comparing the same experiment and geometry using two different meshes: a fine one with 50 million control volumes (CVs) and a coarse one with 10 million CVs (see 3.2.1). The results for both mean velocities and TKE ($=1/2<u_{i}'u_{i}'>$) obtained with the two meshes are depicted in Fig. 13 and 14, respectively.

As can be seen in Fig. 13, the agreement for the mean velocities is very good between both meshes and the results are practically identical. On the other hand, some differences can be appreciated for TKE between the two studied meshes. As observed in the 2D profiles (fig. 14), the fine mesh yields higher levels of TKE than the coarse one. However, both meshes are able to capture the same trend, obtaining the highest levels of TKE in the same positions, which belongs to the shear layer regions. Hence, it is clear that, although first order quantities are well captured by both meshes, special care must be taken if second order quantities must be reproduced accurately, since the latter are very sensitive to mesh resolution. In conclusion, sufficiently fine meshes must be employed in order to correctly capture second order quantities, although first order ones can be obtained accurately with coarser meshes. Obviously, the recommended grid-spacing in this kind of simulations will be the result of a compromise between accuracy and computational cost.

Figure 13: LES2: Comparison of mean velocity contours and profiles at selected stations between coarse (Mesh 1) and fine (Mesh 2) meshes.

Figure 14: LES2: Comparison of turbulent kinetic energy contours and profiles at selected stations between coarse (Mesh 1) and fine (Mesh 2) meshes.

Comparison of different LES subgrid-scale models

In order to assess the influence of the subgrid-scale turbulence model in LES of airflow in the human upper airways, four different turbulence models were compared. Besides the WALE model (Nicoud & Ducros, 1999) previously introduced, three additional turbulence models have been studied. The first one is the model proposed by Smagorinsky (1963), based on the Prandtl mixing length applied to subgrid-scale modelling. In this model the turbulent viscosity is proportional to the strain. The second model is the variational multiscale (VMS) approach applied in WALE. This method was originally formulated for the Smagorinsky model by Hughes et al. (2000). Using this framework the modelling is confined to the effect of a small-scale Reynolds stress, in contrast to classical LES in which the entire SGS stress is modelled. The latter is the model proposed by Verstappen known as the QR eddy-viscosity model (Verstappen et al., 2010; Verstappen, 2011). This model is based on two invariants of the rate-of-strain tensor. The method is computationally very efficient, switches off in the case of back-scattering, laminar and two-dimensional flows, and the subgrid turbulent viscosity is proportional to the cube with respect to the normal direction of the wall.

The results for the in-plane mean velocities and TKE (see section 4 and equations 8-11) are depicted in Fig. 15 and 16, respectively. As can be seen, the four models predict very similar results, except in shear layer regions for TKE levels (just downstream glottis constriction - station D), where some differences can be seen between the different models. However, these deviations are very small and not relevant. Therefore, it can be stated that all LES subgrid-scale models are able to provide results that are in reasonable agreement with the measurements.

Nonetheless, in previous works where different turbulence models in confined flows were compared, it was found that some turbulence models were not able to correctly predict the flow pattern. In the work of Muela et al. (2019), the Smagorinsky model was not able to correctly model the subgrid dissipation in the simulated confined flow, being too dissipative. The main difference to the cases examined in these works is the Reynolds number. In the present study, the Reynolds number in the trachea at 60 L/min is around Re = 4921, while the one of the case presented in Muela et al. (2019) is two orders of magnitude higher ($Re=3.47\cdot 10^{5}$).

In simulations of the airflow in the human upper airways, the Reynolds number is in most cases below Re < 10000 and therefore the flow is either laminar or with low turbulence. Due to that, the influence of the subgrid scales in this kind of flows is small, and the turbulence model is not as relevant as in more turbulent cases. Therefore, the turbulence model in simulations of the airflow in the human upper airways using LES modeling does not play a key role.

Figure 15: Comparison of normalised mean velocity profiles at selected stations -shown in (a)- between PIV and different LES subgrid-scale models.

Figure 16: Comparison of normalised turbulent kinetic energy profiles at selected stations -shown in fig. 14(a)- between PIV and different LES models.

RANS Simulations

Numerical accuracy

---

Contributed by: P. Koullapis^a, J. Muela^b, O. Lehmkuhl^c, F. Lizal^d, J. Jedelsky^d, M. Jicha^d, T. Janke^e, K. Bauer^e, M. Sommerfeld^f, S. C. Kassinos^a — 
^aDepartment of Mechanical and Manufacturing Engineering, University of Cyprus, Nicosia, Cyprus
^bHeat and Mass Transfer Technological Centre, Universitat Politècnica de Catalunya, Terrassa, Spain
^cBarcelona Supercomputing center, Barcelona, Spain
^dFaculty of Mechanical Engineering, Brno University of Technology, Brno, Czech Republic
^eInstitute of Mechanics and Fluid Dynamics, TU Bergakademie Freiberg, Freiberg, Germany
^fInstitute Process Engineering, Otto von Guericke University, Halle (Saale), Germany

Front Page

Description

Test Data

CFD Simulations

Evaluation

Best Practice Advice

© copyright ERCOFTAC 2020