Jourmal of Porous Media,18(11):1059-1074(2015) NUMERICAL SIMULATION OF THERMOMAGNETIC CONVECTION OF AIR IN A POROUS CUBIC ENCLOSURE WITH AN ELECTRIC COIL INCLINED IN GENERAL ORIENTATIONS USING A LTNE MODEL Changwei Jiang,'Er Shi,Xianfeng Zhu,Nan Xie School of Energy and Power Engineering,Changsha University of Science and Technology Changsha 410114,China:and Key Laboratory of Efficient and Clean Energy Utilization, College of Hunan Province,Changsha 410114,China Address all correspondence to Changwei Jiang E-mail:cw jiang@163.com Original Manuscript Submitted:11/9/2014;Final Draft Received:2/25/2015 media can be controlled by the s ence of a gravity field.The Biot-Saoart law is used to calculate the The L ccmfruiamdsoiaa transfer of air can be enhanced or controlled by applying a gradien KEY WORDS:thermomagnetic comection,nmerical simulation,magnetic force,inclined electric coil, porous media 1.INTRODUCTION gated natural convection flows in a vertical annulus filled with a fluid-saturated porous medi um when the inner wall ating Numerical investigati arious thermal boundary conditions was done by ra entific applications,such as in heat exchangers,solar re- makrishna et al.(2013).in which it was found that the hermal vices.ete. ious and heat transfer n porous enclosures (Khalil and Shiv- Choi (2013)theoretically analyzed the onset of convec. akumara,0:Vafai.Pourshaghaghy tion motion in an initially quiescen al and D Ellahi et al.,2013). porous layer transfer in a porous mediumflled right-angled triangular binary fluid enclosure.They found that heat transfer incre ses with Natural convection under the effect of a magnetic field the decreasing aspect ratio.Sankar et al.(2011)investi- is of great importance in many industrial applications. 91-028X/15/S35.by Begell House.Inc. 1059
Journal of Porous Media, 18 (11): 1059–1074 (2015) NUMERICAL SIMULATION OF THERMOMAGNETIC CONVECTION OF AIR IN A POROUS CUBIC ENCLOSURE WITH AN ELECTRIC COIL INCLINED IN GENERAL ORIENTATIONS USING A LTNE MODEL Changwei Jiang,∗ Er Shi, Xianfeng Zhu, & Nan Xie School of Energy and Power Engineering, Changsha University of Science and Technology, Changsha 410114, China; and Key Laboratory of Efficient and Clean Energy Utilization, College of Hunan Province, Changsha 410114, China ∗Address all correspondence to Changwei Jiang E-mail: cw jiang@163.com Original Manuscript Submitted: 11/9/2014; Final Draft Received: 2/25/2015 The natural convection heat transfer of paramagnetic fluid in porous media can be controlled by the gradient magnetic field. Thermomagnetic convection of air in a porous cubic enclosure with an inclined coil around the Y and Z axes is numerically investigated in the presence or absence of a gravity field. The Biot–Savart law is used to calculate the magnetic field. The Darcy model is used to solve the momentum, and the energies for fluid and solid are solved with the local thermal non-equilibrium (LTNE) models. The flow and temperature fields for the thermomagnetic convection of air are presented and the average Nusselt number on the hot wall is calculated and compared. The results show that both the magnetic force and the coil inclination have a significant effect on the flow field and heat transfer in a porous cubic enclosure, and the thermomagnetic convection heat transfer of air can be enhanced or controlled by applying a gradient magnetic field. KEY WORDS: thermomagnetic convection, numerical simulation, magnetic force, inclined electric coil, porous media 1. INTRODUCTION Natural convection in porous media is a subject of prime importance due to its use in recent engineering and scientific applications, such as in heat exchangers, solar receivers, geothermal exploitation, cooling of electronic devices, etc. Various researchers have carried out experimental and numerical investigations to study fluid flow and heat transfer in porous enclosures (Khalili and Shivakumara, 2003; Vafai, 2005; Pourshaghaghy et al., 2007; Bhardwaj and Dalal, 2013; Ellahi et al., 2013). Varol et al. (2006) analyzed the steady-state free-convection heat transfer in a porous medium–filled right-angled triangular enclosure. They found that heat transfer increases with the decreasing aspect ratio. Sankar et al. (2011) investigated natural convection flows in a vertical annulus filled with a fluid-saturated porous medium when the inner wall was subject to discrete heating. Numerical investigation of natural convection within porous square enclosures for various thermal boundary conditions was done by Ramakrishna et al. (2013), in which it was found that the thermal boundary conditions have an important influence on the flow and heat transfer characteristics during natural convection within porous square cavities. Kim and Choi (2013) theoretically analyzed the onset of convection motion in an initially quiescent, horizontal isotropic porous layer. Alloui and Vasseur (2013) investigated natural convection in a shallow porous cavity filled with a binary fluid. Natural convection under the effect of a magnetic field is of great importance in many industrial applications, 1091–028X/15/$35.00 ⃝c 2015 by Begell House, Inc. 1059
1060 Jiang et al NOMENCLATURE B rotation angle around reference rotation angle around the axis (=ui/L)(T) emal diffusivity (m2 s-1) K-1) Da Darcy number (=K/L2) 2 dimensionless mas netic strength parameter (Xob/umgL) n0 nless fuid tem solid-to-fluid heat transfer coefficient =(T-TT-T】 (Wm-2K-1) 8 dimensionless solid temperature (T 、Te切 -1 wm-1- lengthof a cubic enclosure(m) =kd(-e)k1 luid kinematic viscosity(kg m-s-1) ghenc permea nty (H m Prandtl number (=v/) dimensionless position vector fluid density(kg m- position vector (m) solid density (kg m (m3kg-1) Ra reference volume magnetic susceptibility gmasmegnctcsoeptby 3kg-1) hot wall temperature (K) Subscripts solid temperature(K) old erence temperature [=(Th +T)/2](K) nts ( XYZ dimensionless Cartesian coordinates solid 工,2 Cartesian coordinates 0 reference value such as in solar energy technology,crystal growth,and convection heat transfer of Cu-water nanofuid in a cod liquid metal cooling kets used in Tu on rea sure containing a hot inner sinusoid mental studies on therm netic convection in en al.(2012)using the control volume-based finite-element a method. magnetic field. indicate that transfer en nthe absence of ancement ratio d Journal of Porous Media
1060 Jiang et al. NOMENCLATURE B dimensionless magnetic flux b magnetic flux density (T) b0 reference magnetic flux density, (= µmi/L) (T) C = 1+ (1/T0β) cp fluid specific heat at constant pressure (J kg−1 K −1 ) Da Darcy number (= κ/L2 ) g gravitational acceleration (m s−2 ) H dimensionless solid-to-fluid heat transfer coefficient (= hL2/εkf) h solid-to-fluid heat transfer coefficient (W m−2 K −1 ) i electric current in a coil (A) kf fluid thermal conductivity (W m−1 K −1 ) ks solid thermal conductivity (W m−1 K −1 ) L length of a cubic enclosure (m) Num average Nusselt number P dimensionless pressure p pressure (Pa) Pr Prandtl number (= νf/αf) R dimensionless position vector r position vector (m) r radius of the coil (m) Ra Rayleigh number [= gβ (Th − Tc)L 3/αfνf] S dimensionless tangential element of a coil s tangential element of a coil (m) Tc cold wall temperature (K) Tf fluid temperature (K) Th hot wall temperature (K) Ts solid temperature (K) T0 reference temperature [= (Th + Tc)/2] (K) U,V ,W dimensionless velocity components u, v, w velocity components (ms−1 ) X,Y ,Z dimensionless Cartesian coordinates x, y, z Cartesian coordinates xEuler rotation angle around the X axis (◦ ) yEuler rotation angle around the Y axis (◦ ) zEuler rotation angle around the Z axis (◦ ) Greek Symbols αf fluid thermal diffusivity (m2 s −1 ) β thermal expansion coefficient (K−1 ) γ dimensionless magnetic strength parameter (= χ0b 2 0/µmgL) ε porosity θf dimensionless fluid temperature [= (Tf – T0)/(Th – Tc)] θs dimensionless solid temperature [= (Ts – T0)/(Th – Tc)] κ permeability (m2 ) Λ dimensionless thermal conductivity [= εkf /(1 – ε)ks] µf fluid kinematic viscosity (kg m−1 s −1 ) µm magnetic permeability (H m−1 ) µ0 magnetic permeability of vacuum (H m−1 ) νf fluid dynamic viscosity (m2 s −1 ) ρf fluid density (kg m−3 ) ρs solid density (kg m−3 ) χ mass magnetic susceptibility (m3 kg−1 ) χm volume magnetic susceptibility χm0 reference volume magnetic susceptibility χ0 reference mass magnetic susceptibility (m3 kg−1 ) Subscripts c cold f fluid h hot m magnetic s solid 0 reference value such as in solar energy technology, crystal growth, and liquid metal cooling blankets used in fusion reactors. During the past three decades, numerous theoretical and experimental studies on thermomagnetic convection in enclosures have been investigated (Bian et al., 1996; Sarris et al., 2005; Pirmohammadi et al., 2009; Grosan et al., 2009; Sathiyamoorthy and Chamkha, 2010). The natural convection heat transfer of Cu–water nanofluid in a cold outer circular enclosure containing a hot inner sinusoidal circular cylinder in the presence of a horizontal magnetic field was numerically investigated by Sheikholeslami et al. (2012) using the control volume-based finite-element method. The results indicated that in the absence of a magnetic field, the heat transfer enhancement ratio deJournal of Porous Media
Numerical Simulation of Thermomagnetic Convection of Air 1061 alues and found that the magnetic force an be used to on natural convection of non-Newtonian power-law fluids perimentally investigated magnetothermal convection of L2014 ed the fo me dition using the lattice Boltzmann method.They found or oxygen gas can be achieved by the gradient magnetic field.Bednarz et al. (2004,2005. 009)stu ed the ncrease wit r all phas In recent years,with the development of a supercon- alyzed the effect of the inclined angle of the electric co tic induc he location of the electric coil,and the ka number nd h been an interesting research topic investigated by many a paramagnetic field placed in a micro-gravity condition researchers(Qiet al,2001:Tagawa et al,2002).The ef and under a uniform vertical gradient magnetic field in an of the cold side buoyancy force on convection depends on the relative ori- of paramagnetic and diamagnetic fluids.However,only 50 the th ontainer with thermal and field Wang et al.(2007)and Zeng et al.(2007,2009)nu gradients and found that the magnetic buoyancy force onvection in an enclosure out mence of the buoyanc orc gnetic edraG cal wall and cooled from an ing wall ,was loc (Da Ra and mag with a vertical magnetic force numbers,the results of their numerical investig d gradient. Huang et a magpeticforc has a signincant e tion of the diamagnetic field.Gray et al.(2001)obtained plication of a strong magnetic field to a porous me numerical similarity solutions for two-dimensional plums may be found in the fields of medical treatment. met en by the i .200 and allurgy,ma ssing,and com and t and Wakavama (2002)studied the of engineering o erations.Thus.the study of the efTect 0)0) of magnetic force on natural convection in porous me. u et al.(2003)derive dia is impo ant n and enginee the Bor magn ngap pres ral convection of paramagnetic.diamagnetic.and electri ooled from the sidewalls with an electric coil inclined cally conducting fluids n a cubic nclosure with therma ound the Y and Z axes is numeric ally investigated.The and magnene ld gra now and heat magn (2005)studied the natural com ction of nation angle of an electric coil around the y axis paramagnetic and diamagnetic fluids in a cylinder under and Z axis ()Ra number.Da number.and magnetic a gradient magnetic field at different thermal boundary force parameter y. Volume 18,Number 11,2015
Numerical Simulation of Thermomagnetic Convection of Air 1061 creases as the Rayleigh (Ra) number increases, while an opposite trend is observed in the presence of a magnetic field. Kefayati (2014) studied the effect of a magnetic field on natural convection of non-Newtonian power-law fluids in a cavity with a linearly heated wall using the finitedifference lattice Boltzmann method (FDLBM). Mejri et al. (2014) analyzed the effects of a magnetic field on nanofluid flow in a cavity with a sinusoidal boundary condition using the lattice Boltzmann method. They found that the heat transfer and fluid flow decline with the Hartmann number and increase with Ra number for all phase deviations. In recent years, with the development of a superconducting magnet providing strong magnetic induction of 10 T or more, the natural convection of paramagnetic fluids such as oxygen gas or air under a magnetic field has been an interesting research topic investigated by many researchers (Qi et al., 2001; Tagawa et al., 2002). The effect of the magnetic buoyancy force on the convection of paramagnetic fluids was first reported by Braithwaite et al. (1991). They found that the effect of the magnetic buoyancy force on convection depends on the relative orientation of the magnetic force and the temperature gradient. Carruthers and Wolfe (1968) theoretically and experimentally studied the thermal convection of oxygen gas in a rectangular container with thermal and magnetic field gradients and found that the magnetic buoyancy force cancels out the influence of the gravity buoyancy force when the rectangular enclosure, heated from one vertical wall and cooled from an opposing wall, was located in a horizontal magnetic field with a vertical magnetic field gradient. Huang et al. (1998a,b,c); studied the stability of a paramagnetic fluid layer by the action of a nonuniform magnetic field and the thermomagnetic convection of the diamagnetic field. Gray et al. (2001) obtained numerical similarity solutions for two-dimensional plums driven by the interaction of a line heat source and a nonuniform magnetic field. Qi et al. (1999, 2001) and Wang and Wakayama (2002) studied the magnetic control of thermal convection. Kaneda et al. (2002), Shigemitsu et al. (2003), and Akamatsu et al. (2003) derived a model equation for magnetic convection using a method similar to the Boussinesq approximation and studied the natural convection of paramagnetic, diamagnetic, and electrically conducting fluids in a cubic enclosure with thermal and magnetic field gradients at different thermal boundary values. Akamatsu et al. (2004), Filar et al. (2005), and Fornalik et al. (2005) studied the natural convection of paramagnetic and diamagnetic fluids in a cylinder under a gradient magnetic field at different thermal boundary values and found that the magnetic force can be used to control the heat transfer rate of paramagnetic and diamagnetic fluids. Yang et al. (2003, 2004) numerically and experimentally investigated magnetothermal convection of air or oxygen gas in an enclosure by using the gradient magnetic field available from neodymium–iron–boron permanent magnet systems, and pointed out that enhancement or suppression of magnetothermal convection of air or oxygen gas can be achieved by the gradient magnetic field. Bednarz et al. (2004, 2005, 2008, 2009) studied the natural convection of paramagnetic fluids in a cubic enclosure under a magnetic field by an electric coil and analyzed the effect of the inclined angle of the electric coil, the location of the electric coil, and the Ra number on the heat transfer rate of the paramagnetic fluid. Saha (2013) studied thermomagnetic convection and heat transfer of a paramagnetic field placed in a micro-gravity condition and under a uniform vertical gradient magnetic field in an open square cavity with three cold sidewalls. The aforementioned studies were concerned with the effect of the magnetic force on the natural convection of paramagnetic and diamagnetic fluids. However, only a few studies have paid attention to the combined effects of both magnetic and gravitational forces on the natural convection of a paramagnetic fluid in a porous medium. Wang et al. (2007) and Zeng et al. (2007, 2009) numerically investigated natural convection in an enclosure filled with a paramagnetic or diamagnetic fluid–saturated porous medium under a strong magnetic field. Considering the effect of the Darcy (Da), Ra, and magnetic force numbers, the results of their numerical investigations showed that the magnetic force has a significant effect on the flow field and heat transfer in a paramagnetic or diamagnetic fluid–saturated porous medium. The application of a strong magnetic field to a porous medium may be found in the fields of medical treatment, metallurgy, materials processing, and combustion, and there may be plenty of applications in the near future in the field of engineering operations. Thus, the study of the effect of magnetic force on natural convection in porous media is important for both scientific research and engineering applications. In the present study, the thermomagnetic convection of air in a porous cubic enclosure heated and cooled from the sidewalls with an electric coil inclined around the Y and Z axes is numerically investigated. The flow and heat transfer behaviors of thermomagnetic convection are studied in relation to the effect of the inclination angle of an electric coil around the Y axis (yEuler) and Z axis (zEuler), Ra number, Da number, and magnetic force parameter γ. Volume 18, Number 11, 2015
1062 Jiang ct al. 2.PHYSICAL MODEL where is the netie fore is the volu nsideration is sho ity (Hm-1),b is the magnetic flux density (T).pr is the sure,which is kept in a horizontal position.ndon-u fluid density (kgm 1 wall while the other fou U.VU-wU-V where U is the velocity vector (.w):p is the (Pa)is the fuid kinematic viscosity (kgm)gis acceleration (ms and s the perm around thexs(Coil and then around thexis (Coil 2).The gravitational force acts in the minus direc. the following equation: the 0-m+心2+pas (3) 3.MATHEMATICAL FORMULATION 3.1 Governing Equations eratur The hypoh s in thi from Eq.(2)gives pation are assumed to be negligible.According to Braith- (199),the 0rU.VU=-V+uv2U-些U+pX-PmXo】 x Vb2+(Pr-P)g (4) (1) where p=po+p',and p'is the pressure difference due to the perturbe d state (Pa).Beca e pr and x are a functior -Coil 2 (5) g m-m+(2)).-)+… HOT (6 COLD =w+(器).-+ For the paramagnetic fluid,the mass magnetic suscepti- deocgeriondloabotieeampowre bility isinv Coil 1 (7) where m is a constant value:Te is the fluid temperature (K):To =(Ti +T)/2 (K);and subscripts 0,h,and c FIG.1:Schematic diagram of the physical model and co nate system Journal of Porous Medi
1062 Jiang et al. 2. PHYSICAL MODEL A schematic of the system under consideration is shown in Fig. 1. The system consists of a porous cubic enclosure, which is kept in a horizontal position, and a one-turn electric circle coil, which generates a magnetic field. The porous cubic enclosure filled with air is heated isothermally from the left-hand-side vertical wall and cooled isothermally from the opposing wall, while the other four walls are thermally insulated. The electric coil is set centrally and coaxially with the enclosure. Inclination of the coil in general directions is given by the combination of two angles: yEuler and zEuler. The coil is inclined first around the Y axis (Coil 1) and then around the Z axis (Coil 2). The gravitational force acts in the minus Z direction. In the present study, the size of the cubic enclosure (L) is 0.03 m and the radius (r) of the coil is 0.05 m. 3. MATHEMATICAL FORMULATION 3.1 Governing Equations The hypotheses in this model are as follows: the fluid is considered to be steady, incompressible, and Newtonian, and both the viscous heat dissipation and magnetic dissipation are assumed to be negligible. According to Braithwaite et al. (1991), the magnetizing force can be given as follows: fm = χm 2µm ∇b 2 = ρfχ 2µm ∇b 2 (1) FIG. 1: Schematic diagram of the physical model and coordinate system where fm is the magnetic force; χm is the volumetric magnetic susceptibility; µm is the magnetic permeability (Hm−1 ); b is the magnetic flux density (T); ρf is the fluid density (kgm−3 ); and χ is the mass magnetic susceptibility (m3 kg−1 ). The Navier–Stokes equation which includes the magnetic force can be presented as ρfU·∇U=−∇p+µf∇2U− µf κ U+ ρfχ 2µm ∇b 2+ρfg (2) where U is the velocity vector (u, v, w); p is the pressure (Pa); µf is the fluid kinematic viscosity (kg m−1 s −1 ); g is the gravitational acceleration (ms−2 ); and κ is the permeability (m2 ). At the reference state of the isothermal state, there will be no convection. Therefore, Eq. (2) becomes the following equation: 0 = −∇p0 + ρf0χ0 2µm ∇b 2 + ρf0g (3) where p0 is the pressure at the reference temperature (Pa); ρ0 is the fluid density at the reference temperature (kgm−3 ); and χ0 is the mass magnetic susceptibility at the reference temperature (m3 kg−1 ). Subtracting Eq. (3) from Eq. (2) gives ρfU · ∇U = −∇p ′+µf∇2U − µf κ U+ (ρfχ−ρf0χ0) 2µm × ∇b 2 + (ρf − ρf0) g (4) where p = p0+p ′ , and p ′ is the pressure difference due to the perturbed state (Pa). Because ρf and χ are a function of temperature, ρfχ and ρf can be indicated according to the Taylor expansion method, respectively, as follows: ρfχ = (ρfχ) 0 + ( ∂ (ρfχ) ∂Tf ) 0 (Tf − T0) + · · · (5) ρf = ρf0 + ( ∂ρf ∂T ) 0 (Tf − T0) + · · · (6) For the paramagnetic fluid, the mass magnetic susceptibility is inversely proportional to absolute temperature, according to Curie’s law: χ = m Tf (7) where m is a constant value; Tf is the fluid temperature (K); T0 = (Th + Tc)/2 (K); and subscripts 0, h, and c represent the reference value, hot, and cold, respectively. Thus, Eq. (5) can be written as follows: Journal of Porous Media
Numerical Simulation of Thermomagnetic Convection of Air 1063 px-(pxb=(0x-p六),-)+… .Solid phase energy equation =(pBx--T)+… (8) 0=-k(祭+紧+3)+h-间 =--pmxoB(1+3)g-)+… thermal diffusivity (ms- ve is the fluid dynamic vis- cosity(m nlid the U.VU=- 1 mal conductvit WmK the solid-to-fluio x(-To)7b2-B(G-)g heat transfer coefficient (W m-2 (9) K )Tr is the fluid as (K)andT is the soli phase tempera- ts of the pa Equations(10)(15)can be non-dimensionalized as rameters and the superscripts of pressure are omitted,and follows: the governing equations can be written as follows: ·Continuityequation ·Continuity equation ++0 ++%=0 (16 (10) ·Momentum equations ·Momentum equations 完+器+密器原器器+v贺+w器-贺+(祭+别 -冬-0+)-a0w+)--vag (17 2um 密+瑞+-瑞(+贺++w-+(+ +尉)-兰-(+-四+)-贵-9 21 2um Oy (18) 密++密-密+(像±器v+v”w-+震+ )之-吧+器-合单 2um +gB(T:-To) (13) .Fluid phase energy equation .Fluid phase energy equation h(亚+0+)=( v股+v+w-+器+碧 +票+)+a-列 +H(日.-0) (20y (14 .Solid phase energy equation Volume 18,Number 11,2015
Numerical Simulation of Thermomagnetic Convection of Air 1063 ρfχ − (ρfχ) 0 = ( ∂ρf ∂Tf χ − ρf χ Tf ) 0 (Tf − T0) + · · · = ( −ρfβχ−ρf χ Tf ) 0 (Tf −T0)+ · · · (8) = −ρf0χ0β ( 1+ 1 T0β ) (Tf −T0)+ · · · The small higher-order amount in Eq. (8) is omitted and generated into Eq. (4), and Eq. (4) then becomes U · ∇U = − ∇p ′ ρf0 + µf ρf0 ∇2U − ν κ U − χ0β 2µm ( 1+ 1 T0β ) × (Tf − T0) ∇b 2 − β (Tf − T0) g (9) where β is the thermal expansion coefficient (K−1 ). For the sake of convenience, the subscripts of the parameters and the superscripts of pressure are omitted, and the governing equations can be written as follows: • Continuity equation ∂u ∂x + ∂v ∂y + ∂w ∂z = 0 (10) • Momentum equations u ∂u ∂x +v ∂u ∂y +w ∂u ∂z =− 1 ρf ∂p ∂x +νf ( ∂ 2u ∂x 2 + ∂ 2u ∂y 2 + ∂ 2u ∂z 2 ) − νf κ u − ( 1 + 1 T0β ) χ0β (Tf − T0) 2µm ∂ ( b 2 ) ∂x (11) u ∂v ∂x + v ∂v ∂y + w ∂v ∂z = − 1 ρf ∂p ∂y + νf ( ∂ 2v ∂x 2 + ∂ 2v ∂y 2 + ∂ 2v ∂z 2 ) − νf κ v− ( 1+ 1 T0β ) χ0β (Tf −T0) 2µm ∂ ( b 2 ) ∂y (12) u ∂w ∂x + v ∂w ∂y + w ∂w ∂z = − 1 ρf ∂p ∂z + νf ( ∂ 2w ∂x 2 + ∂ 2w ∂y 2 + ∂ 2w ∂z 2 ) − νf κ w − ( 1 + 1 T0β ) χ0β (Tf − T0) 2µm ∂ ( b 2 ) ∂z + gβ (Tf − T0) (13) • Fluid phase energy equation (ρcp) f ( u ∂Tf ∂x + v ∂Tf ∂y + w ∂Tf ∂z ) = εkf ( ∂ 2Tf ∂x 2 + ∂ 2Tf ∂y 2 + ∂ 2Tf ∂z 2 ) + h (Ts − Tf) (14) • Solid phase energy equation 0 = (1−ε) ks ( ∂ 2Ts ∂x 2 + ∂ 2Ts ∂y 2 + ∂ 2Ts ∂z 2 ) +h (Tf −Ts) (15) where x, y, and z are the Cartesian coordinates; u, v, and w are the velocity components (ms−1 ); αf is the fluid thermal diffusivity (ms−1 ); νf is the fluid dynamic viscosity (m2 s −1 ); (cp) f is the fluid specific heat at constant pressure (J kg−1 K −1 ); ε is the porosity; kf is the fluid thermal conductivity (W m−1 K −1 ); ks is the solid thermal conductivity (W m−1 K −1 ); h is the solid-to-fluid heat transfer coefficient (W m−2 K −1 ); Tf is the fluid phase temperature (K); and Ts is the solid phase temperature (K). Equations (10)–(15) can be non-dimensionalized as follows: • Continuity equation ∂U ∂X + ∂V ∂Y + ∂W ∂Z = 0 (16) • Momentum equations U ∂U ∂X + V ∂U ∂Y + W ∂U ∂Z = − ∂P ∂X + Pr( ∂ 2U ∂X2 + ∂ 2U ∂Y 2 + ∂ 2U ∂Z2 ) − Pr Da U − γRaPrθf C 2 ∂ ( B2 ) ∂X (17) U ∂V ∂X + V ∂V ∂Y + W ∂W ∂Z = − ∂P ∂Y + Pr( ∂ 2V ∂X2 + ∂ 2V ∂Y 2 + ∂ 2V ∂Z2 ) − Pr Da V − γRaPrθf C 2 ∂ ( B2 ) ∂Y (18) U ∂W ∂X +V ∂W ∂Y +W ∂W ∂Z =− ∂P ∂Z +Pr( ∂ 2W ∂X2 + ∂ 2W ∂Y 2 + ∂ 2W ∂Z2 ) − Pr Da W −γRaPrθf C 2 ∂ ( B2 ) ∂Z +RaPrθf (19) • Fluid phase energy equation U ∂θf ∂X + V ∂θf ∂Y + W ∂θf ∂Z = ∂ 2θf ∂X2 + ∂ 2θf ∂Y 2 + ∂ 2θf ∂Z2 + H (θs − θf) (20) • Solid phase energy equation Volume 18, Number 11, 2015
1064 Jiang et al. ()(Number Calculation The dimensionless variables and parameters in the above Inorder to con e the total heat equations are defined as wall is defined as follow e V=UL 4 0.5 -I P=PL 1a路c=8 1+ avr 6=五+x 2 0=,C=1+B 3.4 Numerical Procedure and Code Verification v= (22)The go erning equations [Egs.(16(21)1 were dis cretized by the finite-volume method on a uniform grid eaeaR system.Ihe thir and W or convectve kinematics d the components e dimen the di vection and diffusion terms.respectively.The set of dis sionless solid phase temperature:is the hot wall tem- cretized equations for each variable was solved by a perature (K).T is the cold wall temperature(K)Pr is the line-by-line th ure,combining the agonal matrix au numt erence method.The coupn solved by the semi-implicit method for pressure-linked ber,y is the dimensionless magnetic strength parameter equations algorithm.The convergence criterion was that B is the dimen agnetic fux d V:D :(30×30×30.40×40×40.and50 by the Biot-Savart law. B=-∫R (23) nsignificant differen es for the40×4040dst these parameters.Therefore,for all of the computations in this study,.a40×40×40 uniform grid was emp 3.2 Boundary Conditions order to validate the numerical methods and codes The non-slip condition is imposed for the three velocity et al (004)was selected as the benchmark solution for comparison.In Fig.2,our results fo (T)en .9999 09 and yRa =9.06 r=-0.5:(3)vertical adiabatic walls (-0.5 the results in Bednar)It can be seen tha 0.5):00:/aY =00./aY =0;and (4)top and bottom the present results are in very obtained by Bednarz present numerical cod Journal of Porous Media
1064 Jiang et al. ∂ 2θs ∂X2 + ∂ 2θs ∂Y 2 + ∂ 2θs ∂Z2 + ΛH (θf − θs) = 0 (21) The dimensionless variables and parameters in the above equations are defined as X = x L , Y = y L , Z = z L , U = uL αf , V = vL αf , W = wL αf , B = b b0 , θf = Tf −T0 Th−Tc , θs = Ts−T0 Th−Tc , P = pL2 ρfα 2 f , Pr = νf αf , Da = κ L2 , H = hL2 εkf , Λ = εkf (1 − ε) ks , Ra = gβ (Th − Tc)L 3 αfνf , T0 = Th + Tc 2 , b0 = µmi L , C = 1 + 1 T0β , γ = χ0b 2 0 µmgL (22) where X, Y , and Z are the dimensionless Cartesian coordinates; U, V , and W are the dimensionless velocity components; P is the dimensionless pressure; θf is the dimensionless fluid phase temperature; θs is the dimensionless solid phase temperature; Th is the hot wall temperature (K); Tc is the cold wall temperature (K); Pr is the Prandtl number; b0 is the reference magnetic flux density (T); i is the electric current in a coil (A); L is the side of a porous cubic enclosure (m); Ra is the Rayleigh number; γ is the dimensionless magnetic strength parameter; B is the dimensionless magnetic flux density; Da is the Darcy number; H is the dimensionless solid-to-fluid heat transfer coefficient; and Λ is the dimensionless thermal conductivity. The magnetic field is calculated as follows by the Biot–Savart law: B = − 1 4π ∫ R × dS R3 (23) where R is the dimensionless position vector, and S is the dimensionless tangential element of a coil. 3.2 Boundary Conditions The non-slip condition is imposed for the three velocity components on the solid walls. The temperature boundary conditions are as follows: (1) left vertical wall (X = −0.5): θf = θs = 0.5; (2) right vertical wall (X = 0.5): θf = θs = −0.5; (3) vertical adiabatic walls (Y = −0.5, 0.5): ∂θf/∂Y = ∂θs/∂Y = 0; and (4) top and bottom horizontal adiabatic walls (Z = −0.5, 0.5): ∂θf/∂Z = ∂θs/∂Z = 0. 3.3 Nusselt Number Calculation In order to compare the total heat transfer rate, the Nusselt number is used. The average Nusselt number at the hot wall is defined as follows: Num = − Λ 1 + Λ ∫ 0.5 −0.5 ∫ 0.5 −0.5 ∂θf ∂X X=−0.5 dY dZ − 1 1 + Λ ∫ 0.5 −0.5 ∫ 0.5 −0.5 ∂θs ∂X X=−0.5 dY dZ (24) 3.4 Numerical Procedure and Code Verification The governing equations [Eqs. (16)–(21)] were discretized by the finite-volume method on a uniform grid system. The third-order quadratic upstream interpolation for convective kinematics scheme and the second-order central-difference scheme were implemented for the convection and diffusion terms, respectively. The set of discretized equations for each variable was solved by a line-by-line procedure, combining the tridiagonal matrix algorithm with the successive over-relaxation iteration method. The coupling between velocity and pressure was solved by the semi-implicit method for pressure-linked equations algorithm. The convergence criterion was that the maximal residual of all of the governing equations was less than 10−6 . Three grid sizes (30 × 30 × 30, 40 × 40 × 40, and 50 × 50 × 50) were selected for analysis. The average Nusselt numbers for all three grid sizes were monitored at ε = 0.5, Pr = 0.71, H = 10, Λ = 10, Ra = 1 × 105 , Da = 1 × 10−3 , xEuler = 0 ◦ , and γ = 10. The results showed insignificant differences for the 40 × 40 × 40 grids to these parameters. Therefore, for all of the computations in this study, a 40 × 40 × 40 uniform grid was employed. In order to validate the numerical methods and codes of the present study, a recent, similar study by Bednarz et al. (2004) was selected as the benchmark solution for comparison. In Fig. 2, our results for Da = 107 , ε = 0.9999, H = 0, g = 0, and γRa = 9.06 × 105 with different electric inclined angles xEuler are compared with the results in Bednarz et al. (2004). It can be seen that the present results are in very good agreement with those obtained by Bednarz et al. (2004), which validates the present numerical code. Journal of Porous Media
Numerical Simulation of Thermomagnetic Convection of Air 1065 5 30 FIG.2:Comparisons of the average Nusselt number at 0 4.RESULTS AND DISCUSSION and Da ler As indicated by the aforementioned mathematic model. the natural convection under consideration is governed by seven non-dime onal para meters: e.Y.Ra,Pr, Da,F constant at e =0.5.Pr =0.71,H=10.A=10.Ra-1 x10,and Da1x 1 respectively.Therefore.in 0 this st and heat transfer characteristics in the porous enclosure. 4.1 Numerical Results under the Non-Gravitational Condition scribe non-gravity cases.Figures 3-6show the fluid ph aemaSue,widpaseotemast5m2 and Da eft-hand-side vertic al wall and cooled isothermally from ¥0 30° 60° 90° When the coil is in a horizontal position (u FIG.4:Solid phase isothermal surface for different ele ()the magnetic buoyaney force is upward tric inclined anglesandat near the upper part of the hot wall and the lower part of and Da=1x 10-3 under the non-gravitational condition Volume 18,Number 11,2015
Numerical Simulation of Thermomagnetic Convection of Air 1065 FIG. 2: Comparisons of the average Nusselt number at various xEuler between the present study and those of published results (Bednarz et al., 2004) 4. RESULTS AND DISCUSSION As indicated by the aforementioned mathematic model, the natural convection under consideration is governed by seven non-dimensional parameters: ε, γ, Ra, Pr, Da, H, and Λ. In the present study, the porosity (ε), Prandtl number, H, Λ, Rayleigh number, and Darcy number are kept constant at ε = 0.5, Pr = 0.71, H = 10, Λ = 10, Ra = 1 × 105 , and Da = 1 × 10−3 , respectively. Therefore, in this study attention is primarily given to the effects of the magnetic force and the inclination angle on the fluid flow and heat transfer characteristics in the porous enclosure. 4.1 Numerical Results under the Non-Gravitational Condition Here, Ra becomes zero and γ becomes infinity when g = 0; therefore, finite product γRa has to be used to describe non-gravity cases. Figures 3–6 show the fluid phase isothermal surface, solid phase isothermal surface, magnetic force vector, and velocity streamline for γRa = 2.5 × 106 and Da = 1 × 10−3 at various coil inclination angles yEuler and zEuler, respectively. In each graph, the porous cubic enclosure is heated isothermally from the left-hand-side vertical wall and cooled isothermally from the opposing wall, while the other four walls are thermally insulated. When the coil is in a horizontal position (yEuler, zEuler) = (0◦ , 0◦ ), the magnetic buoyancy force is upward near the upper part of the hot wall and the lower part of FIG. 3: Fluid phase isothermal surface for different electric inclined angles yEuler and zEuler at γRa = 2.5 × 106 and Da = 1 × 10−3 under the non-gravitational condition FIG. 4: Solid phase isothermal surface for different electric inclined angles yEuler and zEuler at γRa = 2.5 × 106 and Da = 1 × 10−3 under the non-gravitational condition Volume 18, Number 11, 2015
1066 Jiang ctal it is dow 军重安 ssymmetrical in terms of the middle plane of the porous 0 深国团其 atoward the adiabatic alls at the top and potm aue to 全国国国 which results in four cellular structures with horizontal er a non-gravit 0 深国国厨 =()are the same as those of thenc eonofteeetnceolfiomhoralositonead 0).For round theaxis due to no change in the position of the Maetie ore vecem lectric coil. r the non-gravitat oceearhecOdaldawhecoartowardThro First the cold air moves along the cold wall to the fou 宽魔含 ce upper 30° rection of the ycoordinate.drives the air horizontal flow ear the hot wall.The magnetic force acts in the posi the ion of ad 0° 重置 l and vertical c ounter flows are formed along the top respectively,near the cold with respect to the colplane the same character of flow an be observed along the back-side and bottom adia- 0 30° 60 pects the svmmetrical twisted distribution The magnet FIG.6:Velocity streamline for different electric inclined vector,temperature isosurface,and velocity streamline at Journal o Porous Media
1066 Jiang et al. FIG. 5: Magnetic force vector for different electric inclined angles yEuler and zEuler at γRa = 2.5 × 106 and Da = 1 × 10−3 under the non-gravitational condition FIG. 6: Velocity streamline for different electric inclined angles yEuler and zEuler at γRa = 2.5 × 106 and Da = 1 × 10−3 under the non-gravitational condition the cold wall in the porous cubic enclosure, it is downward near the lower part of the hot wall and the upper part of the cold wall in the porous cubic enclosure, and it is symmetrical in terms of the middle plane of the porous cubic enclosure. Driven by the magnetic buoyancy force, the hot air from the middle plane moves along the left hot wall toward the adiabatic walls at the top and bottom, and then it flows to the right cold wall. When the air is near the right cold wall, due to the attracting magnetic force, the cold air goes downward and upward to the middle plane and returns to the left hot wall along the middle plane, which results in four cellular structures with horizontal symmetry about the middle plane of the enclosure. Because the porous cubic enclosure is under a non-gravity environment, the flow structure and isotherms at (yEuler, zEuler) = (90◦ , 90◦ ) are the same as those of the 90◦ inclination of the electric coil from horizontal position (yEuler, zEuler) = (0◦ , 0◦ ). For yEuler = 0◦ , the flow structure and isotherms do not change with the electric coil inclination around the Z axis due to no change in the position of the electric coil. For (yEuler, zEuler) = (90◦ , 0◦ ), the magnetic buoyancy force near the cold wall draws the cold air toward the coil, while near the hot wall it repels the hot air from the coil. First, the cold air moves along the cold wall to the four corners and then moves toward the coil. The continuity law requires that the hot air being squeezed out to the center of the porous cubic enclosure from the region near the left-hand-side hot wall will then return to the cold wall. For (yEuler, zEuler) = (30◦ , 90◦ ), the magnetic buoyancy force distribution is twisted over the hot and cold walls. In the upper part of porous cubic enclosure, the magnetic buoyancy force, which acts in the reverse direction of the y coordinate, drives the air horizontal flow near the hot wall. The magnetic force acts in the positive Z direction near the junction of the hot and front adiabatic walls; therefore, the air goes upward and moves toward the clod along the top horizontal edge. Horizontal and vertical counter flows are formed along the top and front-side adiabatic walls, respectively, near the cold wall due to the attraction of a stronger magnetic force. Due to the symmetrical distribution of the magnetic force with respect to the coil plane, the same character of flow can be observed along the back-side and bottom adiabatic walls. Hence, four separate roll cells appear along the four adiabatic side walls, and the temperature field respects the symmetrical twisted distribution. The magnetic vector, temperature isosurface, and velocity streamline at (yEuler, zEuler) = (30◦ , 90◦ ) respect the mirror image relationship with those at (yEuler, zEuler) = (60◦ , 90◦ ). When Journal of Porous Media
Numerical Simulation of Thermomagnetic Convection of Air 1067 x 107.From Table I and Fig 7.we also can find that the average Nusselt number does not depend onat did therefore,the local h im mo the mag magnetic field. mlar)=(90°,90)is the same Table 1 lists the average Nusselt numbers for the com a and y n 2g8e0ya之25×102 eocn umbe are in (a =(90 30)when yRa effect of the coil inclination angle on the that th Nuselt number mov t) 4.2 Numerical Results under the Gravitational Condition TABLE 1:Average Nusselt numbers at Da=1x 10-3 for combinations of the two inclination angles and 4.2.1 Effect of the Magnetic Force Number (y) Euler under the non-gravitational conditior Inclination angle oeY ZEuler 0 30° 60°90° The left-hand-side wall is hot and the right-hand-side wall s cold in each graph.Without a magnetic field().the 1.2641.1671.097 1.112 30° 1.264 1.258 1.183 1.097 and then the right-hand side cold wall 60°. 1.264 1386 1362 1.167 000 12641438 1438 1264 as seen in fluid isothermal contours.where it respects a he an Ra=21 lockwise vortex. the porous er enhaaepoe 2.062 1798 1.487 nd t 2.062 202 2 sults show that the flow and temperature field are verv 30 2.062 2.328 2.236 1.798 complex.The counteraction of the magnetic buoyancy 90° 2.0622.4772.4772.062 vRa=50×10 0° 3358 287g 2325 1903 orces enhances the air flow in the lower part of the porous 30o 3358315g2560 600 26门 370 the air 0。 the grav. 3.358 3.935 3.935 3.364 The convection in the porous enclosure is strength- 10 ned with the increase of magnetic force number y,thus 5.288 4.397 3.311 2.279 an 30 5.288 4.724 3.833 3.311 60° 5.288 5454 5.053 4.397 becomes large and the hot fluid along the lef-hand-side 900 5.2885.928 5928 5288 hot wall is repelled downward at the upper-halfheight and Volume 18,Number 11,2015
Numerical Simulation of Thermomagnetic Convection of Air 1067 the coil inclination angle is another angle, the flow and temperature field in porous cubic enclosure are very complex. For the temperature field, with the increase of magnetic force, the temperature difference between the solid and liquid is greatly increased; therefore, the local thermal non-equilibrium model is important in simulating the heat transfer in a porous medium when there is a strong magnetic field. Table 1 lists the average Nusselt numbers for the combination of the two inclination angles at Da = 1 × 10−3 and γRa = 1.0 × 106 , 2.5 × 106 , 5.0 × 106 , and 1.0 × 107 under the non-gravitational condition. Figure 7 shows the effect of the coil inclination angle (yEuler, zEuler) on the average Nusselt number for different magnetic force numbers (γRa) when Da = 1 × 10−3 . Table 1 shows that the heat transfer rate increases with increasing γRa. For example, at (yEuler, zEuler) = (π/6, π/6), the average NusTABLE 1: Average Nusselt numbers at Da = 1 × 10−3 for combinations of the two inclination angles yEuler and zEuler under the non-gravitational condition Inclination angle zEuler yEuler 0 ◦ 30◦ 60◦ 90◦ γRa = 1.0 × 106 0 ◦ 1.264 1.167 1.097 1.112 30◦ 1.264 1.258 1.183 1.097 60◦ 1.264 1.386 1.362 1.167 90◦ 1.264 1.438 1.438 1.264 γRa = 2.5 × 106 0 ◦ 2.062 1.798 1.561 1.487 30◦ 2.062 2.021 1.722 1.561 60◦ 2.062 2.328 2.236 1.798 90◦ 2.062 2.477 2.477 2.062 γRa = 5.0 × 106 0 ◦ 3.358 2.879 2.325 1.903 30◦ 3.358 3.158 2.569 2.325 60◦ 3.358 3.642 3.442 2.879 90◦ 3.358 3.935 3.935 3.364 γRa = 1.0 × 107 0 ◦ 5.288 4.397 3.311 2.279 30◦ 5.288 4.724 3.833 3.311 60◦ 5.288 5.454 5.053 4.397 90◦ 5.288 5.928 5.928 5.288 selt number is 1.258 at γRa = 1.0 × 106 , and it increases to 3.158 at γRa = 5.0 × 106 and to 4.724 at γRa = 1.0 × 107 . From Table 1 and Fig. 7, we also can find that the average Nusselt number does not depend on zEuler at yEuler = 0 ◦ and 180◦ because, in this case, the electric coil inclination around the Z axis affects neither the magnitude nor the direction of the magnetic field. The average Nusselt number at (yEuler, zEuler) = (90◦ , 90◦ ) is the same as that at (yEuler, zEuler) = (0◦ , 0◦ ). The local maximum average Nusselt number appears in (yEuler, zEuler) = (30◦ , 90◦ ) and (60◦ , 90◦ ). The local minimum average Nusselt number appears in (yEuler, zEuler) = (90◦ , 30◦ ) when γRa < 2.5 × 106 ; however, at γRa ≥ 2.5 × 106 the local minimum average Nusselt number moves to (yEuler, zEuler) = (90◦ , 0◦ ). 4.2 Numerical Results under the Gravitational Condition 4.2.1 Effect of the Magnetic Force Number (γ) Figure 8 shows the effect of the γ number on the resultant force vector, fluid phase temperature, solid phase temperature, and velocity vector contours in the vertical cross section at Y = 0.5, xEuler = 0◦ , Ra = 105 , and Da = 10−3 . The left-hand-side wall is hot and the right-hand-side wall is cold in each graph. Without a magnetic field (γ = 0), the flow of air in the porous enclosure is pure gravity convection: the hot air rises up along the left-hand-side hot wall and then descends along the right-hand-side cold wall, as seen in fluid isothermal contours, where it respects a clockwise vortex. When there is a magnetic field, the air in the porous enclosure is then acted upon by both magnetic buoyancy and the gravity buoyancy forces. The results show that the flow and temperature field are very complex. The counteraction of the magnetic buoyancy and gravity buoyancy forces repels the clockwise flow of air in the upper part of the porous enclosure; however, the synergy of the magnetic buoyancy and gravity buoyancy forces enhances the air flow in the lower part of the porous enclosure. When γ < 10, the magnetic buoyancy force is weak; the air flow then completely succumbs to the gravity buoyancy force and respects a clockwise flow structure. The convection in the porous enclosure is strengthened with the increase of magnetic force number γ, thus there is competition between the magnetic buoyancy and gravity buoyancy forces and the air flow appears to be a stratified flow phenomenon. At γ = 50, the magnetic force becomes large and the hot fluid along the left-hand-side hot wall is repelled downward at the upper-half height and Volume 18, Number 11, 2015
1068 Jiang et al Nu. 。 ZRe-1.0x10 Ra2.5×10 02402934039039 7Ra-5.0x10 7R-1.0x10 when ntal osition (O 08)the magnetic buoyancy and gravity buovaney forcesareact ing almost in the same direction in the upper part of the hand-side wall.If y 100,then the air flow state in rous cubic enclosure,and larger magnitudes of th e and the flo ever the ma of two cellular structures with approximately horizontal are acting almost in the opposite direction in the lower symmetry about the middle of the enclosure. orce and 4.2.Effect of Inclined Anglese and he left hot wall toward the top and bottom adiabatic walls Figures9-12 show the fluid phase isothermal surface. from the lower part of electric coil plane,and then moves solid phase isothermal surface.resultant force vector of he magnetic and gravity forces,and velocity streamline e.the Journal of Porous Media
1068 Jiang et al. FIG. 7: Variation of the average Nusselt number distribution for (yEuler, zEuler) at Da = 1 × 10−3 under the nongravitational condition proceeds toward the cold wall where fluid is attracted upward along the cold plate, which forms an anticlockwise flow in the upper-half region. In the lower-half enclosure, usual gravitational convection occurs ascending along the hot left-hand-side wall and descending along the righthand-side wall. If γ > 100, then the air flow state in the porous enclosure is mainly determined by the magnetic buoyancy force, and the flow in the enclosure is that of two cellular structures with approximately horizontal symmetry about the middle of the enclosure. 4.2.2 Effect of Inclined Angles yEuler and zEuler Figures 9–12 show the fluid phase isothermal surface, solid phase isothermal surface, resultant force vector of the magnetic and gravity forces, and velocity streamline for various coil inclination angles yEuler and zEuler when Ra = 1 × 105 , γ = 25, and Da = 1 × 10−3 . When the coil is in a horizontal position (yEuler, zEuler) = (0◦ , 0◦ ), the magnetic buoyancy and gravity buoyancy forces are acting almost in the same direction in the upper part of the porous cubic enclosure, and larger magnitudes of the resultant force and accelerated flow can be expected. However, the magnetic buoyancy and gravity buoyancy forces are acting almost in the opposite direction in the lower part of the porous cubic enclosure, and smaller magnitudes of the resultant force and accelerated flow can be generated. The resultant force vector drives the air along the left hot wall toward the top and bottom adiabatic walls from the lower part of electric coil plane, and then moves to the right cold wall. When the air is near the right cold wall, due to the attracting magnetic force, the cold air Journal of Porous Media