Abstract
The mechanical ventilation of human body is a complex human-computer interaction process. High flow nasal cannula oxygen therapy (HFNC) is a new type of ventilation, which is often measured by lung pressure, respiratory work, and other parameters. The purpose of this paper is to analyse the pressure, flow, and strain rate of upper respiratory tract with different flow and oxygen concentration by using finite element simulation, to guide professionals to adjust the appropriate flow and oxygen concentration parameters of HFNC machine. This paper studies the complex human-computer interaction environment of human respiratory tract and ventilation airflow. The 3D model of respiratory tract established by the conversion of image scanning data was taken as the research object. The flow state of the gases in the respiratory tract was judged by Reynolds equation. After that, RNG K-ε model was applied to the research object, and the simulation diagram of airway pressure, flow rate, strain rate, and trace diagram of flowing particles were obtained under the finite element method. The results explain some clinical phenomena in HFNC and guide people to make better use of mathematical tools to study human-computer complex environment.
1. Introduction
The main clinical manifestations of patients with respiratory diseases are dyspnea, cough, expectoration, and other symptoms, which have serious adverse effects on the physical and mental health and quality of life of patients. At present, many measures and methods for the treatment of respiratory diseases are put forward. Oxygen therapy is the most commonly used treatment of respiratory diseases. Traditional oxygen therapy uses a nasal catheter or mask to inhale oxygen. It can help improve the respiratory function of patients with respiratory diseases, promote the functional indicators of patients to return to normal quickly, and delay and block the development trend of the disease, to improve the prognosis of patients. However, traditional oxygen therapy is difficult to achieve enough humidification degree and temperature, and the oxygen flow rate is limited. High flow nasal cannula oxygen therapy (HFNC) is the latest respiratory support technology, which can meet the needs of airway humidification, proper temperature, and oxygen flow rate at the same time. These advantages make it more comfortable for patients and are rapidly promoted in clinical practice [1]. The oxygen from the electromagnetic flow proportional valve is fully mixed with the air and heated and humidified by flowing through hot water and finally enters the human body with a high flow rate and appropriate temperature and humidity.
HFNC has the following advantages [2]: (1) it can provide more stable and higher inhalation oxygen concentration than the ordinary therapy, which can meet the needs of patients’ spontaneous breathing [3]; (2) high flow airflow can reach or exceed the maximum inspiratory flow of patients’ active inspiration, which can reduce inspiratory resistance and respiratory work [4]; (3) it can warm and humidify the gas to body temperature, which can reduce the heat and water consumption of patients with respiratory distress, and maintain the airway mucociliary function in the best state; (4) high flow airflow washes the dead space of upper airway, reducing anatomical dead space and improving patient ventilation efficiency [5, 6]; (5) high flow airflow provides a certain level of end-expiration positive airway pressure, which can open alveoli and increase lung volume; (6) HFNC does not need a completely closed path and is convenient to eat and communicate [7], while it can increase the comfort of patients [8]. These advantages were also evaluated in vitro under different flow rate, intubation size, and air leakage [9]. HFNC has also been shown to protect against respiratory failure and reintubation after extubation [10].
As a new way of ventilation, it has shown great practicality. Miguel-Montanes pointed out that HFNC significantly improved preoxygenation and reduced the incidence of severe hypoxemia [11]. Messika demonstrated that HFNC is effective for oxygenation in a small number of acute respiratory failure (ARF) through single centre study [12]. Nishimura cites others studies to conclude that HFNC is an effective early treatment for adult respiratory failure [13]. Spoletini suggested that HFNC is not only more comfortable than NIV but also relieves more pain of patients with dyspnea [14]. Manizheh compared the efficacy and safety of HFNC with nasal CPAP for respiratory support in premature infants with respiratory distress syndrome (RDS) [15]. Schlapbach used it in hospital transfer of children with critical disease and found that the treatment rate was significantly higher than that of noninvasive ventilation [16].
The flow field analysis involved in this study is the analysis of nasopharynx and trachea, including nasal cavity, nasopharynx, oropharynx, pharynx, and trachea. These parts are not only the channels of air, but also have the functions of defense, regulating air temperature, and humidity. Due to the complexity of the structure, the simulation results of computational fluid dynamics (CFD) technology have become an important reference for clinical research, which can be used to design and improve the high-flow nasal cannula oxygen therapy. Therefore, ensuring the accuracy of CFD simulation results and making a reasonable analysis of the results became the key to the study.
In this study, the flow field is analysed by FLUENT software, which replaces the physical model to grid model and calculates approximate solution of pressure or velocity in each grid so that the complex problem of the whole flow field state in respiratory tract is turned to a simple grid problem. It is not only highly accurate, but also suitable for the complex shape of respiratory tract. CFD to respiratory tract can not only effectively simulate the air distribution during oxygen therapy but also obtain the velocity and pressure of each point in the trachea. It is of great significance for diagnosing and treating respiratory diseases, adjusting HFNC equipment to the appropriate state accepted by patients, and reducing ventilator injury of patients.
This method is helpful for professionals to understand the airway condition of patients with different flow and oxygen concentration and can provide doctors with certain reference when adjusting the flow and oxygen concentration of HFNC. It can provide sufficient inspiratory flow and oxygen supply level, thus reducing inspiratory resistance and respiratory work, increasing ventilation efficiency, and avoiding airway damage to patients.
2. Methods
The method of flow field analysis is divided into the following steps: establishing geometric model, analysing flow form and flow model, setting material and boundary conditions, setting analysis steps, completing analysis, and postprocessing.
2.1. Geometric Model of Upper Respiratory Tract
In geometric model, two-dimensional (2D) model has simpler modelling, but for complex respiratory system, three-dimensional (3D) model is almost used to describe the structural characteristics of respiratory system in recent years. At present, the following methods are mainly used in the establishment of three-dimensional model: (1) the classic anatomical model is Weibel A model, which is an ideal geometric symmetry model. After that, some scholars put forward the lung model of asymmetric bifurcation region [17] and the asymmetric bifurcation bronchial tree model [18]. (2) The 3D model is established based on the transformation of image scanning data. 3D reconstruction software for computed tomography (CT) or magnetic resonance imaging (MRI) greatly facilitates the generation of 3D models, such as Mimics software. Obviously, the latter modelling method is closer to the real human respiratory tract and has more geometric details. However, due to the influence of noise and image resolution, the important parameters of smaller airway still need to be obtained by anatomical model [19].
The model studied is a 32-year-old female’s upper respiratory tract and trachea. The surface model of respiratory tract wall was established by Mimics software (as shown in Figure 1), and then, SIMENS NX software was used to transform the surface model into a solid model. Then, the solid model of cavity was extracted (as shown in Figure 2) and finally meshed the model in FLUENT software (as shown in Figure 3). In the model of endonasal cavity, the diameter of nasal tube is 5 mm, the diameter of nostril is 10 mm, the axial length of nasal cavity is 73.29 mm, the transverse length of nasal cavity is 30.28 mm, the axial length of throat is 86.14 mm, the axial length of trachea is 109.38 mm, and the diameter of trachea end is 19.5 mm.
During oxygen therapy, the mixed gas flows from the nostril to the lung. The blue section is the air inlet, and the red section is the air outlet. The model has 653723 grids and 194538 nodes.
2.2. Flow Model of Airflow during Ventilation
To obtain the flow form, it is necessary to set the ventilation flow rate higher than the maximum flow rate of nature breathing, which is among 40–60 L/min commonly. In the model, the diameter of nasal tube is assumed to be 5 mm, so the ventilation flow rate is selected as 32 m/s. Simultaneously, there are three kinds of operation conditions: 50 L/min and 60 L/min to describe the conditions of medium flow rate and high flow rate. Therefore, the gas flow rates under the three conditions are 32 m/s, 40 m/s and 48 m/s respectively. The calculation of Reynolds number is as follows:where is the velocity, is the density of air oxygen mixture, is the diameter of gas vent, and is the dynamic viscosity of air oxygen mixture.
The supply gas of the ventilator is a mixture of air and oxygen. In the finite element simulation, the density and dynamic viscosity of the gas mixture should be determined. These two variables are affected by the nature and composition of each molecule in the gas. The proportional valve of the oxygen therapy machine controls the amount of oxygen, and the fan controls the amount of air, which reflects the mixing degree of the gas mixture. The density and dynamic viscosity of the gas mixture can be calculated by the linear relationship of the mixing ratio of air and oxygen. When the mixing ratio of oxygen and air is α (0% < < 100%), the density and dynamic viscosity of mixed gas meet the following requirements:
According to Reynolds number, laminar and turbulent models can be judged. The general fluid motion equation is Navier–Stokes (N–S) equation:
It can be written in rectangular coordinates:where is the density of the fluid, is the pressure, is the velocity vector, , , and are the velocity components of the fluid at point (x, y, and z) at time , is the external force on the fluid per unit volume, and is the dynamic viscosity.
According to the calculation, the Reynolds numbers of 40 L/min, 50 L/min, and 60 L/min at 70% mixing ratio of oxygen and air at the narrowest part in the airway are 11155, 13943.75, and 16732.5 respectively, indicating that the three cases are turbulent flow, and from the physiological point of view, the existence of turbulence can enhance the contact between the airflow and the mucosal boundary layer, and improve the heating and humidifying effects of the air passage on the air.
The solution of N–S equation is very dependent on initial and boundary conditions. For turbulence, the initial condition is to add random turbulence disturbance, which causes the uncertainty and randomness of the obtained turbulent solutions. Therefore, in addition to the general fluid formula, it is necessary to introduce specialized turbulence equations.
In the selection of turbulence model equations, some studies [20] have shown that: in the case of higher gas flow rate, several commonly used turbulence models (LES model, K-ε model, standard k-ω model, SST K-ω model) are compared. Among them, the standard k-ω model has the least difference between the simulation and experimental results. In addition, compared with RNG K-ε model [21], the traditional K-ε model is a high Reynolds number model, but it takes the whole upper respiratory tract as turbulence, which will lead to poor performance. RNG K-ε model [22] is a low Reynolds number turbulence model, which considers the laminar flow and turbulent flow in the simulation process. The RNG model adds a condition to the ε equation, which not only effectively improves the accuracy, but also takes the turbulent vortex into account, which has a high credibility in the numerical simulation of airflow in the respiratory tract. Jeong et al. showed that the RNG K-εe model was more consistent with the experimental values than the traditional K-ε model when simulating the upper respiratory tract of a patient. Therefore, RNG K-εe model equations are used in this study [23].
The RNG K-εe model is shown as follows. K equation is equation of turbulent kinetic energy and ε equation is equation of turbulent dissipation rate:
Among them, , , is the generation of turbulent kinetic energy caused by mean velocity gradient, is the turbulent kinetic energy generated by buoyancy, is the contribution of wave expansion to the total dissipation rate in compressible turbulence, and are the reciprocal of the effective Prandtl numbers of K and ε, and and are the source term. In this study, buoyancy and compressibility of flow are not considered.
The differential equation of turbulent viscosity is as follows:
In (9), the calculation of some parameters is shown as follows:
Expression of turbulent viscosity is as follows:where .
The statistical method of renormalization group theory in RNG theory provides turbulent Prandtl numbers and for formulas 7 and (8), which are set as fixed values in standard K-ε method. Reciprocals of effective Prandtl numbers and are derived from RNG theory:where . In this paper, the Prandtl numbers have little influence on the flow field, but they increase the calculation work, so .
is an additional term in the equation, which is expressed as follows:where and .
3. Materials and Boundary Conditions
The material studied is a mixture of air and oxygen. The mixture density is calculated according to the formula with the oxygen and air mixing ratio mentioned previously. In this study, we analysed the oxygen concentration varying from 21% to 100%, and the oxygen and air mixing ratio is 0% to 100%. The density and dynamic viscosity of the gas mixture can be calculated by these data. In this study, the reduction of total density of mixed gas caused by the introduction of water vapor is not considered because of the very little pressure of saturated vapor pressure of water in the gas from HFNC equipment.
There are four flow surfaces in the nasal cavity and a flow surface at the end of the trachea. These five surfaces are the boundary interfaces. In most cases of inhalation, the nasal cavity does not get additional air from the outside but relies on HFNC ventilation. The mechanical ventilation surfaces are also the flow velocity input surface, and the surfaces from nostril and trachea are set as the pressure output surface and the pressure value is atmospheric pressure. On expiration, the space where the nostrils are not intubated is set as the pressure output surface, and the pressure value is atmospheric pressure. The surfaces of the HFNC ventilation surface and trachea are the flow rate input surface.
3.1. Analysis Steps Setting
The transient solver based on pressure is used in the finite element calculation. The initialization and calculation are started from the inlet. The SIMPLE method is used to solve the coupling of pressure and velocity. The time step is set as 0.1 s, while the number of time steps is 20 and the maximum number of iterations is 30. After calculation, the fluid movement time is 2 s, which is enough for single exhalation or single inhalation.
In the follow-up analysis, simulations are first taken out to compare the different effects of inhalation and exhalation on the airway under the same conditions, and then, the pressure, velocity and path lines are compared to analyse the difference of airway flushing effect under different flow rates and the influence of different oxygen concentration on airway pressure and flow rate.
4. Results
4.1. Comparison of Flow Field in Different Breathing States
Through the comparison of inspiratory and expiratory parts, it can be seen that the parameters of flow field are quite different on inhaling and exhaling. Figure 4 shows the pressure contrast isoline of the airway at 70% oxygen ratio and 40 L/min flow rate. It can be seen that the pressure of 1.5 cm H2O can be reached on inhaling, and the minimum pressure is −1 cm H2O. The highest pressure appears in the nasal cavity, especially at the top of the nasal cavity and the part near nasal concha, while the lowest pressure appears in the trachea. Generally speaking, the pressure from the nasal cavity to the trachea gradually decreases.
During exhalation, the highest pressure (about 2.7 cm H2O) is produced at the end of the trachea, while the lowest pressure (about 1.7 cm H2O) is produced at the place, where the nostrils connect with the outside world. However, the maximum pressure of expiration is still higher than that of inspiration. Simultaneously, the obvious irregular places, oropharynx and pharynx, have no pressure mutation. From the perspective of respiratory work, the pressure difference between the nasal cavity and trachea as well as the unchanged volume of the airway during inhalation provides part of the breathing work, reducing the extra work of the respiratory muscles of the lung and increasing patients’ comfort.
Compared with the pressure, the flow rate shows a significant difference (shown as Figure 5). The maximum velocity (about 10 m/s) appears at the place near nostril. Except the velocity near the gas inlet and junction of larynx and trachea, velocity of all parts presents a low velocity state, most of which are between 0 and 1.5 m/s. During exhalation, the maximum flow rate also appears at gas vent (over 10 m/s). There is also a concentration of the flow velocity at the junction of the larynx and the trachea, and there is a larger velocity area in the trachea. That is because the gas coming out of the lung meets the narrow larynx when exhaling, which causes the velocity increase of this part. The excessive air oxygen mixture gas flows to the atmosphere through the space where the nostrils are not intubated, resulting in the high velocity of the place.
The strain rate can be multiplied by the viscosity to obtain the shear stress, which represents the force on airway wall. The strain rate is larger in the front of nasal cavity, oropharynx, and upper part of trachea, and the maximum value (about 13800) occurs at the junction of larynx and trachea, which is also prone to shear deformation. The maximum strain rate at the nasal cavity is 15000/s, which is generated in the front of the nasal cavity. The maximum stress at the larynx is 13500/s in the transition narrow area between the nasopharynx and oropharynx. The maximum strain rate of the trachea is 13000/s at the upper end of the trachea. In general, both inspiratory and expiratory strains in the front of the nasal cavity and the larynx should be paid more attention. Figure 6 also shows the strain rate of normal breathing on airway. It is found that the strain rate of HFNC at 40 L/min is increased compared with that of normal breathing. Although the maximum strain rate increase is not obvious, it is obvious that the range of larger strain region widens.
As shown in Figure 7, the left figure shows the path lines measured by particle ID during inhalation. In the process of inhalation, high flow gas is injected from nasal tube, and there are obvious gyrations in the vent and inferior nasal meatus, which is more due to the complexity of the air gap in the nasal cavity. Apart from the obvious turbulence at the junction of trachea and larynx, laminar flow features are more obvious in the trachea. The right figure shows the path lines during exhalation. In the process of exhalation, the gas from the trachea affects the airflow state of HFNC, making oxygen and air mixture unable to reach the nasopharynx. In addition, there is an obvious swirling turbulence in the nasal cavity, which finally flows out through the nostril. The gas also presents a similar laminar flow state in the trachea and an obvious turbulence in the nasal cavity. On exhaling, the oxygen and air mixture from the nasal tube forms turbulence in the nasal cavity, which then flows to the atmosphere from the nostril gap. Therefore, a certain amount of oxygen is maintained in the turbulent part, which dilutes carbon dioxide from the lungs and reduces the amount of carbon dioxide brought into the original nasal cavity during the next inhalation to achieve the purpose of flushing the dead zone of nasal cavity.
Figure 8 shows the time map of particles in the airway. When exhaling, the gas from the trachea conflicts with that from HFNC equipment, leading to a large convection in the front of the nasal cavity.
4.2. Comparison of Flow Rate and Oxygen Concentration
Figure 9 shows the maximum pressure comparison in different ventilation flow rates. Under the flow rate of 40 L/min, the pressure of a large area of upper respiratory tract is between 0.5 and 0.8 cm H2O, and the smaller part reaches a pressure of 1.7 cm H2O. The pressure of trachea (about −0.5 to 0.2 cm H2O) is lower than that of upper respiratory tract. At 50 L/min and 60 L/min, the pressure of upper respiratory tract is also higher than that of the lower respiratory tract, whose maximum pressures are 2.0 cm H2O and 2.3 cm H2O and minimum pressures are −0.4 cm H2O and −0.1 cm H2O, respectively. In addition, with the increase of ventilation flow, the maximum pressure and minimum pressure increased in varying degrees. This airway pressure difference increases with the increase of flow rate, which indicates that high flow rate gas can reduce the respiratory work required by respiratory muscles and increase the comfort of patients.
By comparing the airway velocity under different flow rates, the maximum velocity is generated at the contact part between the nasal cavity and the nasal tube, and the velocity is far lower than the ventilation velocity at the position with larger cross-sectional area. At same positions, the change trend of the flow rate is the same. In addition, the second largest velocity is generated at the junction of the larynx and trachea, and the minimum is in the middle and rear part of nasal cavity. Figure 10 shows comparison of maximum velocity at different ventilation flow rates.
According to the pressure comparison of respiratory airway under different oxygen concentrations, it is found that the pressure of nasal cavity is greater than that of throat and trachea. Simultaneously, with the increase of oxygen concentration, the maximum pressure of each part decreases to varying degrees, which can be seen in Figure 11.
On choosing the median plane as the analysis plane to compare the oxygen concentration (as shown in Figure 12), it can be seen that the maximum velocity decreases with the increase of oxygen concentration and the area of high velocity area is also reduced, while the minimum flow rate is basically unchanged. This is because the viscosity of the gas mixture increases with the increase of oxygen concentration, and the friction force of intermolecular motion is large. When the initial velocity is the same, the instantaneous velocity at the same position will decrease with the increase of viscosity. The maximum velocity is produced at the junction of larynx and trachea, and there is no obvious change in other areas, indicating that the difference would be more obvious when the flow diameter is reduced.
5. Discussion
Based on the finite element analysis method, the cavity in the three-dimensional model of respiratory tract established by image scanning data conversion is taken as the research object. It is found that when inhaling, the greatest pressure is concentrated in the nasal cavity, and the highest velocity is concentrated in the vent and junction of the larynx and trachea. The most easily damaged parts by shear stress are the front of the nasal cavity, near the vent, and junction of larynx and trachea. When exhaling, the greatest pressure is concentrated in the trachea, and the highest velocity is concentrated in the nasal vent. The most easily damaged parts by shear stress are the front of the trachea and the upper end of the trachea. In addition, the particle trace states and characteristics of different inlets are obtained by comparison. Then, different flow rates are compared, and it is concluded that even if the change of flow rate is not obvious, the increase of flow rate also leads to the increase of pressure. On comparing different oxygen concentrations, the flow rate decreases with the increase of the oxygen concentration.
In this study there are some assumptions. First, the flow of gas meets the simplified setting. The simplified design is that the fluid is in one-dimensional isentropic flow, and the physical properties are uniform in the inlet and outlet flow cross-section, the flow does not shrink at the nozzle, the gravity is ignored, and the specific heat of the fluid is constant. Second, the fluid-structure interaction between the gas and airway wall is not considered. Without considering the fluid-structure interaction, the elastic deformation caused by gas filling would not be considered, and the influence of actual elastic deformation on the airway would not be considered. The nasal cavity of human body is composed of bone and cartilage and the mucosa and skin covered on its surface. Some scholars regard the soft tissue around the upper respiratory tract as a linear elastic structure and simulate the airflow [24], which is widely used as a preliminary study on the interaction calculation of air flow and soft tissue [25]. However, this study focuses on the flow field state of the inner cavity, fully considering the internal velocity, pressure, strain rate, and path line. This paper also focuses on the possible pressure concentration, strain concentration position, internal pressure range, velocity range, and so on. Finally, the patient will not breathe through the mouth during ventilation.
In the acquisition of research objects, this study used Mimics software to convert the CT scanning data of airway into three-dimensional model, so that the position of ventilation port, nasal cavity, larynx and trachea can be clearly seen in the model. However, the obtained model is STL file, which is suitable for 3D printing directly but cannot be edited, so that it is impossible to extract the lumen directly. Therefore, it is necessary to convert it into a 3D model to the surface solid. SolidWorks software is used and the ideal effect is obtained. After that, the editable solid is obtained in the SIMENS NX software using the commands of combined surface, suture surface, and building solid. Finally, in the geometry module of ANSYS Workbench, the surfaces of inlet and outlet are generated, and the filling command is executed to generate the fluid part, which is the research object.
In terms of conclusion analysis, this study not only uses the overall contrast figure but also uses the longitudinal sectional contrast figure. However, the transverse sectional comparison is not selected because of its less distinctiveness. For example, when studying the velocity comparison of different oxygen concentrations, the maximum velocity at 0.01 m, 0.03 m, and 0.06 m away from the front end of nasal cavity changes slightly, so the middle longitudinal section is selected as the contrast section, which not only shows the whole character of the respiratory tract studied but also reflects the position of the maximum pressure.
The shortcomings of this study are as follows.(1)The fluid-structure interaction of the oxygen air mixture flow on the respiratory tract is not considered, so the research can only focus on the flow field analysis but does not consider the shear damage of the gas flow on the trachea wall.(2)This study is based on the finite element simulation. The mesh fineness can affect the simulation effect. It would be more visual and accurate if more details in the trachea could be analysed using a more refined mesh in the future.(3)The RNG K-ε model used in this paper has some advanced significance, but the complete flow field analysis needs to be verified by experiments. The follow-up study will be based on solving these shortcomings. It is hoped to obtain more accurate and realistic respiratory airway flow field under HFNC.
6. Conclusion
Based on the finite element analysis method, the three-dimensional model of trachea established by the conversion of image scanning data is taken as the research object, and the conclusion that the flow state presents turbulent state under different breathing states, flow rates, and oxygen concentrations is obtained. After that, RNG K-ε model is used to solve the turbulence state of each position. The pressure, flow velocity, and strain rate simulation diagram of the airway under the finite element method are obtained.(1)During inhalation, the greater pressure is concentrated in the nasal cavity, and the high flow velocity is concentrated at the junction of larynx and trachea. The parts easily damaged by shear stress are the front of nasal cavity, near the ventilation port, and junction of larynx and trachea. During exhalation, the greater pressure is concentrated in the trachea, the high flow rate is concentrated in the nasal vent and the front end of the trachea, and the severe shear stress is in the upper part of the trachea.(2)The maximum movement time of particles in inhalation is less than that in exhalation. Oxygenated gases can reduce the dead space by increasing the oxygen content in nasal cavity.(3)Flow rate and oxygen concentration affect the respiratory-related parameters. A large flow of gas can produce a large pressure difference, which can reduce inspiratory resistance and respiratory work. With the increase of oxygen concentration, the maximum flow rate decreases.
In future research, the shear damage of the air flow in the trachea wall will be considered, and a more precise grid will be used to analyse more details in the trachea intuitively and accurately.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was funded by China Postdoctoral Science Foundation and the Open Research Project of the State Key Laboratory of Media Convergence and Communication (no. 2019M660392) and Communication University of China (no. SKLMCC2020KF002).