Thermal convection in a nonlinear non-Newtonian magnetic fluid

We report theoretical and numerical results on thermal convection of a magnetic fluid in a viscoelastic carrier liquid. The viscoelastic properties are described by a general nonlinear viscoelastic model that contains as special cases the standard phenomenological constitutive equations for the stress tensor. In order to explore numerically the system we perform a truncated Galerkin expansion obtaining a generalized Lorenz system with ten modes. We find numerically that the system has stationary, periodic and chaotic regimes. We establish phase diagrams to identify the different dynamical regimes as a function of the Rayleigh number and the viscoelastic material parameters. c © 2015 Elsevier B.V. All rights reserved.


INTRODUCTION
Ferrofluids are magnetic fluids formed by a stable colloidal suspension of magnetic nano particles dispersed in a carrier liquid. Without an applied external magnetic field the orientations of the magnetic moments of the particles are random resulting in a vanishing macroscopic magnetization (magnetic disorder). An external magnetic field, however, easily orients the particle magnetic moments and a large (induced) magnetization is obtained. In the last decades much efforts have been dedicated to the study of the convective mechanisms in ferrofluids. In particular, heat transfer through magnetic fluids has been one of the leading areas of scientific study due to its technological applications [1]. The ferrofluid convection has applications in highpower capacity transformer systems, where the ferrofluid is used as a material in the core as well as a coolant in the transformer. An important application of ferrofluids lies in the biomedicine area where the carrier liquid is blood [2][3][4][5][6] which is known to have also special rheological proprieties [7][8][9]. In addition, when a magnetic field is applied, a ferrofluid can acquire additional rheological properties such as magneto-viscosity, adhesion, and other non-Newtonian behavior [10][11][12][13][14][15][16][17][18][19]. Hence, a detailed study of viscoelastic magnetic fluids is quite important and in order.
Viscoelastic properties of fluids can be described by a constitutive equation, which relates the stress and strain rate tensors. The simplest constitutive equation capable of describing realistically the viscoelastic properties is given by the so-called Oldroyd model [32]. It has been found that, besides the usual stationary convection, also oscillatory states can be obtained at onset. Which type of convection -stationary or oscillatory -appears first, will depend on the values of the rheological parameters. Experimental measurements of oscillatory convection in viscoelastic mixtures were reported by Kolodner [33] in a DNA suspension; and theoretical studies of the convection thresholds for binary Oldroyd mixtures in different types of fluids, can be found in Refs. [34][35][36][37][38]. Recently, studies on stationary and oscillatory convection in Oldroyd magnetic fluids have been done [39][40][41][42]. By a somewhat different approach, a generalized prototype model for non-Newtonian fluids was presented by Pleiner et.al. (PLB) in Ref. [43] to describe the rheology of viscoelastic fluids in terms of the strain tensor. The authors obtained a valid hydrodynamic description of viscoelasticity for arbitrarily large deformations, rotations and flow. Taking the solid limit correctly the structure of the equations is determined completely. In addition, typical viscoelastic models, such as Maxwell, Oldroyd, Giesekus, etc. are contained in the PLB model in certain limits. We use this prototype model and focus on the nonlinear system in the case of two dimensional solutions (roll patterns). We show that the system can exhibit stationary, oscillatory and chaotic regimes depending on the control parameters.

BASIC EQUATIONS
We consider an (infinte) horizontal layer of thickness d of an incompressible magnetic ferrofluid with a viscoelastic carrier liquid in a vertical gravitational field g = gẑ and subject to a positive temperature difference β across the layer (heating from below). An external magnetic field H 0 is assumed along the vertical direction. Within the Boussinesq approximation, the dimensionless equations for the perturbations from the convection-free, heat conducting state can be written as Here {v, θ, φ} are the dimensionless velocity, temperature, and magnetic potential perturbation, respectively. The perturbation equation for U ij reads This implies for the stress tensor in Eq. (2) In Eqs. (1)- (7), the following dimensionless numbers have been introduced, the Rayleigh number, Ra ∼ β, accounting for buoyancy effects due to the heating; the Prandtl number, P , relating viscous and thermal diffusion time scales; the strength of the magnetic force relative to buoyancy, measured by the parameter M 1 ∼ H 2 0 ; the nonlinearity of the magnetization, M 3 − 1 ∼ H 2 0 , a measure for the deviation of the magnetization curve from the linear behavior; the linear and nonlinear relaxation times Γ 1,2 and elastic moduli E 1,2 , and the nonlinear viscosity Z, which are all positive.
On the linear level (Γ 2 = E 2 = Z = 0) it is easy to recover the standard (linear) Oldroyd model containing the Deborah number Γ and the retardation number, Λ. Both descriptions are equivalent with Γ = Γ 1 and Λ = (1 + E 1 Γ 1 ) −1 , revealing however that Λ is restricted to 0 < Λ < 1. The kinematic viscosity ν 1 , used to scale the time in the viscoelastic description, is related to the

TWO DIMENSIONAL SYSTEM
In this section we solve numerically the system of equations using a truncated Galerkin method. For the sake of simplicity, the analysis is limited to twodimensional flows, v = {−∂ z ψ, 0, ∂ x ψ} introducing the stream function ψ. In particular, we assume periodicity with wave number k in the lateral direction, x, describing 2-dimensional convection rolls parallel to the y-axis. By the definition and symmetry of U ij , we only have three components in the 2-dimensional case.
We impose idealized thermal and magnetic boundary conditions [42] and assume rigid boundary conditions for the stream function and the strain tensor components at z = {0, 1}. For the numerical solution we will restrict ourselves to the fundamental mode in the lateral direction, neglecting higher harmonics, while in the z-direction across the layer a multimode description will be used that fulfills the boundary conditions, automatically [44,45].
Substituting the trial functions into the dynamic equations and multiplying , respectively, at r = 2.5. Other fixed parameters as in Fig. 1 these equations by the orthogonal eigenfunctions and integrating over a convection cell, yields a set of ten ordinary differential equations for the time evolution of the amplitudes. They are solved via a standard fourth order Runge-Kutta integration scheme with a fixed time step dt = 0.01 guaranteeing a precision of 10 −8 for the amplitudes. For each set of parameters we evolve the numerical solutions at 10 6 time steps in order to avoid observing purely transient phenomena. This system is a generalization of the Lorenz system, hence we expect that the system can exhibit complex behaviors.
NUMERICAL RESULTS Fig. 1 shows the normalized stream function amplitude, X 1 , as a function of time τ for three different values of Ra. After a transient, X 1 tends to a stationary value, which is zero for r < 1 and, for r > 1, finite and increasing with increasing r. Here, r is the Ra number normalized by the stationary threshold value Ra sc [42]. In Fig. 2 the magnetic field dependence of the saturation amplitude value X 1,sat = X 1 (t → ∞) at r = 2.5 is shown for three different values of M 3 . The amplitude increases with a power law depends on r and the rest of the material parameters. Obviously, a magnetic field has a destabilizing effect increasing the saturation amplitude. Figures 3 and 4 show the system in a periodic state, which is obtained for higher r. In Fig. 3 the time dependence of the normalized stream function, X 1 , is shown in the top frame, while the bottom frame shows the corresponding normalized Fourier power spectrum. There are odd higher harmonic peaks in the spectrum, Ω n = Ω 0 (2n + 1), such that Ω 1 is highest and those for n ≥ 2 are unimportant. A 3D phase portrait of {X 1 , X 3 , X 4 } for this state is shown in Fig. 4. It describes a stable non-symmetric X 2 1 X 2 3 orbit in the Sparrow notation [46]. Figures 5 and 6 show the system in a different periodic state, which is obtained for rather low r values, but for higher elastic moduli, E 1 and E 2 . Again, the time dependence of the normalized stream function, X 1 , is shown in the top frame and its corresponding normalized Fourier power spectrum in the bottom frame. Here, the fundamental peak is the most important one, while the odd higher harmonic peaks are pairwise reduced in height. The appropriate 3D phase portrait of {X 1 , X 3 , X 4 } in Fig. 6 shows a non-symmetric homoclinic orbit. Clearly, this periodic state at high elastic moduli is different from the one obtained at a high r number. Figures 7 and 8 show the system in a chaotic regime at a rather low r value. The aperiodic time behavior is manifest in the upper part of Fig. 7, while the lower part reveals the continuos nature of the corresponding Fourier power spectrum indicating the chaotic nature of this state. Figure  8 shows the appropriate 3D phase portrait having the form of a strange attractor similar to the Lorenz attractor. It should be noted that this chaotic state comes at a much lower r value than the periodic one of Figs  identical). This indicates a rather complicated bifurcation behavior as a function of r with alternating chaotic and periodic states. A comprehensive discussion of the complex bifurcation diagram will be given elsewhere.
Here, we will show the bifurcation diagrams as a function of the linear elastic modulus E 1 and of the nonlinear viscosity Z in Figs. 9 and 10, respectively. These diagrams are obtained by calculating the time series of the amplitudes for many different, random initial conditions. Taking the local maximum values within a certain time slot of, e.g., |X 1 (τ )| reveals the difference between a regular (stationary or oscillatory) and a chaotic behavior. in the former case only one or a few maximum values are found, while in the latter one a continuous-like assembly of such maximum values is present. The bifurcation diagram as a function of the elastic modulus E 1 shows chaos to only exist for intermediate values of E 1 , while for small and large ones regular (stationary and oscillatory) behavior is found. Increasing E 1 further the convection amplitude decreases due to the enhanced stiffness of the fluid. In addition, a large nonlinear viscosity (large Z) suppresses chaos and reduces the convection amplitude in an almost discontinuous manner as is shown in Fig. 10.

SUMMARY
In the present work, Rayleigh-Benard convection in a magnetic viscoelastic liquid is studied. For the viscoelastic properties the PLB model [43] is used. Similar to the Lorenz approach [44], a set of ten coupled linear ordinary differential equations is obtained. In the case of the stationary bifurcation we observe that the absolute value of the stream function amplitude (related to the velocity) increases with the external field and is independent of the viscoelastic properties. Due to the viscoelastic properties the system also has an oscillatory bifurcation. We discuss two dif- ferent periodic states, one at high Rayleigh number and low elastic moduli and vice versa. They are rather different, in particular in their nonlinear, higher harmonic behavior. In addition, chaotic behavior is found. It can already occur at rather small Rayleigh numbers, but period states at higher Rayleigh numbers occur in-between chaotic ones. The existence range for chaos strongly depends on the material properties, in particular on the viscoelastic and magnetic ones. The chaotic regime has been characterized by using bifurcation diagrams and Fourier power spectrum calculations, as well as phase portraits.
ACKNOWLEDGEMENTS D.L. acknowledges the partial financial support from FONDECYT 1120764, Basal Program Center for Development of Nanoscience and Nanotechnology (CE-