# Determination of the Slip Constant in the Beavers-Joseph Experiment for Laminar Fluid Flow through Porous Media Using a Meshless Method

Authors: Magdalena Mierzwiczak; Agnieszka Fraska; Jakub K. Grabski
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In the paper, a numerical simulation of the Beavers-Joseph experiment is presented. The simulation is based on a meshless method which is called the Trefftz method with the special purpose Trefftz functions. The porous medium is modeled as a regular array of fibers. Three kinds of arrays are taken under consideration: triangular, square, and hexagonal. Firstly, the permeability of the fibrous porous medium is determined by consideration of flow between an infinite array of fibers. In the next step the considered region is divided into two parts. The first half of the channel is a porous medium modeled by a regular array of fibers and the second part which is a free fluid flow region. In the paper, the slip constant existing in the Beavers-Joseph boundary condition is determined in terms of the volume fraction of fibers.


---

## Body

## 1. Introduction

The viscous fluid flow in a channel coupled with the flow through an adjacent porous region is often found in several engineering applications including composite materials processing [1], multifilament spinning [2], contaminant transport in porous media [3, 4], porous chemical reactors [5, 6], flow between porous rollers [7], and porous bearings [8]. A set of proper governing equations for the flow through two distinct, neighboring, and interconnected domains has to be established to solve this type of problems numerically. The incompressible flow in the free fluid region is well represented by the Navier-Stokes equations, while the Darcy equation is used to describe the flow in the adjoining porous medium. A problem is considered with a special boundary condition at the interface between these two different regions.Beavers and Joseph (1967) [9] proposed that the interfacial velocities of the freely flowing fluid and the velocity of the fluid flowing through the porous medium can be related by the following interface boundary condition:(1)dudn=αku-q,whereu is the velocity of the fluid flow calculated at boundary (the slip velocity),kis the permeability of the porous medium,nis the normal direction to the boundary, andq is the seepage velocity measured at small distance outside the interface, suggesting the existence of a thin layer just inside the porous medium over which the velocity transition occurs. The dimensionless slip wall coefficient α is an independent of the fluid viscosity and apart from the permeability depends on the structural parameters of the porous medium and it is specific to the geometric features of the interface.For determination of the slip coefficient, Beavers and Joseph described an experiment performed in a parallel-plate channel; one of the bounding walls was made of the porous material while the other one was impermeable. Identical axial pressure gradients were imposed on the channel and the porous medium, thereby giving rise to parallel axial flows.This paper presents a numerical simulation of the Beavers-Joseph’s experiment for the fibrous porous media. The longitudinal laminar flow (the Poiseuille flow) in a parallel-plate conduit is considered. The first half of the examined region is a porous medium and the second one is a region of the free flowing fluid. The porous medium is modeled as a bundle of parallel fibres arranged in a triangular, square, and hexagonal array. The purpose of the present consideration is the determination of the slip constant in the Beavers-Joseph boundary condition for the Newtonian fluid. However, to do this the permeability of the porous medium is required. To determine this, the flow with the pressure gradient in an unbounded porous medium was considered. Essential novelty of this paper is the numerical simulation of the Beavers-Joseph experiment for the Newtonian fluid. Numerical simulations are conducted using the Trefftz method. The influence of the volume fraction of the fibers on the value of the slip constant in the Beavers-Joseph boundary condition is investigated.


## 2. Modelling of the Fibrous Porous Medium

In this paper the porous medium modeled by a regular array of parallel fibers is considered. Square, triangular, and hexagonal arrays of fibers (see Figure1) are analyzed. In Figure 1 one can find the characteristic geometrical features of the regions, which are the diameter of the fibers 2a and the distance between the fibers 2b.Figure 1
A model of the fibrous porous medium arranged in (a) triangular, (b) square, or (c) hexagonal arrays of fibers.
(a) (b) (c)Various authors considered longitudinal flow with respect to the fibers [1, 6, 7, 10–21]. In the present paper the same problem is considered but the solution is obtained using the special purpose Trefftz functions [16, 18]. The method is semianalytical. That means that application of the method gives analytical form of the dimensionless permeability of the porous medium. Only the unknown coefficients of the solution are obtained numerically.Let us introduce the following dimensionless variables:(2)X=xb,Y=yb,R=rb,E=ab,where X,Y are the dimensionless Cartesian coordinates, R,θ are the dimensionless polar coordinates, and E denotes the dimensionless radius of a single fiber.In the case of the unbounded array of fibers and longitudinal laminar flow, the problem is symmetrical with respect to the dashed lines, depicted in Figure1. Thus, it can be considered in a single repeated element of the domain. Furthermore, it can be presented in the same way for each array (triangular, square, or hexagonal) introducing an additional parameter L. The values of the parameter are L=6, L=4, and L=3 for triangular, square, and hexagonal arrays, respectively. The repeated element of the unbounded porous medium with the dimensionless variables is presented in Figure 2.Figure 2
Repeated element of an unbounded porous medium.The volume fraction of the fiber is given by(3)φ=ΩPΩT=π·E2L·tan⁡π/L,where ΩT = 0.5  tan⁡(π/L) is the total area of the repeated element (the fluid and the fiber areas together) and ΩP = 0.5E2∙π/L is the area of the fiber in the repeated element (Figure 2).After determination of the permeability of the porous medium, microstructural flow in the layer of porous medium modeled by a bundle of regular arranged (triangular, square, or hexagonal) arrays of fibers can be considered (Figure3).Figure 3
Porous medium with very high porosity in the presence of free fluid region between two impermeable walls.Since the array of fibers is streaked and periodic in one direction, it is sufficient to consider the problem only in one repeated strip which is depicted in Figure4—for square (Figure 4(a)), triangular (Figure 4(b)), or hexagonal (Figure 4(c)) array. In all cases the repeated strip is divided into smaller elements associated with each of the fibers which are called large finite elements in next parts of this paper. The widths of the porous medium and the free flow area are the same and equal toH. In Figure 4LR denotes the number of rows of fibers in the porous medium.Figure 4
The repeating part of the considered channel with fibers arranged with (a) square, (b) triangular, and (c) hexagonal arrays.
(a) (b) (c)In such an approach in which the slip constant is determined by means of the numerical simulation the permeability of the porous medium is calculated independently. Then the slip constant can be determined in the second numerical simulation. Thus two independent numerical experiments are needed. In both numerical experiments the porous medium is modeled by a bundle of regular fibers arranged in a triangular, square, and hexagonal array. The radius of the fibers is equal toa,and the distance between the fibers is equal to 2b.In the first experiment the considered bundle of fibers is infinite in all directions (Figure 1) because of the symmetry that fluid flow can be considered in the repeating part of the region (Figure 2) and the permeability of the porous medium can be determined. In the second experiment a regular bundle of fibers is located between two parallel walls (Figure 4). However, the porous medium is located only in one half of this channel. Fluid flow is considered only in the repeating element of the considered channel (Figure 5).Figure 5
The large finite elements for which the boundary conditions are satisfied by the special purpose Trefftz functions.


## 3. The Trefftz Method

The method was proposed in 1926 by German mathematician Erich Trefftz [22] as an alternative numerical method to the Ritz method. In the Ritz method the approximate solution is assumed in such a form that satisfies exactly the boundary conditions and the governing equation is satisfying in an approximate way. In the Trefftz method the approximate solution satisfies exactly the governing equation, while the boundary conditions are satisfying approximately. Over the years, the method is known under different names and various modifications of this approach were proposed. Some reviews of different names, modifications, and various applications of the method can be found in [23, 24].A particular case of the Trefftz method is the Trefftz method with the special purpose Trefftz functions in which the trial functions are chosen in such a way that they satisfy not only the governing equations but also some of the boundary conditions. It has been developed in some previous papers [25–30].


## 4. Determination of the Permeability of the Considered Porous Medium

One of the basic methods for determination of the permeability of the porous medium is an experiment in which the flow rate is measured (the filtration velocity at the known pressure gradient). Based on this measurement the permeability can be calculated.Let us consider a steady, fully developed, laminar flow of an incompressible viscous fluid driven by the constant pressure in a system of regular parallel fibers. The flow is longitudinal with respect to fibers which are arranged in regular square (Figure1(a)), triangular (Figure 1(b)), and hexagonal (Figure 1(c)) arrays. The radius of the fibers is equal toa,and the distance between the fibers is equal to 2b. In such a case the equation of motion is reduced to a single partial differential equation expressed in the polar coordinate system r,θ by(4)∂2wr,θ∂r2+1r∂wr,θ∂r+1r2∂2wr,θ∂θ2=1μdpdzinΩF,where w is the velocity component in thez-axis direction [m/s], dp/dz is the constant pressure gradient [Pa/m], μ is the viscosity of the fluid [Pa∙s], and ΩF is the fluid domain. The above governing equation is subject to the following boundary conditions:(5)wr,θ=0forr=a,(6)∂wr,θ∂θ=0forθ=0,πL,(7)∂wr,θ∂x=0forx=1.Let us introduce the following dimensionless axial velocity:(8)WR,θ=-wr,θb2/μdp/dz.Furthermore, the problem can be considered in the repeating part presented in Figure 2. Now, the governing Eq. (4) has the dimensionless form(9)∂2WR,θ∂R2+1R∂WR,θ∂R+1R2∂2WR,θ∂θ2=-1,with the boundary conditions(10)WR,θ=0forR=E,(11)∂WR,θ∂θ=0forθ=0,πL,(12)∂WR,θ∂X=0forX=1.The solution of the boundary value problem described by Eq. (9) subject to the boundary conditions (10)-(12) is a sum given by the following equation:(13)WR,θ=WgR,θ+WpR,θ,where the subscripts g andp denote general and particular solutions, respectively. The latter for Eq. (9) can be expressed as(14)WpR,θ=-14R2-E2.Furthermore, the above particular solution is chosen in such a way that it also satisfies exactly the boundary conditions (10)-(11).The general solution should satisfy the Laplace equation in the following form:(15)∂2WgR,θ∂R2+1R∂WgR,θ∂R+1R2∂2WgR,θ∂θ2=0.In the literature, see, e.g., [21, 24]; it is well known that such a general solution is given by(16)WgR,θ=A0+A1θ+A2θln⁡R+A3ln⁡R+∑n=1∞BnRλn+CnR-λncos⁡λnθ+∑n=1∞DnRλn+EnR-λnsin⁡λnθ,where A0,A1,A2,A3,Bn,Cn,Dn,En, and λn  (n=1,2,3,…) are unknown constants. Thus the whole solution WR,θ can be expressed as(17)WR,θ=-14R2-E2+A0+A1θ+A2θln⁡R+A3ln⁡R+∑n=1∞BnRλn+CnR-λncos⁡λnθ+∑n=1∞DnRλn+EnR-λnsin⁡λnθ.In order to ensure satisfaction of the boundary condition (10) for any θ the constants A1 and A2 should be equal to zero. The logarithmic function is equal to zero when the argument is equal to 1. Because of that the constant A0 can be taken as A0=-A3ln⁡E. Then the above equation can be rewritten in the form(18)WR,θ=-14R2-E2+A3ln⁡RE+∑n=1∞BnRλn+CnR-λncos⁡λnθ+∑n=1∞DnRλn+EnR-λnsin⁡λnθ.The derivative of the above solution with respect to θ takes the following form:(19)dWR,θdθ=-∑n=1∞λnBnRλn+CnR-λnsin⁡λnθ+∑n=1∞λnDnRλn+EnR-λncos⁡λnθ.In order to satisfy the boundary condition (11) simultaneously for θ=0 and θ=π/L the constants Dn and En  (n=1,2,3,…) have to be equal to zero. Furthermore, the constants λn  (n=1,2,3,…) for θ=π/L should satisfy the following relation: (20)λnπL=nπforn=1,2,3,…,because for such angles the sinusoidal function equals zero. Thus (21)λn=nLforn=1,2,3,….In such a case, the approximate solution takes the form (22)WR,θ=-14R2-E2+A3ln⁡RE+∑n=1∞BnRnL+CnR-nLcos⁡nLθ.Furthermore, for R=E the constants Bn and Cn  (n=1,2,3,…) should satisfy the relation(23)BnEnL+CnE-nL=0forn=1,2,3,…It leads to the following formula for the constants Cn  (n=1,2,3,…):(24)Cn=-BnE2nLforn=1,2,3,….Thus, finally the exact solution of Eq. (9) can be expressed using the special purpose Trefftz functions in the following form:(25)WR,θ=-14R2-E2+A3ln⁡RE+∑n=1∞BnRnL-E2nLR-nLcos⁡nLθ.From the infinite series in the above form we choose onlyN– 1 first terms and taking A3=BN, the solution can be rewritten in the form(26)WR,θ=-14R2-E2+∑n=1N-1BnRnL-E2nLR-nLcos⁡nLθ+BNln⁡RE.The solution (26) satisfies exactly the governing equation (9) and the boundary conditions (10)–(11). The unknown coefficients Bn  (n= 1,…,N) are determined by solving the system of linear equations resulting from satisfying of the boundary condition (12) using the boundary collocation technique [21]. Imposing this boundary condition atN collocation points the following set of equations is to solve(27)∑j=1NAkj·Bk=Ck,k=1,2,…,N,where(28)Akj=jLcos⁡θk1-jLcos⁡jL-1θk+Ecos⁡θk2jLcos⁡jL+1θkforj=1,…,N-1,(29)AkN=cos⁡θk2,(30)θk=πLk-1N-1fork=1,2,…,N,(31)Ck=12fork=1,2,…,N.Using the Darcy law given in the following form(32)qz=-kμdpdz,the longitudinal permeability can be related to the average velocity through the repeated element of the fiber system (the filtration velocity). The longitudinal component of the filtration velocity has the form(33)qz=1ΩT∬ΩFwr,θ·dΩF,where ΩF = ΩT – ΩP is the region occupied by the fluid in the repeated element (Figure 2).Using the definition of the dimensionless axial velocity (8) and the fiber volume fraction (3), the longitudinal component of the filtration velocity can be expressed as(34)qz=-b2μFφdpdz,where(35)Fφ=kβ·b2=∬ΩFWX,YdXdYβ·ΩT=2Ltan⁡π/L2∫0π/L∫EsecθWR,θ·RdRdθis the dimensionless component of the permeability tensor in the direction parallel to the fibers. A dimensionless parameter of the porous medium is calculated as β=Ltan⁡π/L and is equal to 2√3 for triangular, 4 for square, and 3√3 for hexagonal arrays of fibers. After an analytical integration in Eq. (35) the dimensionless permeability is a function of the number of collocation pointsN and can be calculated from the following formula:(36)F=E24-112-124cos⁡π/L2-πE48Ltan⁡π/L+2tan⁡π/L∑k=1N-1BkHkE2LkLk-2+GkLk+2+BNln⁡1Ecos⁡π/L-32,where(37)Hk=sin⁡1-Lkπ/L1-Lkcos⁡π/L1-Lk,Gk=sin⁡1+Lkπ/L1+Lkcos⁡π/L1+Lk.


## 5. Determination of the Slip Constant in the Beavers-Joseph Boundary Condition

Beavears and coauthors [9–11] determined the slip coefficient by two measurements of the mass flow rate. In the first case both walls of the considered channel were impermeable, while in the second case one of the walls was a porous medium in which there was filtration flow. The slip coefficient can be calculated from the following formula for the ratio between these two measured flow rates:(38)m˙m˙∗=1+3σ+2ασ1+ασ,where m˙ is the mass rate of flow through channel with one permeable wall, m˙∗ is the mass rate of flow in channel bounded by impermeable walls, σ =h/√k is the permeability parameter, andh is the channel width. It is convenient in numerical simulations to use directly the Beavers-Joseph boundary condition instead of measuring two mass flow rates. Let us consider a plane infinite channel adjoining with a plane infinite porous region in which flow is governed by the Darcy law and the Beavers-Joseph condition can be imposed on the boundary between these two areas in the following form:(39)dudx=αku-kμdpdz,where all quantities are constant. Let us introduce the following dimensionless variables:(40)U=μ·u-b2·dp/dz,σ=bk,where σ is the permeability parameter. It can be calculated based on the quantities presented in the previous section: (41)σ=bk=1β·F.Substituting the dimensionless quantities (40) into Eq. (39) the following dimensionless condition is obtained:(42)dUdX=ασσ2U-1.Assuming that values of the dimensionless slip velocity, the derivative of the dimensionless slip velocity on the boundaries dU/dX, and the dimensionless parameter of the permeability are known, the slip coefficient can be determined from the following formula:(43)α=σdU/dXσ2U-1.The slip velocityU and the derivative of the dimensionless velocity dU/dX on the assumed boundary can be calculated based on the numerical simulation presented in this section. The Trefftz method with the special purpose Trefftz functions is employed in this approach.Let us consider a layer of the porous medium located between two fixed plates (Figure4). In the whole considered domain the fluid flow problem can be described by the dimensionless Poisson equation in the following form:(44)∂2W∂R2+1R∂W∂R+1R2∂2W∂θ2=-1with the no-slip boundary conditionsW = 0 on the immovable wall and on the fibers.The fluid flow problem is solved by means of the Trefftz method. In the fibrous porous area the special purpose Trefftz functions are used to express the numerical solution. For the porous area these functions are associated with the large finite elements (Figure5). For each large finite element the approximate solution is expressed by(45)WR,θ=-14R2-E2+∑k=1N-1BkRLk-E2LkRLkcos⁡Lkθ+BNln⁡RE,which satisfies exactly both the governing Eq. (44) and some of the boundary conditions(46)W=0forR=E,(47)∂W∂θ=0forθ=0,π.For each large finite element the origin of the polar coordinate system R,θ is placed at the center of the fiber; see Figure 5.For the free flow area the approximate solution can be expressed in the Cartesian coordinate system as(48)W∗X,Y=-X2+Y24+∑n=0McnFnX,Y+∑n=1MdnGnX,Y,whereFn(X,Y) andGn(X,Y) are the trial functions(49)FnX,Y=∑k=0n/2-1kXn-2kY2kn-2k!2k!;(50)GnX,Y=∑k=0n-1/2-1kXn-2k-1Y2k+1n-2k-1!2k+1!.For the free flow area the origin of the Cartesian coordinate (X,Y) is placed at the bottom of the symmetry line of the considered channel (at the bottom of the boundary between the porous medium and the free fluid flow region). The unknown coefficientsBk  (k = 1,2,…,N), cn, dn (n = 1,2,…,M), andc0 are determined using the boundary collocation technique by satisfying the boundary conditions (in particular the splitting boundary conditions between the large finite elements) as shown in Figure 6.Figure 6
The repeating part of the considered region with the boundary and contact conditions for (a) square, (b) triangular, and (c) hexagonal arrays.
(a) (b) (c)In the numerical simulation, the large finite element is applied in various configurations with different distribution of the collocation points, as it can be observed in Figure6. These configurations of the large finite elements are shown in Figure 7.Figure 7
The four different configurations of the large finite element with various distribution of the collocation point.
(a) (b) (c) (d)The widths of the porous and free flow area were the same and equal toH = 2∙LRfor square,H = 3∙LRfor hexagonal, andH = 1+LRfor triangular array, whereLR is the number of rows of fibers. For each single edge of the repeated element we have chosenMc = 4 collocation points (Figure 9). For the upper and lower edge of free flow area the number of collocation points is calculated by the formulaH∙Mc/4. The total number of collocation points for the free flow area was equal toNc = 2∙(H∙Mc/4 + 2). The total number of collocation points for the porous area was equal toNc = 2∙LR∙Mc for square,Nc = 2∙ Mc(LR +1) for triangular, andNc = 6∙Mc∙LR for hexagonal array.After solving of the microstructural boundary value problem for the second numerical experiment we can assume that the boundary between the porous and fluid regions is forX = H (in the global coordinate system related to the left immovable wall). We can calculate the average value of the slip velocityU and the average value of derivative dU/dX for this value ofX(at this boundary) by numerical integration. Substituting the obtained values into Eq. (43) and knowing the permeability, the slip constant can be determined.In order to summarize the meshless procedure presented in the paper, the numerical algorithm is shown in Figure8.Figure 8
Numerical algorithm of the presented meshless procedure for determination of the slip constant in the Beavers-Joseph boundary condition.Figure 9
The nondimensional component of the permeability tensor in the direction parallel to the fibers. The numerical results (marks) compared with the analytical solution (line).


## 6. Numerical Results

The numerical calculations were performed for the fibrous porous medium with cylindrical fibers arranged according to the triangular, square and hexagonal arrays. In the first part of the numerical experiment the values of the dimensionless effective permeability and the permeability coefficient as a function of the fibers volume fraction are calculated. The results are presented in Figures9 and 10, respectively. In that case 7 collocation points and the same number of trial functions are used.Figure 10
The permeability parameterσ as a function of the fiber volume fraction φ.The value of the dimensionless permeability as a function of the fibers volume fraction is presented in Figure9. PermeabilityF(φ) for all types of fiber arrays is similar and substantially different for φ > 0.3. The problem was solved forN = 7 collocation points on the boundary. The numerical results (marks) are compared with the analytical presented by Drummond and Tahir [17] (line). Compatibility between the results is great.In the next part the slip constants on the boundary of the considered region are determined for 10 rows. For the hexagonal case as a one row, we considered two parts of fibers and its mean two repeated elements. This is due to the fact that in this case the repeated elements must exist in pairs (Figures6(a), 6(b), and 6(c)).Figure11 shows the slip constant,α, as a function of the volume fraction of fibers, φ. The value of the slip constant is smaller than 1.0 and decreases with increasing the volume fraction of fibers for the high porosity. For the hexagonal and square arrays of the fibers α increases for a large share of fibers φ. For the triangular array the slip constant α decreases also for a low porosity. It is caused by very small values of the permeabilityF for high volume fraction for triangular array of fibers. In such a case the permeability parameter σ is very high. It is visible also in Figure 10. One can notice in Eq. (43) that the slip constant α depends on the permeability parameter σ.Figure 11
The slip constantα as a function of the fiber volume fraction φ for square (4), triangular (6), and hexagonal (3) arrays of fibers forLR = 10 andMc = 4.The effect of the number of rows of the fibersLR on the value of the effective permeability for the triangular, square, and hexagonal arrays was examined. Figures 12, 13, and 14 show the obtained numerical results. The influence of the number of the fibers rows is significant and its effect disappears forLR equal to 12. The influence of the number of the fibers rows is more significant for the porous medium with high porosity (small φ). For small value of the volume fraction of fibers, φ, just four rows of the fibers are sufficient for determining the slip constant, α.Figure 12
Influence of the number of rows at a slip constant for a hexagonal array of fibers.Figure 13
Influence of the number of rows at a slip constant for a square array of fibers.Figure 14
Influence of the number of rows at a slip constant for a triangular array of fibers.


## 7. Conclusions

In this paper the permeability and the slip constant on the boundary of the considered region were determined by numerical simulation of the imaginary physical experiment. The porous medium was modeled by a parallel bundle of straight fibers arranged in a regular triangular or square array. In order to determine the permeability, the flow driven by the pressure gradient in an unbounded porous medium was considered. To determine the slip constant, the shear flow in a flat immovable wall containing both the porous and free flow layers was considered. In both numerical simulations, the Trefftz method with the special purpose Trefftz functions and trial functions was applied.The viscous fluid flow through the permeable porous material belongs to the classic problems of fluid dynamics. In 1856 based on the experimental research, Darcy proposed a formula called today the Darcy’s law. This equation has a theoretical justification for homogeneous materials for sufficiently small Reynolds numbers (laminar flow). The permeability is determined by calculation of a volume flow rate. The flow is considered as a one-dimensional, uniform filtration flow. Obtaining these conditions in a physical experiment is not possible due to the channel walls. If the flow conditions are close to idealization (no channel walls), the Poiseuille flow is realized for small Reynolds numbers at a sufficient distance from the end of the channel. In our paper, this physical experiment is carried out numerically by solving the boundary value problem free from the influence of disturbing factors. In this situation, it can be assumed that the bundle of regular fibers is not limited in a direction perpendicular to the axis of the fibers. Assuming that the flow is steady and fully developed (Poiseuille flow), we provide uniformity of the flow. For this reason, our discussions were carried out for laminar flow and infinite media.In all considered cases the value of the slip constant is lower than one. The slip constant decreases with increasing the volume fraction of fibers for a high porosity. For square and hexagonal array the value of slip constant starts increasing for small porosity. For the same volume fraction of the fibers, both the permeability and the slip constant for the triangular array are smaller than for the square and the hexagonal array. The number of collocation points has inconsiderable effect on the simulation results. To determine the permeability a few collocation points are enough to get an accurate result (Mc = 4). The influence of the number of rows is more visible, especially for high porosity. For all the cases 10 rows of fibers are sufficient to determine the slip constant. The proposed algorithm is easy to implement and accurate and is mesh free.



---

*Source: 1494215-2019-04-08.xml*