Stability Analysis of Gravity Currents of a Power-Law Fluid in a Porous Medium

Abstract

We analyse the linear stability of self-similar shallow, two-dimensional and axisymmetric gravity currents of a viscous power-law non-Newtonian fluid in a porous medium. The flow domain is initially saturated by a fluid lighter than the intruding fluid, whose volume varies with time as . The transition between decelerated and accelerated currents occurs at α = 2 for two-dimensional and at α = 3 for axisymmetric geometry. Stability is investigated analytically for special values of α and numerically in the remaining cases; axisymmetric currents are analysed only for radially varying perturbations. The two-dimensional currents are linearly stable for α < 2 (decelerated currents) with a continuum spectrum of eigenvalues and unstable for α = 2, with a growth rate proportional to the square of the fluid behavior index. The axisymmetric currents are linearly stable for any α < 3 (decelerated currents) with a continuum spectrum of eigenvalues, while for α = 3 no firm conclusion can be drawn. For α > 2 (two-dimensional accelerated currents) and α > 3 (axisymmetric accelerated currents) the linear stability analysis is of limited value since the hypotheses of the model will be violated.

1. Introduction

A thorough comprehension of numerous class of flows in porous media is firmly based on approximations and simplifying assumptions leading to differential problems amenable to analytical or quasi-analytical solutions. A notable example is the shallow water assumption, with negligible vertical velocity and a horizontal length scale much larger than the vertical one [1–4]; the approximation leads in turn to modelling a wide class of propagation problems as one-dimensional ones. An important class of solutions arising in several one-dimensional transient flows driven by gravity is constituted by self-similar solutions, representable in terms of functions of one variable and describing the “intermediate asymptotic” behavior of the system in a region where the solution is no longer dependent on the specific initial and/or boundary conditions adopted [3]. Whenever self-similar solutions can be obtained in analytical or in approximate analytical form, they represent an excellent test to verify model correctness by means of comparison with experiments; they also constitute a reliable benchmark for numerical integration.

Noteworthy examples of self-similar solutions to single-phase gravity driven flows in porous media include flow over a horizontal surface in plane [2] and radial geometry, Lyle et al. [5]; these solutions were extended to incorporate two-layer flow [6], the effect of a sloping bottom [7], the action of impermeable confining boundaries [8], drainage effects [9], and non-Newtonian power-law fluids [10], also in a vertically graded porous medium [11]. The dipole solution, originally derived by Barenblatt [3], was extended by King and Woods [12] to include drainage, and by Mathunjwa and Hogg [13] to incorporate a vertical variation in permeability.

An important step following the derivation of a self-similar solution is the check of its stability. In fluid mechanics, such an analysis can essentially be reduced to the study of the spatiotemporal evolution of disturbances applied to the system. In a broad sense, the evolution of instabilities leads to new solutions of the differential problem describing the flow and provides hints on possible stable solutions. Stability analyses were performed in the context of several filtration problems amenable to self-similar solutions. The linear stability of the Barenblatt-Pattle (B-P) self-similar solutions of the porous medium equation, describing a number of flows including viscous and porous media gravity currents, was investigated by Grundy and McLaughlin [14] for plane and axisymmetric geometry. Their analysis was later extended to nonradially symmetric perturbations by Mathunjwa and Hogg [15]. The linear stability of a class of self-similar solutions to the filtration-absorption equation was demonstrated by Barenblatt et al. [16] and Chertok [17]. Mathunjwa and Hogg [13] showed the linear stability of the dipole self-similar solution of the first kind [3]. With respect to all these analyses the present study is new since it considers non-Newtonian fluids.

The study of gravity currents in porous media was recently extended to non-Newtonian fluids by Di Federico et al. [18, 19], respectively, for two-dimensional and axisymmetric geometry; self-similar solutions in analytical and numerical form were obtained for a power-law fluid spreading in a uniform medium subject to the injection of a volume of intruding fluid increasing with time with the exponent . The motivation for their study was the nonlinear rheological nature displayed by fluids flowing in porous media in a variety of industrial and environmental applications ([20] and references therein).

In this paper, we analyse the stability of the aforementioned solutions. Section 2 presents a linear stability analysis for two-dimensional geometry. First, the special cases (correspondent to an instantaneous release of a finite volume of fluid) and (correspondent to a linear increase in time of injected volume) are discussed, deriving a closed-form expression for the eigenvalues of the associated Sturm-Liouville (SL) problem. Second, the stability of the general case with is demonstrated numerically. In Section 3, a similar analysis is conducted for axisymmetric geometry, considering radially varying perturbations. Stability is investigated analytically for and numerically for . Concluding remarks are outlined in Section 4. Further mathematical details are included in the Appendices.

2. Stability Analysis in Two-Dimensional Geometry

Di Federico et al. [18] (Figure 1) consider the motion of shallow two-dimensional gravity currents of a purely viscous and relatively heavy non-Newtonian fluid of uniform density in a homogeneous porous layer saturated with a lighter fluid of uniform density . The flow is driven by the release from a source at the boundary of a volume , above a horizontal impermeable bottom.

Figure 1 

Sketch of the gravity current of density intruding into a porous medium saturated with fluid of density .

The non-Newtonian fluid in simple shear is described by the rheological Ostwald-DeWaele power-law model: where is the shear stress, the shear rate, the fluid consistency index, and the flow behavior index. When the model describes a shear-thinning behaviour, when , it describes a shear-thickening behaviour, and when , it describes a classical Newtonian behavior.

The flow law for the fluid is a modified Darcy’s law taking into account the nonlinearity of the rheological equation (1) proposed by Bird et al. [21]: where is the generalized pressure, the pressure, the vertical coordinate, the fluid density, the specific gravity, the Darcy velocity, the intrinsic permeability coefficient, and the effective viscosity. The mobility ratio is [22] where denotes the porosity. It is convenient to write (2)-(3) as [23] highlighting the dependence upon permeability .

The Darcy velocity in the direction is given by the one-dimensional counterpart of (4), which reads, for horizontal flow In the shallow water approximation, the pressure gradient driving the flow is thus related to the gradient of the unknown free surface of height by where any capillary forces (surface tension) between the two liquids have been neglected. Introducing (6) into (5), we find that For one-dimensional transient flow, the local continuity condition takes the form Substituting (7) into (8), we obtain the nonlinear partial differential equation for height as: Consider the flux to be such that the total volume per unit width of the dense fluid in the porous medium at any time is given by , where both and are constant. The global mass balance equation then reads The mathematical statement of the problem is completed by the boundary condition: Equations (9)–(11) may be nondimensionalized by setting where the time, space, and velocity scales are thereby generalizing to non-Newtonian case the choice of Vella and Huppert [7], valid for (stability of the special case is discussed later in Section 2.2). Equations (9)–(11) may be written in nondimensional form as where and are dimensionless space and time, is the dimensionless current depth, and represents the length of the gravity current. By introducing the similarity variable , and denoting the value of for by , the similarity solution takes the form with being the solution of subject to , while is given by Another necessary boundary condition is derived for the first derivative at the front end of the current. In the limit a Frobenius series expansion gives [18, 19] and the solution is determined only by the front end conditions.

Once is determined, the gravity current length is given by ; the case acts as a transition between decelerated and accelerated currents. In the cases and , the problem is amenable to an analytical solution. For a generic value of , approximate solutions are available in Frobenius series or in quasi-analytical form; else, a numerical integration of (16) is performed to derive the shape function [18].

To analyze the linear stability of the self-similar solution expressed in functional form by (15), we note that (14a) can be written as The first step is a variable transformation with in order to avoid the difficulties in managing a movable boundary and to reduce the domain to . In the new independent variable , (19) becomes The perturbed solution can be expressed as where is the basic state, is a small parameter, and is a perturbation. Due to perturbation, the boundary is displaced and we deal with bounded perturbed support , where and are the perturbations of the boundaries at the edge of the current. To evaluate the structure of and , we can use the analytical expression of the basic state available for ; that is, . Substituting in one of the boundary conditions and linearizing results in , it follows that the displacement is of and can be expressed as assuming a generic nonsymmetric perturbation; in (22a) and (22b), and are the perturbed edges of the domain in the positive and negative directions, and and are the first coefficients of their series expansion. Substituting the perturbed solution (21) in (20) and comparing powers of , at the equation for the basic state is recovered. At , the following expression in the perturbed variable is derived: The basic state is expressed by (15) and we seek a solution for in (23) by separating the variables according to the following expression: In (24a), (24b), and (24c), the function is the eigenfunction corresponding to the th eigenvalue , while , , and are the coefficients of the series. If the eigenvalues exist and are nonnegative, and the set of eigenfunctions is complete, the perturbation decays in time and the basic state is stable.

The boundary conditions at become in terms of After substitution of (24a), (24b), and (24c) into (23), the perturbation is the solution of the following differential problem, where the prime symbol indicates the derivative with respect to : with and the following boundary conditions: Equation (26) with boundary conditions (28) can be converted in the classical self-adjoint SL form [24]: with The SL structure of the differential problem guarantees that (i) the set of eigenfunctions is complete and (ii) different and orthogonal eigenfunctions correspond to different eigenvalues. In the following, we first evaluate the stability of the solution for the two special cases and , for which a basic solution is available in closed form; we then discuss the stability of the general solution valid for and . We do not explore the case , which is of limited practical interest and subject to the limitations implicit in the shallow water approximation, namely, the condition . The spatial gradient of the current height is readily obtained by (15) as showing that the former condition is time dependent and that the spatial gradient within the domain tends to decay with time if , to remain constant if , and to increase with time if . In the latter case, the basic solution obtained by solving (16)-(17) is valid only for a limited time interval in a portion of the domain, but sooner or later it will not satisfy the shallow water approximation in the entire domain.

Equation (31) can be reformulated introducing the average gradient of the shape function at a given time, equal to ; hence,

Imposing that (a small quantity) at time requires that Figure 2 shows that at short times the limiting condition is even more stringent than the asymptotic condition for .

Figure 2 

Variation of the limiting condition to guarantee that .

2.1. The Case

For , the self-similar solution has the closed-form expression and the SL problem describing the perturbation can be handled with analytical functions. The differential problem for the perturbed variable takes the form (29) with and it is possible to reduce the differential problem to a standard hypergeometric one by using the transformation . The solution in terms of hypergeometric functions is then We further note that replacing by the differential problem admits the solution Consequently, where the second coefficient has been dropped to ensure the regularity in the origin.

The spectrum is continuum with real and positive values [25] and is given by Hence, the eigenvalues are always positive, and the basic solution is stable. The spectrum is linear and the minimum eigenvalue is equal to . Perturbations in shear-thickening fluids gravity currents decay faster than perturbations in shear-thinning fluids gravity currents.

2.2. The Case

If a different scaling is required since the characteristic time scale defined by (13) breaks down and an additional velocity scale arises, namely , in addition to the velocity scale [18]. Defining an arbitrary time scale and spatial length scales and dimensionless variables as ; , , and , the original dimensional problem can be expressed by the following differential equation in dimensionless form: where is the ratio between the two velocity scales in the problem. By introducing the similarity variable , the similarity solution is of the form and the dimensionless length of the gravity current at a given time is given by once is determined. Equation (40) has the following analytical solution [18]: In the linear stability analysis for , the approximations are similar to those already applied for . The details are reported in Appendix B. We look for a solution to (B.1) in the separable forms (24a), (24b), and (24c), with boundary conditions given by (25a), (25b), and (25c).

Substituting (43) into (B.1) gives with With the substitution , (45) can be reduced to a Laguerre form admitting a solution in terms of Laguerre functions with eigenvalues equal to . The point is a regular singularity. For integer eigenvalues, the Laguerre functions become Laguerre polynomials, whose orthogonality requires nonnegative eigenvalues; hence, , and the system is linearly unstable. The less unstable eigenvalue corresponds to shear-thinning fluids.

Note that, , in this case indicating that the value can represent the limiting value between two different behaviours of the sign of the eigenvalues: negative for and positive for .

2.3. The Case of and

Without an analytical expression of the basic state, the eigenvalues and the eigenfunctions cannot be computed in closed form; however, some properties of the problem can be used to evaluate the sign of the eigenvalues. For a regular SL problem, it can be demonstrated ([26], p. 169; see also Appendix A) that if then there exist an infinite number of eigenvalues, which are all real and positive. The theorem can be extended to a singular SL problem in a finite domain with the singularity arising from the conditions , provided that the existence of the eigenvalues and the completeness of the set of eigenfunctions are demonstrated. The singularity at the domain edges renders the boundary conditions redundant but requires the boundedness of at the edges. A sufficient, but not necessary, condition to guarantee that the set of eigenfunctions is complete is that the integral is finite (see, e.g., [27]), with being the associated Green function.

For and , a numerical approach is required to define the functions involved in the Sturm-Liouville problem. Starting from the numerical solution of the basic state, the functions , , and are computed. and are always positive and for . The results are shown in Figure 3 for , 1.0, 1.5, and .

Figure 3 

Function for and , 1.0, and 1.5.

The boundedness condition of is satisfied, the integral in (47) is finite, and the system is linearly stable. This result is also supported by the numerical integration described in Appendix C. The results for (Newtonian fluid) are coincident with a similar analysis due to Grundy and McLaughlin [14].

3. The Axisymmetric Case

The axisymmetric case can be analysed by using the same approximations adopted in the two-dimensional case. The geometry is the one sketched in Figure 1, substituting with . Introducing the dimensionless variables , , , , and , where the time, space, and velocity scales are, for , , , and , and the differential problem in nondimensional form is expressed for by the following equations [19]: with , representing the front position of the gravity current. Equations (48a) and (48b) do not include terms containing the derivatives with respect to azimuthal coordinate , since the present stability analysis is restricted to axisymmetric perturbations, modulated only along the radial coordinate. By means of the similarity variable , the solution to (48a) and (48b) was derived in the form yielding The current is accelerated for and decelerated for .

The analytical solution of the basic state is available only for : For arbitrary (for the special case , see [19]) (49)–(51a) can be solved numerically with (51b) and a second boundary condition, which can be computed developing the asymptotic solution in terms of a Frobenius series expansion near the current tip () obtaining Since results , , in analogy to the plane case, (48a) and (48b) become Expressing the perturbed solution as and substituting in (54), the following equation in the perturbed variable is obtained at : Assuming that the basic state is expressed by (49), we seek its solution by separating the variables with the following boundary conditions: After substitution, the perturbation is the solution of the following differential problem: with which can be converted in the classical SL form [24].

For the special case , the differential equation (57) for the perturbations becomes in SL form: with It is then possible to reduce the differential problem to a standard hypergeometric one by using the transformation . The solution is in terms of hypergeometric functions: while the spectrum is continuum with real and positive values and is given by Hence, the eigenvalues are always positive and the perturbations decay in time for any value of the flow behaviour index . As for the two-dimensional case, the spectrum is linear and the minimum eigenvalue is equal to . Perturbations in shear-thickening fluids gravity currents decay faster than perturbations in shear-thinning fluids gravity currents.

For , a different scaling is requested since the time scale breaks down and the differential problem is similar to (48a) and (48b) except for a coefficient on the left hand side (see [19]), in analogy with (40) in the plane case. Since an analytical solution is not available for and (decelerated currents) no analytical solution for the perturbation can be found, but on the same arguments used in Section 2.3 we found that the conditions for stability (46) are met. In fact the bounded nature of the term at the edges of the domain can be demonstrated and the completeness of the set of eigenfunctions is verifiable. Also a numerical integration conducted with the same methodology adopted for plane geometry (not shown) confirms that the system is stable for . The results for (Newtonian fluid) are coincident with a similar analysis due to Grundy and McLaughlin [14] and to Mathunjwa and Hogg [15].

For the basic state solution is not available and no analytical solution for the perturbation can be found. In this particular case (since it marks the transition from decelerated to accelerated currents) the term is bounded at the edges of the domain but the integral in (47) is not finite. Since the finiteness of the integral is a sufficient, but not necessary, condition to guarantee that the set of eigenfunctions is complete, no firm conclusion can be drawn, even though in analogy with the two-dimensional case a linear instability could be inferred.

For the spatial gradient, given by grows in time; hence, the shallow water approximation will be eventually violated. With the same reasoning adopted for the two-dimensional currents, is the limit of validity of the present linear stability analysis. In addition, recasting (63) in terms of similarity variables indicates that the limits at short times are more stringent than the asymptotic limit for (see Section 2 for a similar analysis in the two-dimensional case).

4. Conclusion

A novel linear stability analysis of self-similar solutions to non-Newtonian power-law gravity currents intruding in a porous layer saturated by a lighter fluid has been performed for two-dimensional and axisymmetric flows; in the latter case, the analysis was limited to radial perturbations. Our investigation leads to the following key results.(1)For two-dimensional flows, closed-form expressions of the stability equation can be obtained in the special cases (instantaneous injection of a finite volume of fluid) and . For and any , the corresponding differential problem is solved in terms of hypergeometric functions with a continuous spectrum of positive eigenvalues; hence, the system is linearly stable for all values of . The spectrum of the eigenvalues is linear and the decay rate for perturbation in shear-thickening fluids gravity currents is larger than the decay rate for perturbations in shear-thinning fluids gravity currents. For the differential problem is solved in terms of Laguerre polynomials with a discrete spectrum of negative eigenvalues; hence, the system is linearly unstable and the eigenvalues are proportional to with a higher growth rate in gravity currents of shear-thickening than in those of shear-thinning fluids. For and (decelerated currents), the existence of a complete set of eigenfunctions is demonstrated and the flow is linearly stable for both shear-thinning and shear-thickening fluids, including the Newtonian case. For (accelerated currents) the present linear stability is of limited value since accelerated currents are characterised by a spatial gradient growing in time; hence, the shallow water approximation will be eventually violated.(2)For axisymmetric flows, an analytical solution is available only for . In this case the differential problem shows a continuum spectrum of real and positive eigenvalues. Hence, the flow is linearly stable. The spectrum of the eigenvalues is linear and decay rate for perturbation in shear-thickening fluids gravity currents is larger than the decay rate for perturbations in shear-thinning fluids gravity currents. For and (decelerated currents), the flow is always linearly stable. For no firm conclusion can be drawn, even though in analogy with the plane case we would expect the arising of instability. For axisymmetric accelerated currents as for the plane case, the present linear stability is of limited values since the limiting approximations of the shallow water approximation will be eventually violated.

Appendices

A. Stability Condition for a Regular Sturm Liouville Problem

Given the Sturm-Liouville problem we multiply each term by and integrate over the domain : Assuming that and in , it follows that if in and , then . The latter condition is always met if the SL problem is singular with .

B. Stability of the Analytical Solution of Plane Case for

Equation (40) can be written as and after introducing the new variable in order to reduce the domain to , (B.1) becomes Consider as in (21) a perturbed solution in the form , where is the basic state, is a small parameter, and is a perturbation. Substituting the perturbed solution in (B.2) and linearizing yields which to order reduces to in the perturbed variable ; the boundary conditions are identical to those listed in (25a), (25b), and (25c) for the general case.

C. Numerical Integration for the Plane Case

In order to check the validity of our conjecture, we perform the numerical integration of (20) at a given time starting from a self-similar solution plus a perturbation. The integration is based on a fourth-order, four-stage explicit Runge-Kutta scheme [28], with centred-space finite difference scheme for the spatial derivatives. The numerical calculations are carried out for , 1.0 and , 1.0, and 1.5, simulating the perturbation with periodic disturbances satisfying the integral boundary condition. Though the numerical solutions conserve volumes to a high degree of accuracy, due to numerical dissipation a volume loss of less than 0.01% is computed after a time interval , starting at .

Figure 4 shows the results of the integration for and and 1.0, respectively, in (a) and (b). The perturbation is a sine function, given by Inspection of (a) and (b) reveals that the initial perturbation tends to decay with time. Similar results are obtained for and (not shown). Stability is thus confirmed.

(a)

(b)

Figure 4 

Time evolution of the perturbation for a shear-thinning fluid and (a) and (b). The dashed line is the perturbation at a given time, and the subsequent curves refer to perturbation at subsequent time, with a time step equal to 0.001. The perturbation experiences an initial growth (not shown) followed by a decay.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.