Linear Stability Analysis of Penetrative Convection via Internal Heating in a Ferrofluid Saturated Porous Layer

Abstract: Penetrative convection due to purely internal heating in a horizontal ferrofluid-saturated porous layer is examined by performing linear stability analysis. Four different types of heat supply functions are considered. The Darcy model is used to incorporate the effect of the porous medium. Numerical solutions are obtained by using the Chebyshev pseudospectral method, and the results are discussed for all three boundary conditions: when both boundaries are impermeable and conducting; when both boundaries are conducting with lower boundary impermeable and free upper boundary; and when both boundaries are impermeable with lower boundary conducting and upper with constant heat flux. The effect of the Langevin parameter, width of ferrofluid layer, permeability parameter, and nonlinearity of the fluid magnetization has been observed at the onset of penetrative convection for waterand ester-based ferrofluids. It is seen that the Langevin parameter, width of ferrofluid layer, and permeability parameter have stabilizing effects on the onset of convection, while the nonlinearity of the fluid magnetization advances the onset of convection.


Introduction
Ferrofluids (i.e., fluids which contain stable colloidal suspensions of one-domain particles of ferromagnetic and ferrimagnetic materials in liquid carriers) belong to a special class of fluid that exhibit both magnetic and fluid properties [1].The particles used in ferrofluids are usually made up of metal materials (ferromagnetic materials) like iron, nickel, cobalt, and also with their oxides (ferrimagnetic materials) like spinel-type ferrites, magnetite (Fe 3 O 4 ), etc. [2].The carrier liquid can be water, ethylene glycol, oils, etc.However, the choice of metal particles and carrier liquid for the preparation of a particular ferrofluid depends on the application for which the ferrofluid is prepared.Ferrofluids can be used as a coolant in thermal management devices (such as in loudspeakers and transformers as a coolant) and/or as a heat transfer medium in energy conversion systems such as heat exchangers and processes including boiling [3].These fluids have been widely used in commercial applications in many areas such as sealing, damping, heat transfer, bearing, and sensing.The first developed and commercialized ferrofluid product was a dynamic process seal, which has applications in areas such as lasers, semiconductor processing, X-ray machines, fiber optics, crystal growing systems, avionics, and heat treating furnaces [4].Ferrofluidic exclusion seals help to protect sensitive environments and critical machinery components, and therefore have important applications in various industrial areas, including robotics, textiles, computer peripherals, and machine tools.In addition to the above-mentioned applications, ferrofluids have been used extensively in bioengineering applications.The area associated with ferrofluid convection in the presence of an applied magnetic field has attracted a great amount of attention from many researchers over the past few decades.The primary study of the generalization of the classical Rayleigh-Bénard convection problem for a ferromagnetic fluid was conducted by Finlayson [5].The author investigated the convective instability of a ferromagnetic fluid for a fluid layer heated from below under the presence of a vertical magnetic field.He analyzed the instability in the presence and absence of gravity.The results were discussed for both shear-free and rigid horizontal boundaries under the framework of linear stability theory.Later, a significant amount of work in this direction was carried out by [6][7][8][9][10].The study of ferromagnetic convection in a porous medium also becomes important due to its significance in the emplacement of geophysically-imageable liquids into particular zones for subsequent imaging, controlled emplacement of liquids, or treatment of chemicals [11].Vaidyanathan et al. [12] studied the thermoconvective instability in a ferromagnetic fluid saturating a porous medium in the presence of a vertical magnetic field.The authors used linear stability theory and discussed the results graphically as well as numerically for free-free boundary condition.They also reported in their study that a porous medium of small permeability tends to stabilize the system.The reader is referred to [13][14][15][16] for more related study.
Penetrative convection is usually defined as a type of motion which may occur whenever convection in a thermally-unstable fluid layer extends from above and below into adjacent stable layers [17].This kind of convection is a widely-occurring phenomenon, and commonly occurs in geophysical and astrophysical flows.Internal heating is one of the most widely employed mechanisms to characterize this phenomenon.Many systems, such as the Earth's mantle, radiating atmospheres, the cores of large main-sequence stars, and nearly any engineered system where chemical and nuclear reactions take place in a fluid environment have internal sources or sinks of buoyancy [18].The presence of an internal heat source or sink can give rise to a situation where one part of a layer is naturally convecting while the other remains stable, and hence allows penetrative convection to occur.Natural convection driven by internal heating can be seen in a large number of physical phenomena (for example, in the field of nuclear energy), and plays an important role in applications in nuclear reactor cores, in post-accident heat removal, etc.This phenomenon also has its importance in the development of a metal waste from spent nuclear fuel, fire and combustion modeling, and the storage of spent nuclear fuel [19].Gasser and Kazimi [20] were the first to study convection in a porous media due to both internal heating and heating from below by using the method of linear stability of small disturbance.Later, the problem of penetrative convection due to purely internal heating in a fluid-saturated porous medium was investigated by Ames and Cobb [21] by applying the methods of linear stability theory and nonlinear energy theory.In the famous book [22], Straughan investigated the onset of convection in a fluid layer with a non-uniform heat source.The author considered four different types of heat supply functions, and discussed the results by obtaining the values of the critical Rayleigh number of energy theory Ra E and the critical Rayleigh number of linear instability theory Ra L .Most recently, Nandal and Mahajan [23] studied the effect of internal heat source on the onset of convection in a fluid layer saturated porous medium.The authors considered four different types of internal heat supply functions, and discussed the results for stress-free and isothermal boundary conditions.Other important studies related to penetrative convection in porous media due to internal heating include [24][25][26][27][28][29][30][31][32][33].
The presence of a uniform heat source (sink) gives rise to a non-uniform temperature gradient, which induces a variation in the magnetic field.These variations of temperature and magnetic field play an important role in controlling ferroconvection in many practical applications, such as in continuous-operation refrigerators.Thus, the study of penetrative convection due to purely internal heating in ferrofluids becomes important.The effect of a uniform heat source on the onset of convection in a magnetic fluid layer under the presence of a transverse applied magnetic field was investigated by Rudraiah and Sekhar [34].Nanjundappa et al. [35] studied the onset of penetrative Benard-Marangoni convection in a horizontal ferromagnetic fluid layer heated from below in the presence of an internal heat source and a uniform vertical magnetic field.In the context of porous medium, the problem of penetrative ferroconvection in a porous layer in the presence of a uniform applied vertical magnetic field was investigated by Nanjundappa et al. [36] by using the internal heating model and the Brinkman-extended Darcy equation.In addition to these studies, some of the related works include [37][38][39].In recent years, the problem of convection in nanofluids induced by internal heating has attracted a great amount of attention from many researchers due to its various industrial and geophysical applications.The effect of an internal heat source in porous medium saturated by nanofluid was first investigated by Yadav et al. [40].There are various studies available in which phenomena associated with the onset of convection in a nanofluid layer induced by purely internal heating have been investigated under different aspects [41][42][43].
In this study, penetrative convection due to purely internal heating in a ferrofluid-saturated porous layer is investigated under the presence of an applied magnetic field.Following Straughan [22], the following four cases are considered for internal heating: (A) when the heat supply function is constant; (B) when the heat supply function is increasing across the layer; (C) when the heat supply function is decreasing across the layer; (D) when the heat supply function heats and cools the layer non-uniformly.The results are discussed graphically for the boundary condition when both boundaries are impermeable and conducting (IMP LU & CON LU ) [44].The effects of Langevin parameter α L , permeability K, the width of ferrofluid layer d, and the nonlinearity of magnetization M 3 are analyzed at the onset of penetrative convection for water-and ester-based ferrofluids.To the best of our knowledge, this problem has not yet been examined.

Formulation
Consider an infinite horizontal layer of incompressible ferrofluid saturating a low permeability porous medium.The fluid is assumed to occupy the layer z ∈ [0, d] with the gravity g acting along the vertical axis (z-axis) in the negative direction.A magnetic field, H = H ext 0 k also acts along the vertical axis from outside the layer, as shown in Figure 1.

Incompressible ferrofluid saturating
a porous layer The convection is driven due to the presence of an internal heat supply function of strength Q.The presence of the internal heat generation in the system may be due to the absorption of external radiation, as well as the heat generation due to a viscous dissipation [38].The following four different types of heat supply functions are considered: where Q 0 is a constant.The current forms for Q(z) are chosen in such a way that the non-dimensional form of Q(z), say Q * (z * ), is such that In case A, Q is constant.For Q 0 > 0, Q increases across the layer in case B and decreases across the layer in case C. In case D, heat supply function Q heats and cools the layer in a non-uniform way [22].
Following [5,36,44] and using the Boussinesq approximation for an incompressible ferrofluid, the equations of continuity, momentum, temperature, and Maxwell (in the magnetostatic limit) can be written as (1) (ρc where u is the filter velocity, ρ f is the fluid density, is the porosity, t is the time, p is the pressure, µ is the viscosity, K is the medium permeability, µ 0 is the magnetic permeability of vacuum, M is the magnetization, α is the thermal volumetric expansion coefficient, T 0 is the reference temperature, (ρc) m is the effective heat capacity of porous medium, T is the ferrofluid temperature, c f is the ferrofluid specific heat, k 1 is the ferrofluid thermal conductivity, B is the magnetic induction, and µ 0 is the magnetic permeability of the vacuum.We write u = (ϑ, v, w).At equilibrium, the magnetization is aligned with the stationary magnetic field, and is a function of magnetic field and temperature.It is assumed to be governed by Langevin formula such as [5] where , here α L , M s , m and k B are the Langevin parameter, magnetic saturation, magnetic moment of a single particle, and Boltzmann constant, respectively.In order to obtain the steady state solution, following Finlayson [5], the magnetization equation M eq is firstly linearized in the following form χ 2 is the chord magnetic susceptibility, H 0 is the uniform magnetic field of the ferrofluid layer when placed in an external magnetic field H = H ext 0 k, χ is tangent magnetic susceptibility, and K m = χH 0 T 0 represents a function of temperature and magnetic field.Using the above expression for M eq , Equation ( 5) now becomes For a different value of α L , the tangent magnetic susceptibility χ and chord magnetic susceptibility χ 2 can be calculated using Langevin Formula (5) as [45] ).
The temperature is assumed to take the same values at the boundaries; therefore, the boundary conditions can be expressed as For the magnetic boundary conditions, the normal component of magnetic induction and the tangential component of the magnetic field are assumed to be continuous across the boundary.
We introduce the following dimensionless variables: Thus, Equations ( 1)-( 6) take the form (after dropping the superscript *) where In Equations ( 9)-( 12), the following non-dimensional parameters have been introduced: Here P r is the Prandtl number, D a is the Darcy number, V a is the Vadasz number, Ra is the internal Rayleigh number, and M 1 , M 2 are the magnetic parameters.

The Steady Solution
The set of Equations ( 9)-( 12) possesses a steady-state solution of the form Using Equation ( 14), the set of Equations ( 9)-( 12) becomes Solving Equations ( 16)-( 18) with the boundary conditions (13), the basic state solutions can be obtained as follows: where

The Linear Stability Problem
The system is now perturbed by considering a small perturbation of amplitude (0 < 1) to the steady state Equation ( 14).This gives After substituting the perturbed variables into the set of Equations ( 9)- (12), linearizing them about the steady state, dropping primes and collecting the O( ) provides Here represents the measure of nonlinearity of magnetization.Since ∇ × H = 0, there exists a potential function ψ such that H = ∇ψ.In the above set of Equations ( 20)-( 22), Equation (20) has been obtained by taking the vertical component of (∇ × ∇×) of the linearized momentum equation.The Chebyshev pseudospectral-QZ method is now applied to solve the set of Equations ( 20)- (22).In this procedure, to match the domain of Chebyshev pseudospectral-QZ method, the present domain is re-adjusted from [0, 1] to [−1, 1] by applying a coordinate transformation from z to 2z − 1.In order to perform the normal mode analysis, all the perturbation quantities w, θ, ψ are assumed in the form {w, θ, ψ} = {w(z), θ(z), ψ(z)} exp{σt + i(k x x + k y y)}, where k x and k y are the wave number in x and y directions, respectively.Substituting (23) into the set of Equations ( 20)- (22) gives where D = d dz and k = k 2 x + k 2 y is the wave number.The boundary conditions become

Method of Solution
In order to solve the above system of equations, the method and algorithm of Kaloni and Lou [46] are closely followed.Equations ( 24)-( 26) with the boundary conditions ( 27) are solved by using the Chebyshev pseudospectral method [47].For a given value of the strength of the internal heat source Q, Langevin parameter α L , wave number k, H 0 , d, K, and other physical parameters, the system of equations is solved by using QZ algorithm and EIG function in MATLAB to obtain the leading eigenvalue σ = σ r + iσ i .Here the leading eigenvalue is assumed to be the one for which the real part σ r is largest among the whole set of eigenvalues.Thereafter, by adjusting the value of Q with the use of secant method, that particular value of Q is calculated corresponding to which the real part σ r of the leading eigenvalue σ = σ r + iσ i approaches zero.This procedure generates only a single point in the neutral stability curve.In order to obtain the desired neutral stability curve, this procedure is repeated for different values of the wave number k [48].The critical value of the strength of the internal heat source Q c with the critical wave number k c can be characterized as The function FMINBND of MATLAB-which is a combination of golden section and parabolic method-is used to minimize Equation (28).
To check the validity of our method, we have solved the same problem in the absence of a magnetic field for Case A by considering IMP LU & CON LU boundary condition.From Table 1 it is clear that results obtained using the above-mentioned method of solution agree well with the results given by Nouri-Borujerdi et al. [49] and Nield and Kuznetsov [32].

Numerical Results and Discussion
In this section, the numerical results are presented for water-and ester-based ferrofluids.The value of the permeability parameter K and porosity are taken as 2.0 × 10 −7 m 2 and 0.35, respectively.The values of the other physical quantities are taken from Kaloni and Lou [46] (Table 1, p. 7) and Rosensweig [45] (Table 2.4, p. 71).The results are discussed for a 1 mm (0.0001 m) thick layer of ferrofluid.
In order to see the effect of magnetic field, neutral curves associated with three different values of the Langevin parameter α L are plotted in Figure 2 for four different cases viz.Case A, Case B, Case C, Case D. The curves are drawn for water-and ester-based ferrofluids by considering the IMP LU & CON LU boundary condition.It can be seen from the figures that as the value of α L increases, the value of the internal Rayleigh number Ra c also increases.The forces arise due to the presence of an internal heat source dominating the magnetic forces that exist because of the presence of an applied magnetic field, and therefore as the strength of the magnetic field increases, disturbances in the ferrofluid slow down, leading to the larger value of Ra c .Thus, the Langevin parameter delays the onset of convection in all four cases.It is also observed that the ester-based ferrofluids are more resilient to convection than the water-based ferrofluids.In Figure 3, we have plotted the neutral curves for three different values of the width of ferrofluid layer d.As d increases, the value of the internal Rayleigh number Ra c increases, and the neutral curves shift upward.A higher value of d suppresses the random motion of particles in a ferrofluid, which results in the higher value of Ra c in all four cases.Thus, the parameter d has a stabilizing effect on the system.Figure 4 shows the neutral curves for three different values of the permeability parameter K for water-and ester-based ferrofluid.The plots indicate that the value of the internal Rayleigh number Ra c increases as the value of K increases.Thus, the permeability parameter K delays the onset of convection.It can also be seen from the figure that the behavior of the parameter K remains the same in all four cases.In Figure 5, the effect of the nonlinearity of the fluid magnetization is examined at the onset of penetrative convection in all four cases.The neutral curves show that the value of the internal Rayleigh number Ra c decreases as the value of M 3 (the measure of nonlinearity of magnetization) increases.It is also observed from the figure that as the value of M 3 increases, the value of critical wave number k c decreases.Thus, the parameter M 3 advances the onset of convection and also shows its influence via expansion of the convection cell size.This behavior is exactly the same as reported earlier by Nanjundappa et al. [38] in the case of uniform heating.We would like to mention here that we have solved the same problem for the combination of other two boundary conditions viz. when both boundaries are conducting with lower boundary impermeable and free upper boundary (IMP L , CON LU & FRE U ); and when both boundaries are impermeable with lower boundary conducting and upper with constant heat flux (IMP LU , CON L & CHF U ) [44].Since the results of these boundary conditions are qualitatively the same, all the results are not included in this work.The boundary conditions in mathematical terms are taken as The subscripts L and U refer to lower and upper boundaries, respectively.Here L L , L U are Biot numbers taking limit value zero for constant heat flux across boundary and infinite at conducting boundary; the coefficients K L and K U take discrete values, zero for constant pressure at the surface, and infinite at the impermeable surface.
Table 2 represents the values of the critical wave number k c and the critical internal Rayleigh number Ra c , with variation in the value of α L for three different boundary conditions.Our first observation from the table is that the value of critical internal Rayleigh number Ra c is higher in the case of IMP LU & CON LU boundaries, and is least in the case of IMP LU , CON L & CHF U boundaries.Thus, the system is most stable when both boundaries are impermeable and conducting, and least stable when both boundaries are impermeable with lower boundary conducting and upper with constant heat flux.We have also noticed that in Case A to Case D, the value of the Ra c first decreases as the value of α L increases from 1 to 2 and then it starts increasing with further increase in the value of α L .This kind of behavior was reported earlier by Kaloni and Lou [46] for the magnetic fluids in a non-porous medium.It is also observed from the table that for any given value of α L , the value of the critical internal Rayleigh number Ra c is larger in case B, where the heat supply function is increasing, and Ra c is least in case D, where the heat supply function Q heats and cools the layer in a non-uniform way.A similar behavior was reported by Straughan [22] (p.92).

Conclusions
Penetrative convection is studied for a thin layer of ferrofluid saturating a porous medium under the presence of an applied magnetic field and internal heating.Four different types of heat supply functions are considered in this study.We used the linear stability theory to examine the onset of penetrative convection, and the resulting eigenvalue problem is solved by employing the Chebyshev pseudospectral method.The following conclusions can be drawn from this study: 1.The effect of d, α L , K is to stabilize the system, while the parameter M 3 has a destabilizing effect on the system.2. The system is most stable for IMP LU & CON LU boundary condition, while it is least stable for IMP LU , CON L & CHF U boundary condition.3. The water-based ferrofluids are less stable than the ester-based ferrofluids.4. The value of Ra c is higher in the case when the heat supply function is increasing, while it is the least in the case when the heat supply function heats and cools the layer in a non-uniform way.

Figure 2 .
Figure 2. Neutral curves for different values of the Langevin parameter α L for (a) water-based ferrofluid and (b) ester-based ferrofluid for Case A to Case D. The fixed parameter values are d = 0.001 m, = 0.35, and K = 2.0 × 10 −7 .

Figure 3 .
Figure 3. Neutral curves for different values of the width of ferrofluid layer d for (a) water-based ferrofluid and (b) ester-based ferrofluid for Case A to Case D. The fixed parameter values are α L = 2, = 0.35, and K = 2.0 × 10 −7 .

Figure 4 .
Figure 4. Neutral curves for different values of the permeability parameter K for (a) water-based ferrofluid and (b) ester-based ferrofluid for Case A to Case D. The fixed parameter values are α L = 2, = 0.35, and d = 0.001 m.

Figure 5 .
Figure 5. Neutral curves for different values of the nonlinearity of magnetization parameter M 3 for (a) water-based ferrofluid and (b) ester-based ferrofluid for Case A to Case D. The fixed parameter values are α L = 2, = 0.35, K = 2.0× 10 −7 , and d = 0.001 m.

Table 1 .
Comparison of critical Rayleigh number Ra c and the critical wave number k c for IMP LU & CON LU boundary condition.

Table 2 .
The values of the critical internal Rayleigh number Ra c and the critical wave number k c for water-and ester-based ferrofluids.

IMP LU & CON LU IMP L , CON LU & FRE U IMP LU , CON L & CHF U
T 0 Temperature at the lower and upper surfaces (K)