# On the Asymptotic Boundary Condition at the Free-Fluid/Porous-Medium Interface in Periciliary Layer due to the Ciliary Movement

Authors: Sudaporn Poopra; Kanognudge Wuttanachamsri
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The purpose of this study is to study the movement of the periciliary (PCL) fluid due to the ciliary locomotion. In this research, because the bundle of cilia is considered instead of a single cilium, Stokes–Brinkman equations in a macroscopic scale are employed to find the velocity of the PCL fluid. When the cilia are perpendicular to the horizontal plane, the PCL consists of only the cilia. The inclination of cilia (cilia make an angleθθ<90∘ to the horizontal plane) results in two different domains in the PCL, the regions comprising and not comprising cilia. The main objective of this study is to determine the appropriate boundary conditions of the velocity between these two regions where the PCL fluid is moved by self-propelled cilia rather than the pressure gradient. A matched asymptotic expansion method is applied to the Stokes–Brinkman equations to determine the constraints. Two boundary conditions at the interface are obtained and the velocity of the PCL fluid at the top of PCL can be used to be the boundary conditions at the bottom of the mucus layer to determine the velocity of mucus. This model can help engineers to build devices to treat patients who have problem with the respiratory system. Applications include modeling fluid flow through filters such as engine filter and rice fields.


---

## Body

## 1. Introduction

Breathing is one of the most important bodily functions that keep us alive. When we breathe air in through the nose it is not just oxygen but also small substances containing particles and pathogens that pass through lungs. However, the respiratory system has filters to trap and remove the strange pathogens, thus allowing us to breathe without irritation. The filters lining the bronchus within the respiratory system are called cilia, hair-like structure, scattered throughout the ciliated cells, which are vital components of the immune system. Beside cilia, there are goblet cells scattering among the ciliated cells in the respiratory system. One of the primary functions of the goblet cells is to secrete mucus to trap the particles or microorganisms that pass through our respiratory system. Then mucus forms a layer at the tips of cilia. The cilia aid in transporting these foreign particles out of the body by sweeping the mucus upwards towards the throat. The tips of the cilia make contact with the mucus layer on the forward stroke and bend sideways and backwards on the reverse stroke. Essentially, this causes the mucus to be propelled only in the forward direction on the forward stroke [1]. The excreted mucus is ejected through the vocal cords and into the pharynx. The primary innate defense mechanism is called mucociliary clearance [2]. Under optimal lung conditions, cilia beat at 12 to 15 Hz in coordinated waves propelling mucus at 4 to 20 mm/min [1].Figure1 shows a small portion of the trachea cross section. Inside the segment there exist three main layers, air, mucus, and the periciliary layers (PCL). The PCL is a moist layer composed of cilia covering the ciliated cells. The fluid in this layer is called PCL fluid. Below the PCL is epithelial cells. Goblet cells secret mucus to trap the strange particles entering the body to protect the lungs from contaminants found in the epithelium. Among the goblet cells, there exist ciliated cells, where cilia are found on the top and wave in a rhythmic or pulsating motion and use that motion to keep passage ways free of mucus or foreign particles.Figure 1
PCL, mucus, and air layers in the respiratory tracts.In this work, we determine the velocity of the PCL fluid due to the movement of cilia. We consider that the movement of cilia results in two different patterns in the PCL as shown in Figure2. The left of Figure 2 shows the cilia perpendicular to the horizontal plane. The PCL has only one domain, which contains only cilia. When the cilia move forward and form an angle θ,θ<90∘ with the horizontal plane, the PCL can be divided into two domains, the layer containing cilia and the layer above the cilia as shown on the right of Figure 2. Since the PCL layer consists of the PCL fluid, assisting to treat the cilia to operate normally, and the solid phases, cilia, this part in PCL is considered as a porous medium. The layer above the cilia in the PCL is called free-fluid domain, as shown on the right of Figure 2. The PCL domain is denoted by ω=ω1∪ω2 where ω1 is the free-fluid region above the porous layer ω2 when the cilia form an angle θ,θ<90∘ with the horizontal plane. With the perpendicular case, we have only ω2. The variable yStoke is the height of the domain ω2. Notice that yStoke changes according to the angle θ. In this study, the first aim is to figure out the velocity of the PCL fluid in both domains to determine the boundary conditions between the domains ω1 and ω2.Figure 2
From left to right, the cilia are perpendicular and form an angleθ,θ<90∘ with the horizontal plane, respectively.Mathematical models used to calculate the velocity in the free-fluid region and the porous medium are numerous [3–8]. The equation applied in the free-fluid region was described by the Stokes equation or Navier–Stokes equations, while flow in the porous-medium region was usually modeled by using Darcy’s Law or Brinkman equations [6, 8–13]. In general, the Brinkman model is used for the porous media with high porosity [14]. Examples of flow in free-fluid region are diverse. For instance, Abbas et al. [15] studied the effect of the flow for an incompressible liquid with nonlinear differential equations. For flow in a porous medium, Ado-Elkhair et al. [16] presented the influence of the magnetic field and Darcy medium on oscillatory wavy walled (Peristalsis) flow due to the cilia motion by using a Adomian decomposition method (ADM) with a system of nonlinear partial differential equations. Zaher et al. [17] studied the effect of the ciliated walls in the uterine tube to strengthen the sperm to reach the egg. They used Darcy model with an oscillating wall and non-Newtonian magnetized fluid. An Adomian analysis method was employed to solve the equation. The illustrations of simultaneous flow in both free-fluid and porous regions are various; for example, Valdes-Parada and Lasseux [10] provided a new approach for modeling flow in porous and free-fluid domains based on Beaver and Joseph (BJ) configuration by using Darcy’s Law and Navier–Stokes equation. Carraro et al. [12] figured out a high-precision simulation of the slow viscous flow over a porous bed by applying a finite element method to Stokes equation in free-fluid domain and Darcy’s Law in porous bed. Naqvi and Bottaro [13] studied boundary conditions between an isotropic porous-medium and a free-fluid region by using the Stokes-like equation and homogenization theory. Du and Zu [18] proposed a local and parallel finite element method to the mixed Navier-Stokes/Darcy model with the Beaver–Joseph interface condition for an incompressible fluid flow. Wuttanachamsri [6] used one-dimensional Stokes–Brinkman equations to find the shape of free interface between a porous-medium and a free-fluid region.The two-layer configuration that uses different equations in different domains requires proper boundary conditions at the interface between the free-fluid and porous-medium regions. The difficulty in determining boundary conditions results from the fact that the orders of the corresponding differential operators are often different [14]. The boundary conditions at the interface between the porous-medium and the free-fluid regions had been studied by several authors [10, 19–23]. For example, Beavers and Joseph [22] proposed a slip boundary condition based on repeated experiments with a Poiseuille flow over a permeable porous block. Sahraoui and Kaviany [23] determined an appropriate hydrodynamic boundary condition at the interface between a porous-medium and a plain domain using the pressure correction method. Valdes-Parada and Lasseux [10] studied three different boundary conditions: creeping flow under no-slip conditions, inertial flow no-slip conditions, and creeping flow with slip conditions with Navier-Stokes–Darcy equations. The solutions were calculated by using Direct Numerical Simulation (DNS).Although the boundary conditions had been investigated extensively by many researchers, the studies of boundary conditions at the interface between the free-fluid region and the porous medium in the human lungs are limited. This study aims to determine the boundary conditions at the interface between the free-fluid region and the porous medium in PCL by using the method of matched asymptotic expansion, which has been used widely by considerable researchers. For example, Chandesris and Jamet [9] studied the Poiseuille flow over a permeable block by using the method of matched asymptotic expansion and determined the boundary conditions at the top of the permeable block. Poopra and Wuttanachamsri [11] calculated the velocity of the PCL fluid with the BJ boundary condition when cilia formed angles with the horizontal plane by using the asymptotic expansion method.The results in the reviewed literature, however, are not appropriated for our research problem. In most of the studies, solid phases in the porous domain are stationary. In our study, on the contrary, cilia in the porous domain are mobile. Although Poopra and Wuttanachamsri [11] studied the PCL with self-propelled solid phases, they did not seek for matching conditions to match the solutions at the transition zone at the free-fluid/porous-medium interface and they did not provide a proper boundary condition at the interface, which was practical to the problem. Thus, this research has been conducted to figure out the boundary conditions suiting for the problem.We first provide our mathematical model and boundary conditions at the top and bottom of our domain in dimensionless forms in Section2. Second, the velocity of the PCL fluid is calculated by using the asymptotic expansion method in Sections 3.1 and 3.2. Third, because unknown constants occurred during finding the asymptotic solutions, matching conditions presented in Section 3.3 are needed to match the outer and inner results presented in Sections 3.1 and 3.2. Fourth, the zeroth-order solution after applying the matching conditions to the inner and outer results is shown in Section 3.4. Next, to verify our solutions, comparisons of our results with experimental data are described in Section 4. Then, the boundary conditions at the free-fluid/porous-medium interface for this problem are illustrated in Section 5. Finally, the conclusion is drawn in the last section.


## 2. Mathematical Model and Boundary Conditions

We provide the governing equations, the Stokes–Brinkman equations, and their normalized models, in a one-dimensional domain and boundary conditions in this section. Because we consider the movement of the carpet of cilia (not just the movement of a single cilium), the equations used are in a macroscopic scale obtained from the Hybrid Mixture Theory (HMT), which is an upscaling method of changing a microscale equation to be a macroscale equation [24].


### 2.1. Stokes–Brinkman Equations

In this section, we provide the Stokes–Brinkman equations derived in [7] in a macroscopic scale in the one-dimensional domain ω. That is,(1)2μεld2εlvldy2−μk11−1εlvl+ρg=dPdx−μk11−1εlvs,where μ is the dynamic viscosity, k11−1 is the inverse of the permeability tensor, εl is the porosity, vl and vS are the velocities of the liquid and solid phases, respectively, p is pressure, g is gravity, and ρ is the fluid density. We assume that the fluid flows along the x-axis. With the macroscopic scale, we project the average variables on the y-axis. That is, other variables depend on only the y-direction and in the y-direction the pressure gradient is zero and dP/dx is a constant. The Stokes equation is equation (1) without the terms consisting of the permeability k11 and εl=1. That is,(2)2μd2vldy2+ρg=dPdx,and the Brinkman equation is equation (1). Therefore (1) is called Stokes–Brinkman equations. The Brinkman equation used in ω2 has been normalized in [11], which is(3)2ℓ1d2U+dy+2−ℓ2U++ℓ3=dP+dx+−By+,and the normalized Stokes equation is(4)2d2U+dy+2+ℓ3=dP+dx+,where ℓ1=1/εl,ℓ2=h2/K0K+, ℓ3=h2ρg+g0/μU0, h is the characteristic length, U0 is a volumetric average velocity in the porous medium, K0 is characteristic permeability, g0 is characteristic gravity, and By+=h2εlv+/K0K+ where the normalized variables(5)y+=yh,U+=εlvlU0,v+=vsU0,x+=xh,K+=k11K0,P+=PhU0μandg+=gg0.For the Stokes equation, the variableU+=vl/U0 because εl=1. In this study, we assume that the velocity of cilia is linear, that is,(6)By+=h2εlk11a1y++a0,where a0 and a1 are constants. Next, we provide the boundary conditions at the top and the bottom of the domains used in this research.


### 2.2. Boundary Conditions

Because cilia are immobilized at the bottom of the domain, we assume that the PCL fluid has a no-slip condition aty=0. That is,(7)U+=0aty+=0.For the boundary condition at the top of the domain, we assume that the velocity is bounded by a constant. That is,(8)limy+⟶1U+y+=D0,where D0 is a constant as illustrated in Figure 3.Figure 3
Boundary conditions used in this study.For the special cases of the boundary condition, Vanaki et al. [25] used(9)U+=0aty+=0,∂U+∂y=0aty+=1,and U+ is periodic on the left and right sides of the boundaries for the study of the movement of the PCL fluid in case of diseased and healthy conditions. Liron and Meyer [26] studied the PCL flow at y+=1 by using the boundary condition ∂U+/∂y=0. In this study, we use bounded boundary at y+=1 so that one may apply this study to other related problems.We now have the normalized Stokes–Brinkman equations and boundary conditions. Next, we find the solutionU+ by using the method of matched asymptotic expansion in Section 3.


## 2.1. Stokes–Brinkman Equations

In this section, we provide the Stokes–Brinkman equations derived in [7] in a macroscopic scale in the one-dimensional domain ω. That is,(1)2μεld2εlvldy2−μk11−1εlvl+ρg=dPdx−μk11−1εlvs,where μ is the dynamic viscosity, k11−1 is the inverse of the permeability tensor, εl is the porosity, vl and vS are the velocities of the liquid and solid phases, respectively, p is pressure, g is gravity, and ρ is the fluid density. We assume that the fluid flows along the x-axis. With the macroscopic scale, we project the average variables on the y-axis. That is, other variables depend on only the y-direction and in the y-direction the pressure gradient is zero and dP/dx is a constant. The Stokes equation is equation (1) without the terms consisting of the permeability k11 and εl=1. That is,(2)2μd2vldy2+ρg=dPdx,and the Brinkman equation is equation (1). Therefore (1) is called Stokes–Brinkman equations. The Brinkman equation used in ω2 has been normalized in [11], which is(3)2ℓ1d2U+dy+2−ℓ2U++ℓ3=dP+dx+−By+,and the normalized Stokes equation is(4)2d2U+dy+2+ℓ3=dP+dx+,where ℓ1=1/εl,ℓ2=h2/K0K+, ℓ3=h2ρg+g0/μU0, h is the characteristic length, U0 is a volumetric average velocity in the porous medium, K0 is characteristic permeability, g0 is characteristic gravity, and By+=h2εlv+/K0K+ where the normalized variables(5)y+=yh,U+=εlvlU0,v+=vsU0,x+=xh,K+=k11K0,P+=PhU0μandg+=gg0.For the Stokes equation, the variableU+=vl/U0 because εl=1. In this study, we assume that the velocity of cilia is linear, that is,(6)By+=h2εlk11a1y++a0,where a0 and a1 are constants. Next, we provide the boundary conditions at the top and the bottom of the domains used in this research.


## 2.2. Boundary Conditions

Because cilia are immobilized at the bottom of the domain, we assume that the PCL fluid has a no-slip condition aty=0. That is,(7)U+=0aty+=0.For the boundary condition at the top of the domain, we assume that the velocity is bounded by a constant. That is,(8)limy+⟶1U+y+=D0,where D0 is a constant as illustrated in Figure 3.Figure 3
Boundary conditions used in this study.For the special cases of the boundary condition, Vanaki et al. [25] used(9)U+=0aty+=0,∂U+∂y=0aty+=1,and U+ is periodic on the left and right sides of the boundaries for the study of the movement of the PCL fluid in case of diseased and healthy conditions. Liron and Meyer [26] studied the PCL flow at y+=1 by using the boundary condition ∂U+/∂y=0. In this study, we use bounded boundary at y+=1 so that one may apply this study to other related problems.We now have the normalized Stokes–Brinkman equations and boundary conditions. Next, we find the solutionU+ by using the method of matched asymptotic expansion in Section 3.


## 3. Asymptotic Solutions

In this section we employ the method of asymptotic expansion to find the solutions in both domainsω1 and ω2. Because there exists a transition layer at the free-fluid/porous-medium interface, the approximation of this part of the domain is called inner solution. The approximation of other parts of the domain ω is named outer solution. Then the solutions are combined through matching conditions provided in Section 3.3. The inner and outer solutions, derived in [11], are briefly introduced in Sections 3.1 and 3.2. The zeroth-order results after applying the matching conditions are provided in Section 3.4.


### 3.1. Outer Solutions

In this section, we briefly provide the outer solutions found in [11] in two different domains, ω1 and ω2, by using the method of asymptotic expansion. The solutions in domain ω2 and ω1 are obtained by substituting(10)U+=U+0+εU+1+ε2U+2+Oε3,into (3) and (4), respectively. Thus, we have the zeroth-order and first-order solutions for the outer region ω2 with the boundary condition (7) are(11)U+0=w1eℓ2/2ℓ1y+−w1e−ℓ2/2ℓ1y++−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0e−ℓ2/2ℓ1y++h2εlℓ2k11a1y++ℓ3ℓ2−1ℓ2dP+dx++h2εlℓ2k11a0,(12)U+1=w5eℓ2/2ℓ1y+−w5e−ℓ2/2ℓ1y+,respectively, where w1 and w5 are constants, and the zeroth-order and first-order solutions for the outer region ω1 with the boundary condition (8) are(13)U+0=1/2−ℓ3/2y+22+w3y++D0−1/2−ℓ3/22−w3,(14)U+1=w7y++D0−w7,respectively, where w3 and w7 are constants. Next, we provide the solution in the inner region.


### 3.2. Inner Solutions

In this section, we provide the inner solutions in the transition zone between the porous medium and the free-fluid domain. In the inner region, for any variablef, we use f¯ instead of a physical variable f. For example, y+ in the inner region is y¯+ and U+ in the inner region is U¯+. Define(15)y¯+=y+−yStokeεandU¯+y¯+=U+y+−yStokeε.Applying (15) to (3) with the properties, obtained from (15),(16)dy¯+dy+=1εandd2U+dy+2=1ε2d2U¯+dy¯+2,and using the same process of finding the outer solution, we have the zeroth-, first-, and second-order solutions as follows:(17)U¯+0=m1y¯++m2,(18)U¯+1=m3y¯++m4,(19)U¯+2=ℓ212ℓ1m1y¯+3+ℓ2m22ℓ1−ℓ32ℓ1+12ℓ1dP+dx+−h2εl2ℓ1k11a1yStoke−h2εl2ℓ1k11a0y¯+22+m5y¯++m6,where m1,m2,m3, and m4 are constants. The outer and inner solutions will be matched by using the matching conditions provided in the next section.


### 3.3. Matching Conditions

In this section, we derive matching conditions to match the variables defined in the outer and inner regions. The matching condition used in this work is developed from Van Dyke’s rule [9, 27]. With the outer and inner variables,(20)fy+,ε=f0y++εf1y++οε2.(21)f¯y¯+,ε=f¯0y¯++εf¯1y¯++οε2,and applying Van Dyke’s rule matching principle, we have(22)limy¯+⟶±∞f¯y¯+=limy+⟶yStoke±fy+.We find the matching condition off0 by substituting (20) and (21) into (22). Therefore,(23)limy+⟶yStoke±f0y++εf1y++Oε2=limy¯+⟶±∞f¯0y¯++εf¯1y¯++Oε2.We suppose thatlimy+⟶yStoke±fiy+ and limy¯+⟶±∞fiy+ exist, i=0,1,2,... Then, considering the zeroth-order term of ε, we have (24)limy¯+⟶±∞f¯0y¯+=limy+⟶yStoke±f0y+.Next, we find another matching condition by applying Taylor’s series expansion to (20) around the point yStoke+ or yStoke−, where the superscript plus and minus signs mean the right and left sides of yStoke and we assume that fiy+,i=0,1,2,3,… are continuous at y=yStoke so that(25)fiyStoke+=limy+⟶yStroke±fiy+,i=0,1,2,3,….Then, we have(26)fy+,ε=limy+⟶yStoke±f0y++y+−yStoke±limy+⟶yStoke±df0y+dy++y+−yStoke±22limy+⟶yStoke±d2f0y+dy++Οy+−yStoke±3+εlimy+⟶yStoke±f1y++y+−yStoke±limy+⟶yStroke±df1y+dy++Οy+−yStoke±2+Οε2.Substituting (15) into (26), we obtain(27)fy+,ε=limy+⟶yStoke±f0y++εy¯+limy+⟶yStoke±df0y+dy++εy¯+22limy+⟶yStoke±d2f0y+dy++Οεy¯+3+εlimy+⟶yStoke±f1y++εy¯+limy+⟶yStoke±df1y+dy++Οεy¯+2+Οε2.Differentiating (27) with respect to y¯+, taking y+⟶yStoke± on both sides, and omitting the big O, then, we obtain(28)limy+⟶yStoke±dfy+,εdy¯+=εlimy+⟶yStoke±df0y+dy++ε2y¯+limy+⟶yStoke±d2f0y+dy++ε2limy+⟶yStoke±df1y+dy+,where in this work we assume that the function f is continuous enough so that the differentiation can be switched with the limit. To find the matching condition for the first derivative of the functions f and f¯, we differentiate (21) with respect to y¯+ and taking limit y¯+⟶±∞ on both sides, we obtain(29)limy¯+⟶±∞df¯y¯+,εdy¯+=limy¯+⟶±∞df¯0y¯+dy¯++εdf¯1y¯+dy¯++ε2df¯2y¯+dy¯+.Since we assume that the differentiation can be switched with the limit, from (22), we have(30)limy+⟶yStoke±dfy+,εdy¯+=limy¯+⟶±∞df¯y¯+,εdy¯+.This fact demonstrates that (28) and (29) are equal. Then, we get(31)εlimy+⟶yStoke±df0y+dy++ε2y¯+limy+⟶yStoke±d2f0y+dy++ε2limy+⟶yStoke±df1y+dy+=limy¯+⟶±∞df¯0y¯+dy¯++εdf¯1y¯+dy¯++ε2df¯2y¯+dy¯+.Collecting the same order ofε in (31), we obtain(32)ε0:limy¯+⟶±∞df¯0y¯+dy¯+=0,(33)ε1:limy¯+⟶±∞df¯1y¯+dy¯+=limy+⟶yStoke±df0y+dy+,(34)ε2:limy¯+⟶±∞df¯2y¯+dy¯+=y¯+limy+⟶yStoke±d2f0y+dy++limy+⟶yStoke±df1y+dy+.Next, we determine the zeroth-order matching solutions for the couple Stokes–Brinkman equations.


### 3.4. The Zeroth-Order Solution

In this section, we couple the zeroth-order solutions of the governing equations (3), (4) by using the matching conditions provided in Section 3.3. Before using the matching conditions, we differentiate the zeroth-order and first-order inner solutions (17) and (18) with respect to y¯+ and taking the limit y¯+⟶+∞ on both sides, we obtain(35)limy¯+⟶+∞dU¯+0dy¯+=m1,(36)limy¯+⟶+∞dU¯+1dy¯+=m3,respectively. According to the matching condition (32), we have m1=0. Then, from (17), we obtain(37)U¯+0=m2.Next, we provide the matching process of the outer solutions in domainsω2 and ω1. Differentiating the zeroth-order outer solutions (11) and (13) for the outer regions ω2 and ω1, respectively, with respect to y+, and taking the limit y+⟶yStoke± on both sides, we obtain(38)limy+⟶yStoke+dU+0dy+=w1ℓ22ℓ1eℓ2/2ℓ1yStoke++w1ℓ22ℓ1e−ℓ2/2ℓ1yStoke+−ℓ22ℓ1−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0e−ℓ2/2ℓ1yStoke++h2εlℓ2k11a1,(39)limy+⟶yStoke+dU+0dy+=12−ℓ32yStoke++w3,respectively. Substituting (36) and (38) into (33), we obtain(40)m3=w1ℓ22ℓ1eℓ2/2ℓ1yStoke++w1ℓ22ℓ1e−ℓ2/2ℓ1yStoke+−ℓ22ℓ1−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0e−ℓ2/2ℓ1yStoke++h2εlℓ2k11a1,in region ω2. Similarly, substituting (36) and (39) into (33), we have(41)m3=12−ℓ32yStoke++w3,in region ω1. Since (40) is equal to (41), we then have(42)w3=−12−ℓ32yStoke++w1ℓ22ℓ1eℓ2/2ℓ1yStoke++ℓ22ℓ1e−ℓ2/2ℓ1yStoke+−ℓ22ℓ1−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0e−ℓ2/2ℓ1yStoke++h2εlℓ2k11a1.Define(43)z1=ℓ22ℓ1,z2=12−ℓ32,z3=12ℓ1dP+dx+−ℓ32ℓ1,z4=−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0andz5=h2εlℓ2k11a1,which are known constants. Applying (43) to (42) and (11), this yields(44)w3=−z2yStoke++w1z1ez1yStoke++z1e−z1yStoke+−z1z4e−z1yStoke++z5,(45)U+0p=w1ez1y+−w1e−z1y++z4e−z1y++z5y+−z4;0<y+<yStoke,respectively, where U+0p means the zeroth-order solution in the porous-medium domain ω2. Substituting (43) and (44) into (13), we have(46)U+0f=z2y+22+D0−z22+−z2yStoke++w1z1ez1yStoke++z1e−z1yStoke+−z1z4e−z1yStoke++z5y+−1;yStoke<y+<1,where U+0f means the zeroth-order solution in the free-fluid domain ω1. Now, we obtain U+0p and U+0f in terms of w1.To illustrate our results, the solution (45) is plotted in Figure 4, where the angles between the cilia and the horizontal plane considered in this work are 40∘,50∘,60∘,70∘,80∘, and 90∘, respectively. In this research, we assume that the cilia have the highest velocity at θ=90∘ and stop beating at θ=40∘. The values of the variables used in (45) are presented in Tables 1 and 2. The values in Table 1 are used for all angles θ and the values of the variables in Table 2 are different for each angle θ, which are obtained from [8, 28], where they have obtained the data from a biological experiment [29]. We concentrate on only one viscosity because we intend to find the solutions that are close to the real problem and the solutions do not change with the different viscosities. The results for the different angles are interesting because we determine the PCL velocities due to the movement of cilia, which means the angle changes. The Earth gravity, solid velocity, and permeability variables in (5) are fixed according to our designed problem. Pressure can be modified but for this problem the most effect on the velocity of the PCL fluid is the solid phases, rather than pressure gradient. The variables x and y in (5) are independent variables, which give the position of the fluid velocity on the domain.Figure 4
The asymptotic solutionU+0, equation (45).Table 1
Values of the variables applied to the solution (45) to plot Figure 4.
VariableValueUnitVariableValueUnita001U01μm/sh7.5μmK01μm2ρ992.2×10−15g/μm3dP+/dx+1⌈1⌉μ3×10−6g/μm⋅sw10⌈1⌉g9.81×106μm/s2Table 2
Values of the variables used to plot Figure4 for each angle θ.
Variable50∘60∘70∘80∘90∘yStoke0.76720.86380.93530.98181.0000K+0.00120.00150.00160.00170.0018εl0.67170.70990.73310.74390.7487a11.21.41.61.82Figure4 shows the velocity of the PCL fluid with different angles θ in domain ω2. The highest velocity occurred at θ=90∘ as expected. This is because, in this study, we assume that the maximum velocity of cilia occurs at θ=90∘. Notice that the velocities decrease when the angles between the cilia and the horizontal plane decrease.In the next section, we compare our solution with an experimental data where the velocity of the solid phases in our solution is set to be zero.


## 3.1. Outer Solutions

In this section, we briefly provide the outer solutions found in [11] in two different domains, ω1 and ω2, by using the method of asymptotic expansion. The solutions in domain ω2 and ω1 are obtained by substituting(10)U+=U+0+εU+1+ε2U+2+Oε3,into (3) and (4), respectively. Thus, we have the zeroth-order and first-order solutions for the outer region ω2 with the boundary condition (7) are(11)U+0=w1eℓ2/2ℓ1y+−w1e−ℓ2/2ℓ1y++−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0e−ℓ2/2ℓ1y++h2εlℓ2k11a1y++ℓ3ℓ2−1ℓ2dP+dx++h2εlℓ2k11a0,(12)U+1=w5eℓ2/2ℓ1y+−w5e−ℓ2/2ℓ1y+,respectively, where w1 and w5 are constants, and the zeroth-order and first-order solutions for the outer region ω1 with the boundary condition (8) are(13)U+0=1/2−ℓ3/2y+22+w3y++D0−1/2−ℓ3/22−w3,(14)U+1=w7y++D0−w7,respectively, where w3 and w7 are constants. Next, we provide the solution in the inner region.


## 3.2. Inner Solutions

In this section, we provide the inner solutions in the transition zone between the porous medium and the free-fluid domain. In the inner region, for any variablef, we use f¯ instead of a physical variable f. For example, y+ in the inner region is y¯+ and U+ in the inner region is U¯+. Define(15)y¯+=y+−yStokeεandU¯+y¯+=U+y+−yStokeε.Applying (15) to (3) with the properties, obtained from (15),(16)dy¯+dy+=1εandd2U+dy+2=1ε2d2U¯+dy¯+2,and using the same process of finding the outer solution, we have the zeroth-, first-, and second-order solutions as follows:(17)U¯+0=m1y¯++m2,(18)U¯+1=m3y¯++m4,(19)U¯+2=ℓ212ℓ1m1y¯+3+ℓ2m22ℓ1−ℓ32ℓ1+12ℓ1dP+dx+−h2εl2ℓ1k11a1yStoke−h2εl2ℓ1k11a0y¯+22+m5y¯++m6,where m1,m2,m3, and m4 are constants. The outer and inner solutions will be matched by using the matching conditions provided in the next section.


## 3.3. Matching Conditions

In this section, we derive matching conditions to match the variables defined in the outer and inner regions. The matching condition used in this work is developed from Van Dyke’s rule [9, 27]. With the outer and inner variables,(20)fy+,ε=f0y++εf1y++οε2.(21)f¯y¯+,ε=f¯0y¯++εf¯1y¯++οε2,and applying Van Dyke’s rule matching principle, we have(22)limy¯+⟶±∞f¯y¯+=limy+⟶yStoke±fy+.We find the matching condition off0 by substituting (20) and (21) into (22). Therefore,(23)limy+⟶yStoke±f0y++εf1y++Oε2=limy¯+⟶±∞f¯0y¯++εf¯1y¯++Oε2.We suppose thatlimy+⟶yStoke±fiy+ and limy¯+⟶±∞fiy+ exist, i=0,1,2,... Then, considering the zeroth-order term of ε, we have (24)limy¯+⟶±∞f¯0y¯+=limy+⟶yStoke±f0y+.Next, we find another matching condition by applying Taylor’s series expansion to (20) around the point yStoke+ or yStoke−, where the superscript plus and minus signs mean the right and left sides of yStoke and we assume that fiy+,i=0,1,2,3,… are continuous at y=yStoke so that(25)fiyStoke+=limy+⟶yStroke±fiy+,i=0,1,2,3,….Then, we have(26)fy+,ε=limy+⟶yStoke±f0y++y+−yStoke±limy+⟶yStoke±df0y+dy++y+−yStoke±22limy+⟶yStoke±d2f0y+dy++Οy+−yStoke±3+εlimy+⟶yStoke±f1y++y+−yStoke±limy+⟶yStroke±df1y+dy++Οy+−yStoke±2+Οε2.Substituting (15) into (26), we obtain(27)fy+,ε=limy+⟶yStoke±f0y++εy¯+limy+⟶yStoke±df0y+dy++εy¯+22limy+⟶yStoke±d2f0y+dy++Οεy¯+3+εlimy+⟶yStoke±f1y++εy¯+limy+⟶yStoke±df1y+dy++Οεy¯+2+Οε2.Differentiating (27) with respect to y¯+, taking y+⟶yStoke± on both sides, and omitting the big O, then, we obtain(28)limy+⟶yStoke±dfy+,εdy¯+=εlimy+⟶yStoke±df0y+dy++ε2y¯+limy+⟶yStoke±d2f0y+dy++ε2limy+⟶yStoke±df1y+dy+,where in this work we assume that the function f is continuous enough so that the differentiation can be switched with the limit. To find the matching condition for the first derivative of the functions f and f¯, we differentiate (21) with respect to y¯+ and taking limit y¯+⟶±∞ on both sides, we obtain(29)limy¯+⟶±∞df¯y¯+,εdy¯+=limy¯+⟶±∞df¯0y¯+dy¯++εdf¯1y¯+dy¯++ε2df¯2y¯+dy¯+.Since we assume that the differentiation can be switched with the limit, from (22), we have(30)limy+⟶yStoke±dfy+,εdy¯+=limy¯+⟶±∞df¯y¯+,εdy¯+.This fact demonstrates that (28) and (29) are equal. Then, we get(31)εlimy+⟶yStoke±df0y+dy++ε2y¯+limy+⟶yStoke±d2f0y+dy++ε2limy+⟶yStoke±df1y+dy+=limy¯+⟶±∞df¯0y¯+dy¯++εdf¯1y¯+dy¯++ε2df¯2y¯+dy¯+.Collecting the same order ofε in (31), we obtain(32)ε0:limy¯+⟶±∞df¯0y¯+dy¯+=0,(33)ε1:limy¯+⟶±∞df¯1y¯+dy¯+=limy+⟶yStoke±df0y+dy+,(34)ε2:limy¯+⟶±∞df¯2y¯+dy¯+=y¯+limy+⟶yStoke±d2f0y+dy++limy+⟶yStoke±df1y+dy+.Next, we determine the zeroth-order matching solutions for the couple Stokes–Brinkman equations.


## 3.4. The Zeroth-Order Solution

In this section, we couple the zeroth-order solutions of the governing equations (3), (4) by using the matching conditions provided in Section 3.3. Before using the matching conditions, we differentiate the zeroth-order and first-order inner solutions (17) and (18) with respect to y¯+ and taking the limit y¯+⟶+∞ on both sides, we obtain(35)limy¯+⟶+∞dU¯+0dy¯+=m1,(36)limy¯+⟶+∞dU¯+1dy¯+=m3,respectively. According to the matching condition (32), we have m1=0. Then, from (17), we obtain(37)U¯+0=m2.Next, we provide the matching process of the outer solutions in domainsω2 and ω1. Differentiating the zeroth-order outer solutions (11) and (13) for the outer regions ω2 and ω1, respectively, with respect to y+, and taking the limit y+⟶yStoke± on both sides, we obtain(38)limy+⟶yStoke+dU+0dy+=w1ℓ22ℓ1eℓ2/2ℓ1yStoke++w1ℓ22ℓ1e−ℓ2/2ℓ1yStoke+−ℓ22ℓ1−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0e−ℓ2/2ℓ1yStoke++h2εlℓ2k11a1,(39)limy+⟶yStoke+dU+0dy+=12−ℓ32yStoke++w3,respectively. Substituting (36) and (38) into (33), we obtain(40)m3=w1ℓ22ℓ1eℓ2/2ℓ1yStoke++w1ℓ22ℓ1e−ℓ2/2ℓ1yStoke+−ℓ22ℓ1−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0e−ℓ2/2ℓ1yStoke++h2εlℓ2k11a1,in region ω2. Similarly, substituting (36) and (39) into (33), we have(41)m3=12−ℓ32yStoke++w3,in region ω1. Since (40) is equal to (41), we then have(42)w3=−12−ℓ32yStoke++w1ℓ22ℓ1eℓ2/2ℓ1yStoke++ℓ22ℓ1e−ℓ2/2ℓ1yStoke+−ℓ22ℓ1−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0e−ℓ2/2ℓ1yStoke++h2εlℓ2k11a1.Define(43)z1=ℓ22ℓ1,z2=12−ℓ32,z3=12ℓ1dP+dx+−ℓ32ℓ1,z4=−ℓ3ℓ2+1ℓ2dP+dx+−h2εlℓ2k11a0andz5=h2εlℓ2k11a1,which are known constants. Applying (43) to (42) and (11), this yields(44)w3=−z2yStoke++w1z1ez1yStoke++z1e−z1yStoke+−z1z4e−z1yStoke++z5,(45)U+0p=w1ez1y+−w1e−z1y++z4e−z1y++z5y+−z4;0<y+<yStoke,respectively, where U+0p means the zeroth-order solution in the porous-medium domain ω2. Substituting (43) and (44) into (13), we have(46)U+0f=z2y+22+D0−z22+−z2yStoke++w1z1ez1yStoke++z1e−z1yStoke+−z1z4e−z1yStoke++z5y+−1;yStoke<y+<1,where U+0f means the zeroth-order solution in the free-fluid domain ω1. Now, we obtain U+0p and U+0f in terms of w1.To illustrate our results, the solution (45) is plotted in Figure 4, where the angles between the cilia and the horizontal plane considered in this work are 40∘,50∘,60∘,70∘,80∘, and 90∘, respectively. In this research, we assume that the cilia have the highest velocity at θ=90∘ and stop beating at θ=40∘. The values of the variables used in (45) are presented in Tables 1 and 2. The values in Table 1 are used for all angles θ and the values of the variables in Table 2 are different for each angle θ, which are obtained from [8, 28], where they have obtained the data from a biological experiment [29]. We concentrate on only one viscosity because we intend to find the solutions that are close to the real problem and the solutions do not change with the different viscosities. The results for the different angles are interesting because we determine the PCL velocities due to the movement of cilia, which means the angle changes. The Earth gravity, solid velocity, and permeability variables in (5) are fixed according to our designed problem. Pressure can be modified but for this problem the most effect on the velocity of the PCL fluid is the solid phases, rather than pressure gradient. The variables x and y in (5) are independent variables, which give the position of the fluid velocity on the domain.Figure 4
The asymptotic solutionU+0, equation (45).Table 1
Values of the variables applied to the solution (45) to plot Figure 4.
VariableValueUnitVariableValueUnita001U01μm/sh7.5μmK01μm2ρ992.2×10−15g/μm3dP+/dx+1⌈1⌉μ3×10−6g/μm⋅sw10⌈1⌉g9.81×106μm/s2Table 2
Values of the variables used to plot Figure4 for each angle θ.
Variable50∘60∘70∘80∘90∘yStoke0.76720.86380.93530.98181.0000K+0.00120.00150.00160.00170.0018εl0.67170.70990.73310.74390.7487a11.21.41.61.82Figure4 shows the velocity of the PCL fluid with different angles θ in domain ω2. The highest velocity occurred at θ=90∘ as expected. This is because, in this study, we assume that the maximum velocity of cilia occurs at θ=90∘. Notice that the velocities decrease when the angles between the cilia and the horizontal plane decrease.In the next section, we compare our solution with an experimental data where the velocity of the solid phases in our solution is set to be zero.


## 4. Comparison with an Experimental Study

In this section, we compare the mass flow rateΦ of our result with the study of Beavers and Joseph [22], where(47)Φ=Mf−MpMp,where Mf is the mass flow rate in the free-fluid region and Mp is the mass flow rate in the porous medium. In their study, BJ carried out the experiment to measure the flow rate of a long porous block and a small uniform gap immediately above the block, while in this study, we analytically find the mass flow rate of the fluid flow in the PCL in human lung by using the asymptotic expansion method. The porous medium in the experimental study is stationary porous blocks, while the porous medium in this study is the ciliary layer, which is a motion medium. Therefore, we compare our solution with the experimental data at only one porosity and the solid velocity is set to be zero. Since Beavers and Joseph provide the data, the mass flow rate, only on the free-fluid domain, as the result, only the solution on the free-fluid domain is compared.By applying the zeroth-order solutions (45), (46) provided in Section 3.4 to (47), the mass flow rate obtained from both equations is(48)Φ=1/21/2−σ2K0ρg+g0/2μU0σ2K0+D0−1/41−σ2K0ρg+g0/μU0+−1/2−σ2K0ρg+g0/2μU0yStoke++w1σ1/2K+eσ1/2K+yStoke++σ1/2K+e−σ1/2K+yStoke++σ1/2K+ρg+g0K0K+/μU0+h2/σ2K0e−σ1/2K+yStoke+σK0−1w1eσ2εlK0/2K+−w1e−σ2εlK0/2K+−ρg+g0K0K+/μU0+h2/σ2K0e−σ2εlK0/2K++ρg+g0K0K+/μU0+h2/σ2K0−1,where σ=y+/K0 and yStoke≤y+≤1. To compare the result with the experimental data, we use the same values of the variables as in the experiment. The dynamic viscosity and the density of water are 8.9×10−4Pa⋅s and 997kg/m3, respectively [27]. The gravity g+g0 is 9.8m/s2. The length of the free-fluid domain ranges from 3.37046×10−4 to 1, where yStoke is assumed to be 3.37046×10−4. The porosity is 0.78 and the characteristic permeability K0 is 7.1×10−9. The values of other variables used in (48) are K+=k11/K0,k11=0.4,w1=−2.5,U0=1,D0=1 and h=7.5, and Φ is normalized by 104.Figure5 shows the comparison of the mass flow rate in this study (solid line) through the free-fluid domains as a function of σ, with the mass flow rate in the experimental study of Beavers and Joseph [22]. As can be seen from Figure 5, the solution of our study (solid line) is in a good agreement with the experimental data of Beavers and Joseph [22].Figure 5
The mass flow rate in this study and the experimental data [22].


## 5. Boundary Conditions at the Free-Fluid/Porous-Medium Interface

In this section we find the boundary condition at the free-fluid/porous-medium interface by using the two solutions in the two domainsω1 and ω2. Based on solutions (45) and (46), we can now obtain the boundary condition at the free-fluid/porous-medium interface, which is(49)U+0f=z2y+22+D0−z22+−z2yStoke+−z1z4e−z1yStoke++z5y+−1+1ez1y+−e−z1y+U+0p−z4e−z1y+−1−z5y+z1ez1yStoke++z1e−z1yStoke+y+−1.Next, we provide another boundary condition at the free-fluid/porous-medium interface.Differentiating the zeroth-order solutions (45) and (46) with respect to y+ both sides, we have(50)dU+0Pdy+=w1z1ez1y++w1z1e−z1y+−z1z4e−z1y++z5,(51)dU+0fdy+=z2y+−z2yStoke++w1z1ez1yStoke++z1e−z1yStoke+−z1z4e−z1yStoke++z5,respectively. Therefore we obtain another boundary condition at the free-fluid/porous-medium interface, which is(52)dU+0fdy+=z2y+−yStoke++ez1yStoke++e−z1yStoke+ez1y++e−z1y+dU+0Pdy++z1z4e−z1y+−z5−z1z4e−z1yStoke++z5.Thus, we obtain two boundary conditions at the free-fluid/porous-medium interface, which are proper for the fluid flow in the PCL due to the movement of cilia.


## 6. Summary and Conclusions

The movement of the cilia in the PCL layer plays an important role in the clearance of mucus and in preventing foreign materials from entering the lower respiratory tract. Determining the velocity of the PCL fluid flow in motile cilia is necessary for a better understanding of the clearance mechanism. The velocity of PCL fluid can be determined using an ordinary differential equation (ODE), which requires boundary conditions. This research quantifies the boundary conditions at the interface between the porous-medium and free-fluid regions using Stokes–Brinkman equations. To achieve our aim:(i)
We apply the Stokes equation to the free-fluid domain and the Brinkman equation to the porous-medium region.(ii)
We also assumed two boundary conditions: The first boundary condition is no-slip boundary condition, applied at the bottom of the porous-medium domain; the second boundary condition is that the velocity of the PCL fluid is bounded wheny tends to a constant, applied at the top of the free-fluid region.(iii)
The velocity of the PCL fluid due to the movement of cilia is obtained by using the matched asymptotic expansion.(iv)
The results are illustrated in Figure4 with different angles, θ=40∘,50∘,60∘,70∘,80∘,90∘, that the cilia make with the horizontal plane. The variables employed to plot the graph of the solutions are shown in Tables 1 and 2.(v)
We find that the velocity of the PCL fluid reaches its maximum value whenθ=90∘ and decreases when θ decreases, which satisfies our assumption that the highest velocity of cilia occurs when the cilia are perpendicular to the horizontal plane and the cilia velocities decrease with decreasing angle θ.(vi)
We also found that the solution in this study is consistent with the experimental data studied by Beavers and Joseph, which is shown in Figure5.As mentioned above, not only are there two domains, porous domain and free-fluid domain, but also there is the interface between the free-fluid domain and porous medium. The boundary conditions at the interface between the porous-medium and the free-fluid regions are determined from the zeroth-order solutions provided in Section3.4.(i)
The first boundary condition is the relationship between the two zeroth-ordered solutions in the two domains while the second boundary condition is the derivative of those two solutions.(ii)
It should be noted that, in this study, the solid phases in the porous domain are self-propelled motion and we provide proper boundary conditions for the problem. On the contrary, in previous studies, they are stationary. Although they are moved, available literature has not provided appropriate boundary conditions for the fluid flow in PCL.(iii)
The limitation of this study is that the asymptotic expansion method applied in this study is good for a linear equation, but it is not easy to find the solutions for nonlinear equations.(iv)
The disadvantage of this study is that we calculate the zeroth-order approximation solution, where the other orders are cut off.(v)
In future work, the more realistic solution such as every angle is considered in one calculation rather than one angle. Applications of the boundary conditions to some other problems are, for example, flows in rice field [30–34].



---

*Source: 1390486-2022-03-08.xml*