#### 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.

#### 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:where* u* is the velocity of the fluid flow calculated at boundary (the slip velocity),* k *is the permeability of the porous medium,* n *is the normal direction to the boundary, and* q* 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 Figure 1) are analyzed. In Figure 1 one can find the characteristic geometrical features of the regions, which are the diameter of the fibers 2*a* and the distance between the fibers 2*b*.

**(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:where are the dimensionless Cartesian coordinates, are the dimensionless polar coordinates, and 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 Figure 1. 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 . The values of the parameter are , , and for triangular, square, and hexagonal arrays, respectively. The repeated element of the unbounded porous medium with the dimensionless variables is presented in Figure 2.

The volume fraction of the fiber is given bywhere Ω_{T} = 0.5 is the total area of the repeated element (the fluid and the fiber areas together) and Ω_{P} = 0.5* E*^{2}* ∙π*/

*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 (Figure 3).

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 Figure 4—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 to* H*. In Figure 4* LR* denotes the number of rows of fibers in the porous medium.

**(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 to* a, *and the distance between the fibers is equal to 2*b. *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).

#### 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 (Figure 1(a)), triangular (Figure 1(b)), and hexagonal (Figure 1(c)) arrays. The radius of the fibers is equal to* a, *and the distance between the fibers is equal to 2*b*. In such a case the equation of motion is reduced to a single partial differential equation expressed in the polar coordinate system bywhere is the velocity component in the* z*-axis direction [m/s], is the constant pressure gradient [Pa/m], is the viscosity of the fluid [Pa* ∙*s], and is the fluid domain. The above governing equation is subject to the following boundary conditions:Let us introduce the following dimensionless axial velocity:Furthermore, the problem can be considered in the repeating part presented in Figure 2. Now, the governing Eq. (4) has the dimensionless formwith the boundary conditionsThe 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:where the subscripts and

*p*denote general and particular solutions, respectively. The latter for Eq. (9) can be expressed asFurthermore, 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:In the literature, see, e.g., [21, 24]; it is well known that such a general solution is given bywhere , and () are unknown constants. Thus the whole solution can be expressed asIn order to ensure satisfaction of the boundary condition (10) for any the constants and should be equal to zero. The logarithmic function is equal to zero when the argument is equal to 1. Because of that the constant can be taken as . Then the above equation can be rewritten in the formThe derivative of the above solution with respect to takes the following form:In order to satisfy the boundary condition (11) simultaneously for and the constants and () have to be equal to zero. Furthermore, the constants () for should satisfy the following relation: because for such angles the sinusoidal function equals zero. Thus In such a case, the approximate solution takes the form Furthermore, for the constants and () should satisfy the relationIt leads to the following formula for the constants ():Thus, finally the exact solution of Eq. (9) can be expressed using the special purpose Trefftz functions in the following form:From the infinite series in the above form we choose only* N *– 1 first terms and taking , the solution can be rewritten in the formThe solution (26) satisfies exactly the governing equation (9) and the boundary conditions (10)–(11). The unknown coefficients (*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 at* N* collocation points the following set of equations is to solvewhereUsing the Darcy law given in the following formthe 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 formwhere Ω_{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 aswhereis the dimensionless component of the permeability tensor in the direction parallel to the fibers. A dimensionless parameter of the porous medium is calculated as 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 points* N* and can be calculated from the following formula:where

#### 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:where is the mass rate of flow through channel with one permeable wall, is the mass rate of flow in channel bounded by impermeable walls, *σ* =* h*/√*k* is the permeability parameter, and* h* 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:where all quantities are constant. Let us introduce the following dimensionless variables:where is the permeability parameter. It can be calculated based on the quantities presented in the previous section: Substituting the dimensionless quantities (40) into Eq. (39) the following dimensionless condition is obtained:Assuming that values of the dimensionless slip velocity, the derivative of the dimensionless slip velocity on the boundaries d*U*/d*X*, and the dimensionless parameter of the permeability are known, the slip coefficient can be determined from the following formula:The slip velocity* U* and the derivative of the dimensionless velocity d*U*/d*X* 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 (Figure 4). In the whole considered domain the fluid flow problem can be described by the dimensionless Poisson equation in the following form:with the no-slip boundary conditions* W* = 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 (Figure 5). For each large finite element the approximate solution is expressed bywhich satisfies exactly both the governing Eq. (44) and some of the boundary conditionsFor each large finite element the origin of the polar coordinate system 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 aswhere* F*_{n}(*X,Y*) and* G*_{n}(*X,Y*) are the trial functionsFor 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 coefficients* B*_{k} (*k* = 1,2,…,*N*), , (*n* = 1,2,…,*M*), and* c*_{0} 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.

**(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 Figure 6. These configurations of the large finite elements are shown in Figure 7.

**(a)**

**(b)**

**(c)**

**(d)**

The widths of the porous and free flow area were the same and equal to* H* = 2* ∙LR *for square,

*H*= 3

*for hexagonal, and*

*∙*LR*H*= 1+

*LR*for triangular array, where

*LR*is the number of rows of fibers. For each single edge of the repeated element we have chosen

*Mc*= 4 collocation points (Figure 9). For the upper and lower edge of free flow area the number of collocation points is calculated by the formula

*H*/4. The total number of collocation points for the free flow area was equal to

*∙*Mc*Nc*= 2

*(*

*∙**H*/4 + 2). The total number of collocation points for the porous area was equal to

*∙*Mc*Nc*= 2

*for square,*

*∙*LR*∙*Mc*Nc*= 2

*(*

*∙*Mc*LR*+1) for triangular, and

*Nc*= 6

*for hexagonal array.*

*∙*Mc*∙*LRAfter 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 for* X = H* (in the global coordinate system related to the left immovable wall). We can calculate the average value of the slip velocity* U* and the average value of derivative d*U*/d*X* for this value of* X *(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 Figure 8.

#### 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 Figures 9 and 10, respectively. In that case 7 collocation points and the same number of trial functions are used.

The value of the dimensionless permeability as a function of the fibers volume fraction is presented in Figure 9. Permeability* F*(*φ*) for all types of fiber arrays is similar and substantially different for *φ* > 0.3. The problem was solved for* N* = 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 (Figures 6(a), 6(b), and 6(c)).

Figure 11 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 permeability

*F*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 .

The effect of the number of rows of the fibers* LR* 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 for* LR* 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, *α*.

#### 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.

#### Data Availability

All data can be accessed in the Results section of this article.

#### Conflicts of Interest

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

#### Acknowledgments

The work was supported by the grant 2015/19/D/ST8/02753 funded by the National Science Center, Poland. During realization of the work Dr. Jakub K. Grabski was supported by a scholarship funded by Foundation for Polish Science (FNP).