## Abstract

The nonlinear stability of a rectangular porous channel saturated by a fluid is here investigated. The aspect ratio of the channel is assumed to be variable. The channel walls are considered impermeable and adiabatic except for the horizontal top which is assumed to be isothermal. The viscous dissipation is acting inside the channel as internal heat generator. A basic throughflow is imposed, and the nonlinear convective stability is investigated by means of the generalized integral transform technique. The neutral stability curve is compared with the one obtained by the linear stability analysis already present in the literature. The growth rate analysis of different unstable modes is performed. The Nusselt number is investigated for several supercritical configurations in order to better understand how the system behaves when conditions far away from neutral stability are considered. The patterns of the neutrally stable convective cells are also reported. Nonlinear simulations support the results obtained by means of the linear stability analysis, confirming that viscous dissipation alone is indeed capable of inducing mixed convection. Low Gebhart or high Péclet numbers lead to a transient overheating of the originally motionless fluid before it settles in its convective steady state.

## 1. Introduction

The models of seepage flow in porous media include the effects of viscosity on the momentum transfer. Under sufficiently intense pressure gradients, and if the average thermal conductivity of the saturated porous medium is sufficiently small, then the fluid viscosity has also an important effect on heat transfer. By assuming that the seepage flow can be modelled by Darcy’s law, the term in the local energy balance equation that yields viscous heating is proportional to the square of the local velocity [1–4]. The thermal effects of viscosity, viscous dissipation and internal heating have been widely investigated with reference to either fluids [5–10] or fluid-saturated porous media [2,11–13].

Viscous dissipation acts as a heat source and, as such, is capable of building up a temperature gradient within a fluid-saturated porous medium. Under suitable boundary conditions, this temperature gradient may yield an unstable temperature stratification and, possibly, it may produce a thermal instability of the saturated porous medium. This possibility has been first investigated for fluids [5,14–18] and, recently, it has been extended to fluid-saturated porous media [19–22]. A review of the literature on this topic can be found in Barletta [23]. Different mechanisms can drive thermal instabilities, such as the temperature dependence of viscosity, buoyancy force and surface tension gradients. In this study, we will focus our attention on the instability driven by the buoyancy force. Thus, the instability is essentially of the Rayleigh–Bénard type, but the non-uniform temperature gradient is induced by the inner viscous heating and not by the temperature difference between the boundaries.

A characteristic of the existing literature on the thermal instability induced by viscous dissipation is that the studies reported are relative to linear analyses. With reference to porous media, linear stability analyses of flows with viscous dissipation were carried out by Rees *et al.* [11], Barletta *et al.* [19], Barletta & Rees [24], Barletta *et al.* [25], Barletta & Nield [26] and Alves *et al.* [27] The linear stability analysis provides an efficient tool for the investigation of the threshold conditions for the onset of instability. Exceptions emerge in those cases where the transition to instability depends on the amplitude of the perturbations acting on the basic flow. When this is the case, linearization of the governing equations for the perturbation fields is not a reliable approximation, and one needs to tackle the nonlinear problem to correctly predict the conditions for the onset of instability. Subcritical instabilities are one example, where critical parameters are perturbation-amplitude dependent. Another domain where linearization fails to provide reliable results is the strongly supercritical regime where the parametric conditions are far away from neutral stability. The energy method [28,29] proved to be an effective tool to study the possible emergence of subcritical instability. On the other hand, the deeply supercritical regime can be either tackled, at least for a limited parametric range, by a weakly nonlinear stability analysis (see, for instance, chapter 7 of [30]) or, more generally, by a direct numerical solution of the fully nonlinear governing equations. The latter approach is the one that will be employed in this paper. In order to do so, the numerical method adopted here to solve the nonlinear governing equations is the generalized integral transform technique (GITT) [31]. This technique is particularly useful for such studies because it does not introduce dissipative or dispersive errors in its treatment of the governing equation spatial resolution [32]. The former might lead to artificial steady states, i.e. stability of otherwise unstable asymptotic solutions, whereas the latter might introduce incorrect phase speeds or even spurious high-frequency modes.

The aim of this paper is to reconsider the results of the linear stability analysis carried out by Barletta *et al.* [19] by studying the nonlinear governing equations for the perturbations instead of their linearized formulation. Barletta *et al.* [19] analysed the linear stability of Darcy’s flow with viscous dissipation in a horizontal porous layer where the lower plane boundary is adiabatic and the upper plane boundary is isothermal. Here, this analysis will be revisited in a nonlinear perspective by assuming a lateral confinement of the porous layer by adiabatic vertical sidewalls. The GITT will be adopted to transform the governing disturbance equations into ordinary differential equations expressing the evolution in time. Numerical integration of these equations yields the nonlinear evolution in time of the disturbances. The comparison between linear and nonlinear evolutions shows that no subcritical instability emerges. Furthermore, the behaviour in the supercritical regime will be investigated. The flow patterns at nonlinear saturation will be determined and the Nusselt number at the upper isothermal boundary will be evaluated. As far as the authors are aware, this is the first nonlinear demonstration of viscous-dissipation-induced thermal instabilities in the literature.

## 2. Mathematical model

The throughflow in a horizontal fluid-saturated porous channel is here studied with respect to the onset of convective instability. The channel is assumed to be rectangular with arbitrary aspect ratio and impermeable walls. From the thermal viewpoint, the boundaries are all adiabatic except for the upper horizontal wall that is assumed to be isothermal. The contribution of viscous dissipation is taken into account in such a way that it is the only possible source of convective instability. Viscous friction can indeed generate a non-trivial basic temperature profile that becomes unstable under particular conditions. The instability is triggered by the buoyancy force that is modelled through the Oberbeck–Boussinesq approximation applied to Darcy’s law. The governing equations are written by assuming that

— the Oberbeck–Boussinesq approximation is a valid model of the thermal buoyancy in the fluid;

— the momentum balance within the saturated porous medium can be expressed by Darcy’s law;

— local thermal equilibrium between the fluid and the porous material holds true; and

— viscous dissipation cannot be neglected.

Thus, the governing equations in their dimensionless form, together with their boundary conditions, are given by
** u**=(

*u*,

*v*,

*w*) is the velocity vector,

*T*is the temperature,

*e*_{y}is the unit vector in the

*y*-direction,

*t*is the time,

**=(**

*x**x*,

*y*,

*z*) is the Cartesian coordinate vector,

*Ge*the Gebhart number and

*s*is the aspect ratio. The rescaling used to obtain equations (2.1) is

*σ*is the ratio between the volumetric heat capacity of the porous medium and the volumetric heat capacity of the fluid,

*H*is the channel height,

*α*is the average thermal diffusivity,

*T*

_{0}is the reference temperature,

*ν*is the kinematic viscosity,

*g*is the gravitational acceleration modulus,

*β*is the thermal expansion coefficient of the fluid,

*K*is the permeability of the porous medium,

*c*is the fluid specific heat and

*L*is the channel width. A sketch of the system is drawn in figure 1.

## 3. Basic state

Equations (2.1) admit a stationary solution, which is here referred to as ‘basic state’. This basic velocity field is fully developed and characterized by a constant value equal to the Péclet number, *Pe*. The basic temperature field depends on the vertical *y*-coordinate, on the Gebhart number and on the Péclet number. They are given by
*Pe*=*U*_{0}*L*/*α*, *U*_{0} is the prescribed dimensional velocity of the basic throughflow, the subscript *b* refers to the basic state and
*Pe*≠0 and *Ge*≠0.

## 4. Perturbed governing equations

In order to investigate the nonlinear stability of the basic state defined in equation (3.1), one perturbs this state with disturbances of arbitrary magnitude, namely
*Θ*_{0}(*x*,*y*,*z*)=*Θ*(*x*,*y*,*z*,0), namely
*x*,*y*)-plane, the system here investigated does not display any sensitivity to the *z*-direction. Hence, all the contributions in this direction, from this point on, are neglected. This allows the introduction of the streamfunction *U*=∂*Ψ*/∂*y* and *V* =−∂*Ψ*/∂*x*. Thus, equations (4.3) reduce to

## 5. Generalized integral transforms technique

The GITT is a hybrid analytical–numerical method that is here employed to perform the nonlinear stability analysis. Its application starts by defining an eigenfunction expansion for the different fields in all spatial variables. This expansion is based on auxiliary eigenvalue problems. One then proceeds to integral transform the governing equations, eliminating their spatial dependence and generating an infinite series of ordinary differential equations in time for the transformed fields. Because these eigenvalue problems belong to the Sturm–Liuville class, they give rise to a convergent series solution. Hence, this series can be properly truncated in order to guarantee a user prescribed error tolerance. This system is thus solved numerically, and the transformed fields obtained are used to recover the original temperature and streamfunction fields.

### (a) The auxiliary eigenvalue problems

The auxiliary eigenvalue problems for the streamfunction disturbance fields in the *x*- and *y*-directions are defined as

On the other hand, the auxiliary eigenvalue problems for the temperature disturbance fields in the *x*- and *y*-directions are defined as

The inverse-transform pairs for the streamfunction disturbances are

### (b) Integral transform procedure

Once the auxiliary eigenvalue problems have been defined, it is possible to integral transform governing equations (4.5). One first operates this transformation along the *x*-direction, i.e. multiply equation (4.5a) by *x*≤*s* to obtain
*x*-direction, this procedure is repeated in the *y*-direction. In order to do so, equation (5.11) is multiplied by *y*≤1 to yield

## 6. Numerical analysis

Once all the governing equations and initial conditions are transformed, the truncated version of the system of equations (5.14), (5.16) and (5.18) is solved numerically by means of the software Mathematica. In particular, the tool *NDSolve* is employed. The transformed disturbance fields *Θ* and streamfunction *Ψ* are recovered using the inverse transformations in equations (5.5)–(5.8). In order to verify if a given parameter set leads to unstable conditions, we make sure that the velocity of the steady convective state at different random positions inside the domain has reached a non-zero value even though a motionless initial state was imposed. For this work, *t*_{max}=10^{3} was considered a large enough time value to verify if the perturbed system converged to a steady convective state. If we define *tol*=10^{−5}, then we assumed to reach the steady state when ∂*Θ*/∂*t*<*tol*. For all the simulations performed here, we reached the nonlinear saturation for 10^{1}<*t*<10^{2}.

A key parameter in the present analysis is named here *v*_{T}, which is the dimensionless threshold velocity used to identify the onset of instability. In other words, the velocity field departs from its initial conduction state when its magnitude is higher than the value provided by this parameter. It distinguishes the incoherent motion owing to Brownian diffusion from the coherent motion owing to convection. We identify this threshold with the mass diffusion velocity *v*_{D} of a fluid in a porous medium. It should be stressed that we are dealing with a secondary flow superposed to a basic flow. In order to understand the base flow convective instability, the mass diffusion velocity associated with the disturbance is subtracted from the sought overall solution. We thus consider *v*_{D} as a constant value independent of the magnitude of the basic flow. A typical value for the effective diffusivity coefficient for solutes in liquid saturated porous media is *D*∼10^{−9} m^{2} s^{−1} [33]. The respective dimensional velocity is then given by *v*_{D}=*D*/*H*. On the other hand, because the thermal diffusivity of a fluid is of order ∼10^{−7} m^{2} s^{−1}, on average, equation (2.2) can be used to define the dimensionless velocity *v*_{T} as
*N* was imposed for all truncated series in this paper. Because they are convergent series, *N* has to be chosen high enough to guarantee a user prescribed error tolerance. However, in order to increase the accuracy of the results without increasing the computational cost, a reordering scheme was applied to each summation series pair in equations (5.14) and (5.16). The idea underneath this reordering is that terms with indexes such as (*i*,*j*)=(2,2) are more important than terms with indexes such as (*i*,*j*)=(1,15) and, hence, should be added first in the summation process. In doing so, convergence is achieved using a significantly smaller number of terms in the truncated equations. The criterion employed here to reorder these series *a priori* is based on the classical pure thermal diffusion eigenvalues from the auxiliary eigenvalue problems in equations (5.1) and (5.3). A comparison among results obtained with different number of equations *N* is shown in figure 2 for *N*=5,10,20,40,80. It is possible to note that *N*=5 is not sufficient to draw the lower neutral stability branch. With *N*=10, the first part of the lower branch can be drawn, but not the other parts. The whole lower branch can only be drawn with at least *N*=40, where the results obtained with *N*=80 do not add any additional information to this figure.

Finally, it is important to note that a sensitivity test of the results with respect to different shapes and intensities of the initial condition *Θ*_{0}(*x*,*y*) was performed. It revealed that the asymptotic stability analysis results do not depend on the choice of initial condition. Minor changes in disturbance amplitude and shape have only a minor effect on the unsteady solution behaviour. Here follows a list of the approximations employed to perform the nonlinear stability analysis

— the tolerance is set to

*tol*=10^{−5};— the truncation threshold is set to

*N*=40; and— the threshold velocity used to identify the onset of instability is set to

*v*_{T}=10^{−2}.

## 7. Results

### (a) Comparison with linear stability analysis

Figure 3 displays the threshold values of the governing parameter *Λ* for the onset of convective instability as a function of the aspect ratio *s*. Different neutral stability curves and the benchmark linear stability analysis from Barletta & Storesletten [34] (solid lines) are drawn in figure 3 for comparison purposes. The vertical dashed lines are taken from this linear analysis as well to define the separation between two different convection cell patterns at the onset of instability. In order to compare the results from this linear stability analysis with the present nonlinear simulations, figure 3 includes both limiting cases considered in the nonlinear analysis: *Ge*→0 (dotted lines) and *Ge*=1 (dashed lines). These two limits span the entire range of values that the Gebhart number can assume in a physically sensible set-up. In most practical cases, the Gebhart number is very small. The condition *Ge*=1 is far beyond any conceivable application, except for some geophysical cases. In fact, to obtain such a high value of *Ge* with typical working fluids or gases, the thickness reference value *H* would be of several kilometres. Moreover, the case *Ge*→0 indicates a negligible nonlinear contribution of the viscous dissipation, whereas *Ge*=1 indicates an extremely high heating owing to the nonlinear viscous dissipation. It is important to point out that one should not consider the limit *Ge*→0 as necessarily equivalent to the limit *Λ*→0. Equation (3.2) suggests an alternative scenario where *Ge*→0 with a large enough value of the Péclet number, *Pe*≫1, in such a way that *Λ*∼*O*(1). If one considers a channel thicknesses of 10^{−1} *m*, engine oil as working fluid, and a flow velocity compatible with Darcy’s law, 10^{−2} m s^{−1}, the flow would be characterized by *Pe*∼10^{4}. The limiting case *Ge*→0 with *Pe*≫1 is physically sensible and it lets the system display an important heat production even when small values of the Gebhart number are considered. In such a case, the very high basic flow velocity compensates for the vanishingly small Gebhart number. More details on the conditions when viscous dissipation starts playing a role can be found in Barletta *et al.* [19]. Figure 3 shows a nearly perfect agreement between the linear stability analysis and the nonlinear results obtained here using the GITT. Furthermore, figure 3 highlights the fact that the threshold for the onset of convective instability does not depend on *Ge* alone, but it rather depends on *Λ*. As mentioned in §6, the onset of instability is independent of the disturbance amplitude or shape. Hence, it is possible to conclude that subcritical instabilities are not present. This behaviour is consistent with the non-dissipative nature of the GITT, which does not force the solution to converge towards an unstable steady state even when the initial conditions are very close to it [32]. Furthermore, it is important to note that an absolute instability is observed in the nonlinear simulations, because perturbations grow exponentially in time beyond threshold values of *Λ*. Nevertheless, this onset of absolute instability must coincide with the onset of convective instability from linear theory for this model. This occurs, because growth in the base flow direction has been precluded owing to the removal of this coordinate in equations (4.5).

Figure 4 displays the growth rate *r* of supercritical configurations with *Λ*=65 as a function of the aspect ratio *s*. Once again, linear stability results from Barletta & Storesletten [34] (solid line) are used to validate the present nonlinear simulation results for *Ge*→0 (plus symbol) and *Ge*=1 (times symbol). This slightly supercritical configuration, i.e. *Λ*=65, was chosen to guarantee that the linear stability results are accurate enough. Nevertheless, a nearly perfect match among all three datasets can be observed in this figure, highlighting the fact that growth rates also depend weakly on *Ge* for a fixed *Λ*. Furthermore, figure 4 shows four local maxima that identify the aspect ratios that yield the fastest-growing disturbances for *Λ*=65, i.e. the preferred modes of convective instability that will be observed in an unforced experiment subject to random noise at each respective parametric configuration. The wavenumbers of these disturbances are related to the aspect ratio of the channel, because the horizontal bounds force the system to be perturbed by waves characterized by wavelengths which are multiples of the aspect ratio.

### (b) Nonlinear stability analysis

Now that the nonlinear code has been validated against benchmark linear stability results, it is possible to use it with confidence to investigate the destabilising role played by viscous dissipation away from marginal stability conditions. Figure 5 displays the behaviour of the average Nusselt number *Nu* as a function of time *t* at the horizontal upper wall, *y*=1, mathematically defined as

for *s*=1.285, which is the aspect ratio that yields the critical *Λ* for a single cell pattern. Figure 5*a* displays curves for different values of *Λ* and *Ge*=1, whereas figure 5*b* displays curves for different values of *Ge* and *Λ*=65. They show two asymptotic trends, a linear exponential disturbance growth at early times and a stationary asymptotic condition at late times. The former is predicted by a linear stability analysis, but the latter can only be predicted by a nonlinear analysis. It shows the steady cell pattern reached by the system in the specific supercritical configuration considered. One can note in figure 5 that *Λ* and *Ge*. Furthermore, nonlinear saturation occurs earlier as *Λ* increases, which is expected because the linear growth rate *r* increases strongly with *Λ*. On the other hand, a vanishing *Ge* leads to a vanishing *Nu* at saturation. This can be explained if we consider that *Ge*→0 means a vanishing nonlinear internal heat generation and, because there is no other heat source but viscous dissipation, the relative value of average Nusselt number would also be vanishing. Finally, figure 5 also shows that *Nu* is not always a monotonically increasing function of time, i.e. nonlinear saturation does not always lead directly to a stationary condition. Such a novel behaviour creates a maximum between both asymptotic time limits. It is clear in figure 5*b* that such maxima occur for low *Ge*, but figure 5*a* indicates that they also occur at high *Pe*. This is implicit in the fact that *Λ* can only be increased with a fixed *Ge* by increasing *Pe*. The maxima behaviour as functions of *Λ* are summed up in the following correlation from figure 5*a*

Figure 6 displays the behaviour of the local steady Nusselt number, *Nu*_{x}, relative to the horizontal upper wall, *y*=1. All results shown from this point on were generated from steady data, *Nu*_{x} is rescaled by *Λ* in figure 6, so that all resulting curves have the same order of magnitude. They are drawn for different values of *Λ* with *Ge*=1 and *s*=1.285, again chosen to yield the critical *Λ* for a single cell pattern. Figure 6 highlights the strong dependence of the local Nusselt number shape and magnitude on *Λ*. At slightly supercritical *Λ* values, such as *Λ*=65, *Nu*_{x} amplitude is small and its shape is essentially sinusoidal. As *Λ* increases to higher values, such as 100 and 200, the maximum amplitude increases proportionally to *Λ*^{0.377}, whereas the sinusoidal shape gets increasingly distorted. An important characteristic of this flow already observed in figures 2–4 is the strong cell pattern dependence on *s*. On the other hand, figure 6 also shows the existence of a cell pattern dependence on *Λ*, although relatively weaker. This can be identified by the *Nu*_{x} wavelength halving as *Λ* increases from 200 to 400, which means a pattern change from one to two cells. Hence, the discontinuous changes in both *Nu*_{x} and *Nu* help us to identify the parametric configurations which display a transition between two different cell patterns. Figure 7 displays the behaviour of *Nu* over the horizontal boundary at *y*=1 for *Ge*=1. Figure 7*a* shows this behaviour as a function of *s* for *Λ*=65, 75 and 85, whereas figure 7*b* shows this behaviour as a function of *Λ* for *s*=1.755, 3.000, 4.200 and 5.420. The values of *s* in figure 7*b* are chosen to illustrate consecutive cell pattern changes from one to five. In order to evaluate *Nu*, *Ge*=1 was chosen because it guarantees a non negligible heat production inside the domain. This additional understanding of the supercritical region highlighted in figures 6 and 7 allows us to complete figure 3, including the cell pattern transitions on the stability plane (*s*,*Λ*). Such a plane is shown in figure 8, where the lines *Λ*_{1}, *Λ*_{2}, *Λ*_{3} and *Λ*_{4} divide the plane in four supercritical regions. Each region identifies convective modes characterized by a particular number of cells, which are illustrated by the number of convective cell symbols (concentric circles) drawn inside every region. The lines which divide the four regions are described by the following equations
*s* as in figures 9 and 10 but for *Λ*=70, which brings us into the supercritical region. One may note that the isolines are here distorted with respect to figure 9 because of the effect of convection. The analogous to figure 10 for the supercritical region is not shown because it is qualitatively the same, because the cell configuration will change only at much higher values of *Λ*, as it is shown in figure 8.

## 8. Conclusion

The effect of viscous dissipation on the flow of a Newtonian fluid in a confined rectangular porous channel modelled by Darcy’s law has been investigated. Viscous dissipation was included in the energy equation as the only potential source of thermal instabilities. The study has been focused on two-dimensional modes, so that only longitudinal modes have been considered. Previous linear stability analyses published in the literature have shown that this source term is capable of inducing an unstable thermal stratification in the basic flow which can lead to the onset of convective instability in the absence of any other effect. In this paper, a fully nonlinear stability analysis of this flow has been performed to verify the results of the linear stability and probe into the supercritical region. A generalized integral transform technique using orthogonal basis functions was employed to eliminate the spatial dependence of the model equations, generating a system of ordinary differential equations in time for the nonlinear disturbances which has been solved using the software Mathematica. As far as the authors are aware, this is the first attempt in the literature to perform a nonlinear stability analysis of viscous-dissipation-induced thermal instabilities.

The major findings of this nonlinear study are summarized below:

— nonlinear simulations support all linear stability findings, confirming for the first time that viscous dissipation alone is indeed capable of inducing mixed convection;

— low Gebhart

*Ge*or high Péclet*Pe*numbers lead to a transient overheating of the originally motionless fluid before it settles in its convective steady state;— the number of convective cells formed depends weakly on the parameter

*Λ*=*Ge**Pe*^{2}, which controls the strength of viscous dissipation as a source term for instabilities; and— the maximum amplitude of the average Nusselt number

*Nu*increases with*Λ*and its originally sinusoidal shape in the marginally stable region is increasingly distorted as this control parameter increases.

## Authors' contributions

All authors conceived the mathematical model, interpreted the computational results and wrote the paper. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

This research was supported by Alma Mater Studiorum Università di Bologna, Fundação Cultural Ítalo-Brasileira (FIBRA), FAPERJ through grant no. E-26/103.254/2011 and CNPq through grant no. 481072/2012-8.

- Received January 14, 2016.
- Accepted April 19, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.