Verification of the boundary condition at the porous medium – fluid interface

We study the hydrodynamic stability of liquid flowing down over the inclined layer of uniform porous medium and compare the results obtained in the different frameworks. The flow in porous medium is described by the Brinkman model and by the Darcy model with corresponding boundary conditions at the interface between the homogeneous fluid and porous medium. It is shown the critical Reynolds number is calculated for the Darcy model much lower than in case of the Brinkman model, while the flow velocity is the same in the both models. It is a pure mathematical effect, which can be used to verify the models and to determine the empirical coefficients in the boundary conditions from an experimental study of flow instability.


Introduction
The fluid flow over the saturated porous medium enthralls the fluid inside it.The drag force in a porous medium leads to the rapid decrease of the tangential velocity, and it permits the flow instability development.The coupled flows have been studied since 1950s [1][2][3][4][5][6][7].The problem has many applications in the environmental science and technology: contaminant transport in the rivers and atmosphere [8,9], solidification of the liquid metals [10][11][12], heat and mass transfer between the contacting media [13], optimization of the flows in the fuel cells [5], etc.
There is also a fundamental interest.The problem of the conditions at the interface between the coupled flows is not solved until the end.Many authors offer the different kinds of the boundary condition being applied to the different filtration models [1,3,10,[14][15][16][17][18].The porous media models by Darcy and Brinkman are the most frequently used in known studies [1].
The Darcy model usually uses the interface boundary conditions by Beavers and Joseph [3].These conditions were obtained from the experimental data assuming the Darcy law determines the filtration.The flow velocity is discontinuous in that case.An alternative for the Darcy model is the approximate conditions by le Bars and Worster [10].It is more precise than the first, because the velocity stays continuous in that model.This approach postulates the transition layer exists between the flows.In spite of the transition zone is located in the porous layer, the flow velocity inside it is assumed to have a parabolic profile which continues the overlying fluid velocity.
Another description of the filtration is given by the Brinkman model [19].It was introduced as empirical, but the further studies show that the Brinkman equations can be obtained from the Navier-Stokes equations by averaging at the microscopic scale [4,14,15,20].Moreover, it gives the boundary conditions at the interface between the porous a e-mail: kbtsiberkin@psu.rumedium and homogeneous fluid.Ochoa-Tapia and Whitaker obtained one of the most used condition for the Brinkman model [14,15].Despite it is obtained by the strict formal way, these conditions include other empirical parameter sets a tangential stress jump at the interface.Other conditions are proposed in the elder and the latest studies [16][17][18], but these variants are not well known yet.
It is shown the flow structure and the velocity magnitude obtained in the different models may be matched by tuning of the empirical parameters in the interface boundary conditions [13].The velocities are the same in the overlying layer and in the porous medium.The difference exists only in the thin transition zone between the flows.Thus, the models cannot be distinguished by the stationary flow measurements.However, we have found the difference of the calculated flow structure makes a great influence on the results of the theoretical analysis of flow stability.The critical Reynolds number and wave number of the flow in the framework of the Darcy model are significantly different from results for the Brinkman model.We suppose the comparison of the theoretical and experimental data about the flow stability can verify the filtration model and interface boundary conditions, can determine the empirical parameters of these models, and may clarify their applicability area.
In this study, we present a comparative, theoretical analysis of the coupled flow stability using the Darcy model with the interface boundary conditions by Beavers and Joseph, and the Brinkman model with the conditions by Ochoa-Tapia and Whitaker.In addition, we analyse an influence of the heavy admixture dissolved in the porous layer on flow stability.
The paper is organized as follows.After the Introduction, Section 2 describes the problem geometry and flow models.The stationary solutions are presented in Section 3, while the results of numerical linear stability analysis and comparison between the models are given in Section 4. Section 5 contains a brief description of the linear stabil- 2 Two-layer system

Problem geometry
We consider the layer of incompressible homogeneous fluid is flowing down over the inclined uniform porous medium, which is saturated with the same fluid (Figure 1).The overlying fluid is bounded by a free rigid surface from above, and the porous layer has an impermeable rigid wall from below.
The porous layer thickness is h p , and the full thickness of the system is h.The layers are inclined at angle θ to the horizon.The permeability and porosity of the porous medium are K and ϕ, respectively.The fluid has the density ρ and kinematic viscosity ν.
The fluid moves only because of the gravity force, and there is no additional longitudinal pumping.Besides, the system is isothermal, and the fluid does not contain any dissolved or weighted admixture excluding the problem considered in Section 5.

Governing equations
The overlying fluid layer is described by the Navier-Stokes equations: where v is flow velocity, p is pressure, Re is the Reynolds number, θ is the slope angle of the layers, and γ is a unit vector contrary the gravity acceleration g.We present all the equations in dimensionless form.The variable scales and parameters are listed in Section 2.2.3.
We use the Darcy model and the Brinkman model to describe the filtration [1,19].The Darcy equation sets the fluid velocity is proportional the pressure gradient: where v p is fluid velocity in pores, ϕ is medium porosity, and q replaces the Darcy number and determines the medium permeability (see Section 2.2.3).We save the inertial terms in the filtration equations.They are usually negligible, but we consider the medium with high permeability and porosity, thus the inertial terms are significant because of the pore-scale Reynolds number is large [1].It equals 10 2 ÷ 10 3 in the water plants, and 10 1 ÷ 10 2 in the foametal layer.The time-derivative is in the equation because of flow instability is a rapid process.The two-layer flow instability is oscillatory, and the fast travelling waves occur in the system.Thus, we save the inertial terms to describe these factors correctly.
The Brinkman equations are as follows: It includes an additional effective viscous term.The Brinkman model works well if the medium has the large porosity and permeability.This model was proposed as an empirical [19], but it may be also obtained theoretically by the averaging of the Navier-Stokes equation at the pore-scale level [4,14,15,20].

Boundary conditions
The filtration models include the equations of a different order and require the different boundary conditions at the walls and at the interface between homogeneous flow and porous medium.
The top surface z = 1 is rigid and free, and the conditions are as follows: It includes the undeformability condition and absence of the tangential viscous stress.The bottom wall z = 0 is rigid and impermeable.The condition for the Darcy model sets only the wall is impermeable: while the Brinkman model requires no-slip condition: At the same time, many different kinds of the interface boundary conditions were developed.These conditions usually include the empirical coefficients connected with the properties of the porous medium, and one may vary the parameters to equate the stationary flow velocity [13].A unified approach to describe the coupled flows is still not developed [1,13,16].
We apply the two most frequently used kinds of the interface boundary conditions at z = d.The conditions by Beavers and Joseph developed empirically for the Darcy model are as follows [3]: EPJ Web of Conferences where α is the parameter should be found from the experiment.The conditions by Beavers and Joseph determine the jump of the tangential flow velocity.The normal velocity component and pressure are continuous.The Brinkman model uses the conditions by Ochoa-Tapia and Whitaker obtained by the direct averaging of the Navier-Stokes equations at the pore-scale level [14,15]: Here, β is another empirical parameter.These conditions set the continuity of the flow velocity and normal stresses, while the tangential viscous stress has a jump at the interface.

Dimensionless parameters
We use the following units of the variables: Namely, the length is measured in the units of the whole system thickness.The time-scale is a characteristic viscous time.The velocity unit is doubled maximum velocity of the flow over the inclined plane.The pressure scale is determined by the longitudinal hydrostatic gradient.
The dimensionless equations and boundary conditions include the following parameters: where Re is the Reynolds number, d is a relative thickness of the porous layer and q is an equivalent of the Darcy number Da.Here, the Reynolds number is larger than it is usually determined for the channel flow.To compare the results of the current study with the channel flow stability [21], one should multiply Re by (1 − d) 3 .

Stationary flow
The equations ( 1)-( 6) permit the stationary solution for the plane-parallel flow.The velocity profiles calculated in the Darcy as well as Brinkman model frameworks are presented in Figure 2 and Figure 3, and the structure of flow near the "fluid-porous medium" interface is magnified.The velocity magnitudes are the same in the both layers.Nevertheless, the boundary conditions (10) and (11) produce the difference between the solutions.The Brinkman flow has a continuous profile at the interface (Figure 2).It reaches the filtration velocity determined by the Darcy law within the boundary layer with thickness is of order q −1 , which becomes wider in the high permeable medium.If q is low, the velocity profile tends to the parabolic Poiseuille-like flow in the whole system.Besides, the stationary velocity profile has a break at the interface is connected to the tangential stress jump.
On the other hand, the Darcy model gives a constant value of the velocity in the porous layer, thus the tangential velocity is discontinuous at the interface (Figure 3).That makes a significant influence on the results of the theoretical flow stability analysis.

Linear stability analysis 4.1 Shooting method
In this paragraph, we present a brief review of the method applied to the linear stability analysis.First, we introduce the small perturbations of the stationary velocity and pressure.Substitution where |v| U and p P, and linearization of the full equations and boundary conditions (1)-(11) produce the system of linear ordinary differential equations for the perturbations.Assuming they are longitudinal waves growing exponentially, i.e. v, p ∼ exp (λt + ikx), (16) and expressing the pressure through the velocity, we obtain the linear system of ordinary differential equations for the amplitude functions depending on z.It includes the Orr-Sommerfeld equation for the homogeneous fluid and its' equivalent for the porous medium [21].The equation system is of the sixth order if the Darcy model applied, and  eighth order in the case of the Brinkman model.It has no an analytic solution, and stability analysis is provided by the numerical shooting method [22][23][24].
We solve the equations in the both fluid layers numerically by the adaptive Runge-Kutta-Merson 4th order method, and try to satisfy the interface boundary conditions.It is possible only with some specific Re, k and other parameters.Thus, the semi-automatic numerical algorithm (namely, two-dimensional secant method) varies the Re, k and λ = iω until the boundary conditions are not satisfied with some precision.The method has some additional features, such as forced orthogonalization of the solution vectors must be linearly independent [24].The described method is well proven by the studied of many different problems of hydrodynamic and convective instability.

Results
The numerical linear stability analysis shows the huge difference between the flow models we use.We set d = 0.50, α = 1 and β = 0, and q = 100 and 150.The relative difference of stationary flow velocity calculated in the both frameworks is about 0.5%.Nevertheless, the critical Reynolds number for the Darcy model is almost one hundred times lower than for the Brinkman model (Figure 4).Table 1 contains the values of the wave number and critical Reynolds number calculated for the both models at the different q and fixed d, α and β.
The structure of the neutral curves is also different.The long-wave minimum of the curve dominates in the Brink-man model, while the short-wave perturbations are critical in the Darcy model.The great difference is a result of the flow structure changes at the interface between the overlying fluid and the porous medium.The velocity jump in the Darcy model determines the instability mechanism closes the Kelvin-Helmholtz.It also makes an influence on the Brinkman model, but it is mitigated here.To compare the stability threshold with the well-known data for the channel flow (Re ≈ 5700) [21], the critical Reynolds number may be scaled to the overlying fluid thickness (see Section 2.2.3).It gives the critical Re for the Brinkman model is comparable to the known value for the channel flow, while it is much lower for the Darcy model.
One should remember the obtained difference is the pure theoretical result, and it shows the way to verify the filtration models and boundary conditions.Comparing the critical Re and k with experimental data, one may check the reliability of the mathematical models and determine the empirical parameter in the boundary conditions.
The bimodality of the neutral curves is a specific feature of the multilayer systems [25,26].It means there are two physical mechanisms of the instability.The first generates the large-scale vortices cover the upper layer and have a relatively large longitudinal wavelength, while the second produces the short wave vortices localized near the interface.The detailed study of the two-layer system instability, its' main peculiarities and mechanisms of instability developlment are presented in papers [27,28].

Admixture transport model
In addition, we consider an influence of the heavy dissolved admixture in a porous layer on the flow stability.The contaminant is initially located in the porous medium.It forms the diffusional boundary layer near the interface, and the flow enthralls the admixture even at the stable state.
We use the Brinkman model ( 1), ( 2), ( 5) and ( 6) with the boundary conditions (11) to describe the fluid motion with an extended gravitational force term: where β C characterises the fluid density dependence on concentration: The dimensionless admixture transport equations are as follows: Here, the Schmidt number Sc = ν/D is introduced, where D is a molecular diffusion coefficient.The top and bottom boundaries of the system are impermeable for the contaminant: The admixture diffuses from the porous layer into the external flow and the diffusional boundary layers occur at the interface.Their dimensionless widths δ and δ p are as follows: We do not consider the light admixture.In that case, a convective instability may occur at the time moment when the concentration gradient in the diffusional boundary layers becomes large enough, even if there is no longitudinal flow.

Velocity profile and flow stability
A typical diffusion time is several thousand times larger than the characteristic viscous time.The instability development is the most rapid process in the two-layer system.Thus, we may provide the linear stability analysis using quasi-static concentration profile and corresponding velocity profile as the basic state.The analytical structure of the exact solution for the plane-parallel flow is huge, therefore we apply a piecewise approximation of the concentration profile and the basic profile of the velocity corresponding to it (Figure 5).Unfortunately, the approximate solution is also huge, and we do not present it here.The velocities of the pure fluid and flow with admixture have no large difference, therefore the critical Reynolds numbers close each other too (Figure 6).The contaminant presence makes a little flow acceleration in the porous medium and slight destabilization of the basic state.
The main effect of the heavy admixture on the flow stability is change of the neutral curve structure.The additional minimum occurs on the curves (Figure 7).It corresponds to a damping of the perturbation with some specific wavelength.The critical wavelength of that minimum obtained from the numerical data rises by the t 1/2 -law like the diffusional boundary layer width (23) until the both minima of the curve do not merge.Furthermore, we study an influence of the heavy admixture dissolved in the porous layer on flow stability.The critical Reynolds number and wave number do not change significantly, while the new minimum occurs at the neutral curve.Its' wavelength changes in time like the diffusional boundary layer width, and it means the additional perturbation damping occurs within the diffusional boundary layer near the interface.The heavy admixture does not produce a new instability mechanism, and only effect on the pure hydrodynamics mechanisms described previously.
Contrary to that, the light admixture causes a concentrational convection even if there is no flow.Here, the instability occurs when the concentration gradient becomes large enough.These results prove it is possible to use the heavy visualizing admixtures to study the stationary planeparallel flow stability and structure of the critical motions without large error.
We suppose the empirical parameters of the models can be determined from the experimental study of the flow stability.The flow structure in the transition zone cannot be defined with good precision.Therefore, the stability study can be more effective than the velocity measurements.The study of instability development does not require any information about the flow within the porous medium.The wave number and critical Reynolds number can be found from the flow structure only in a homogeneous fluid.
Another possible method to determine the parameters of the models is a measurement of the viscous friction force acting on porous medium.It is determined by the velocity magnitude in the main flow and by its' gradients in the transition zone.Thus, the calculating of the viscous force in the both models we consider and the result comparison with the experiments can be effective to obtain the empirical parameter values.Nevertheless, it seems to be not so effective to specify the filtration model because it is based on the stationary flow magnitude.

DOI: 10 Figure 1 .
Figure 1.Geometry of the two-layer system and schematic structure of the flow.

Figure 2 .
Figure 2. Profiles of the stationary plane-parallel flow in the Brinkman model.

Figure 3 .
Figure 3. Profiles of the stationary plane-parallel flow in the Darcy model.

Figure 5 .
Figure 5. Quasistatic velocity profile of the plane-parallel flow at q = 30 and q = 100.

Figure 6 .
Figure 6.Comparison of the neutral curves of the flow without and with the heavy admixture dissolved in the porous layer, at d = 0.35, and q = 100 and q = 130.

Figure 7 .
Figure 7. Neutral curves of the flow with heavy admixture dissolved in the porous layer at the different time moments, q = 100 and d = 0.35.

Table 1 .
Critical Re and k for the both flow models at different permeability of the porous layer, d = 0.50, α = 1 and β = 0.