Lattice Boltzmann Simulation of the Air Flow Through the Children Respiratory Tract

. Aerosol particles like dust involved in the air can pass through the respiratory tract and can be deposited there. The long-term deposition of ﬁne solid particles in the lower airways can stimulate lung diseases. In contrary, for the healing of certain lung diseases, the delivering of medicament’s to the lower airways is important. Computational Fluid Dynamics (CFD) enables to simulate the ﬂow ﬁeld and the particles deposition in the lungs. By this method, it is possible to identify the probability of reaching particles into the lower airways through the tracheo-bronchial tree. In this paper, the ﬂow ﬁeld in the child lungs (5-year-old) in a constant sedentary breathing regime was simulated using the Lattice Boltzmann Method (LBM) with the LES Smagorinsky turbulence model. To verify the calculated velocities by LBM in OpenLB software, the ﬁnite volume RANS simulation in Star-CCM + was performed. The results showed the good agreement in the upper part of the airways. Some discrepancies were found in the lower part. Nevertheless, LBM and OpenLB were proven to be a capable tool to solve the air ﬂow through the respiratory tract and will be used for particle deposition.


Introduction
The last decades the research of aerosol particles and its deposition in the human lungs was carried out across the scientific community.Inhalation of these particles can cause diseases, e.g.lung cancer or pulmonary diseases.For instance, Raaschou-Nielsen et al. [1] revealed a correlation between long-term exposure to ambient air pollution and lung cancer incidence in European populations.They analysed data from 17 cohort studies in Europe countries and showed a statistical significant correlation between 10 µm particulate matter (PM 10 ) and lung cancer.On the other hand, S.P.Newmann in [2] discusses the importance of pulmonary drug delivery into the lungs which is successfully used for the treatment of diseases, i.e., asthma or chronic obstructive pulmonary disease (COPD).One possibility to reduce the effect of the behavioural mechanical barrier for preventing the inflow of the drugs into the lungs or airways is the usage of inhalator devices.To investigate the inhalation efficiency or aerosol spreading inside the respiratory tract the particle tracking method can be used.In the simulations it means to calculate the velocity field across the respiratory system and then simulate the particle trajectories.By these simulations, it is possible to evaluate phenomena, which are hard to measure on living humans.Therefore, an increased interest in the flow conditions inside the respiratory tract leads to a number of studies using computational fluid dynamics (CFD) methods.
To perform the CFD simulations of the airways, the 3D geometry model must be prepared.Earlier, simplified models were used (idealized symmetrical Weibel model [3], more sophisticated Horsfield model [4]) resulting in high error rates.In the last decades screening techniques, such as Computer Tomography (CT) and medical imaging methods, can be performed to obtain realistic 3D models of a respiratory tract of real human.These simulations on realistic geometry provide further information about local particle deposition inside the respiratory system.Breathing is a naturally transient process and could be distinguished into inspiratory and respiratory cycles and both are locally laminar and turbulent.Elcner et.al. [5] showed a good agreement between the velocity values obtained by RANS simulation and experimental data obtained by Phase Doppler Anemometry (PDA) at appropriate locations in the tracheobronchial airways.They investigated the temporal evolution during the respiration cycle and at distinct points in the geometry the findings resulted in numerical discrepancies.The particle deposition was not carried out in this study.With the improvement of the computational technologies, direct numerical simulation (DNS) gives even more accurate results.Elghobashi et.all ( [6]) applied the Lattice Boltzmann method for DNS of the nose cavity and upper part of the trachea at the maximal inspiration mouth velocity with a tidal volume of 30 l/min.The geometry domain contained 208 millions voxels.Lintermann and Schröder [7] simulated the deposition of particles in the upper tracheobronchial tract.Two particle tracking simulation runs by the Lagrangian approach were computed with various particle radii.The analysis of the particle deposition showed that 32% of the heavier and larger particles deposit deeper in the airways, while for the particle mixture only 0.69% deposit in the primary, lobar, and segmental bronchi.This investigation reveals that small and light particles penetrate into the lung regions below the sixth generation.Nevertheless, the laryngeal tract was not assembled to the geometry which could influence the inflow conditions and thus impact the resulting deposition.Furthermore, they consider just the steady state solution and the time-dependent dynamic of the inspiration was not taken into account.
This paper describes the first step of our research, i.e., verification of the LBM in OpenLB to simulate velocities in the airways properly with the finite volume method using Star-CCM+.Particularly, these velocity profiles could be used to simulate the particle deposition (in the following periods of our study).
At first, the numerical method LBM is briefly introduced, then the simulation parameters and boundary and initial conditions are listed.In the results, the air velocities in predetermined points are compared between both numerical methods.Finally, relevant statements and ideas are discussed and the main findings are highlighted in the conclusion.

Methods
The Lattice Boltzmann method is a mesoscopic method based on the molecular dynamics.It is a statistical approach considering the probability of the position and velocities of particles in the domain.The governing equation is the Boltzmann equation [8] describing the time change of the particle probability distribution function f .The discretisation process derived by Chapman and Enskog [9] lead to a discretised set of Lattice-Boltzmann equations.It was proved the equivalence with the well-known macroscopic Navier-Stokes equations of fluid dynamics.Each time step of the equation is divided into a collision step and a streaming step.The whole algorithm is schematically demonstrated in the Figure 1.This approach is from the family of Direct Numerical Simulations (DNS) approaches which are very time expensive.Therefore, especially for turbulent flows, a combination with a Large Eddy Simulation (LES) model is effi-cient.The most common choice of a subgrid scale model is the Smagorinsky model [11].In this model, only eddies larger than a defined grid length are simulated while the smaller ones are being modeled.Therefore, a modeled artificial turbulent viscosity is added to the kinematic viscosity.The filtering length is determined by the grid size length multiplied by the Smagorinsky constant C f [12].
The resulting velocity field can be used for the simulation of particle motion (not part of this study).The aerosols can be represented by spherical solid particles and their motion is tracked by the one-way Lagrangian approach, i.e., the movement of the particles does not influence the air flow.If the value of particle Reynolds number is much less than the acting force on the particles, approximation by the Stokes force [7] is allowed.

Simulation setup
The 3D geometry of the airways used in this study was created by scaling of the adult airway model which was produced at Brno University of Technology and used for computational and experimental purposes [5].The original model involves the mouth and nose cavity, nosopharynx, larynx, and the tracheobronchial tree up to the 7th branching generation.For the scaling properties, the information from [13] and [14] were used.In this paper, the trachea is determined as the distance between the carina of the first bifurcation and the glottis.The distance in our model is 136 mm and was reduced for the 5-yearold child to 73 mm.Furthermore, the tracheobronchial tree was simplified from seven to two branching generations (see Fig. 2).Respiration is considered as a combined breathing, i. e., both through the oral and nasal airways, extended by a inhalator mask to investigate the efficiency of medical inhalation.
M. Roy and C. Courtey [15] measured the human respiration parameters for distinct age groups.According to them, the tidal volume for a 5-year child in a sedentary regime is 10.5 l/min with a respiration cycle period 1.607 s.The data presented in [15] lead to specific flow rates through the inlet and outlets of each trachea branch in the state of averaged flow where constant piston profiles were considered.
The airway walls are assigned to a no-slip interpolated boundary condition for curved walls derived by Bouzidi et al. [16].At the inhalator inlet a Dirichlet boundary condition with a constant velocity profile of 2.203 m s −1 .This corresponds to the Reynolds number 1190 which is still in the laminar regime.However, deeper in the human airways a turbulent regime with vortical structure is developed.For the stabilisation of the flow we assigned a sinusoidal range up of the velocity from zero to the maximal value for a time period of 0.25 s.After that the velocity at the inlet boundary was hold at the maximal velocity which expects the averaged value of the sedentary respiration.The simulated fluid is air at normal atmospheric state, e.g., pressure of 101 325 Pa and constant temperature 20 °C, density of 1.3 kg m −3 and a kinematic viscosity of 1.67 mm 2 s −1 .At the starting time, the airways are occupied by the described air with zero velocity as the initial condition.The created computational mesh was uniform with a grid length of 176 µm and time step 1.764 71 × 10 −6 s.This results in 5.7 millions voxels in the whole geometrical domain.The stability is governed by the relaxation time τ and the collision operator depends on this parameter.The value is given by the kinematic viscosity and time and spatial resolution.In this simulation the considered relaxation parameter results in 0.50314 which leads to a stable solution.Another important stability issue affecting the simulation is the initial condition.Dramatic velocity gradients in the domain reduces the stability and therefore a starting period at the velocity boundary is added.This period lasts 0.25 s in which the velocity magnitude grows sinusoidal from zero to the maximal value.The total simulation time was 1.5 s.We suppose that the stabilisation of the flow takes 1 s and then the last 0.5 s was determined for the resulting averaged velocity to outfilter the velocity disturbances produced by the turbulent character of the flow.The simulation is triggered with an applied LES subgrid scale shear-improved Smagorinsky model employed by Jafari et al. [11].The Smagorinsky parameter was assigned to a value of 0.1.This significantly reduces the computational resources and still computes velocity field with a sufficient accuracy to perform particle tracking.
The simulations were obtained by OpenLB software, release 1.4-0 [17], [18].This software has been developed at Karlsruhe Institute of Technology (KIT) by the Lattice Boltzmann Research Group (LBRG) of Dr. Krause.The C++ code of OpenLB software enables an efficient parallelisation of the solving problems, which divides the computation into more threads to reduce the simulation time significantly.The simulations were performed on the CPU cluster Barbora/Karolina in IT4Innovations in the project DD-21-19.For postprocessing and visualisation, the Gnuplot and Paraview software were used.For verification of the computed velocity field by LBM, another simulation was performed by the steady solution Reynolds Averaged Navier-Stokes (RANS) finite volume method (FVM) -Menter's modification of Wilcox's model -SST kω model [19] -in the multiphysics CFD software Star-CCM+.As the output, the velocity magnitudes were used for comparison with the LBM.

Results
In this chapter the comparison of the results for both method (FVM vs LBM) is presented.It was selected 67 distinct point probes, where the velocity magnitudes were compared.The resulting plots of LBM simulations were post-processed in Paraview.To study the influence of the whole geometry (e.g., wall treatment, mixing of the fluid at the branching area of the trachea or in the mouth and nose cavity), probe points were selected across the domain.The points were located on 6 cross-sectional planes through the domain: one in the mouth inlet, one each in the rare section of the nose and mouth chamber, one in the middle of the trachea and in the left and right bronchus after the first bifurcation.In every plane, a couple of points on two perpendicular lines were selected including the wall vicinity to study the effect of the wall boundary condition.A point set was evaluated at the tracheal axis, too.The positions of the cross sections used for comparison are depicted in Fig. 3.For the more accurate location of the point probes, the spatial directions are denoted by letters showing the given direction.For that purpose we used the coronal plane view of the airways from the back.
According to the applied system of coordinates, we used the left direction abbreviated by L corresponding to the X+ direction, the right direction R to the X-direction, the top side T by Z+, the bottom side by Z-, the front side F by Y-and the back side B by Y+.For an easier understanding of the position of the cross sections, these direc- tions are written down to the velocity fields in the following figures.
In Fig. 4, 5 and 6 we present the averaged velocity magnitudes on the desired cross sections.The point positions are shown on the cross section where the velocity field is depicted, too.The comparisons of each group of point probes are visualised in the charts below the pictures.
In Fig. 4 we can see the velocity field in the middle distance from the inlet tube and the outflow of the mouth and nose chamber as marked in Fig. 3. Due to the inlet Reynolds number of 1190, which corresponds to the laminar regime in the pipe, a symetrical flat profile evolves from the inflow piston profile at the inflow boundary which indicates a typical entrance region of a pipe with a further development of a fully parabolic profile.The velocity in the inlet channel shows a good agreement between both methods with a middle relative error rate of 8.4% with higher error rate near the walls.Moving to the cross sections in the mouth and nose cavity, we can observe slightly higher velocities in the mouth cavity.There are 14 probes determined through the domain in the mouth cavity and 7 in the nose cavity.The velocities with the Lattice Boltzmann method in the nasal cavity are lower than by the FVM method.The middle error rate in the mouth cavity is 15.4% and in the nose 37.4%.
In Fig. 5 the velocity field at a cross section in the middle part of the trachea and on the left and right bronchus after the main bifurcation is displayed.Due to the narrowing of the airways, the velocities are higher up to 8 m s −1 than in the the upper airways.We can notice differences in the right and left bronchus cross section in the order of multiple of 10% relative, especially the velocity field in the left bronchus is predicted higher than from the FVM.Contrary, the velocity field in the tracheal cross section shows a good agreement in all points.Slightly greater difference is noticeable in the point number 5 which is in the wall vicinity.The relative error of all points in this cross section is 14.6 %.
In Fig. 6 we show the results at the axis of the trachea.The points were predetermined approximately in the middle of a cross section where every point is located.Here we can clearly distinguish a growing difference between the simulation values with a movement deeper into the airways towards the bifurcation.The highest discrepancy is in the point number 2 in the relative value of 10.6 % which is already acceptable.However, the relative error of this point sequence is 3.9 % which is showing a good agreement indeed.
Analysing the relative errors in the considered probes, some statements can be postulated.First, the values of the velocities in the upper airways above the bifurcation show a better agreement than under the bifurcation.We suppose some inaccuracy produced by the pressure boundary condition applied at the outlets.These are not aligned parallel  to the Cartesian axis and therefore cause some discrepancies acting on the bulk, too.Second, the flow through the nose cavity in the investigated probes predicted by the LBM indicates lower values than by the RANS FVM and therefore a lower mass flux through this part of geometry is simulated by LBM.Third, the relative errors near the walls are higher than in the core of the flow.This could be caused by the considered boundary wall condition.Eventually, a finer grid could also improve the results and could contribute to a more accurate distribution of the flow into the mouth and nose cavities.The middle relative error of all points in the domain is about 21.4 %.
To discuss the quality of the comparison between both methods, one could look at both approaches in more detail.Large eddy simulation ignores the smallest length scale by using the low pass filtering (LES filter) on both time and length scale.The RANS simulation is based on the Reynolds decomposition into its time average and fluctuating quantities.Therefore, the turbulences which are capable by the RANS approach are of a higher length scale than in the case of LES.Although the time averaging in the LES approach was obtained, some discrepancies could have been produced by the differences of both methods.To improve the confidence of the validation, a LES simulation in FVM should be performed and a comparison with experimental measurements, too, which is planned as a next step of our study.Another possibility could also be the effect of the time averaging period and a longer time for the relaxation of the fluid flow.Longer time averaging period and a longer range-up time could increase the accuracy, too.Despite of several discrepancies, a relatively good agreement between both methods is verified.After some improvements, the simulated velocity field ought to be satisfactory for considering particle deposition simulation.

Conclusion
In this paper, a Lattice Boltzmann simulation of constant breathing of a 5 year-old child with a respiratory mask in the sedentary regime was obtained.The LBM results showed a good agreement with the FVM method, especially in the upper part of the airways.Some discrepancies were found in the lower part.Nevertheless, LBM and OpenLB were proven to be a capable tool to solve the air flow through the respiratory tract and it will be incorporated in the aerosol particles tracking to investigate the deposition.Furthermore, more realistic breathing regimes with inspiratory and expiratory cycles will be performed for 10-month infants.In this case, we expect a more detailed computational mesh to retain smaller-scale flow phenomena.

Fig. 2 .
Fig. 2. Front view and side view of the rescaled STL geometry used in this simulation.

Fig. 3 .
Fig. 3. Position of the investigated cross sections in the geometry matched by black spheres and direction description.

Fig. 4 .
Fig. 4. Comparison of velocity fields at the investigated positions in the inlet and in the mouth and nose cavity.

EPJFig. 5 .
Fig. 5. Comparison of velocity fields at the investigated positions in the tracheal cross section and in the right and left bronchus.

Fig. 6 .
Fig. 6.Comparison of velocity fields at the investigated positions on the axis of trachea.