CFD Simulations AC7-01
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 speciﬁed 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 ﬁner 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 , 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 – Airﬂow
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 ﬂow in the realistic airway geometries. Previous studies have shown that this model performs well in transitional ﬂows in the human airways (Radhakrishnan & Kassinos, 2009; Koullapis et al., 2016). The airﬂow is described by the ﬁltered 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.
Table 5: Characteristics of Meshes 1 – 3. is the height of the wall-adjacent cells. the average expansion ratio of the prism layers and , the average cell volume in the domain. |
The governing equations are discretized using a ﬁnite volume method and solved using
OpenFOAM, an open-source CFD code. In this framework, unstructured boundary ﬁtted
meshes are used with a collocated cell-centred variable arrangement. The ﬁnite 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 signiﬁcant 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 ﬂow 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 ﬂowrates 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 proﬁle is imposed. For the higher ﬂowrate of 60 L/min, the Reynolds number at the inlet, based on the inﬂow bulk velocity and inlet tube diameter, is which lies in the turbulent regime. In order to generate turbulent inﬂow 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 ﬁeld generated in a precursor pipe ﬂow LES. During the simulation of the airway geometry, the velocity ﬁeld from the mid-plane of the pipe domain was mapped to the inlet boundary. Scaling of the velocities was applied to enforce the speciﬁed bulk ﬂow rate. In this manner, turbulent ﬂow is sustained in the extended pipe section, and a turbulent velocity proﬁle enters the mouth inlet. The volumetric ﬂowrates 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 deposition 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 ﬂuid velocity interpolated at the position of the particle, and is the particle response time, deﬁned as:
with being the particle diameter, the dynamic ﬂuid viscosity and the particle Reynolds number. is the Cunningham correction factor, deﬁned 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 coefﬁcient, , 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 coefﬁcient, 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. Reﬂection 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 ﬂow-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 eﬂicient 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 ﬂow ﬁeld 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 ﬂow 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 eﬂicient 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 ﬂow 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 ﬁnal destination. Thus, the maximum time step used to track a particle is the one deﬁned for the continuous phase simulation . At the particle new position, which is either on a face that has been crossed or the ﬁnal 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 ﬂuid 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 ﬂow equations while is computed once in the beginning of the simulation for each cell and ﬁtted 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 ﬁnds it, linear interpolation is performed based on the 4 vertices of the tetrahedron (Baker, 2003). This strategy leads to an interpolated ﬂuid 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 ﬁnal destination at the end of the ﬂow time step.
Numerical accuracy
The sensitivity of the calculations (airﬂow ﬁelds and particle deposition) on the mesh resolution was tested at the lower ﬂowrate 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 ﬂowrate of 15 L/min even the coarser mesh could be used. The main ﬂow 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 ﬂowrate increases but larger velocity ﬂuctuations 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
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 ﬁgure 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 ﬂowrates 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 ﬂowrates.
The numbering of the segments is shown in
figure4 (deposition in the 10 larger outlet segments is excluded). Deposition reported in
ﬁgure 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 signiﬁcant 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 conﬁguration 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 ﬁve 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 reﬁning the mesh close to the wall. Speciﬁcally, 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.
Solution strategy and boundary conditions — Airﬂow
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 ﬁelds, 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 ﬂow 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.
Table 6: Characteristics of Meshes 1 – 2. is the height of the wall-adjacent cells, the average expansion ratio of the prism layers and , 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 reﬁnements near the
wall for a proper solution. A mapped boundary condition is used at the inlet in order to
obtain a developed inlet ﬂow 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: Conﬁguration of the boundary conditions for the gas calculations considering a gas ﬂow 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:
where are the coordinates of the particle position, is the particle linear velocity, is the mass of the particle and 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 ﬂow ﬁeld:
- The time required for a particle to cross a control volume
- The particle response time
- The integral time scale of turbulence
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:
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 ﬂuid-particle systems and can be expressed by Sommerfeld et al. (2008)
where and are the fluid and particle density, respectively, is the particle diameter and is the instantaneous ﬂuid velocity. The drag coefficient is given as a function of the particle Reynolds number :
For small particles the drag coefﬁcient may be remarkably reduced due to the slip effect (often also refereed to rarefaction effect). Hence, this reduction of the drag coefﬁcient may be accounted for by a correction function, the so called Cunningham correlation Cu (valid for: ):
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 coeﬂicient):
where is the standard deviation of the mean relative velocity between gas molecules, is the mean free path of the gas molecules, is the ﬂuid absolute viscosity and is the system pressure.
Especially in the airway systems of lungs with wide range of ﬂow 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):
where is the ﬂuid rotation, is the particle Reynolds number of the shear flow and the lift coefficient is expressed as being the ratio of the extended lift force to the Saffman force (Mei, 1992) as given below (Saffman, 1965):
where .
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 ﬂow. They observed that Brownian motion is the dominant mechanism for dispersion of particles smaller than 0.1. For particles larger than 0.5, the effect of Brownian diffusion was negligibly small.
where is the Boltzmann constant, is the ﬂuid temperature at particle position, is a random number with equal probability in the range between zero and one and is the particle time-step. The gravitational force including buoyancy effects acting on the particle can be expressed as:
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 ﬂuid ﬂuctuating 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 ﬂuid velocity ﬂuctuation experienced by the particle (Sommerfeld et al., 1993):
where the superscripts denote the time step and the subscripts the spatial component. is the Lagrangian time step and is the spatial separation between the virtual fluid element and the particle during the time . Here, turbulence is considered to be isotropic so that represents the rms value of the fluid velocity ﬂuctuation and denote Gaussian distributed random numbers with zero mean and unit variance. The ﬁrst term on the right hand side of the equation represents the correlated part, the second term the random contribution to the velocity ﬂuctuation, which depend on the degree of correlation and hence the turbulent Stokes number of the particles. The correlation functions have Lagrangian and Eulerian components:
The Lagrangian correlation describes the instantaneous velocity ﬂuctuation along the way of a virtual fluid element and depends on the Lagrangian integral time scale:
On the other hand the Eulerian correlation reﬂects the deviation of the particle trajectory from the path of the virtual fluid element, the so-called crossing trajectory effect:
where and are the longitudinal and transverse two-point correlation functions (Sommerfeld et al., 1993; von Karman & Horwarth, 1938). The required integral time and the turbulent length scale of turbulence are estimated with the turbulent kinetic energy and the dissipation rate :
From numerous studies related to the dispersion of ﬁne particles, it is known that in anisotropic turbulence such ﬁne 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 modiﬁes turbulence properties depending on the particle-wall distance :
where is the corrector factor proposed by Matida et al. (2004) and is the dimensionless distance to the wall.
The particle phase consists of spheres having a density of 914 (di-ethylhexyl sebacate) and different sizes according to the experiments executed by the Brno University of Technology team (i.e. 0.5—10). 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 ﬂuctuation 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. In figure 18, the deposition pattern of particles with size 4.3 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 ﬂow 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 reﬁnement in the near wall region for a better solution of these properties. For all further calculations, only the reﬁned 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. |
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 ﬂowrates examined. Figure 19(b)—(d) show deposition fractions per segment at the three examined ﬂowrates. 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^{a}, F. Lizal^{b}, J. Jedelsky^{b}, L. Nicolaou^{c}, K. Bauer^{d}, O. Sgrott^{e}, M. Jicha^{b}, M. Sommerfeld^{e}, S. C. Kassinos^{a} —
^{a}Department of Mechanical and Manufacturing Engineering, University of Cyprus, Nicosia, Cyprus
^{b}Faculty of Mechanical Engineering, Brno University of Technology, Brno, Czech Republic
^{c}Division of Pulmonary and Critical Care, School of Medicine, Johns Hopkins University, Baltimore, USA
^{d}Institute of Mechanics and Fluid Dynamics, TU Bergakademie Freiberg, Freiberg, Germany
^{e}Institute Process Engineering, Otto von Guericke University, Halle (Saale), Germany
© copyright ERCOFTAC 2019