CFD Simulation of Airflow Dynamics During Cough Based on CT-Scanned Respiratory Airway Geometries

The airflow dynamics observed during a cough process in a CT-scanned respiratory airway model were numerically analyzed using the computational fluid dynamics (CFD) method. The model and methodology were validated by a comparison with published experimental results. The influence of the cough peak flow rate on airflow dynamics and flow distribution was studied. The maximum velocity, wall pressure, and wall shear stress increased linearly as the cough peak flow increased. However, the cough peak flow rate had little influence on the flow distribution of the left and right main bronchi during the cough process. This article focuses on the mathematical and numerical modelling for human cough process in bioengineering.


Introduction
Coughing is an important human defense reflex that protects the respiratory system from infections and improves the clearance of secretions [1][2][3][4].It is usually caused by inflammation and respiratory infections of the respiratory airways or larynx.A normal cough process has four steps.It begins with an inspiration to expand the lung volume.Then, the glottis closes and the pectoral and abdominal muscles contract to generate pressure in the chest.Finally, the glottis opens which releases a sudden burst of gas [5].The cough peak flow rate (CPFR) is used to assess the effectiveness of a cough [6,7].
The airflow dynamics of a cough can provide useful information about the size distributions of the virus droplets, the quantity of virus in the droplets, predictions of virus and drug particle transmission and delivery patterns [8][9][10][11][12][13].Several studies on coughs have been carried out.Zhu used high speed photography technology to visualize the cough process [14].Singh and Mahajan determined the relationships between the cough peak flow rate, the cough expired volume and the peak velocity time which are considered characteristics of coughing and age, height, weight and gender of human subjects [15,16].However, none of these studies involved the airflow dynamics of coughing.Gupta analyzed the flow dynamics and characterization of a cough and proposed a flow rate fitting model based on the gamma-probability-distribution functions [17].However, he did not study the airflow dynamics within the respiratory airways.
Considering the difficulties in simulating the airflow dynamics within the highly asymmetric branching structure of the respiratory airways, computational fluid dynamics (CFD) offers a method of analyzing the airflow based on a physically realistic model.The accuracy of CFD simulations in the respiratory system has been proven through experiments [18][19][20].Along with the development of computerized tomography (CT) and image reconstruction processing software, it has become easy to generate a real respiratory airways model for CFD simulations.Numerous studies based on the CT-scanned lung model and image reconstruction processing software have been carried out to investigate the airflow in respiratory airways.Zhang, Ertbruggen, and other researchers have concentrated on droplet dispersion and particle deposition [21][22][23][24].Rochefort studied the air velocity in the respiratory airways through an experimental study and CFD simulation [20].Paz used the Eulerian wall film model in CFD and a gas-liquid two-phase interaction simulation to simulate the coughing clearance process for the upper airways [25,26].Qi conducted an airflow transient dynamics simulation of the human airway tree [27].Luo and Liu modelled the bifurcating flow in a human lung airway [28,29].
In this paper, we analyzed the airflow dynamics during a cough process based on a CT-scanned respiratory airway model.Additionally, the influence of the cough peak flow rate on airflow dynamics and flow distribution were also analyzed.

Governing Equation
Considering the temperature change was very small in this study, the energy transmission was neglected.Only the continuity equation (Equation ( 1)) and momentum equation (Equation ( 2)) were considered and solved. ∇ where ρ and µ are the density and viscosity of air, respectively.
→ v is the air velocity vector, p is the air pressure and → g is the gravity vector.

Turbulent Model
Due to the high air velocity in narrow airways, there is a turbulent flow phenomenon (Re = ρvd/µ) during a cough process in the airways.The applicability of the k-ω model and the shear stress transport (SST) sub-model to a respiratory airflow simulation has been confirmed by numerous studies [18,30].Therefore, the k-ω model and SST sub-model were selected to simulate the turbulent flow in this study.

Boundary Condition
The air was assumed to be a Newtonian fluid with a density of 1.185 kg/m 3 and a viscosity of 1.81 × 10 5 Pa•s in 20 • C and a standard atmosphere.The mass flow inlet was chosen and the mass flow rate was coded through the Fluent user defined function using the model proposed by Gupta which is presented in Equations ( 3) and (4) [17].
where M is the non-dimensional flow rate, τ is the non-dimensional time, CEV is the cough expired volume and PVT is the peak velocity time, which is the time when the flow rate reaches its peak value CEV(l) = 0.138CPFR(l/s) + 0.2983, for male; CEV(l) = 0.0204CPFR(l/s) − 0.043, for female; PVT(ms) = 1.360CPFR(l/s) + 65.860, for male; PVT(ms) = 3.152CPFR(l/s) + 64.631, for female. ( A cough peak flow rate of 6 L/s was selected to generate the cough flow rate.The airflow variation with time is presented in Figure 1.Four points, (a) 0.02 s, (b) 0.08 s, (c) 0.2 s and (d) 0.5 s were selected to analyze the airflow dynamics.
A cough peak flow rate of 6 L/s was selected to generate the cough flow rate.The airflow variation with time is presented in Figure 1.Four points, (a) 0.02 s, (b) 0.08 s, (c) 0.2 s and (d) 0.5 s were selected to analyze the airflow dynamics.The pressure outlet was set at the mouth exit and the pressure was set to be 100 kPa (absolute pressure).The stationary wall with no slip was set to all wall domains.Gravity was also considered in the simulation.

Numerical Modelling
The CT images used in this study were taken from a 34-year-old Chinese male without any history of respiratory disease.The 3D geometry of the lung airways was reconstructed using imageprocessing software (Mimics) and the generated real lung model is shown in Figure 2. Considering the medical statistical regularity in five-generation airways, we extracted the upper-fifth-order airways from the mouth exit to the bronchi, to conduct the study.The pressure outlet was set at the mouth exit and the pressure was set to be 100 kPa (absolute pressure).The stationary wall with no slip was set to all wall domains.Gravity was also considered in the simulation.

Numerical Modelling
The CT images used in this study were taken from a 34-year-old Chinese male without any history of respiratory disease.The 3D geometry of the lung airways was reconstructed using image-processing software (Mimics) and the generated real lung model is shown in Figure 2. Considering the medical statistical regularity in five-generation airways, we extracted the upper-fifth-order airways from the mouth exit to the bronchi, to conduct the study.
A cough peak flow rate of 6 L/s was selected to generate the cough flow rate.The airflow variation with time is presented in Figure 1.Four points, (a) 0.02 s, (b) 0.08 s, (c) 0.2 s and (d) 0.5 s were selected to analyze the airflow dynamics.The pressure outlet was set at the mouth exit and the pressure was set to be 100 kPa (absolute pressure).The stationary wall with no slip was set to all wall domains.Gravity was also considered in the simulation.

Numerical Modelling
The CT images used in this study were taken from a 34-year-old Chinese male without any history of respiratory disease.The 3D geometry of the lung airways was reconstructed using imageprocessing software (Mimics) and the generated real lung model is shown in Figure 2. Considering the medical statistical regularity in five-generation airways, we extracted the upper-fifth-order airways from the mouth exit to the bronchi, to conduct the study.Three models with an increasing number of mesh which are coarse mesh (190941), basic mesh (444204) and fine mesh (2563271 cells) were generated and analyzed.The air velocities along the y direction on the outlet surface were compared and presented in Figure 3.The results showed a limited difference among three models and suggested that the basic mesh could be used in further analysis.
Symmetry 2018, 10, x FOR PEER REVIEW 4 of 11 Three models with an increasing number of mesh which are coarse mesh (190941), basic mesh (444204) and fine mesh (2563271 cells) were generated and analyzed.The air velocities along the y direction on the outlet surface were compared and presented in Figure 3.The results showed a limited difference among three models and suggested that the basic mesh could be used in further analysis.The transient pressure-based solver and the Semi-Implicit Method (SIMPLE) were adopted for the simulation.The Quadratic Upwind Interpolation (QUICK) scheme was adopted for spatial discretization.The flow rate time step size rather than the numerical time step was set to be 0.0001 s.The hybrid-mesh technique was used and the number of cells, faces and nodes were 444,204, 1031,575 and 167,000, respectively.The minimum, maximum and mean quality values, which represented the mesh quality were 0.21, 0.99 and 0.79, respectively.The maximum skewness was 0.73.The CFD solver Fluent Inc. ANSYS 17.0 was adopted for the simulation.

Validation
In order to validate the model and computational methodology, our simulation results were compared with the results measured by Rochefort using magnetic resonance (MR) gas velocimetry under the same boundary conditions [20].The comparisons of the cross-sectional area and volume flow rate are shown in Table 1.The location of the extracted face, X direction and Y direction are presented in Figure 4a.The velocity comparison between our CFD simulation and Rochefort's results at the right main bronchus is presented in Figure 4b,c.In addition, the CFD simulation results were consistent with the MRI results which validates the model and computational methodology used in this study.

Right Main Bronchus Current CFD Simulation Measurement (Rochefort (2007))
A (cm 2 ) 4.75 5.10 Q (L/s) 0.139 0.140 ± 0.011 The transient pressure-based solver and the Semi-Implicit Method (SIMPLE) were adopted for the simulation.The Quadratic Upwind Interpolation (QUICK) scheme was adopted for spatial discretization.The flow rate time step size rather than the numerical time step was set to be 0.0001 s.The hybrid-mesh technique was used and the number of cells, faces and nodes were 444,204, 1031,575 and 167,000, respectively.The minimum, maximum and mean quality values, which represented the mesh quality were 0.21, 0.99 and 0.79, respectively.The maximum skewness was 0.73.The CFD solver Fluent Inc. ANSYS 17.0 was adopted for the simulation.

Validation
In order to validate the model and computational methodology, our simulation results were compared with the results measured by Rochefort using magnetic resonance (MR) gas velocimetry under the same boundary conditions [20].The comparisons of the cross-sectional area and volume flow rate are shown in Table 1.The location of the extracted face, X direction and Y direction are presented in Figure 4a.The velocity comparison between our CFD simulation and Rochefort's results at the right main bronchus is presented in Figure 4b,c.In addition, the CFD simulation results were consistent with the MRI results which validates the model and computational methodology used in this study.

Air Velocity, Wall Pressure, and Wall Shear Stress of the Cough Process
The airflow velocity, wall pressure, and wall shear stress of the four selected points are presented in Figures 5-7.
Figure 5 shows that the airflow velocity streamlines at four time points during the cough process.The highest velocity reached was 35.8 m/s at 0.08 s.By comparing Figure 5a-d, it can be observed that the streamlined patterns at the bronchi were similar.However, there was an obvious difference at the mouth exit where the airflow was more chaotic due to an increase in the airflow velocity, which may have resulted from the airway's turning structure.

Air Velocity, Wall Pressure, and Wall Shear Stress of the Cough Process
The airflow velocity, wall pressure, and wall shear stress of the four selected points are presented in Figures 5-7.The spatial distribution of wall pressure is shown in Figure 6.The wall pressure decreased gradually from the end of the bronchi to the mouth exit at atmospheric pressure.The highest wall pressure reached was 1.012 × 10 5 Pa (absolute pressure) which appeared at the end of the bronchi.Due to the high airflow velocity, the wall pressure at some places, which have been labeled by a red circle, was below the atmospheric pressure.The spatial distribution of wall pressure is shown in Figure 6.The wall pressure decreased gradually from the end of the bronchi to the mouth exit at atmospheric pressure.The highest wall pressure reached was 1.012 × 10 5 Pa (absolute pressure) which appeared at the end of the bronchi.Due to the high airflow velocity, the wall pressure at some places, which have been labeled by a red circle, was below the atmospheric pressure.The spatial distribution of the wall shear stress is presented in Figure 7.We can see that the wall shear stress corresponded to the airflow velocity during the coughing process.The wall shear stress at the bronchi part was generally higher than that at the main airways because of the narrow structure.Additionally, the maximum was 18.69 Pa at the terminal bronchi.

Local Flow Properties
Four cross sections, captured from the trachea, joint, left main bronchus and right main bronchus, were selected to analyze the local flow properties.They were named f0, f1, f2 and f3, respectively.The velocity magnitudes and vectors at these four sections during a cough process are presented in Figure 8. From Figure 8 we can observe that the velocity magnitude at f0 was higher than that at f2 and f3 because of the airflow confluence.Since the area of f1 is larger than f0, the velocity magnitude at f0 was also higher than that at f1.A second flow was observed at four cross sections, as shown in Figure 8b,c, when the airflow was higher than others.The secondary flow direction was from the inner to the outer of the wall.Figure 5 shows that the airflow velocity streamlines at four time points during the cough process.The highest velocity reached was 35.8 m/s at 0.08 s.By comparing Figure 5a-d, it can be observed that the streamlined patterns at the bronchi were similar.However, there was an obvious difference at the mouth exit where the airflow was more chaotic due to an increase in the airflow velocity, which may have resulted from the airway's turning structure.
The spatial distribution of wall pressure is shown in Figure 6.The wall pressure decreased gradually from the end of the bronchi to the mouth exit at atmospheric pressure.The highest wall pressure reached was 1.012 × 10 5 Pa (absolute pressure) which appeared at the end of the bronchi.
Due to the high airflow velocity, the wall pressure at some places, which have been labeled by a red circle, was below the atmospheric pressure.
The spatial distribution of the wall shear stress is presented in Figure 7.We can see that the wall shear stress corresponded to the airflow velocity during the coughing process.The wall shear stress at the bronchi part was generally higher than that at the main airways because of the narrow structure.Additionally, the maximum was 18.69 Pa at the terminal bronchi.

Local Flow Properties
Four cross sections, captured from the trachea, joint, left main bronchus and right main bronchus, were selected to analyze the local flow properties.They were named f0, f1, f2 and f3, respectively.The velocity magnitudes and vectors at these four sections during a cough process are presented in Figure 8. From Figure 8 we can observe that the velocity magnitude at f0 was higher than that at f2 and f3 because of the airflow confluence.Since the area of f1 is larger than f0, the velocity magnitude at f0 was also higher than that at f1.A second flow was observed at four cross sections, as shown in Figure 8b,c, when the airflow was higher than others.The secondary flow direction was from the inner to the outer of the wall.

Local Flow Properties
Four cross sections, captured from the trachea, joint, left main bronchus and right main bronchus, were selected to analyze the local flow properties.They were named f0, f1, f2 and f3, respectively.The velocity magnitudes and vectors at these four sections during a cough process are presented in Figure 8. From Figure 8 we can observe that the velocity magnitude at f0 was higher than that at f2 and f3 because of the airflow confluence.Since the area of f1 is larger than f0, the velocity magnitude at f0 was also higher than that at f1.A second flow was observed at four cross sections, as shown in Figure 8b,c, when the airflow was higher than others.The secondary flow direction was from the inner to the outer of the wall.

Influence of the Cough Peak Flow Rate on Airflow Dynamics
Eight different cough peak flow rates from 3 to 10 L/s were involved in the simulation study.The relationship between airflow dynamics and the cough peak flow rate is given in Figure 9.We found a positive correlation between the maximum velocity, wall pressure, wall shear stress and cough peak flow rate.As the cough peak flow rate increased, the maximum velocity, wall pressure and wall shear stress increased linearly.The first-order linear fitted curves for the maximum velocity, wall pressure and wall shear stress are presented in     2. From the results, we found that the cough peak flow rate had little influence on the flow distribution of the left and right main bronchi.This indicates that the airflow distribution to the left and right main bronchi is almost the same during a cough process.

Figure 2 .
Figure 2. The 3D geometry of the lung airways used in this study.

Figure 2 .
Figure 2. The 3D geometry of the lung airways used in this study.Figure 2. The 3D geometry of the lung airways used in this study.

Figure 2 .
Figure 2. The 3D geometry of the lung airways used in this study.Figure 2. The 3D geometry of the lung airways used in this study.

Figure 3 .
Figure 3.The comparison of air velocity on the outlet surface of three mesh models.

Figure 3 .
Figure 3.The comparison of air velocity on the outlet surface of three mesh models.

Figure 4 .
Figure 4. Comparison of the velocity profile for the right main bronchus.(a) The locations of the extracted face, X axis, and Y axis.(b) The velocity along the X direction.(c) The velocity along the Y direction.

Figure 4 .
Figure 4. Comparison of the velocity profile for the right main bronchus.(a) The locations of the extracted face, X axis, and Y axis.(b) The velocity along the X direction.(c) The velocity along the Y direction.

Figure 8 .
Figure 8. Local flow properties of the four cross sections at four time points (a) 0.02 s; (b) 0.08 s; (c) 0.2 s; (d) 0.5 s.

Figure 9 .
When the cough peak flow rate reached 10 L/s, the maximum velocity, wall pressure and wall shear stress reached 63.93 m/s, 1.0326 × 10 5 Pa, and 35.91 Pa, respectively.
correlation between the maximum velocity, wall pressure, wall shear stress and cough peak flow rate.As the cough peak flow rate increased, the maximum velocity, wall pressure and wall shear stress increased linearly.The first-order linear fitted curves for the maximum velocity, wall pressure and wall shear stress are presented in Figure9.When the cough peak flow rate reached 10 L/s, the maximum velocity, wall pressure and wall shear stress reached 63.93 m/s, 1.0326 × 10 5 Pa, and 35.91 Pa, respectively.

Figure 9 .
Figure 9.The variation in airflow dynamics with the cough peak flow rate (a) maximum velocity; (b) maximum wall pressure; (c) maximum wall shear stress.

Figure 9 .
Figure 9.The variation in airflow dynamics with the cough peak flow rate (a) maximum velocity; (b) maximum wall pressure; (c) maximum wall shear stress.The flow distributions of the left and right main bronchi under the different cough peak flow rate conditions are shown in Table2.From the results, we found that the cough peak flow rate had little influence on the flow distribution of the left and right main bronchi.This indicates that the airflow distribution to the left and right main bronchi is almost the same during a cough process.
Symmetry 2018, 10, x FOR PEER REVIEW 3 of 11 in a single cough process.The CEV and PVT can be determined by CPFR through regression analyses, as shown in Equation (4).

Table 1 .
Comparisons of cross-sectional area and volume flow rate.

Table 1 .
Comparisons of cross-sectional area and volume flow rate.

Table 2 .
Flow distribution at different cough peak flow rates.