Thermal convection in a Brinkman–Darcy–Kelvin–Voigt fluid with a generalized Maxwell–Cattaneo law

We investigate thoroughly a model for thermal convection of a class of viscoelastic fluids in a porous medium of Brinkman–Darcy type. The saturating fluids are of Kelvin–Voigt nature. The equations governing the temperature field arise from Maxwell–Cattaneo theory, although we include Guyer–Krumhansl terms, and we investigate the possibility of employing an objective derivative for the heat flux. The critical Rayleigh number for linear instability is calculated for both stationary and oscillatory convection. In addition a nonlinear stability analysis is carried out exactly.


Introduction
The Brinkman-Darcy equations have been used extensively to analyse flow of a viscous incompressible fluid in a porous medium which is not too dense in the sense that the porosity is not too small.Thermal convection in saturated porous media, which is the subject of this paper, has been intensively studied with the saturating fluid described by Brinkman-Darcy theory, see e.g.Rees [1], Postelnicu and Rees [2], Capone and Gianfrani [3], Capone et al. [4], and the references therein.The range of permeability values where the Darcy-Brinkman equations are valid is critically reviewed in Gentile and Straughan [5].
Recent research has also focussed on flow of viscoelastic fluids in porous media, where the fluid stress displays history dependent behaviour, see e.g.Amendola and Fabrizio [6], Antontsev and Rodrigues [7].Since flow of viscoelastic fluids in porous media has immense importance practically in, for example, the oil industry, in vascular dynamics Cavallini et al. [8], we here analyse thermal convection in a model for viscoelastic flow in a porous medium by employing the Kelvin-Voigt equations in a Brinkman-Darcy porous medium.
Kelvin-Voigt fluids have been analysed in great depth from the analytical viewpoint of existence, attractors, structural stability, see e.g.Oskolkov [9], Kalantarov and Titi [10], Damazio et al. [11], and the references therein.Thermal convection in a layer of a Kelvin-Voigt fluid has been investigated recently by Straughan [12].In addition [13,14] has demonstrated continuous dependence of the solution for the equations of Brinkman-Darcy-Kelvin-Voigt theory for the isothermal, improperly posed, backward in time problem, and in Straughan [15] for the Kelvin-Voigt equations forward in time.
It is important to emphasize that Maxwell-Cattaneo theory for heat flow is a very rich area in modern research.In addition to heat flow, the idea behind Maxwell-Cattaneo theory where a relaxation effect is included for an appropriate flux is being employed in many diverse research areas.Such models do not generally suffer from the effect of the Fourier law of infinite speed of heat propagation, and allow heat to travel with a finite wavespeed.The technique behind Maxwell-Cattaneo theory whereby the flux equation is transformed from a constitutive equation to an evolution equation has been applied for example, in nanoscale mechanics, Sellitto et al. [16], Jou et al. [17], in sound waves in porous media, Jordan et al. [18], in temperature wave motion, Jordan and Lambers [19], Carillo and Jordan [20], Christov [21,22], in convection in stellar atmospheres, Herrera and Falcon [23], Herrera [24], in collapse in stellar structures, Govender [25], in destruction of tumours, Andres and Pinnau [26], in drug delivery models, Ferreira and Oliveira [27], in describing how vegetation forms into patterns on landscapes, Consolo et al. [28,29], in chemotaxis, Barbera and Valenti [30], in pollution, Barbera et al. [31], and other examples may be found in the book by Straughan [32].
Notable works on thermal convection using a Maxwell-Cattaneo model, some of which include magnetohydrodynamic effects are by Bissell [33], Eltayeb et al. [34], Hughes et al. [35], and Papanicolaou et al. [36].Other treatments of the Cattaneo technique to include shock evolution, bifurcation theory, are by Jordan and Lambers [19].
The present work derives a novel model for thermal convection using Brinkman-Darcy-Kelvin-Voigt theory but we allow the heat flux to be of Maxwell-Cattaneo type.The novelties involve employing Guyer-Krumhansl theory in conjunction with Maxwell-Cattaneo theory, cf.Jou et al. [17], Van et al. [37]; splitting the heat flux into an instantaneous Fourier part and a memory part depending on the temperature gradient, using an extra flux idea of Mariano [38]; and allowing the heat flux equation to be objective, which necessarily introduces other terms, cf.Morro [39,40].Incorporation of Guyer-Krumhansl theory is pertinent since Van [37], Fulop et al. [41], suggest that this theory may be more relevant than Fourier theory, even at room temperatures, see also Berezovski [42,43], Capriz et al. [44], Carlomagno et al. [45], Cimmelli [46], Fama et al. [47], Rogolino and Cimmelli [48].
A complete stability analysis for the Bénard problem for the Brinkman-Darcy-Kelvin-Voigt equations is presented here.In addition to a detailed linear instability analysis for both stationary and oscillatory convection we include a fully nonlinear energy stability analysis.Such analyses for systems incorporating Maxwell-Cattaneo effects have previously proved difficult, cf.Straughan [49, p. 192].

Basic equations
Let v i (x, t), p(x, t), T (x, t) and Q i (x, t) denote the velocity, pressure, temperature, and heat flux at position x and time t.The Brinkman-Darcy-Kelvin-Voigt equations consist of a balance of momentum equation, a balance of mass equation a balance of energy equation and a Maxwell-Cattaneo like equation for the heat flux of form In these equations λ is the Kelvin-Voigt coefficient, ρ is the constant density of the fluid, ν is the kinematic viscosity, α is the thermal expansion coefficient, g is gravity, k = (0, 0, 1), μ 1 is the Darcy coefficient which is essentially the kinematic viscosity divided by the permeability of the porous medium, ζ is a positive constant, τ is the relaxation time, κ is the thermal diffusivity, and ξ1 , ξ2 are the Guyer-Krumhansl coefficients.Throughout this article we employ standard indicial notation in conjunction with the Einstein summation convention, so for example, where v ≡ (u, v, w) and x ≡ (x, y, z).For a nonlinear example The term denotes the Laplacian in R 3 .
Equations ( 1)-( 4) are defined on the horizontal layer {(x, y) ∈ R 2 } × {z ∈ (0, d)} for t > 0. The boundary conditions are that Equation ( 3) is suggested by Payne and Song ( [50], pp.181-189), and I believe it may be justified from work of Mariano [38].We interpret the ζ T term to arise as an extra flux induced by the microstructure of the porous skeleton.If we assume v i ≡ 0 in Eqs. ( 3) and (4) then one may eliminate Q i and one finds the temperature satisfies the equation where ξ = ξ1 + ξ2 .Thus, the temperature field satisfies a second order in time differential equation but this equation is not hyperbolic.As the constants ξ, τ, ζ, κ are positive the ∂ T /∂t terms lead to strong damping.Alternatively, we may think of the heat flux being in two parts as in Mariano [38], so the total heat flux is Q i + F i where Q i is given by a Cattaneo-like evolution equation such as (4), whereas is a Fourier term.If we do not include the Guyer-Krumhansl terms then the system becomes (when v i ≡ 0) Equation (7) 3 may be integrated employing an integrating factor, and assuming fading memory of T ,i we may show Thus, the total heat flux is which is a viscoelastic-like term for the temperature gradient plus a term emphasizing the instantaneous temperature gradient.Such a representation is common in viscoelasticity, cf.Boltzmann [51], Miller [52], or for a modern appreciation of the work of Boltzmann see Markowitz [53].Eliminating Q i and F i from (7) leads in this case to the equation We choose a Guyer-Krumhansl theory in (4), cf.Jou et al. [17], Straughan [32], Van et al. [37], Fulop et al. [41]; recent interesting and inspiring articles dealing with Guyer-Krumhansl type theories, their generalizations to other areas in Continuum Mechanics, and novel applications are Berezovski [42,43], Capriz et al. [44], Carlomagno et al. [45], Cimmelli [46], Fama [47], Rogolino and Cimmelli [48].Indeed, Van et al. [37] is an interesting article containing experimental and theoretical work which suggests Guyer-Krumhansl theory may be more accurate than Fourier or traditional Maxwell-Cattaneo theory, even for temperatures typical of laboratory experiments.
In Eq. ( 4) the term DQ i /Dt can possibly have several forms.One such form is the material derivative and this derivative is shown to lead to a Galilean invariant theory by Christov and Jordan [54].On the other hand Morro [39,40] notes that ( 9) is not objective and therefore one should employ a suitable objective derivative.A general objective derivative is given by Morro [39,40] as where for our purposes γ is a constant, and It may be enlightening to employ (10) in a complete analysis of thermal convection, but we here follow Christov [55] and employ a Lie derivative for which γ = −1.This derivative was employed by Ciarletta and Straughan [56] and Tibullo and Zampoli [57], and by Bissell [33], Eltayeb et al. [34], Hughes et al. [35].The resulting theory employing (10) with γ = −1 is often referred to as Cattaneo-Christov theory.
It is of interest to note that Morro [39] describes various objective derivatives and he observes that the form of derivative in (10) with γ = −1 was suggested by Truesdell [58].Furthermore, Straughan ([32], pp.[22][23][24] observes that employing a relaxation effect like that in (4) was suggested by Graffi [59] in the context of electromagnetism.Thus, when a relaxation theory is employed with an objective derivative (10) where γ = −1, such a theory could be associated with the names of Graffi and Truesdell.
In this work we now employ (4) with DQ i /Dt given by (10) with γ = −1.For clarity we write Eq. (4) in this case as The steady conduction solution to Eqs. ( 1), ( 2), ( 3), ( 11) and ( 5) in whose stability we are interested has form where β is the temperature gradient, The steady pressure is a quadratic function of z determined from (1).

Thermal convection
To analyse the instability/stability of the steady solution ( 12) we define perturbation variables u i , θ, q i , π by These expressions are substituted into Eqs.( 1), ( 2), ( 3) and ( 4) and are nondimensionalized with the scalings where Sg is a parameter introduced in Papanicolaou et al. [36].The Rayleigh number Ra is introduced as In this way we arrive at the non-dimensional perturbation equations, where the *s are dropped, where (13) hold on R 2 × (0, 1) × {t > 0}, with w ≡ u 3 .
Care must be taken with the boundary conditions when Guyer-Krumhansl terms are present.The non-dimensional perturbation boundary conditions we employ are or together with horizontal periodicity of the solution, cf.Straughan ([49], p. 51).The period cell of the solution is denoted by V .When we use (15) one additionally requires V q 3 dx = 0, to exclude constant vertical heat flux perturbations.

Linear instability
To develop a linear instability analysis for ( 13) and ( 14) we drop the quadratic terms and introduce a time dependence like and we then have an eigenvalue problem for the growth rate σ .To solve this eigenvalue problem for fixed boundary conditions results in a heavy numerical computation.One must take curlcurl of (13) 1 and then one is left with solving three fourth order equations for u, v, w (≡ (u 1 , u 2 , u 3 )), three second order equations for q 1 , q 2 , q 3 , and one second order equation for θ .We thus need twenty boundary conditions.Sixteen of these follow from ( 5) and ( 2) and are w = 0, ∂w ∂z = 0, θ = 0, q 1 = 0, To determine the remaining four boundary conditions we take curlcurl of (13) 1 and evaluate the result for components 1 and 2 on the boundaries z = 0, 1.In this way one derives the further four conditions as and Since this is the first analysis of the convection model presented here we solve Eq. ( 13) subject to stress free surface boundary conditions, as is done in e.g.Eltayeb et al. [34], Hughes et al. [35].
The condition n j σ i j = 0 on z = 0, 1, reduces to and then with the continuity equation one sees that w ,zz = 0, on z = 0, 1.
For linear instability we then reduce Eq. ( 13) to where For stress free boundary conditions we may show from ( 16) that w, θ and may be represented by sin series of form w = ∞ n=1 w n sin nπ z, with a similar form for θ, .The stationary convection boundary is found by evaluating a determinant derived from ( 16) to be where = π 2 + a 2 , with a being the wavenumber.(Actually, = n 2 π 2 + a 2 , but it may be shown n = 1 yields the minimum.)The critical values of the Rayleigh number for stationary convection, Ra stat = R 2 stat , are then found by minimizing numerically R 2 in (17) in a 2 .
To find the oscillatory convection boundary we put σ = iω, ω ∈ R, in (16), and then solve a determinant equation which yields a cubic in σ .The real and imaginary parts of this equation are resolved and the oscillatory convection Rayleigh number, R 2 osc , may be shown to follow from where and The critical values of the Rayleigh number for oscillatory convection may be found by minimizing R 2 osc from (18) in a 2 .

Remark 1
It is worth pointing out that for linear instability via stationary or oscillatory convection the results obtained here for the Lie derivative where γ = −1 in ( 10) are exactly the same as would be obtained if we had employed a material derivative for Q i in (4).Of course, the nonlinear theory is different in both cases.The point here is that if we had not employed a Truesdell derivative in (10), i.e. had we selected a value of γ = −1, then even in the linearized theory the results are different from when one simply employs the material derivative.

Nonlinear stability
Let • and (, •, ) denote the norm and inner product on L 2 (V ).Further, let • q be the norm on L q (V ) where we omit the index when q = 2.
To investigate nonlinear stability of (12), we multiply (13) 1 by u i , (13) 3 by θ , (13) 4 by q i , and integrate each over a period cell V .After some integrations by parts and use of the boundary conditions we may obtain the identities, d dt and and By addition of these equations and integration by parts on the (θ ,i , q i ) term we may arrive at the energy equation where the energy function is given by the dissipation, D E , is the production term has form and where the cubic nonlinearity has form We henceforth discard the ξ 2 term and work with the following energy inequality [which follows from ( 22)], where the dissipation function, D, is given by where is the space of admissible solutions.To obtain nonlinear stability from ( 27) we obtain and we suppose now R < R E .To handle the cubic nonlinear term we use the Cauchy-Schwarz inequality to find and then we employ the Sobolev inequality q 4 ≤ c where c 1 is a constant depending on V .Then from (30) one deduces where Provided E 1/2 (0) < b/c one may employ a continuity argument to show, cf.Straughan [49, pp. 15-16], where Hence, we have demonstrated nonlinear stability in the E measure provided To find R E we obtain the Euler-Lagrange equations from (29)as where is a Lagrange multiplier.We solve Eq. ( 36) for stress free boundary conditions.In this case we reduce Eq.(36) to 123 The term q 3,3 is eliminated from (37) 1 and the remaining equations may be reduced to the relation, The critical Rayleigh number of nonlinear stability is determined by minimizing (38) in a 2 .

Remark 2
We have discarded the ξ 2 term in the analysis leading to (27).This term may be retained but the resulting analysis of the subsequent Euler-Lagrange system is substantially more involved and leads to very little improvement in the stability threshold.Since this is the first analuysis of this system we prefer to employ the simpler argument and not obscure the physical highlights with technical details.

Remark 3
If we had employed a material derivative for Q i in (4) rather than the form (11) then the energy stability results would be for global stability instead of the conditional result here, see (35).

Unconditional nonlinear stability
It is of interest to observe that we may obtain a result of unconditional nonlinear stability (i.e. for all initial data) from Eq. ( 13).To achieve this we replace (13) 4 by differentiating that equation by ∂/∂ x i to obtain the evolution equation for as One now works with the energy identities (19), (20), and the equation obtained by multiplying (39) by and integrating over V (where we are assuming q 3,3 = 0 on z = 0, 1, ), namely We now add (19), (20) and (40) to see that where the energy function, F, is now while whereas, Define now where H is the space of admissible solutions and then from (41) one may find Inequality ( 47) is integrated after use of Poincaré's inequality on D F to see that for a constant c 3 depending on b 3 .Unconditional nonlinear stability ensues from inequality (48).The limit for R U is R U = 1 and this yields the energy stability threshold.Henceforth, we take R U = 1.The Euler-Lagrange equations arising from (45) with R U = 1 are and from these equations we deduce, for two stress free surfaces, The unconditional energy stability limit may be found by minimizing R

Remark 4
One may develop another energy estimate by employing (19), (40), then multiplying (13) 3 by − θ and integrating over V .This avoids the and θ quadratic terms in (41), and would lead to a sharper Rayleigh number threshold than (49).In this case one derives an "energy" equation of form where and the cubic nonlinearity has form The drawback with this is that the nonlinear stability so obtained is conditional.
We observe that and then from (38), where ( 17) is recognized.Thus, we know as it should be.Firstly, from Tables 1, 2 and 3 we observe that for Sg in the range 10 −2 -10 −4 oscillatory convection will not be observed.It is possible to witness oscillatory convection, but Sg has to be much larger, perhaps of order 100, and then a 2 osc is very small.Consolo [28,29] do suggest τ values large enough for Sg in the range O( 102 ), and hence oscillatory convection may not be negligible in all cases.We should point out that the Kelvin-Voigt parameter is only instrumental in the oscillatory convection case, although it is essential in the nonlinear theory.
From Tables 1, 2 and 3 we note that the variation in Ra stat , Ra en , and Ra unc is significant as ζ increases.Indeed, Ra stat , Ra en and Ra unc increase substantially with increasing ζ and this means the system is more stable and convective motion commences less easily.On the other hand, as ξ 1 increases Ra stat decreases, but Ra en and Ra unc increase, but the variation is less than that observed with ζ varying.As ξ 1 increases the nonlinear values, Ra en , Ra unc , become much closer to Ra stat , as may be seen in Tables 1, 2 and 3, and as shown in Fig. 1.

Conclusions
We have introduced a new model for thermal convection of a Kelvin-Voigt viscoelastic fluid in a Darcy-Brinkman porous material with the temperature field being governed by a procedure suggested by [50], and incorporating Maxwell-Cattaneo-Guyer-Krumhansl theory of heat flow.Attention is given to which derivative is employed for the heat flux since this may affect the instability values for convective motion substantially.
Linear instability values and nonlinear stability thresholds are calculated for the case of two stress free surfaces, and in addition to a thorough analysis of stationary and oscillatory convection, an exact nonlinear analysis is included.We have not seen such a nonlinear analysis for a model employing Maxwell-Cattaneo theory previously.

Remark 5
We have dealt here with the Brinkman-Darcy-Kelvin-Voigt system of equations.The mathematical structure of the Navier-Stokes-Voigt equations with a linear friction term is exactly the same as Eqs.( 1), ( 2), ( 3) and (4), cf.DiPlinio et al. [60], provided we employ the thermal structure as introduced in Payne and Song [50] and as used here.DiPlinio et al. [60] analyse an analogous isothermal system, although they also include an extra viscoelastic memory term.In the pure fluid case, the −λ u i,t term represents the Kelvin-Voigt part in the Navier-Stokes-Voigt equations and the −ξ u i piece is the Rayleigh friction contribution.

Fig. 1
Fig. 1 Critical values of Rayleigh numbers for stationary convection, Ra stat , energy stability, Ra energy , and unconditional energy stability, Ra uncond , against ξ 1 .The other parameters have values Pr = 6, ζ = 1, λ = 0.1, Sg = 0.01, ξ = 0.1, ξ 2 = 0.1 2 unc in a 2 .It is straightforward to show that R 2 stat ≥ R 2 unc and some numerical results for the unconditional critical Rayleigh number, Ra unc , are given in Sect.7. The results for the unconditional critical Rayleigh number are less than the linear instability threshold, although they do represent a bound for global nonlinear stability.As ζ increases Ra unc becomes much closer to Ra stat , the linear stationary convection threshold.