The Velocity of PCL Fluid in Human Lungs with Beaver and Joseph Boundary Condition by Using Asymptotic Expansion Method

Abstract

Humans breathe air into the respiratory system through the trachea, but with all the pollutants in our environment (both outside and inside), the air we breathe may not be clean. When that is so, the respiratory system secretes mucus to trap dirt that is inhaled through the nostrils. The respiratory tract contains hair-like structures in the epithelial tissue, called cilia: These wave back and forth to help expel particles of dust, dirt, mucus, and contaminants from the body. Cilia are found in this layer (a porous medium) and the fluid in this layer is called the periciliary layer (PCL). This study aims to determine the velocity of the PCL fluid flow in motile cilia. Usually, fluids move due to pressure changes, but in this study, the velocity of solids or of the cilia moves the PCL fluid. Stokes-Brinkman equations are used to determine the velocity of PCL fluid flow when cilia form an angle with the horizontal plane. The Beavers and Joseph boundary condition is applied in this study. The asymptotic expansion method is adapted in order to determine the velocity of PCL from the movement of the cilia.

1. Introduction

The human body contains numerous cilia, which are omnipresent inside and outside of the body. Cilia are hair-like organelles that provide locomotion to liquids throughout the body. This study focuses on cilia inside the body, particularly in the respiratory system. The human respiratory system is illustrated in Figure 1, with labeling of the nose, trachea, and lungs. We breathe in and out all the time, and as we inhale, some dust and pathogens will enter the body. However, the human body has a defense mechanism to protect against these foreign particles. In a closer look at the immune system, a diagram representing a cross-section of the trachea is displayed in Figure 2, and Figure 3 shows a close-up view of the trachea.

Figure 3 shows the platform of the innate immune system of the throat cells, which consists of three layers: The air layer, the mucus layer, and the periciliary layer (PCL). Cilia are found in the PCL layer. The function of the cilia is to filter out dust and other inhaled foreign particles. Cilia work as a part of the immune system, protecting the body against pathogens in the air. Goblet cells (vital component cells of the immune system) secrete mucus to trap the inhaled particles and cilia help to transport these particles out of the body by cilia-generated flow. To calculate the mucus velocity, the velocity of PCL fluid is measured to determine the bottom boundary of the mucus layer.

Numerical studies on the PCL and mucus layers have been carried out by several researchers. For example, Machemer [1] observed the ciliary system and found that frequency of peristomal cilia decreases with increasing viscosity, which leads to an increase of the average wavelength from 10.7 at 1 cP to 14.3 at 40 cP. Smith et al. [2] developed a new model for the nodal flow, utilizing the regularized stokeslet method. Fulford and Blake [3] studied a two-layer Newtonian fluid model for muco-ciliary transportation in the lung, and discovered that the viscosity of the upper mucous layer is much greater than the viscosity of the lower periciliary layer. Jayathilake et al. [4] developed a three-dimensional numerical model to simulate human pulmonary cilia motion in the PCL using the immersed boundary method combined with the projection method. Their numerical results indicated that the phase differences of cilia—in both stream-wise and span-wise directions—resulted in the maximum PCL velocity in the stream-wise direction.

Another approach on the study of the PCL and mucus layers is experimental study, which has been examined by a number of researchers. For example, Serafini and Michaelson [5] determined ciliary length and the percentage of ciliated cells from six mongrel dogs and ten humans. Matsui et al. [6] studied the movements of mucus and PCL liquid in airway surfaces, using conventional and confocal microscopy of fluorescent microspheres photoactivated fluorescent dyes and well-differentiated human tracheobronchial epithelial cell cultures that exhibited spontaneous, radial mucociliary transport. Their findings showed that the entire PCL liquid was transported at approximately the same rate as mucus, 39 ± 24.7 and 39.8 ± 4.2 µm/s. Moreover, they found that the removal of the mucus layer reduced PCL transport by >80%, 4.8 ± 0.6 µm/s, which is close to the value predicted by theoretical analyses of the ciliary beat cycle. Therefore, it has been suggested that the movement of PCL liquid depends on the transport of mucus.

Although the study of PCL and mucus layers, both numerical and experimental studies, has been investigated extensively, these studies covered only movement, length, and direction of cilia; the research on velocity of cilia in the PCL layer is limited. In general, the asymptotic expansion method is used to examine the PCL fluid flow order to determine the velocity of PCL fluid. Cilia in the PCL layer are shown in Figure 4; cilia making a 90° angle with the horizontal plane are displayed in Figure 4a, while cilia at an angle of less than 90° are illustrated in Figure 4b. The domain ${\mathsf{\Omega }}_1$ is the PCL with an absence of cilia. The second domain, ${\mathsf{\Omega }}_2$ is the PCL containing cilia.

As the domain ${\mathsf{\Omega }}_2$ contains both liquid and solid phases, it is considered to be a porous medium. For the mathematical model used in the free fluid and the porous medium domains, some studies have used Navier-Stokes in the free fluid layer and Darcy in the porous medium [7], or Stokes equation in the free fluid layer and the Brinkman equation in the porous medium [8,9]. The Beavers-Joseph condition has been widely used in many numerical studies as well as to compare solutions, especially between a porous medium and a free fluid domain. For example, Francisco et al. [10] derived boundary conditions that complete the statement of a two-domain approach for a one-dimensional momentum transport in a system containing a porous medium and free fluid under a constant pressure gradient. These boundary conditions involve a jump in both the velocity and the viscous stress and are thus expressed in terms of jump coefficients.

In this paper, we use the Stokes equation in the free fluid region and the Brinkman equation in the porous medium in order to analytical study the Stoke–Brinkman equations. The asymptotic expansion method has been widely used in a number of studies [11,12,13]. Chandesris and Jamet [14] found the boundary condition at the interface between the porous medium and the free fluid based on the Poiseuille flow over a permeable block. The problem was solved using the matched asymptotic expansion method, therefore, both a heterogeneous transition zone and a homogeneous zone between the two homogeneous regions (porous medium and free fluid) were considered. The two homogeneous regions were described using volume averaged transport equation, with the assumption that this equation still holds in the heterogeneous transition zone by considering variable porosity and permeability.

However, to the authors’ knowledge, no study to date has used asymptotic expansion related to the velocity of the cilia in the model. The model used in the present study involves the velocity of solid helping to move the PCL fluid. In this study, we provide the velocity of the PCL fluid by using asymptotic expansion due to the movement of cilia where the Beaver–Joseph boundary condition is employed at the free fluid/porous medium interface.

In Section 2, we write n-dimensional Stokes–Brinkman equations in one dimension by using the indicial notation. The dimensionization of the governing equations is derived in Section 3. The solutions of the Stokes–Brinkman equations are computed with the asymptotic expansion method in Section 4. In Section 5, we determine the relation of the constants. We describe result and discussion in Section 6, and in Section 7 we draw our conclusions.

2. Mathematical Model and Boundary Conditions

In this section, we introduce our governing equations: the n-dimensional Stokes–Brinkman equations and the boundary conditions used in the PCL.

2.1. Stokes–Brinkman Equations

The n-dimensional Stokes–Brinkman equations were derived using an upscaling method, which is a method that helps to change the equations from a microscopic scale to a macroscopic scale. Then, we simplified the equations to one dimension by using an indicial notation. The macroscale Stokes–Brinkman equations [12] and the continuity equation [12] in n dimension are

$$\mu {\mathrm{k}}^{-1}\cdot \left({\epsilon }^l{\mathrm{v}}^l\right)+\nabla p-\frac{\mu }{{\epsilon }^l}\Delta \left({\epsilon }^l{\mathrm{v}}^l\right)=\rho \mathrm{g}+\mu {\mathrm{k}}^{-1}\cdot {\epsilon }^l{\mathrm{v}}^s+\frac{\mu }{{\epsilon }^l}\nabla f,$$

(1)

$$\nabla \cdot \left({\epsilon }^l{\mathrm{v}}^l\right)=f,$$

(2)

respectively, where the function $f=\frac{-{\epsilon }^{\cdot l}}{\left(1-{\epsilon }^l\right)}+\nabla \cdot {\epsilon }^l{\mathrm{v}}^s;$ $\mu$ is a dynamic viscosity; $k^{-1}$ is the inverse of the permeability tensor; ${\epsilon }^l$ is the porosity; $v^l$ and $v^S$ are the velocities of the liquid and solid phases, respectively; $p$ is the pressure; $\rho$ is the fluid density; $g$ is gravity; and ${\epsilon }^{\cdot l}$ is the material time derivative of the porosity with respect to the solid phase, ${\epsilon }^{\cdot l}=\frac{\partial {\epsilon }^l}{\partial t}+{\mathrm{v}}^s\cdot \nabla {\epsilon }^l$. Equation (1) is the Brinkman equation. Without the first and fifth terms in the equation, it is a Stokes equation. Thus, in general, Equation (1) is called a Stokes–Brinkman equations. Applying the indicial notation to the Stokes–Brinkman equations we obtain

$$\mu {\mathrm{k}}_{\mathrm{ij}}^{-1}\left({\epsilon }^l{v}^l{}_j\right)+p_{,i}-\frac{\mu }{{\epsilon }^l}{\left({\epsilon }^l{v}_{\mathrm{i}}{}^l\right)}_{,jj}=\rho {\mathrm{g}}_{\mathrm{i}}+\mu {\mathrm{k}}_{\mathrm{ij}}{}^{-1}{\epsilon }^l{v}_j{}^s+\frac{\mu }{{\epsilon }^l}{f}_{,i},i=1,2,3,\dots ,n$$

(3)

$${\epsilon }^l{v}_{i,i}^l=f,i=1,2,3,\dots ,n$$

(4)

where the repeat index indicates the summation. Substituting Equation (4) into the last term of Equation (3), for the one-dimensional domain, we have

$$\mu {\mathrm{k}}_{11}^{-1}{\epsilon }^l{v}_1{}^l+p_{,1}-\frac{\mu }{{\epsilon }^l}{\left({\epsilon }^l{v}_1{}^l\right)}_{,11}=\rho {\mathrm{g}}_1+\mu {\mathrm{k}}_{11}^{-1}\left({\epsilon }^l{v}_1{}^s\right)+\frac{\mu }{{\epsilon }^l}{\left({\epsilon }^l{v}_{1,1}{}^l\right)}_{,1}$$

(5)

Note that ${\left({\epsilon }^l{v}_{1,1}{}^l\right)}_{,1}={\left({\epsilon }^l{v}_1{}^l\right)}_{,11}$. Simplifying Equation (5) yields

$$\frac{2\mu }{{\epsilon }^l}{\left({\epsilon }^l{v}_1{}^l\right)}_{,11}-\mu {\mathrm{k}}_{11}^{-1}\left({\epsilon }^l{v}_1{}^l\right)+\rho {\mathrm{g}}_1=p_{,1}-\mu {\mathrm{k}}_{11}^{-1}\left({\epsilon }^l{v}_1{}^s\right)$$

(6)

In this study, we assume that the solid velocity depends on only the $Y$ direction; the fluid flow along the $x$-axis; and the pressure depends on the $X$ direction. As a result, Equation (6) becomes

$$\frac{2\mu }{{\epsilon }^l}\frac{d^2\left({\epsilon }^l{v}_1{}^l\right)}{d{y}^2}-\mu {\mathrm{k}}_{11}^{-1}\left({\epsilon }^l{v}_1{}^l\right)+\rho {\mathrm{g}}_1=\frac{dp}{dx}-\mu {\mathrm{k}}_{11}^{-1}\left({\epsilon }^l{v}_1{}^s\right)$$

(7)

Equation (7) is a Brinkman equation, which is used in the domain containing cilia.

In region ${\mathsf{\Omega }}_1$, we use a Stokes equation,

$$\frac{2\mu }{{\epsilon }^l}\frac{d^2\left({\epsilon }^l{v}_1{}^l\right)}{d{y}^2}+\rho {\mathrm{g}}_1=\frac{dp}{dx}-\mu {\mathrm{k}}_{11}^{-1}\left({\epsilon }^l{v}_1{}^s\right)$$

(8)

The continuity equation applied in this region is

$$\frac{\partial \left({\epsilon }^l{v}^l\right)}{\partial y}=0$$

(9)

For the free fluid region, we use a Stokes equation in domain ${\mathsf{\Omega }}_1$. The equation is actually the Brinkman equation, without the permeability term.

Notice that we have two different and adjacent domains. Therefore, a boundary used at the interface is special. In this work we applied the Beavers–Joseph boundary condition at the interface between the free fluid region and the porous medium domain.

2.2. Boundary Conditions

No free fluid exists in the PCL when the cilia are perpendicular to the horizontal plane. Figure 5a shows that the PCL has only domain ${\mathsf{\Omega }}_2$. Figure 5b shows our domains ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$; the boundary condition between ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2$ is at $y_{Stoke}$ at the tips of the cilia. Note that $y_{Stoke}$ depends on the angle $\theta$. We used 3 boundary conditions in this study $u=0$ at $y=0$, $u=G$ at $y=h$ and $y_{Stoke}$ Beavers–Joseph boundary condition at $y_{Stoke}$.

Let $u={\epsilon }^l{v}_1{}^l$. We now have the system of equations with the unknown $u$ and boundary conditions;

$$u=0\mathrm{at}y=0$$

(10)

$$u=G\mathrm{at}y=h$$

(11)

$${\frac{du}{dy}|}_{y=y^{+}{}_{stoke}}-\frac{1}{{\phi }_p}{\frac{du}{dy}|}_{y=y^{-}{}_{stoke}}=\frac{\beta }{\sqrt{K}}{u|}_{y=y_{stoke}}\mathrm{at}y=y_{Stoke}$$

(12)

where $\beta$ is a dimensionless parameter of the order of one; ${\epsilon }^l$ is porosity of the porous medium; $y_{Stoke}$ is length of porous medium; and $G$ is a function depending on $X$ and the angle $\theta$.

3. Dimensionless Stokes–Brinkman Equations

In this section, we normalize the Stokes–Brinkman equations by introducing the following new variables.

$$y^{+}=\frac{y}{h},U^{+\left(0\right)}=\frac{u}{U_0},v^{+}=\frac{v_1^s}{U_0},x^{+}=\frac{x}{h},K^{+}=\frac{k_{11}}{K_0},P^{+}=p\frac{h}{U_0\mu }\mathrm{and}{g}^{+}=\frac{g}{g_0}$$

(13)

where $h$ is the characteristic length; $U_0$ is volumetric average velocity in the porous medium; $K_0$ is characteristic permeability; and $g_0$ is characteristic gravity.

3.1. Dimensionless of Brinkman Equation

We normalize the Brinkman equation, of Equation (7) by substituting Equation (13) for Equation (7). We then have

$$\frac{2}{{\epsilon }^l}\frac{d^2{U}^{+\left(0\right)}}{d{y}^{+}{}^2}-\frac{h^2}{K_0{K}^{+}}{U}^{+\left(0\right)}+\frac{h^2\rho {\mathrm{g}}^{+}{\mathrm{g}}_0}{\mu U_0}=\frac{d{P}^{+}}{d{x}^{+}}-\frac{h^2{\epsilon }^l{v}^{+}}{K_0{K}^{+}}$$

(14)

Equation (14) can be rewritten as

$$2M_1\frac{d^2{U}^{+\left(0\right)}}{d{y}^{+}{}^2}-M_2{U}^{+\left(0\right)}+M_3=\frac{d{P}^{+}}{d{x}^{+}}-H\left(y^{+}\right),$$

(15)

where $M_1=\frac{1}{{\epsilon }^l},M_2=\frac{h^2}{K_0{K}^{+}},M_3=\frac{h^2\rho {\mathrm{g}}^{+}{\mathrm{g}}_0}{\mu U_0}$are constants and

$$H\left(y^{+}\right)=\frac{h^2{\epsilon }^l{v}^{+}}{K_0{K}^{+}}\mathrm{is}\mathrm{a}\mathrm{function}\mathrm{of}{y}^{+}.$$

(16)

3.2. Dimensionless Stokes Equation

We normalize the Stokes Equation (8) by substituting Equation (13) for Equation (8). We then obtain

$$\frac{2}{{\epsilon }^l}\frac{d^2{U}^{+\left(0\right)}}{d{\left(y^{+}\right)}^2}=\frac{d{P}^{+}}{d{x}^{+}}-\frac{h^2\rho g^{+}{g}_0}{\mu U_0}$$

(17)

Equation (17) can be rewritten as

$$2M_1\frac{d^2{U}^{+}}{d{y}^{+}{}^2}=\frac{d{P}^{+}}{d{x}^{+}}-M_3$$

(18)

where $M_1$ and $M_3$ are the constants defined in Equation (16).

3.3. Dimensionless Boundary Conditions

In this section we normalize the boundary conditions at the free fluid/porous medium interface used in our domain.

Substituting Equation (13) for Equation (12) we have

$${\frac{d{U}^{+\left(0\right)}}{d{y}^{+}}|}_{y^{+}=\frac{y_{stoke}}{h}}-\frac{1}{{\phi }_p}{\frac{d{U}^{+\left(0\right)}}{d{y}^{+}}|}_{y^{+}=\frac{y_{stoke}}{h}}=\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}{U^{+\left(0\right)}|}_{y^{+}=\frac{y_{stoke}}{h}.}$$

(19)

Next, we normalize the bottom and top boundary conditions by substituting Equation (13) for Equation (10) and Equation (11), giving us

$$U^{+\left(0\right)}=0\mathrm{at}{y}^{+}=0.$$

(20)

$$U^{+\left(0\right)}=\frac{G}{U_0}\mathrm{at}{y}^{+}=1.$$

(21)

respectively.

We have now normalized the Stokes–Brinkman equations and the boundary conditions. Consequently, the system of equations used in the next sections and the boundary conditions are

$$2M_1\frac{d^2{U}^{+\left(0\right)}}{d{y}^{+}{}^2}-M_2{U}^{+\left(0\right)}+M_3=\frac{d{P}^{+}}{d{x}^{+}}-H\left(y^{+}\right)\mathrm{in}{\mathsf{\Omega }}_2$$

(22)

$$2M_1\frac{d^2{U}^{+}}{d{y}^{+}{}^2}=\frac{d{P}^{+}}{d{x}^{+}}-M_3\mathrm{in}{\mathsf{\Omega }}_1$$

(23)

$$U^{+\left(0\right)}=0\mathrm{at}{y}^{+}=0.$$

(24)

$$U^{+\left(0\right)}=\frac{G}{U_0}\mathrm{at}{y}^{+}=1.$$

(25)

$${\frac{d{U}^{+\left(0\right)}}{d{y}^{+}}|}_{y^{+}=\frac{y_{stoke}}{h}}-\frac{1}{{\phi }_p}{\frac{d{U}^{+\left(0\right)}}{d{y}^{+}}|}_{y^{+}=\frac{y_{stoke}}{h}}=\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}{U^{+\left(0\right)}|}_{y^{+}=\frac{y_{stoke}}{h}.}\mathrm{at}{y}^{+}=\frac{y_{Stoke}}{h}$$

(26)

4. Asymptotic Expansion Method of the Stokes–Brinkman Equations

As mentioned in Section 2, we have two different domains where the domain ${\mathsf{\Omega }}_1$, is the region where $y>y_{Stoke}$ and the region ${\mathsf{\Omega }}_2$ is the layer where $y<y_{Stoke}$. In this section, we find the analytical solution of the governing equations with the asymptotic expansion. Section 4.1 shows the result of the Brinkman equation and Section 4.2 shows the result of the Stokes equation.

4.1. Asymptotic Expansion Method of the Brinkman Equation

Using the method of asymptotic expansion, we assume that

$$U^{+}=U^{+\left(0\right)}+\epsilon U^{+\left(1\right)}+{\epsilon }^2{U}^{+\left(2\right)}+\dots$$

(27)

In the case $y^{+}<\frac{y_{Stoke}}{h}$, we substitute Equation (27) for Equation (22), giving us

$$2M_1\frac{d^2(U^{+\left(0\right)}+\epsilon U^{+\left(1\right)}+o\left({\epsilon }^2\right)}{d{y}^{+}{}^2}-M_2\left(U^{+\left(0\right)}+\epsilon U^{+\left(1\right)}+o\left({\epsilon }^2\right)\right)+M_3=\frac{d{P}^{+}}{d{x}^{+}}-H\left(y^{+}\right)$$

(28)

Considering the zeroth-order of ${\epsilon }^0$; we obtain

$$2M_1\frac{d^2{U}^{+\left(0\right)}}{d{y}^{+}{}^2}-M_2{U}^{+\left(0\right)}+M_3=\frac{d{P}^{+}}{d{x}^{+}}-H\left(y^{+}\right)$$

(29)

Since $H\left(y^{+}\right)$, the velocity of cilia, is a known source term in this work, we assume that it is a linear function. That is

$$H\left(y^{+}\right)=\frac{h^2{\epsilon }^l}{k_{11}}\left(c_1{y}^{+}+c_0\right)$$

where $c_0$ and $c_1$ are constants. Thus, Equation (29) becomes

$$2M_1\frac{d^2{U}^{+\left(0\right)}}{d{y}^{+}{}^2}-M_2{U}^{+\left(0\right)}=-M_3+\frac{d{P}^{+}}{d{x}^{+}}-\frac{h^2{\epsilon }^l}{k_{11}}\left(c_1{y}^{+}+c_0\right)$$

(30)

First, we find the general solution of the homogeneous Equation (30). Solving the homogeneous part of the ordinary differential equation

$$2M_1\frac{d^2{U}^{+\left(0\right)}}{d{y}^{+}{}^2}-M_2{U}^{+\left(0\right)}=0.$$

(31)

We obtain the general solution $U_c{}^{+\left(0\right)}=w_1{e}^{\sqrt{\frac{h^2{\epsilon }^l}{2K_0{K}^{+}}}{y}^{+}}+w_2{e}^{-\sqrt{\frac{h^2{\epsilon }^l}{2K_0{K}^{+}}}{y}^{+}}$. To find the particular solution of Equation (30), we use the method of undetermined coefficients. Therefore, the solution of the differential equation is

$$U^{+\left(0\right)}=w_1{e}^{\sqrt{\frac{M_2}{2M_1}}{y}^{+}}+w_2{e}^{-\sqrt{\frac{M_2}{2M_1}}{y}^{+}}+\frac{h^2{\epsilon }^l}{M_2{k}_{11}}{c}_1{y}^{+}+\frac{M_3}{M_2}-\frac{d{P}^{+}}{M_{2}d{x}^{+}}+\frac{h^2{\epsilon }^l}{M_2{k}_{11}}{c}_0.$$

(32)

Let

$$J_1=\sqrt{\frac{M_2}{2M_1}},J_2=\frac{d{P}^{+}}{M_{2}d{x}^{+}},J_3=\frac{h^2{\epsilon }^l}{M_2{k}_{11}},J_4=\frac{M_3}{M_2},J_5=\frac{d{P}^{+}}{M_{1}d{x}^{+}},\mathrm{and}{J}_6=\frac{M_3}{M_1}$$

(33)

all of which are constants. Therefore Equation (32) can be rewritten as

$$U^{+\left(0\right)}=w_1{e}^{J_1{y}^{+}}+w_2{e}^{-J_1{y}^{+}}+J_3{c}_1{y}^{+}+J_4-J_2+J_3{c}_0.$$

(34)

Applying Equation (24) to Equation (34), we have

$$w_2=J_2-J_4-J_3{c}_0-w_1$$

(35)

which yields

$$U^{+\left(0\right)}=w_1{e}^{J_1{y}^{+}}-w_1{e}^{-J_1{y}^{+}}+\left(J_2-J_4-J_3{c}_0\right)e^{-J_1{y}^{+}}+J_3{c}_1{y}^{+}+J_4-J_2+J_3{c}_0.$$

(37)

We now have the solution of the Brinkman equation.

4.2. Asymptotic Expansion Method of Stokes Equation

Similar to the previous section, we calculate the velocity of the PCL fluid in domain ${\mathsf{\Omega }}_1$ using the asymptotic expansion method with the Stokes Equation (23). We assume that

$$U^{+}=U^{+\left(0\right)}+\epsilon U^{+\left(1\right)}+{\epsilon }^2{U}^{+\left(2\right)}+\dots$$

(37)

In the case $y^{+}>\frac{y_{Stoke}}{h}$, we substitute Equation (37) for Equation (23). Therefore

$$2M_1\frac{d^2(U^{+\left(0\right)}+\epsilon U^{+\left(1\right)}+o\left({\epsilon }^2\right)}{d{y}^{+2}}=\frac{d{P}^{+}}{d{x}^{+}}-M_3.$$

(38)

Considering the zeroth-order term of $\epsilon$, ${\epsilon }^0$, we obtain

$$\begin{array}{c}\frac{d^2{U}^{+\left(0\right)}}{d{y}^{+2}}=\frac{1}{2M_1}\frac{d{P}^{+}}{d{x}^{+}}-\frac{M_3}{2M_1}\\ \mathrm{or}\frac{d^2{U}^{+\left(0\right)}}{d{y}^{+2}}=\frac{1}{2}{J}_5-\frac{1}{2}{J}_6.\end{array}$$

(39)

As we assume that $\frac{d{P}^{+}}{d{x}^{+}}$ is a constant, by integrating Equation (39) twice, we have

$$U^{+\left(0\right)}=\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)y^{+}{}^2}{2}+w_3{y}^{+}+w_4$$

(40)

Applying the boundary condition (25) to Equation (40), we obtain

$$\begin{array}{c}\frac{G}{U_0}=\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)}{2}+w_3+w_4.\mathrm{Hence}{w}_4=\frac{G}{U_0}-\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)}{2}-w_3.\\ \mathrm{Then}{U}^{+\left(0\right)}=\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)y^{+}{}^2}{2}+w_3{y}^{+}+\frac{G}{U_0}-\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)}{2}-w_3.\end{array}$$

(41)

We can find the solutions at and ${\mathsf{\Omega }}_2$, yet the constants $w_1$ and $w_3$ of the Stokes–Brinkman equation are unknown.

5. The Relation between the Constants

In this section, we determine the relation of the constants $w_1$ and $w_3$ with Beavers–Joseph boundary condition. Substituting Equation (36) and Equation (41) for Equation (26), we obtain

$$\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)\frac{y_{Stoke}}{h}+w_3-\frac{1}{{\phi }_p}{w}_1{J}_1{e}^{J_1\frac{y_{Stoke}}{h}}-\frac{1}{{\phi }_p}{w}_1{J}_1{e}^{-J_1\frac{y_{Stoke}}{h}}+\frac{1}{{\phi }_p}{J}_1\left(J_2-J_4-J_3{c}_0\right)e^{-J_1\frac{y_{Stoke}}{h}}-\frac{1}{{\phi }_p}{J}_3{c}_1=\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}{U^{+}|}_{y^{+}=\frac{y_{stoke}}{h},}$$

(42)

where ${U^{+}|}_{y^{+}=\frac{y_{stoke}}{h}}$ (on the right-hand side of Equation (42)) can be obtained from the PCL velocity in ${\mathsf{\Omega }}_1$${\mathsf{\Omega }}_2$, or both. However, the velocity $\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}{U^{+}|}_{y^{+}=\frac{y_{stoke}}{h}}$ is a slip velocity [13,14]. Therefore, in this work, we consider the case that ${U^{+}|}_{y^{+}=\frac{y_{stoke}}{h}}={U^{+}|}_{y^{+}=\frac{y_{stoke}{}^{+}}{h}}$.

Substituting Equation (41) for the right-hand side of Equation (42) and simplifying, we have

$$w_3=\frac{1}{\left(1-\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{y_{Stoke}}{h}+\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\right)}\left[\begin{array}{l}{w}_1\left(\frac{1}{{\phi }_p}{J}_1{e}^{J_1\frac{y_{Stoke}}{h}}+\frac{1}{{\phi }_p}{J}_1{e}^{-J_1\frac{y_{Stoke}}{h}}\right)+\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right){\left(\frac{y_{Stoke}}{h}\right)}^2}{2}+\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{G}{U_0}\\ -\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)}{2}-\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)\frac{y_{Stoke}}{h}-\frac{1}{{\phi }_p}{J}_1\left(J_2-J_4-J_3{c}_0\right)e^{-J_1\frac{y_{Stoke}}{h}}+\frac{1}{{\phi }_p}{J}_3{c}_1\end{array}\right].$$

(43)

We now have the relation between the constants $w_1$ and $w_3$. Inserting Equation (42) into Equation (41), we have the solution of the Stokes equation depending on $w_1$. That is

$$\begin{array}{l}{U}^{+\left(0\right)}=\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)y^{+}{}^2}{2}+\left\{\frac{1}{\left(1-\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{y_{Stoke}}{h}+\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\right)}\left[\begin{array}{l}{w}_1\left(\frac{1}{{\phi }_p}{J}_1{e}^{J_1\frac{y_{Stoke}}{h}}+\frac{1}{{\phi }_p}{J}_1{e}^{-J_1\frac{y_{Stoke}}{h}}\right)+\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right){\left(\frac{y_{Stoke}}{h}\right)}^2}{2}+\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{G}{U_0}\\ -\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)}{2}-\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)\frac{y_{Stoke}}{h}-\frac{1}{{\phi }_p}{J}_1\left(J_2-J_4-J_3{c}_0\right)e^{-J_1\frac{y_{Stoke}}{h}}+\frac{1}{{\phi }_p}{J}_3{c}_1\end{array}\right]\right\}{y}^{+}\\ +\frac{G}{U_0}-\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)}{2}-\left\{\frac{1}{\left(1-\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{y_{Stoke}}{h}+\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\right)}\left[\begin{array}{l}{w}_1\left(\frac{1}{{\phi }_p}{J}_1{e}^{J_1\frac{y_{Stoke}}{h}}+\frac{1}{{\phi }_p}{J}_1{e}^{-J_1\frac{y_{Stoke}}{h}}\right)+\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right){\left(\frac{y_{Stoke}}{h}\right)}^2}{2}+\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{G}{U_0}\\ -\frac{h\beta }{\sqrt{K_p^{+}{K}_0}}\frac{\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)}{2}-\left(\frac{1}{2}{J}_5-\frac{1}{2}{J}_6\right)\frac{y_{Stoke}}{h}-\frac{1}{{\phi }_p}{J}_1\left(J_2-J_4-J_3{c}_0\right)e^{-J_1\frac{y_{Stoke}}{h}}+\frac{1}{{\phi }_p}{J}_3{c}_1\end{array}\right]\right\}.\end{array}$$

(44)

We now obtain the general solution of the Stokes–Brinkman equations, depending only the constant $w_1$ with the asymptotic expansion method, in which the Beavers–Joseph condition is employed at the free fluid/ porous medium interface.

6. Results and Discussion

The solutions of Stokes–Brinkman equations, Equation (36) and Equation (44), calculated in Section 4 and Section 5 are plotted on a graph and discussed in this section. To plot the graph, we assume that the graphs $\theta$ between the cilia and the horizontal plane are ${40}^{\circ },{50}^{\circ },{60}^{\circ },{70}^{\circ },{80}^{\circ },{90}^{\circ }$ where the cilia have the highest velocity at $\theta ={90}^{\circ }$ during the forward stroke and stop beating at $\theta ={40}^{\circ }$. So, at $\theta ={40}^{\circ }$ we do not need to calculate the velocity of the PCL. The variables used in Equation (44) are given in

Table 1. Variable of the experimental data in the solution of the Stoke–Brinkman equations.

Variable	Value	Unit
h	7.5	$\left[\mu m\right]$
$\rho$	$992.2\times {10}^{-15}$	$\left[N/m^2\right]$
$\mu$	$3\times {10}^{-6}$	$\left[M{L}^{-1}{T}^{-2}\right]$
$\mathrm{g}$	$9.81\times {10}^6$	$\left[\mu m\right]$
$U_0$	1.00	$\left[\mu m\right]$
$\beta$	1	$\left[1\right]$
$v^{+}$	1	$\left[1\right]$
$\frac{d{P}^{+}}{d{x}^{+}}$	1	$\left[1\right]$
$w_1$	0	$\left[1\right]$
$G$	1	$\left[1\right]$

and

Table 2. Values of variables used in Equation (44) used to plot Figure for each angle $\theta$.

Variable	${50}^{\circ }$	${60}^{\circ }$	${70}^{\circ }$	${80}^{\circ }$	${90}^{\circ }$
$y_{Stoke}$	0.7672	0.8638	0.9353	0.9818	1.0000
$K_{P^{+}}$	0.0012	0.0015	0.0016	0.0017	0.0018
${\epsilon }^l$	0.6717	0.7099	0.7331	0.7439	0.7487
$v^s$	0.0023	0.0024	0.0044	0.0049	0.0050

. The values of $y_{Stoke}$ and $K_p{}^{+}$ are obtained from [4], and $v^s$ is the velocities of the solid, cilia, where the maximum velocity of motile cilia value is assumed at $\theta ={90}^{\circ }$ in descending order of degree, terminating at $\theta ={40}^{\circ }$, $h$ is the cilium height, $U_0$ is one cycle of ciliary beat, and ${\phi }_p$ and $\beta$ were obtained from [2]. We assume that the rate of pressure changes $\frac{d{P}^{+}}{d{x}^{+}}$ equals 1. The boundary condition at the top of the free fluid domain equals 1. The variable $w_1$ in Equation (44) equals 0 because in this study we aim to verify if the solid velocities applied to the equation are valid, the solutions correlate with the variable velocity employed, where the maximum velocity is assumed at $\theta ={90}^{\circ }$ and decreases with decreasing angles.

These solutions are the velocity of PCL fluid in the layer containing cilia with different angle $\theta$ $\left({50}^{\circ },{60}^{\circ },{70}^{\circ },{80}^{\circ },{90}^{\circ }\right),$ respectively. It can be noted that the velocity decrease with decreasing angles and bends sharply at the seam of the porous medium and the free fluid region. The asymptotic solution $U^{+(0)}$ of the Stoke–Brinkman Equation (44) is shown in Figure 6 (the code is given in Appendix here).

The asymptotic expansion method is the tool for finding analytical approximate solutions to complicated practical problems, for example, the problem in ordinary differential equations in terms of regular perturbation and singular perturbation. We construct different asymptotic solutions inside and outside the region of rapid change and match them together to determine a global solution. Other methods used for finding exact and approximate solutions for linear and nonlinear partial differential equations is the homotopy perturbation method which is only a special case of the homotopy analysis method. The basic ideas of the homotopy perturbation method is a coupling of the traditional perturbation method and homotopy in topology deforms continuously the problem in hand to a simple one which can be easily solved. The advantage of this method is that it can be applied to various nonlinear problems, and the disadvantage is that we should suitably choose an initial guess, or infinite iterations are required [15]. The asymptotic expansion method and the homotopy perturbation method are principally based on a Taylor series with respect to an embedding parameter.

Another method used to describe the samples with porous media, is the fractal derivative model. This method is developed from fractal calculus for calculating the solution of problem phenomena in porous media or hierarchical structures [16]. For instance, in a study by Wang, Shi, He, and Li [17], to find an optimal hair length of a polar bear for thermal protection, which helps maintain a normal body temperature, by using a fractal derivative model of one-dimensional heat conduction along the hair. The basic ideas of the fractal calculus begin from Fourier’s law of thermal conduction, which can be expressed as $q=k\frac{\Delta T}{\Delta x},$ where $q$ is the heat flux, $k$ is the material’s conductivity, and $\frac{\Delta T}{\Delta x}$ is the temperature gradient, but the asymptotic expansion method is adopted to determine the velocity of ciliary motion that moves foreign particles out of the body.

While the asymptotic expansion method can be used to study the interface between the porous media region and the free fluid region, the homotopy perturbation method and the fractal derivative model are only adopted for the porous media domain. These three approaches are the tools for finding analytical approximate solutions to the problem in ordinary differential equations or partial differential equations. Another difference is that the asymptotic expansion method can be applied to solve linear ordinary differential equations and perturbation equations. On the other hand, the homotopy perturbation method and the fractal derivative model are used to derive the solutions to nonlinear ordinary differential equations or partial differential equations.

7. Conclusions

In this article, we studied the fluid flow in the periciliary layer (PCL) of the respiratory tract using the of asymptotic expansion method to determine the velocity of PCL fluid. The PCL is divided into two domains: the free fluid and the porous medium domains. The free fluid domain is the upper layer of the PCL from cilia are absent. The porous medium domain is the lower layer of PCL and contains cilia. The Stokes equation is employed in the free fluid domain and the Brinkman equation is applied in the porous medium. The Beavers–Joseph boundary condition is employed at the interface between the porous medium and the free fluid region. We assume that the velocity of the PCL fluid at the bottom of the porous medium domain (at $y=0$) is zero and the velocity of the PCL fluid in the porous medium depends on angle $\theta$. The boundary condition at the top of the free fluid domain is assumed to be a function of $G$, where $G$ depends on $X$ and the angle $\theta$. The explicit formula of velocity of PCL fluid $U^{+\left(0\right)}$ is obtained by applying asymptotic expansion to the Stokes–Brinkman equations.