《高等传热学》课程教学资源（文献资料）Numerical simulation of thermomagnetic convection of air in a porous cubic enclosure with an electric coil inclined in general orithentations using a LTNE model

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 sci￾entific applications, such as in heat exchangers, solar re￾ceivers, geothermal exploitation, cooling of electronic de￾vices, etc. Various researchers have carried out experi￾mental and numerical investigations to study fluid flow and heat transfer in porous enclosures (Khalili and Shiv￾akumara, 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) investi￾gated 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 Ra￾makrishna 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 nat￾ural convection within porous square cavities. Kim and Choi (2013) theoretically analyzed the onset of convec￾tion motion in an initially quiescent, horizontal isotropic porous layer. Alloui and Vasseur (2013) investigated nat￾ural 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. Dur￾ing the past three decades, numerous theoretical and ex￾perimental studies on thermomagnetic convection in en￾closures 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 de￾Journal 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 finite￾difference 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 con￾dition using the lattice Boltzmann method. They found that the heat transfer and fluid flow decline with the Hart￾mann number and increase with Ra number for all phase deviations. In recent years, with the development of a supercon￾ducting magnet providing strong magnetic induction of 10 T or more, the natural convection of paramagnetic flu￾ids 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 ef￾fect 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 ori￾entation of the magnetic force and the temperature gradi￾ent. Carruthers and Wolfe (1968) theoretically and experi￾mentally 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 verti￾cal 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 sta￾bility of a paramagnetic fluid layer by the action of a non￾uniform magnetic field and the thermomagnetic convec￾tion 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 non￾uniform 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 simi￾lar to the Boussinesq approximation and studied the natu￾ral convection of paramagnetic, diamagnetic, and electri￾cally conducting fluids in a cubic enclosure with thermal and magnetic field gradients at different thermal bound￾ary 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 diamag￾netic fluids. Yang et al. (2003, 2004) numerically and ex￾perimentally investigated magnetothermal convection of air or oxygen gas in an enclosure by using the gradi￾ent magnetic field available from neodymium–iron–boron permanent magnet systems, and pointed out that enhance￾ment 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 en￾closure under a magnetic field by an electric coil and an￾alyzed 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) nu￾merically investigated natural convection in an enclosure filled with a paramagnetic or diamagnetic fluid–saturated porous medium under a strong magnetic field. Consid￾ering the effect of the Darcy (Da), Ra, and magnetic force numbers, the results of their numerical investiga￾tions showed that the magnetic force has a significant ef￾fect on the flow field and heat transfer in a paramagnetic or diamagnetic fluid–saturated porous medium. The ap￾plication of a strong magnetic field to a porous medium may be found in the fields of medical treatment, met￾allurgy, 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 me￾dia is important for both scientific research and engineer￾ing 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 con￾vection are studied in relation to the effect of the incli￾nation 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 enclo￾sure, 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 isother￾mally 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 cen￾trally 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 direc￾tion. 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 dissi￾pation are assumed to be negligible. According to Braith￾waite 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 co￾ordinate system where fm is the magnetic force; χm is the volumetric magnetic susceptibility; µm is the magnetic permeabil￾ity (Hm−1 ); b is the magnetic flux density (T); ρf is the fluid density (kgm−3 ); and χ is the mass magnetic sus￾ceptibility (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 perme￾ability (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 suscepti￾bility 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 pa￾rameters 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 vis￾cosity (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 ther￾mal 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 tempera￾ture (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 co￾ordinates; U, V , and W are the dimensionless velocity components; P is the dimensionless pressure; θf is the dimensionless fluid phase temperature; θs is the dimen￾sionless solid phase temperature; Th is the hot wall tem￾perature (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 num￾ber; γ 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 bound￾ary 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 dis￾cretized 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 con￾vection and diffusion terms, respectively. The set of dis￾cretized 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 Nus￾selt 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 emp￾loyed. 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 dif￾ferent 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 num￾ber, 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 de￾scribe non-gravity cases. Figures 3–6 show the fluid phase isothermal surface, solid phase isothermal surface, mag￾netic 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 elec￾tric 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 elec￾tric 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 in￾clined 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 down￾ward 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. Be￾cause 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◦ incli￾nation 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 cen￾ter 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 buoy￾ancy 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 di￾rection of the y coordinate, drives the air horizontal flow near the hot wall. The magnetic force acts in the posi￾tive Z direction near the junction of the hot and front adi￾abatic walls; therefore, the air goes upward and moves toward the clod along the top horizontal edge. Horizon￾tal 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 adia￾batic walls. Hence, four separate roll cells appear along the four adiabatic side walls, and the temperature field re￾spects the symmetrical twisted distribution. The magnetic vector, temperature isosurface, and velocity streamline at (yEuler, zEuler) = (30◦ , 90◦ ) respect the mirror image rela￾tionship 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 com￾plex. For the temperature field, with the increase of mag￾netic force, the temperature difference between the solid and liquid is greatly increased; therefore, the local ther￾mal 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 com￾bination 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 num￾bers (γRa) when Da = 1 × 10−3 . Table 1 shows that the heat transfer rate increases with increasing γRa. For ex￾ample, at (yEuler, zEuler) = (π/6, π/6), the average Nus￾TABLE 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 mag￾nitude 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 min￾imum 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 temper￾ature, 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 convec￾tion: 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 mag￾netic buoyancy and the gravity buoyancy forces. The re￾sults 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 grav￾ity buoyancy force and respects a clockwise flow struc￾ture. The convection in the porous enclosure is strength￾ened 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 non￾gravitational condition proceeds toward the cold wall where fluid is attracted up￾ward 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 right￾hand-side wall. If γ > 100, then the air flow state in the porous enclosure is mainly determined by the mag￾netic 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 act￾ing almost in the same direction in the upper part of the porous cubic enclosure, and larger magnitudes of the re￾sultant force and accelerated flow can be expected. How￾ever, 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 magni￾tudes 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