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