A Numerical Simulation of Air Flow in the Human Respiratory System Based on Lung Model

Md. Kamrul Hasan^1*, Mahtab U. Ahmmed², Md. Samsul Arefin²

¹Department of Mathematics, Khulna University of Engineering and Technology, Khulna, Bangladesh.
²Department of Mathematics, Jahangirnagar University, Dhaka, Bangladesh.

DOI: 10.4236/jamp.2023.118142                

Abstract

The lung is an important organ that takes part in the gas exchange process. In the study of gas transport and exchange in the human respiratory system, the complicated process of advection and diffusion (AD) in airways of human lungs is considered. The basis of a lumped parameter model or a transport equation is modeled during the inspiration process, when oxygen enters into the human lung channel. The quantitative measurements of oxygen are detached and the model equation is solved numerically by explicit finite difference schemes. Numerical simulations were made for natural breathing conditions or normal breathing conditions. The respiratory flow results for the resting conditions are found strongly dependent on the AD effect with some contribution of the unsteadiness effect. The contour of the flow rate region is labeled and AD effects are compared with the variation of small intervals of time for a constant velocity when breathing is interrupted for a negligible moment.

Keywords

Lumped Model, Lumped Model Channel, Mass Flow Rate, Ideal Law of Gas, 2D Advection Diffusion Equation, Finite Difference Scheme

1. Introduction

A complete cycle of human respiratory system consists of two processes, one is inspiration and other is expiration and the respiration in the human lung channel is a consequence of oscillatory flow (advection) and stagnation in transition (diffusion) in nature. During inspiration while we are in rest, the inhaled oxygen gas has a partial pressure about and temperature about after reaching trachea in the flow system [1] . Branching of breathing airways of the human lung was introduced by Weibel (1963) which leads to an idea of fluid flow system in our lung model. In the fluid flow system, the motion is governed by flow rate with three basic parameters: inductance, resistance and capacitance [2] . The advection diffusion of inhaled gas has actually been observed in the 18th generation (G18) channel. Along a circular channel laminar, turbulent and oscillatory dispersion is investigated by H.K. Chang [3] . To study airflow distribution a lumped parameter model has been developed by Elad et al. [4] . The transportation of fluid flow problem which describes advection diffusion process can be modeled using a two dimensional advection diffusion equation. Heat, mass, velocity, energy, concentration profile etc. can be elaborated by ADE [5] . Numerical techniques are used in solving ADE, finite difference method (FDM) is one of them. The groundwater flow problem is solved using the Crank-Nicolson finite difference scheme in second order time and space by Lingyu Li, Zhe Yin [6] . A mathematical model is developed in the form of two dimensional advection diffusion equation for calcium profile and solved analytically by B. K. Jha et al. [7] . A. Fedi et al. [8] analytically solved two dimensional advection diffusion equation in porous media where a laterally bounded domain is considered. In many engineering problems such as pollutant transport in rivers and streams, thermal pollution in river systems, the ADE is used as a model equation [9] . The stability analysis for different finite difference schemes is discussed by ( [10] [11] [12] [13] [14] ) for the two dimensional ADE. Most of the work as stated above has been performed for open channels. But along closed channels the fluid (liquid or gas) may flow downstream by advection or spread out by diffusion. The two-dimensional modified equation method has been successfully used with weighted discretization to develop several new explicit methods and explain the effect of two-dimensional advection-diffusion equation. In particular it gives an idea of the wave propagation characteristics of the methods [15] . High frequency oscillatory ventilation (HFOV) is a type of mechanical ventilation that helps in breathing of lung affected people. HFOV controls the constant oscillatory flow along human lung channels and during HFOV the influence of asymmetric airway compliances on redistribution of gas particles has been investigated numerically by Hirahara et al. [16] . The airways network has quite small dimension and the smaller airways deep down into the lungs are inaccessible. A computational fluid dynamics (CFD) modelling approach is used to study the unsteady respiratory airflow dynamics within a human lung. Numerical simulation was made for two breathing conditions: 1) resting or normal breathing condition and 2) maximal exercise condition. The respiratory flow results for the both conditions are found strongly dependent on the convective effect and the viscous effect with some contribution of the unsteadiness effect R.K. Calay et al. [17] . Their study was to use the proposed computational model to investigate the dynamic ability of the nasal cavity to heat and moisten the inhaled air in the studied area, similar to the nose so that it is possible to comprehensively study the structural components ( [18] [19] ).

In this paper we propose a transport equation of lung model channel involving flow rate (mass flow rate) on a quantitative analysis of inhaled O₂ gas. The transport equation is solved numerically. The stability condition of the scheme is studied and flow rate region is presented by contour line. Comparison of advection diffusion effect for different radial and longitudinal position with same initial and boundary conditions are also shown graphically. The goal to understanding the airflow mechanisms can be approached using numerical modelling of the respiratory airflow network.

2. Respiratory Branches and Function of Lungs

The lung is an important organ that takes part in the gas exchange process. In the pulmonary hyperinflation process, the muscles of the thoracic cage extend the lung tissues to increase the volume and then air is breathed into it. In the artificial ventilation process air is sucked into the human lung by increased outside pressure. Human lung consists of a structure like a branching tree to distribute the breathing air in different parts of consecutive generations of channels. The area of each cross section of respiratory bronchiole is decreased in the next generation and for that the air flow is dropped off to reach the lowest velocity in last generations.

In the generations (G17-G19) of respiratory bronchiole there are the bronchioles are small air passages found at the end of the bronchiole tree in the lungs. These bronchioles are responsible for conducting air from the bronchi to the alveoli of the lungs, where gas exchange takes place. The respiratory bronchioles also have a role in oxygen uptake and carbon dioxide elimination. The walls of the respiratory bronchioles are lined with small clusters of air sacs, known as alveoli, which are responsible for gas exchange with the blood vessels. The thin walls of the alveoli allow for rapid gas exchange, enabling the respiratory bronchioles to efficiently oxygenate the blood and remove carbon dioxide. Apart from its role in gas exchange, the respiratory bronchioles also help to humidify and filter the air that enters the lungs. They do this by having specialized epithelial cells that secrete mucus, which helps to trap foreign particles and prevent them from entering the lungs. Overall, the respiratory bronchioles are an essential component of the respiratory system, playing a vital role in ensuring that oxygen is delivered to the body’s tissues and carbon dioxide is efficiently eliminated from the body. According to respirational function and its anatomical configuration, the lung consists of two regions called air transport and gas diffusion. The airways from trachea to terminal bronchioles are fractionalized repeatedly and gas is delivered without any gas exchange.

According to the measurement of lung tree and its anatomy, the human airway of lung is shown in Figure 1. It is drawn by performing the calculation of Weibel’s real data analysis for adult lung.

Figure 1. Pattern of airway branching (Weibel, 1963).

3. The Mathematical Model

To simulate the air flow in the G18 channel, the basic equations for the conservation of mass, lumped parameter model, equation of mass flow rate, ideal gas law, temperature and density of oxygen are determined as follows:

$\frac{\partial q}{\partial t}+\left(\nabla \cdot q\right)q=-\frac{1}{\rho }\nabla p+\vartheta {\nabla }^{2}q$ (1)

In absence of driving force for incompressible fluid flow along G18 channel a lumped parameter model can be derived as:

$L\frac{\text{d}q}{\text{d}t}+Rq+\frac{1}{C}{\displaystyle \int q\text{d}t}=0$ (2)

where,

Inductance: $L=8.75\times {10}^3\text{Pa}\cdot {\text{s}}^2\cdot {\text{m}}^{-3}$ .

Resistance, $R=1.44\times {10}^7\text{Pa}\cdot \text{s}\cdot {\text{m}}^{-3}$ .

Compliance, $C=4.44\times {10}^{-12}{\text{m}}^3\cdot {\text{Pa}}^{-1}$ and $q=q\left(t\right)$ be the flow rate (mass flow rate) of the fluid.

$\frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=\frac{L{u}^2}{R}\left(\frac{{\partial }^{2}q}{\partial x^2}+\frac{{\partial }^{2}q}{\partial y^2}\right)$ (3)

$\frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=D\left(\frac{{\partial }^{2}q}{\partial x^2}+\frac{{\partial }^{2}q}{\partial y^2}\right)$ (4)

This is a two dimensional advection diffusion equation (ADE) for lung model channel. Where, D is diffusion coefficient, $q\left(x,y,t\right)$ is the flux or velocity vector, t is the real time, T is the temperature, $\rho ,\vartheta$ and u are the density, kinematic viscosity and convective term respectively. For numerical study, second-order spatial discretization schemes for length and diameter were used. Thus, the real problem will be solved in a discrete form by applying the finite difference method.

The mass flow rate within the G18 bronchial tube by the formula:

$q=\rho \times A\times v$ .

where,

$\rho$ : Density of the oxygen at 37˚C.

v: Velocity of the fluid in the tube.

A: Cross sectional area of the tube.

The ideal gas law is written as:

$PV=nRT$ (5)

where P, V and T are the pressure, volume and absolute temperature of gas; n is the number of moles of gas; and $R=0.0821\text{Latm}\cdot {\text{K}}^{-1}\cdot {\text{mol}}^{-1}$ is the ideal gas constants. Also $n=m/M$ where m is the weight of gas; and M is the molecular weight of gas [20] . Inserting $n=m/M$ into Equation (5) we have a formula for calculating density and which is given by

$\rho =\frac{PM}{RT}$ (6)

Since the molecular weight of oxygen O₂ is 32 gm then at 37˚C temperature and 13 kPa (0.128) pressure the density of the O₂ gas is =1.6 × 10⁻¹⁰ Kgm·m⁻³.

However, the 18th generation tube is very short in length, about 1.2 mm and the diameter is very narrow, about 0.5 mm [Weibel’s Table 1]. Then the cross-sectional area of the tube is $A=\pi \times r^2$ , where r is the radius of the tube. Thus the mass flow rate within the tube for velocity 20 mm·s⁻¹ is $q=\rho \times \pi \times r^2\times v=6.3\times {10}^{-10}\text{Kg}/\text{s}$ .

4. Numerical Technique for Model Equation

Without driving force and compliance effect for the unsteady incompressible flow along a rigid channel is an two dimensional ADE for physical domain of channel length $0\le x\le l$ and diameter $0\le y\le d$ (creating the surface that we get by cut off the channel along its length about the center of the channel), is

$\begin{array}{l}{q}_t+u{q}_x=D\left(q_{xx}+q_{yy}\right)\\ q\left(x,y,0\right)=q_0(x,y)=A,\text{}0\le x\le l,0\le y\le d\\ q{\left(x,y,t\right)}_{|x=l}=B,\text{}0\le y\le d,0<t\le T\\ q{\left(x,y,t\right)}_{|y=0}=C,\text{}0\le x\le l,0<t\le T\\ q{\left(x,y,t\right)}_{|y=d}=D,\text{}0\le x\le l,0<t\le T\end{array}\}$ (7)

To develop numerical scheme by finite difference method (FDM), we discretize xy-plane with mesh size $\Delta x\times \Delta y$ . The spatial and temporal coordinate at the grid point $q\left(t^n,x_i,y_j\right)$ is defined as:

$\begin{array}{l}{x}_i=x_0+i\Delta x;i=0,1,2,\cdots ,M\\ y_j=y_0+j\Delta y;j=0,1,2,\cdots ,N\\ t^n=t_0+n\Delta t;t=0,1,2,\cdots ,K\end{array}$

Then the approximate solution at grid points $q\left(t^n,x_i,y_j\right)$ is $q_{i,j}^n\in R^n$ so that $q_{i,j}^n\approx q\left(t^n,x_i,y_j\right)$ .

The forward time difference formula

$\frac{\partial q}{\partial t}\approx \frac{q_{i,j}^{n+1}-q_{i,j}^n}{\Delta t}$ (8)

The forward space difference formula

$\frac{\partial q}{\partial x}\approx \frac{q_{i+1,j}^n-q_{i,j}^n}{\Delta x}$ (9)

The second order centered space difference formula for second order partial derivative

$\frac{{\partial }^{2}q}{\partial x^2}\approx \frac{q_{i+1,j}^n-2q_{i,j}^n+q_{i-1,j}^n}{{\left(\Delta x\right)}^2}$ (10)

$\frac{{\partial }^{2}q}{\partial y^2}\approx \frac{q_{i,j+1}^n-2q_{i,j}^n+q_{i,j-1}^n}{{\left(\Delta y\right)}^2}$ (11)

Substituting (8), (9), (10) and (11) into (7) and then rearranging to the time level, we get

$q_{i,j}^{n+1}=\left(1+\alpha -2\beta -2\gamma \right)q_{i,j}^n+\left(\beta -\alpha \right)q_{i+1,j}^n+\beta q_{i-1,j}^n+\gamma q_{i,j+1}^n+\gamma q_{i,j-1}^n$ (12)

where,

$\alpha =\frac{u\Delta t}{\Delta x},\beta =\frac{D\Delta t}{{\left(\Delta x\right)}^2},\gamma =\frac{D\Delta t}{{\left(\Delta y\right)}^2}$

We perform the stability analysis as a convex combination technique and also implement the schemes for various parameters. The stability condition is $2\beta +2\gamma -1\le \alpha \le \beta$ , $0\le \beta \le 1$ and $0\le \gamma \le 1$ .

5. Problem Description

An incompressible fluid flow along human lung model channel (G18 channel of Table 1) of length $l=1.2\text{mm}$ , diameter $d=0.5\text{mm}$ for a constant velocity $u=20\text{mm}/\text{s}$ . The diffusion rate is $D=\frac{L{u}^2}{R}$ where u is the constant velocity

and L and R are the lab experiment values. The computed time step and grid width are also taken as our convenience. We divide the domain 150 times for length (x-axis) and 100 times for diameter (y-axis), that is (Δx × Δy) = (150 × 100). Our assumption is that the flow rate is initially the same within the G18 channel and if the velocity is zero (breathing interrupted) just at the start of the G19 channel then what actually happened to the flow rate due to inhalation. The computed time step and spatial grid width are $\Delta t=5\times {10}^{-6}$ , $\Delta x=0.008$ and $\Delta y=0.005$ .

Table 1. Approximation quantification of the human bronchial system (Weibel’s model).

6. Numerical Solution

The air flow in the G18 channel of human lung of a person plays an important role in many physiological functions of the respiratory system, such as temperature and moisturizing the flow of air and others. In this section, the proposed model is used to predict air flow and related transport phenomena in human lung. Normal respiration was chosen as a reference base, and then the effect of changes in mass flow rate was investigated for variation of velocity. This study serves as the basis for a better understanding of transport phenomena in the G18 channel, which is the function of the lung. The studied area was identical to the second test problem of Figure 1. The temperature of the oxygen gas is 37˚C, diameter and length is taken from (Weibel’s Model) Table 1. From the obtained results, it can be noted that the global behavior of the air flow was not changed, but, however, the maximum longitudinal velocity increased.

The diffusion effect for different values of time is shown in Figure 2. It is observed that flow is higher for higher value of time and diffusivity. Different values of effective time (0.04, 0.05, and 0.06) are taken and the contour in Figure 2 shows the flow phenomena. In Figure 2, we found the diffusion effect along the diameter of the 18^th generation bronchial tube for different time, where the labeled maximum flow rate region is about 6.0 × 10⁻¹⁰ Kg·s⁻¹ for the oxygen gas at 0.06 seconds after interruption of breathing. When breathing time increased more than 0.06sec the flow rate increased.

In Figure 3, we present a diffusion profile with different values of times for a fixed velocity along the diameter and length. We observe that the AD effect for 0.04, 0.05, 0.06 seconds are almost same up to half of the length of the channel and also observe the diffusion; as the time goes, the flow rate is diffusing more.

The diffusion effect along the diameter and length of the G18 bronchial tube are shown in Figure 4. We observe that flow rate are increasing as the diffusion increase within the constraint channel about 6.3 × 10⁻¹⁰ Kg·s⁻¹, 6.6 × 10⁻¹⁰ Kg·s⁻¹, 6.9 × 10⁻¹⁰ Kg·s⁻¹ for the velocities 20 mm·s⁻¹, 21 mm·s⁻¹ and 22 mm·s⁻¹ respectively, at the start of the channel and then advection diffusion occurred parallay up to the end of the channel according to velocity variation at 0.06 seconds.

Figure 2. Contour plot of the flow rate region for different values of time: (a) t = 0.04, (b) t = 0.05 and (c) t = 0.06 seconds where, u = 20 mm/s.

Figure 3. Diffusion effect of flow rate along diameter and length for different values of time with fixed velocity is 20 mm/s.

Figure 4. Diffusion effect of flow rate along diameter and length for different values of velocity with fixed time 0.06 sec.

As the pressure gradient is increasing, the flow rate is gradually increasing and developing to parabolic. It can be seen that no flow occurs if the pressure gradient is zero. In Figure 2 it is noted that other velocity is kept fixed. Figure 4 shows the variations of flow rate for different values of diffusion coefficient. The diffusion indicates how much fluid can flow through the G18 bronchial tube of the lung channel. It can be seen from Figure 4 that when the diffusion is increased the flow rate is also increased in the channel. As a result the flow of oxygen will be increased. So it can be said that the flow rate is an increasing function of diffusion for variation of velocity.

7. Conclusion

This study is a comprehensive work on modeling the processes of oxygen transfer in a 2D ADE model of the human lung. The purpose of this study was to use the proposed computational model to investigate the dynamic ability of the human lung to inhale air in the studied area. Transport functions depend on accurate predictions of the airflow nature. Our model channel is very narrow and short in length. The AD effect is investigated with the variation of time and velocity at a fixed velocity and a fixed time respectively and shown graphically for different radial and longitudinal positions respectively. The maximum flow rate region is also labeled which is shaped as a nice parabolic pattern. With the variation of a small interval of time, the AD effect remains the same almost half of the length of the tube, and then diffuses as time increases. Along the diameter of the tube the flow rate diffuses more with the increase of time. Therefore, the numerical model is tested in an analytically accurate computational model by comparing simulated velocity profiles with values from the measurement.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Gerard, T.J. and Nicholas, A.P. (1987) Principles of Anatomy and Physiology. Fifth Edition, Harper & Row, New York.
[2]	Zamir, M. (2005) The Physics of Coronary Blood Flow. Springer, Berlin.
[3]	Chang, H.K. (1984) Mechanisms of Gas Transport during Ventilation by High Frequency Oscillation. Journal of Applied Physiology, 56, 553-563. https://doi.org/10.1152/jappl.1984.56.3.553
[4]	Elad, D., Shochat, A. and Shiner, R.J. (1998) Computational Model of Oscillatory Airflow in a Bronchial Bifurcation. Respiration Physiology, 112, 95-111. https://doi.org/10.1016/S0034-5687(98)00005-X
[5]	Noye, B.J. (1990) Numerical Solution of Partial Differentials Equations. Lecture Notes.
[6]	Li, L.Y. and Yin, Z. (2017) Numerical Simulation of Groundwater Population Problems Based on Convection Diffusion Equation. American Journal of Computational Mathematics, 7, 350-370. https://doi.org/10.4236/ajcm.2017.73025
[7]	Jha, B.K., Adlakha, N. and Mehta, M.N. (2012) Analytic Solution of Two Dimensional Advection Diffusion Equation Arising in Cytosolic Calcium Concentration Distribution. International Mathematical Forum, 7, 135-144.
[8]	Fedi, A., Massabo, M., Paladino, O. and Cianci, R. (2010) A New Analytic Solution for the 2D Advection-Dispersion Equation in Semi-Infinite and Laterally Bounded Domain. Applied Mathematical Sciences, 4, 3733-3747.
[9]	Chatwin, P.C. and Allen, C.M. (1985) Mathematical Models of Dispersion in Rivers and Estuaries. Annual Review of Fluid Mechanics, 17, 119-149. https://doi.org/10.1146/annurev.fl.17.010185.001003
[10]	Chauhry, M.H., Cass, D.E. and Edinger, J.E. (1983) Modeling of Unsteady Flow Water Temperatures. Journal of Hydraulic Engineering, 109, 657-669. https://doi.org/10.1061/(ASCE)0733-9429(1983)109:5(657)
[11]	Kalita, J.C., Dalal, D.C. and Dass, A.K. (2002) A Class of Higher Order Compact Schemes for the Unsteady Two-Dimensional Convection—Diffusion Equation with Variable Convection Coefficients. International Journal for Numerical Methods in Fluids, 38, 1111-1131. https://doi.org/10.1002/fld.263
[12]	Ekebjaerg, L. and Justesen, P. (1991) An Explicit Scheme for Advection-Diffusion Modelling in Two Dimensions. Computer Methods in Applied Mechanics and Engineering, 88, 287-297. https://doi.org/10.1016/0045-7825(91)90091-J
[13]	Dehghan, M. (2004) Numerical Solution of the Three Dimensional Advection Diffusion Equation. Applied Mathematics and Computation, 150, 5-19. https://doi.org/10.1016/S0096-3003(03)00193-0
[14]	Greenside, H.S. and Cross, M.C. (1985) Stability Analysis of Two-Dimensional Models of Three-Dimensional Convection. Physical Review A, 31, 2492-2501. https://doi.org/10.1103/PhysRevA.31.2492
[15]	Noye, B.J. and Tan, H.H. (1989) Finite Difference Methods for Solving the Two Dimensional Advection Diffusion Equation. International Journal for Numerical Methods in Fluids, 9, 75-98. https://doi.org/10.1002/fld.1650090107
[16]	Hirahara, H., Iwazaki, K., Ahmmed, M.U. and Nakamura, M. (2011) Numerical Analysis of Air Flow in Dichotomous Respiratory Channel with Asymmetric Compliance under HFOV Condition. Journal of Fluid Science and Technology, 6, 932-948. https://doi.org/10.1299/jfst.6.932
[17]	Calay, R.K., Kurujareon, J. and Holdø, A.E. (2002) Numerical Simulation of Respiratory Flow Patterns within Human Lung. Respiratory Physiology & Neurobiology, 130, 201-221. https://doi.org/10.1016/S0034-5687(01)00337-1
[18]	Tsega, E.G. (2018) Computational Fluid Dynamics Modeling of Respiratory Airflow in Tracheobronchial Airways of Infant, Child and Adult. Computational and Mathematical Methods in Medicine, 2018, Article ID: 9603451. https://doi.org/10.1155/2018/9603451
[19]	Issakhov, A., Zhandaulet, Y., Abylkassymova, A. and Issakhov, A. (2021) A Numerical Simulation of Air Flow in the Human Respiratory System for Various Environmental Conditions. Theoretical Biology and Medical Modeling, 18, Article No. 2. https://doi.org/10.1186/s12976-020-00133-8
[20]	Castka, J.F., Metcalfe, H.C., Davis, R.E. and Williams, J.E. (2002) Modern Chemistry. Holt, Rinehart and Winsto, Austin.