Simulation of Thermomagnetic Convection in a Cavity Using the Lattice Boltzmann Model

Abstract

Thermomagnetic convection in a differentially heated square cavity with an infinitely long third dimension is numerically simulated using the single relaxation time lattice Boltzmann method (LBM). This problem is of considerable interest when dealing with cooling of microelectronic devices, in situations where natural convection does not meet the cooling requirements, and forced convection is not viable due to the difficulties associated with pumping a ferrofluid. Therefore, circulation is achieved by imposing a magnetic field, which is created and controlled by placing a dipole at the bottom of the enclosure. The magnitude of the magnetic force is controlled by changing the electrical current through the dipole. In this study, the effects of combined natural convection and magnetic convection, which is commonly known as “thermomagnetic convection,” are analysed in terms of the flow modes and heat transfer characteristics of a magnetic fluid.

1. Introduction

Natural convection in enclosures is of interest in many engineering applications, such as the cooling of electronic equipment and solar energy collection, and it has been extensively studied experimentally and numerically by several researchers [1–4]. The enhancement of convective heat transfer in these devices is an active area of research. For small length-scale applications, such as the cooling of microelectronic devices, thermomagnetic convection has been identified as a viable approach for augmenting and controlling the convective heat transfer. This involves controlling the motion of colloidal suspensions containing magnetic nanoparticles, known as ferrofluids, using external magnetic fields [5, 6].

The study of flow and heat transfer characteristics of natural convection of a magnetic fluid has received considerable attention lately. Finlayson [7] in 1970 explained how an external magnetic field imposed on a ferrofluid with varying magnetic susceptibility, for example, due to a temperature gradient, results in a nonuniform magnetic body force, which leads to thermomagnetic convection. This phenomenon is complex in nature, and so a thorough understanding of the relation between an applied magnetic field and the resulting heat transfer is necessary for the proper design and control of thermomagnetic devices [8].

The numerical studies of heat transfer with magnetic fluids reported in the literature are relatively scarce and invariably use as governing equations the Navier-Stokes and energy equations. In what follows, relevant work in the field is succinctly reviewed. Yamaguchi et al. [9] experimentally and numerically studied the heat transfer characteristics and flow behavior for natural convection of a magnetic fluid in a rectangular cavity under an imposed uniform vertical magnetic field. Al-Najem et al. [10] numerically calculated the flow and temperature fields under uniform magnetic field in a tilted square cavity with isothermal vertical walls and adiabatic horizontal walls at a low Prandtl number of 0.71 and demonstrated that the suppression effect of the magnetic field on convection currents and heat transfer is more significant for low inclination angles and high Rayleigh numbers. Krakov and Nikiforov [11] determined the influence of the angle between the direction of the temperature gradient and that of a uniform magnetic field on the convection structure and intensity of heat flux in a square cavity.

A few other researchers have considered spatially nonuniform magnetic fields for the experimental or numerical investigations but did not fully describe the variations in those fields [12–15]. For example, Sawada et al. [12] and Kikura et al. [13] experimentally investigated the influence of a varying magnetic field in a cubical enclosure and in concentric horizontal annuli, respectively, but the field gradient was only measured in a single direction. Moreover, the orientation of the magnetic field was not described. Hadavand et al. [16] simulated thermomagnetic convection in a cavity which was in close proximity to a line dipole. The imposed nonuniform magnetic field was considered only in one direction and the effect of temperature on the magnetic susceptibility was ignored. Ganguly et al. [8] performed a study to simulate thermomagnetic convection in a square enclosure by introducing a line dipole using a Navier-Stokes solver. The effect of temperature on the magnetic susceptibility was considered; because they used a very strong nonuniform magnetic field, the buoyancy was neglected in their simulation. They found that the average Nusselt number on the wall increases with increasing the magnetic dipole strength and temperature but decreases by increasing the fluid viscosity.

The magnetic control has enabled numerous developments dealing with technical and medical applications, and these suspensions are already being considered the next generation heat transfer fluids as they offer the possibility of achieving heat transfer rates much higher than those of conventional heat transfer fluids and fluids containing microsized metallic particles [17].

In the present work, the objective is to numerically study the effects of combined natural convection and magnetic convection in a square enclosure using the lattice Boltzmann method. The square enclosure is under the influence of an imposed nonuniform two-dimensional magnetic field that satisfies Maxwell’s equations. The variation of the ferrofluid susceptibility with the temperature difference is considered in this study.

2. Physical Model

In the absence of an external magnetic field, the nanoparticles are oriented randomly in the ferrofluid; once the external magnetic field is applied, the nanoparticles align with the applied magnetic field. The generated force between the magnetic field and the homogenously distributed magnetic nanoparticles enables the manipulation of the ferrofluid by adjusting the applied external magnetic field. This uniformly distributed magnetic force manifests itself as a body force analogous to the gravity field. These suspensions, in general, exhibit normal liquid behavior coupled with super paramagnetic properties. This leads to the possibility of controlling the properties and the flow of these liquids with relatively moderate magnetic field strengths. Each solid magnetic nanoparticle in the ferrofluid is covered with a dispersant material in order to keep the particles separated, and, in this way, preventing agglomeration of the particles due to gravity, coagulation, or an external magnetic field.

Figure 1 presents a schematic of the enclosure used for the numerical simulation, which has height D and its vertical walls are at constant, but different, temperatures. The third dimension is considered to be infinite; therefore, the flow in this configuration can be assumed to be two-dimensional. The left and right vertical walls are at constant high temperature, T_h, and constant low temperature, T_c, respectively. The upper and lower walls are adiabatic. A line dipole, which provides the external magnetic field, is placed adjacent to the lower adiabatic wall at a distance d. The magnetic field is governed by the Maxwell relations in static form [18], as follows:∇⋅𝐵=0,∇×𝐻=0, (2.1) where 𝐵 is the magnetic field inside the ferrofluid due to the line dipole which can be expressed as follows [8]:𝐵=𝜇₀1+𝜒𝑚sin𝜃𝑟²𝑒_𝑟−sin𝜃𝑟²𝑒_𝜃, (2.2) where 𝐻 and 𝐵 are related through the following equation:1𝐻=𝜇₀1+𝜒𝐵, (2.3) where 𝜇₀ is the magnetic permeability of free space, 𝜇₀=4𝜋×10⁻⁷(N/A²); m denotes the magnetic dipole moment of the electromagnet coil (having two parallel counter-conductors separated by a distance 𝑏 and carrying current 𝐼) per unit length (𝑚=𝐼𝑏/2𝜋). The magnetic fluid susceptibility, 𝜒, varies with temperature as follows:𝜒𝜒=⁰1+𝛽𝑇−𝑇₀, (2.4) where 𝑇₀ is the reference temperature and 𝜒₀ is the magnetic field susceptibility at reference temperature.

Figure 1 

Schematic of the cavity and position of the line dipole magnet.

The magnetization of the ferrofluid occurs when it becomes polarized in the presence of an external magnetic field. Consequently, attractive forces acting on each particle are produced due to the interaction of the ferrofluid polarization and the external magnetic field. The attractive force acting on the particles can then be treated as a body force acting on the ferrofluid, which is dealt with as a homogenous medium. The attractive force on the ferrofluid per unit volume is 𝐹=(𝑀⋅∇)𝐵, (2.5) where 𝑀 is the magnetization vector and is related to 𝐻 by 𝑀=𝜒𝐻. The magnetic field strength diminishes with increasing distance from the dipole. The Kelvin body force, see (2.5), can be written as follows:𝐹=𝜇₀𝜒()11+𝜒𝐻⋅∇𝐻=2𝜇₀𝜒1+𝜒∇(𝐻⋅𝐻)+𝜇₀𝜒𝐻𝐻⋅∇𝜒. (2.6)

The first term on the right-hand side of (2.6) has a same nature as the pressure term in the momentum equation. The second term becomes important when there is a spatial gradient in the fluid susceptibility. In the absence of a temperature gradient, this term is zero. The Kelvin body force is simplified in another study [19] to the following format by defining an effective pressure as 𝑃^∗𝜇=𝑃−⁰𝜒₀2𝐻², (2.7) and it yields 1𝐹=2𝜇₀𝜒₀1−𝛽𝑇−𝑇₀∇(𝐻⋅𝐻)+𝜇₀𝜒²₀𝛽(𝐻⋅∇𝑇)𝐻. (2.8)

In this work, the single relaxation time lattice Boltzmann method is used to simulate the thermomagnetic convection in a two-dimensional square enclosure.

3. Numerical Model

The lattice Boltzmann method (LBM) is a discrete particle-based mesoscopic approach for viscous fluid flow simulation. The LBM has been shown to be an efficient tool for flow simulation in complex geometries, when compared with conventional fluid dynamics approaches [20]. In this method, finite volumes of fluid are presented as discrete particles and the fluid dynamics is obtained based on a probability distribution function of the particles and a discrete set of prescribed velocities.

The thermal lattice Boltzmann model, which is employed in the present work, is a double-distribution-function thermal model. Two distinct distribution functions for the flow and temperature fields are 𝑓 and 𝑔, respectively, where 𝑓 is the density (or fluid particles) distribution function and 𝑔 is the temperature distribution function. This present thermal lattice Boltzmann model uses the BGK approximation [21] and it is based on the work of He et al. [22] and Shi et al. [23]:𝑓_𝑖𝑥+𝑐_𝑖𝛿𝑡,𝑡+𝛿𝑡−𝑓_𝑖(𝑥,𝑡)=−𝛿𝑡𝜏_𝜈𝑓_𝑖(𝑥,𝑡)−𝑓^eq_𝑖(𝑥,𝑡)+𝑓_𝑔𝑔+𝐹,_𝑖𝑥+𝑐_𝑖𝛿𝑡,𝑡+𝛿𝑡−𝑔_𝑖(𝑥,𝑡)=−𝛿𝑡𝜏_𝑐𝑔_𝑖(𝑥,𝑡)−𝑔^eq_𝑖(𝑥,𝑡), (3.1) where 𝛿𝑥 and 𝛿𝑡 are the lattice grid spacing and the lattice time step, respectively. These two are related to each other through the streaming speed, 𝑐, as 𝑐=𝛿𝑥/𝛿𝑡. However, in most of the LBM simulations on a uniform lattice, including those of the present work, 𝛿𝑥 and 𝛿𝑡 for simplicity are considered to be equal to one, and 𝑓_𝑔 is the buoyant body force term, which is formulated by using the Boussinesq approximation [24], namely:𝑓_𝑔=3𝜔_𝑖𝛽𝜌𝑔_𝑦𝑇−𝑇_ave𝑐_𝑖𝑦, (3.2) where β is the thermal expansion coefficient, 𝑔_𝑦 is the acceleration of gravity acting in the 𝑦-direction of the lattice links, and 𝑐_𝑖𝑦 is the y-component of c_i. The terms 𝜌 and T are the local density and temperature and are calculated at each lattice site. 𝑇_ave is the average temperature of hot and cold wall. The Kelvin body force, 𝐹, is the force that a magnetic fluid experiences in nonuniform magnetic field that can be calculated based on (2.5).

The most common two-dimensional lattice is the D2Q9 lattice, which is used in this study, where “2” denotes the number of space dimensions and “9” refers to the discrete set of nine velocities as proposed by Qian et al. [25]. The D2Q9 lattice and its labelling for the lattice directions are shown in Figure 2. For the D2Q9 lattice, the discrete velocity set is formulated as follows:𝑐_𝑖()𝜋=𝑐cos𝑖−12()𝜋,sin𝑖−12𝑐,𝑖=1,2,3,4,_𝑖=√()𝜋2𝑐cos2𝑖−14()𝜋,sin2𝑖−14𝑐,𝑖=5,6,7,8,₉=0. (3.3)

Figure 2 

Topology of the D2Q9 lattice and labelling of the directions.

The corresponding equilibrium distribution functions for the density and temperature distribution functions are 𝑓^eq_𝑖 and 𝑔^eq_𝑖, respectively. These two parameters are in the form of a second-order truncated expansion of the Maxwell-Boltzmann equilibrium function:𝑓^eq_𝑖=𝜌𝜔_𝑖11+𝑐²_𝑠𝑐_𝑖1⋅𝑢+2𝑐⁴_𝑠𝑐_𝑖⋅𝑢²−12𝑐²_𝑠,𝑔𝑢⋅𝑢_eq𝑖=𝜔_𝑖𝑇11+𝑐²_𝑠𝑐_𝑖1⋅𝑢+2𝑐⁴_𝑠𝑐_𝑖⋅𝑢²−12𝑐²_𝑠𝑢⋅𝑢, (3.4) where c_s is the lattice sound speed and is equal to 1/√3 for a two-dimensional lattice with 9 velocities (D2Q9). 𝜔_𝑖 is the weighting factor and is defined as follows for the D2Q9 lattice:𝜔_𝑖=491,𝑖=9,91,𝑖=1,2,3,4,36,𝑖=5,6,7,8. (3.5)

The macroscopic density (𝜌), pressure (𝑃), velocity (𝑢), and temperature (𝑇) on each lattice site are calculated as follows:𝜌=_𝑖𝑓_𝑖,𝑃=𝜌𝑐²_𝑠,𝜌𝑢_𝑖=_𝑖𝑐_𝑖,𝑇=_𝑖𝑔_𝑖. (3.6)

The kinematic viscosity and the thermal diffusivity are given by𝜈=𝜏_𝜈−12𝑐²_𝑠,𝛼=𝜏_𝑐−12𝑐²_𝑠. (3.7)

The no-slip solid-wall boundary condition is set for the vertical and horizontal walls of the enclosure by using the bounce back method. As mentioned before, the top and bottom walls are adiabatic, and the vertical walls are isothermal, with the left wall at higher temperature than that of the right wall.

4. Results and Discussion

The physical domain considered in this study is a two-dimensional square cavity. The parameters D and d are depicted in Figure 1, where D is set to 2 mm and d is equal to 𝐷/2. The density and dynamic viscosity of the ferrofluid at 300 K were set equal to 1180 kg/m³ and 1.0×10⁻³ kg/ms, respectively. The Prandtl number, Pr, and fluid compressibility, 𝛽, are taken as 5.5 and 5.6×10⁻⁴K⁻¹, respectively.

A uniform lattice system is employed for the simulations. The numerical solution should be independent of the lattice size, and to this purpose, three different lattice sizes of 81×81, 101×101, and 121×121 were examined; results with a lattice density of 101×101 lattice units were found to be adequate to establish lattice-independent solutions for the range of the parameters used in the present study. In order to validate the thermal LBM predictions, laminar natural convection of air, Pr=0.71, in a cavity was simulated for Rayleigh numbers of 10⁴ and 10⁵. The predicted natural convection results with the 101×101 grid system for Ra=10⁴ and 10⁵, respectively, were in excellent agreement with the benchmark steady flow patterns and temperature fields reported by de Vahl Davis [26]. In Table 1, the predicted average Nusselt numbers on the hot wall are compared against the values reported by de Vahl Davis [26] for the same conditions and, as it can be observed, are in excellent agreement.

Table 1 

Comparison of the predicted average Nusselt number for two different Rayleigh numbers against the values reported in [26].

The average Nusselt number on the side walls is determined from the following relation:1Nu=−Δ𝑇^𝐷1𝑑𝑇𝑑𝑥. (4.1)

The two dimensionless parameters which characterize the natural convective flow, Rayleigh number (Ra) and the Prandtl number (Pr), are defined as follows:𝑔Ra=^𝑦𝛽Δ𝑇𝐷³,𝜈𝜈𝛼Pr=𝛼, (4.2) where g_y, β, ΔT, and D are the gravitational acceleration, working fluid thermal expansion coefficient, temperature difference between the two vertical walls, and the characteristic length (equal to the width of the cavity), respectively.

Presenting the magnetic dipole will disturb the flow motion in the cavity; consequently, the temperature distribution will be affected along with the heat transfer rate at the walls of the cavity. The isotherms for three different values of m equal to 0.0, 0.05, and 0.1 Am, are shown in Figures 3(a)–3(c) for Δ𝑇=30K (Ra=10⁴). The influence of the induced magnetic field on the temperature distributions is apparent. If the magnetic field magnitude is relatively weak, the buoyancy effect is dominant and the isotherms and streamlines are similar to those of pure natural convection in the absence of a magnetic field. For a stronger magnetic field the isotherm stratification in the core will change, as shown in Figure 3(c); this is due to the suppression of the free convective flow by the magnetic field. By increasing the magnetic field strength, the effect of the natural convection is eventually eliminated and the flow motion of the ferrofluid is governed by the strength of the magnetic field.

(a)

(b)

(c)

Figure 3 

Isotherms for the thermomagnetic flow in a square cavity for Δ𝑇=30K, Ra=104: (a) 𝑚=0.0, (b) 𝑚=0.05, and (c) 𝑚=0.1Am.

Figures 4(a)–4(c) report the isotherms at Δ𝑇=60K (Ra=2×10⁴) for 𝑚=0.0, 0.05, and 0.1 Am. From these figures, it can be observed that for large values of Δ𝑇, the effect of increasing the magnetic strength on the temperature distribution in the enclosure is more prominent compared to that in the case of a smaller Δ𝑇 (Ra=10⁴, Figure 3). This is because the magnetic susceptibility depends on temperature and influences the motion of the ferrofluid in the cavity. In these circumstances, the colder fluid close to the lower temperature wall moves down towards the region of larger field strength due to its higher magnetic susceptibility, while the warmer fluid with lower susceptibility is displaced away from the line dipole and is pushed upward along the hot wall.

(a)

(b)

(c)

Figure 4 

Isotherms for the thermomagnetic flow in a square cavity for Δ𝑇=60K, Ra=2×104: (a) 𝑚=0.0, (b) 𝑚=0.05, and (c) 𝑚=0.1Am.

For sufficiently large values of 𝑚, the natural convection is suppressed and the magnetic field effect is the dominant factor. As a result, the isotherms will be almost independent of Ra number for a relatively high magnetic field strength. The streamlines for the thermomagnetic flow in the square cavity for Δ𝑇=30K and Δ𝑇=60K and for different values of 𝑚, 0.0, 0.05, and 0.1 Am, are presented in Figures 5(a)–5(c) and 6(a)–6(c), respectively.

(a)

(b)

(c)

Figure 5 

Streamlines for the thermomagnetic flow in a square cavity for Δ𝑇=30K, Ra=104: (a) 𝑚=0.0, (b) 𝑚=0.05, and (c) 𝑚=0.1Am.

(a)

(b)

(c)

Figure 6 

Streamlines for the thermomagnetic flow in a square cavity Δ𝑇=60K, Ra=2×104: (a) 𝑚=0.0, (b) 𝑚=0.05, and (c) 𝑚=0.1Am.

Figure 7 shows the predicted values of the average Nusselt number on the hot wall for a range of temperature difference, 15≤Δ𝑇≤90, 5×10³≤Ra≤3×10⁴, and 0≤𝑚≤0.2. It can be observed that the average Nusselt number increases with both magnetic field strength and temperature difference. The average Nusselt number increases more than double for 𝑚=0.2Am compared to the Nusselt number for free convection only.

Figure 7 

Average Nusslet number for different temperature difference and magnetic field strength.

To show the effect of magnetic field strength on heat transfer enhancement, the buoyancy force was turned off and the average Nusselt numbers on the hot wall were compared with those from the combination of buoyancy and magnetic force. In Figure 8, the average Nusselt numbers are presented for different magnetic field strengths and a range of temperature difference. For magnetic field moment, 𝑚, lower than 0.1 Am the buoyancy force has a noticeable effect on the Nusselt number, but with increasing magnetic field strength the effect of the buoyancy force is reduced and the Nusselt number can be calculated by considering only the magnetic force for 𝑚 larger than 0.15 Am.

Figure 8 

Average Nusslet number based on the magnetic force and combination of magnetic and buoyancy force. (The solid lines are the Nuave based on the combination of the magnetic and buoyancy force. The marked lines are calculated based on the magnetic force only.)

5. Conclusion

Simulations for thermomagnetic convection were conducted for a range of Rayleigh numbers and magnetic field strengths in the laminar regime. The effects of these parameters on the heat transfer characteristics were analyzed; the numerical results showed that with the application of an external magnetic field, the temperature and velocity fields were significantly modified.

For weak magnetic fields, both buoyancy force and magnetic force have effects on the heat transfer. However, for sufficiently high magnetic field strength, the magnetic convection is suppressed for all the examined values of the Rayleigh number. This study shows that even for relatively weak magnetic fields, the increase in heat transfer for small scale devices is considerable.

Nomenclature

𝑐𝑖:	Particle discrete velocity set
𝑐𝑠:	Speed of sound
𝑑:	Distance between dipole and enclosure
𝐷:	Enclosure height
𝑓:	Distribution function for the flow field
𝑔:	Distribution function for the temperature field
𝑔𝑦:	Acceleration of gravity in the 𝑦-direction
𝐻:	Defined in  (2.3)
𝑚:	Magnetic field moment
𝑀:	Magnetization
Nu:	Nusselt number
Pr:	Prandtl number
Ra:	Rayleigh number
𝑡:	Time
𝑇:	Temperature
𝑥, 𝑦:	Cartesian coordinate
𝑢:	Velocity component.

Greek Letters

𝛼:	Thermal diffusivity
𝛽:	Fluid compressibility
𝜒:	Magnetic susceptibility
𝜇:	Viscosity
𝜇0:	Magnetic permeability
𝜌:	Density
𝜏:	Lattice relaxation time
𝜈:	Fluid kinematic viscosity
𝜔𝑖:	Lattice weighting factor.

Acknowledgments

The authors acknowledge the support received from Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant 12875 (ACMS) and Foundation for Science and Technology-Portugal (FCT) Grant PTDC/EME-MFE/105031/2008 (ACMS). M. Hadavand would like to thank Dr. R. Ganguly for helpful discussions.