Numerical investigation of natural convection in a rectangular cavity under different directions of uniform magnetic field

International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Numerical investigation of natural convection in a rectangular cavity under different directions of uniform magnetic field P.X. Yu, J.X. Qiu, Q. Qin, Zhen F. Tian ⇑ Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, PR China

a r t i c l e

i n f o

Article history: Received 20 July 2012 Received in revised form 3 December 2012 Accepted 22 August 2013

Keywords: Magnetohydrodynamic (MHD) natural convection Streamfunction-velocity formulation Inclined uniform magnetic field Finite difference Compact scheme

a b s t r a c t Natural convection flows of an electrically conducting fluid under a uniform magnetic field at different angles h with respect to horizontal plane are investigated numerically in rectangular cavities. A new compact finite difference algorithm, involving a second-order compact scheme for the streamfunction-velocity (w  u) form of Navier–Stokes equations and a fourth-order compact scheme for the energy equation, is employed to solve the steady-state laminar magnetohydrodynamic (MHD) natural convection problems. Numerical simulations are carried out in a wide range of Rayleigh number (Ra) and Hartmann number (Ha) at the Prandtl number Pr = 0.025. The computed results show that the heat transfer is not only determined by the strength of the magnetic field, but also influenced by the inclination angle. Especially, when the aspect ratio (A) is less or more than 1, it is found that the inclination angle plays a great role on flow and heat transfer. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction In the recent decades, laminar natural convection of electrically conducting fluids in the presence of magnetic field has been one of the major interesting research subjects due to its widely engineering applications [1,2]. For example, it can control growing crystals because a magnetic field can help to suppress convective flows. The phenomena has also been successfully used in solidification processes to weaken the buoyancy driven flows. Since the velocity and the temperature governing equations are coupled due to the buoyancy force, and the Lorentz force driven by the magnetic field is connected with the velocity, the study of the MHD natural convection is complicated. The MHD natural convection has been investigated by many researchers using experimental, analytical and numerical methods. Okada and Ozoe [3] carried out the experiments using molten gallium with Prandtl number (Pr) of 0.024 in a cubical enclosure heated from one side wall and cooled from the opposite wall with all other walls are insulated in the presence of different magnetic field direction, and found that the external magnetic field in the vertical direction was more effective than the magnetic field applied parallel to the heated vertical wall. An analytical solution was proposed by Garandet et al. [4] to be used to model the effect of a transverse magnetic field on natural convection in a two ⇑ Corresponding author. Permanent address: No. 220 Handan Road, Shanghai 200433. Tel.: +86 21 55664597. E-mail addresses: [email protected], [email protected] (Z.F. Tian). 0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.08.087

dimensional (2D) cavity. Vasseur et al. [5] obtained analytical and numerical solutions for an inclined tall cavity under transverse magnetic field for a porous medium. Besides these, many numerical results have been reported by means of various numerical methods. Ozoe and Maruo [6] investigated numerically the natural convection of an electrically conducting fluid (Pr = 0.054) in the presence of magnetic field based on 2D streamfunction–vorticity formulation of Navier–Stokes equations and obtained the Nusselt number correlation in terms of Rayleigh (Ra), Prandtl and Hartmann (Ha) numbers. After that, Ozoe and Okada [7] employed a three dimensional vector potential-vorticity formulation finite difference method to study the effect of magnetic direction. Rudraiah et al. [8] developed a finite difference scheme consisting of modified ADI (Alternating Direction Implicit) method and SLOR (Successive Line Over Relaxation) method to solve the streamfunction–vorticity formulation of the problem. In their study, numerical predictions were obtained at Pr = 0.733, and they found the average Nusselt number decreases with an increase of Hartmann number. Ben Hadid et al. [9,10] used an ADI method combining a second-order central differentiation and a fourth-order compact Hermitian method based on the streamfunction–vorticity formulation to validate their analytical solutions for the central region and for the turning flow region under some assumptions in 2D cases and also developed a projection method to solve the 3D cases. Al-Najem et al. [11] used a control volume algorithm consisting of ADI procedure to numerically investigate the laminar natural convection in tilted enclosure with transverse magnetic field. Pirmohammadi [12] used SIMPLER algorithm for laminar natural

1132

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

convection flows in the presence of a longitudinal magnetic field at Pr = 0.733 and extended the method to a tilted square enclosure [13]. Recently, Lo [14] proposed a high-resolution method by the differential quadrature method based on the 2D velocity–vorticity form of Navier–Stokes equations to simulate the effect of a transverse magnetic field on buoyancy-driven magnetohydrodynamic flow in an enclosure and concluded that the heat transfer rate is at its maximum for higher Pr and in the absence of MHD effects. In addition, some scholars extended the convectional problem in the different boundary conditions. For example, Venkatachalappa and Subbaraya [15] have presented numerical results for a fluid of Pr = 0.733 in a rectangular enclosure under the influence of a vertical magnetic field with uniform heat flux from the side wall. Kandaswamy [16] numerically simulated the MHD convection in a square cavity with partially thermally active vertical walls by the control volume method with QUICK scheme. Sathiyamoorthy [17] numerically studied the cavity with the uniformly heated bottom wall and the linearly left vertical wall whereas the right vertical wall is linearly heated or cooled while the top wall kept thermally insulated. However, there are just few literatures to discuss effects of the magnetic direction in detail besides the work by Okada and Ozoe [3]. In the previous studies, the numerical methods mentioned were based on the streamfunction–vorticity form or velocity–vorticity form of the 2D incompressible Navier–Stokes equations. Although these methods have been proved effectively, one has to develop numerical boundary conditions for the vorticity because these two formulations lack the simple physical boundary conditions for the vorticity at the no-slip boundaries. In order to avoid this drawback, the compact streamfunction–velocity or the streamfunction formulation of the 2D incompressible fluid flows has been utilized by some authors for solving Navier–Stokes equations [18– 22]. The advantage of the compact streamfunction–velocity or the streamfunction formulation mainly embodies in that the boundary conditions of streamfunction and velocities are generally known, so that the numerical difference schemes proposed are very effective and efficient [18–22]. Combined with a fourth-order accurate compact scheme for the temperature equation, the compact scheme based on this formulation has been utilized in the natural convection problem and its stable and effective numerical performance has also been showed [23]. However, to the best of authors knowledge, no report has been found so far using the streamfunction–velocity formulation-based compact difference method to solve the MHD convection problems. The objective of this study is to develop a streamfunction– velocity formulation-based second-order compact difference method for solving the MHD convection problem and to examine systematically the characteristics of flow and heat transfer under an inclined uniform magnetic field with an electrically conducting fluid (Pr = 0.025) in rectangular cavities using the compact scheme based on the streamfunction–velocity formulation. The remainder of this paper is organized as follows. Section 2 gives a nondimensional streamfunction–velocity formulation of the governing equations for the MHD natural convection problem. In Section 3, the numerical method based on the streamfunction–velocity formulation is proposed and verified with other numerical results. The effects of three parameters involving Hartmann number, inclination angle and aspect ratio on the flow and heat transfer performance are discussed in Section 4. At last, the whole work is summarized in Section 5.

2. Model description and governing equations As shown in Fig. 1, the problem being solved is the flow of an incompressible electrical conducting fluid of Pr = 0.025 in a 2D

Fig. 1. Schematic diagram of the cavity along with the boundary conditions

rectangular cavity, whose height and width are given by H and L respectively. In the present study, the two horizontal walls are assumed to be insulated, while the vertical walls are maintained at constant temperature. The left wall with temperature Th is hotter than the right wall with temperature Tc. The magnetic Reynolds number (Rm) which is defined as the ratio of the inducted magnetic b0 Þ and total external magnetic field ~ B is assumed to be small field ð~ (Rm  1), hereby, the induced magnetic field generated by electromagnetic induction phenomenon can be neglected. Thus, Ohm’s law is utilized for calculation of the electric current density ~ J (i.e. ~ Eþ~ V ~ BÞ, where r is the electric conductivity. Moreover, J = rð~ base on the hypothesis of weak electric field ~ E, the two dimensional governing equations of the problem can be described as

@u @ v þ ¼ 0; @x @y

ð1Þ

! @u @u 1 @p @2u @2u r 2 u þv ¼ þm þ  B2 ðusin h  v cos h sin hÞ; @x @y @x2 @y2 q @x q ð2Þ u

@v @v 1 @p @2v @2v þv ¼ þm þ @x @y q @y @x2 @y2 

! þ gbðT  T 0 Þ

r 2 B ðv cos2 h  u cos h sin hÞ; q

! @T @T @2u @2u u þv ¼j þ ; @x @y @x2 @y2

ð3Þ

ð4Þ

where u and v are the velocity components in x- and y-directions, p, T and B are the pressure, temperature and the uniform external magnetic field respectively, h stands for the inclined angle which is formed by a uniform magnetic field and the positive direction of the horizontal axis, g is the gravity acceleration. The kinematic viscosity m, the thermal diffusivity j and the coefficient of the volumetric thermal expansion of fluid b are determined by the physical properties of the fluid at the reference density q and the reference temperature T0 (Tc is taken as the reference temperature in this paper). The velocity and temperature boundary conditions are expressed as

u ¼ v ¼ 0;

T ¼ Th

at x ¼ 0; ð0 6 y 6 HÞ

u ¼ v ¼ 0;

T ¼ Tc

at x ¼ L; ð0 6 y 6 HÞ

u ¼ v ¼ 0;

T y ¼ 0 at y ¼ 0; ð0 6 x 6 LÞ

u ¼ v ¼ 0;

T y ¼ 0 at y ¼ H: ð0 6 x 6 LÞ

ð5Þ

1133

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

In order to make the above model more clearly describe the intrinsic characteristics of the practical problems, the dimensionless process would be employed. The following dimensionless parameters are defined

ðx; yÞ ðu; v ÞL ; ðu ; v  Þ ¼ ; L j 2 T  Tc pL ; p ¼ : T ¼ Th  Tc qj2

ðx ; y Þ ¼

ð6Þ

Here x⁄, y⁄, u⁄, v⁄, T⁄, p⁄ are dimensionless variables. For the sake of convenience, the asterisks will be omitted in the next mathematical formulation. Thus, Eqs. (1)–(4) are transformed into the following dimensionless forms by using the above dimensionless variables

@u @ v þ ¼ 0; @x @y

ð7Þ !

u

@u @u @p @2u @2u 2 þv ¼  þ Pr þ  Ha2 u sin h þ Ha2 v cos h sin h ; @x @y @x @x2 @y2

Because the Eq. (13) is only determined by the streamfunction

w and its first derivatives @w ð¼ v Þ and @w ð¼ uÞ, it is termed as pure @x @y streamfunction or streamfunction–velocity (w  u) formulation of Navier–Stokes equation. Considering the physical meanings, the dimensionless boundary conditions applied to the computational regions are given by

w ¼ u ¼ v ¼ 0;

T ¼ 1 at x ¼ 0; ð0 6 y 6 AÞ

w ¼ u ¼ v ¼ 0;

T ¼ 0 at x ¼ 1; ð0 6 y 6 AÞ

w ¼ u ¼ v ¼ 0;

T y ¼ 0 at y ¼ 0; ð0 6 x 6 1Þ

w ¼ u ¼ v ¼ 0;

T y ¼ 0 at y ¼ 1: ð0 6 x 6 1Þ

Here, A ¼ stands for the aspect ratio. The heat transfer characteristics across the cavity are often depicted by the Nusselt number. In this problem, we mainly pay attention to the Nusselt number on the hot wall, which is computed by

 @T  NuðyÞ ¼   @x x¼0

Nu0 ¼

Z

ð9Þ

NuðyÞdy:

ð19Þ

Another important output concerned for the flow characteristics is the total kinetic energy, which can reflect the movement intensity and defined by

E¼

Z

A 0

@T @T @ 2 T @ 2 T þv ¼ þ ; @x @y @x2 @y2

A

0

! @2v @2v 2 2 2 ; þ Pr þ  Ha v cos h þ Ha u cos h sin h @x2 @y2

u

ð18Þ

and the average Nusselt number is defined as

ð8Þ @v @v @p u þv ¼  þ RaPrT @y @x @y

ð17Þ

H L

Z

1

ðu2 þ v 2 Þdx dy:

ð20Þ

0

ð10Þ 3. Numerical method and validation

where

Ra ¼

gbðT h  T c ÞL3

mj m Pr ¼ ; jsffiffiffiffiffiffiffiffiffiffiffiffiffi rB2 L2 Ha ¼ mq

For convenience, some standard finite difference operators at mesh point (x, y) are given by

;

ð11Þ

are three dimensionless parameters. The above equations can be simplified through eliminating the pressure p with Eqs. (8) and (9) and introducing the streamfunction according to Eq. (7)

u¼

@w ; @y

v ¼

@w : @x

ð12Þ

Then, Eqs. (7)–(10) are equivalent to

@4w @4w @4w þ2 2 2þ 4 4 @x @x @y @y " ! !# 1 @2v @2v @2u @2u v ¼ u þ þ Pr @x2 @y2 @x2 @y2     @T @v @u @u @ v 2 þ Ra sin h cos h ; þ Ha2  cos2 h þ sin h þ  @x @y @x @y @x ð13Þ u¼

@w ; @y

v ¼ u

@w ; @x

@T @T @ 2 T @ 2 T þv ¼ þ : @x @y @x2 @y2

ð14Þ

ð15Þ

ð16Þ

d2x dy / ¼

/5 þ /6  /7  /8  2ð/2  /4 Þ

dx d2y / ¼

/5  /6  /7 þ /8  2ð/1  /3 Þ

3

2h

3

2h

d2x / ¼

/1  2/0 þ /3

d2y / ¼

/2  2/0 þ /4

h

h

2

2

;

ð21Þ

;

ð22Þ

;

dx / ¼

/1  /3 ; 2h

ð23Þ

;

dy / ¼

/2  /4 ; 2h

ð24Þ

where h denotes the space increment in x- and y-directions, and 0, 1, 2, 3, 4, 5, 6, 7 and 8 represent grid points (xi, yj), (xi+1, yj), (xi, yj+1), (xi1, yj), (xi, yj1), (xi+1, yj+1), (xi1, yj+1), (xi1, yj1), and (xi+1, yj1), respectively (see Fig. 2). In [23], Yu and Tian have developed a second order accurate compact scheme to solve the natural convection problem based on the w  u formulation. It is noted that comparing with the conventional natural convection problem in the absence of magnetic field [23], Eq. (13) is equivalent to add a magnetic field Fb in the right hand term

    @v @u @u @ v 2 sin h cos h F b ¼ Ha2  cos2 h þ sin h þ  @y @x @y @x " #   2 2 @ w @ w @u @ v 2 2 ; sin h cos h cos h þ sin h þ  ¼ Ha2 @x2 @y2 @x @y

ð25Þ

which is a combination of the streamfuction and velocities and can be discretized as

1134

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

2 •

6 •

A1 ¼ 32  2hð4u0 þ 3u1 þ u2  u3 þ u4 Þ

• 5

2

þ h ð4u20 þ u0 u1  u0 u3 þ u2 v 0  u4 v 0 Þ; A2 ¼ 32  2hð4v 0 þ v 1 þ 3v 2 þ v 3  v 4 Þ 2

þ h ð4v 20 þ u0 v 1 þ v 0 v 2  u0 v 3  v 0 v 4 Þ;

0 •

3 •

• 1

A3 ¼ 32 þ 2hð4u0  u1 þ u2 þ 3u3 þ u4 Þ 2

þ h ð4u20  u0 u1 þ u0 u3  u2 v 0 þ u4 v 0 Þ; A4 ¼ 32 þ 2hð4v 0 þ v 1  v 2 þ v 3 þ 3v 4 Þ

4 •

7 •

2

þ h ð4v 20  u0 v 1  v 0 v 2 þ u0 v 3 þ v 0 v 4 Þ;

• 8

2

A5 ¼ 8 þ hð4u0  u2 þ u4  4v 0  v 1 þ v 3 Þ þ 2h u0 v 0 ;

Fig. 2. Computational stencil.

h i 2 2 F b ¼ Ha2 cos2 hd2x w þ sin hd2y w þ sin h cos hðdx u  dy v Þ þ Oðh Þ: ð26Þ

2

A6 ¼ 8 þ hð4u0 þ u2  u4  4v 0 þ v 1  v 3 Þ  2h u0 v 0 ; 2

A7 ¼ 8 þ hð4u0  u2 þ u4 þ 4v 0  v 1 þ v 3 Þ þ 2h u0 v 0 ;

Hence, following [23], Eq. (13) can be approximated as 2

4  X 4 48w0  12 wj ¼ 6hðv 1  v 3  u2 þ u4 Þ þ h dx d2y v  d2x dy u j¼1

þ

4 4 X X v 0 uj  u0 v j

1 2 h Pr

j¼1

A8 ¼ 8 þ hð4u0 þ u2  u4 þ 4v 0 þ v 1  v 3 Þ  2h u0 v 0 :

! 4

þ h ðRaT x0 þ F b Þ:

j¼1

ð27Þ Substituting (26) into Eq. (27) and omitting the truncation error O(h6), yields 4  X 4 48w0  12 wj ¼ 6hðv 1  v 3  u2 þ u4 Þ þ h dx d2y v  d2x dy u j¼1 4 4 X X v 0 uj  u0 v j

1 2 þ h Pr 4

2

h

j¼1 2

þ h Ha cos

ð20  2hðv 2  v 0 ÞÞT 0 ¼ ð8  2hðv 2  v 0 ÞÞT 2 þ ð4  hðu2  u0 ÞÞT 1 þ ð4 þ hðu2  u0 ÞÞT 3 þ 2ðT 5 þ T 6 Þ;

!

2

þ sin i þ sin h cos hðdx u  dy v Þ :

þ h RaT x0

ð20  2hðv 0  v 4 ÞÞT 0 ¼ ð8  2hðv 0  v 4 ÞÞT 4 þ ð4 þ hðu0  u4 ÞÞT 1 þ ð4  hðu0  u4 ÞÞT 3 þ 2ðT 7 þ T 8 Þ:

hd2y w ð28Þ

Combining the approximation (28) rearranged with the fourth-order compact scheme suggested for the temperature and velocities in [23], a stable and efficient streamfunction–velocity formulation-based second order accurate compact method is obtained for solving the MHD convection problem, which is 2

2

ð48 þ 2Ha2 h Þw0  ð12 þ Ha2 h cos2 hÞðw1 þ w3 Þ 2 2

1 4 1 T3  T1 ðT x Þ1 þ ðT x Þ0 þ ðT x Þ3 ¼ : 6 6 6 2h

4

ð29Þ 1 4 1 w w v1 þ v0 þ v3 ¼ 3 1 ; 6 6 6 2h

ð30Þ

1 4 1 w  w4 u2 þ u0 þ u4 ¼ 2 ; 6 6 6 2h

ð31Þ

8 X Aj T j ¼ 0;

ð32Þ

j¼0

where

A0 ¼ 160  8hðu1 þ u3  v 2 þ v 4 Þ  8h ðu20 þ v 20 Þ;

ð35Þ

On the left and the right boundaries, the value of (Tx)0 can be calculated by

1 ð17T 1;j þ 9T 2;j þ 9T 3;j  T 4;j Þ; 6h

 ð12 þ Ha h sin hÞðw2 þ w4 Þ ¼ 6hðv 1  v 3  u2 þ u4 Þ þ Ha2 h sin h cos hðdx u  dy v Þ ! 4 4  1 X X 4 2 4 2 2 þ h dx dy v  dx dy u þ h v 0 uj  u0 v j þ h RaðT x Þ0 ; Pr j¼1 j¼1

ð34Þ

In the inner domain, the value of (Tx)0 in the Eq. (29) can be computed by

ðT x Þ1;j þ 3ðT x Þ2;j ¼

2

2

ð33Þ

4

j¼1

hd2x w

It is easily found that the above second-order compact approximation (29) for the streamfunction–velocity formulation at grid point (xi, yj) only involves the five values of w in compact stencil and the associated algebraic systems is linear and diagonally dominant, which means the scheme (29) can be performed efficiently. For the top and bottom walls, the following fourth-order discretizations for the Neumann boundary conditions are used [23]:

ðT x ÞN;j þ 3ðT x ÞN1;j ¼

1 ð17T N;j  9T N1;j  9T N2;j þ T N3;j Þ 6h

ð36Þ

ð37Þ

respectively. It is worthy pointing out that the algebraic systems arising from the newly proposed finite difference approximations (29) is still symmetric and positive, which allows algorithms like Jacobi, Gauss–Seidel, SOR and conjugate gradient (CG) to be used. In this paper, the multigrid method is employed to accelerate calculation. The entire strategy is summarized in the following algorithm: 1. Begin with wn, un, vn and Tn. 2. Obtain wn+1 in the inner region [h, 1  h]  [h, 1  h] by solving the discretized streamfunction Eq. (29). 3. Solve the velocity un+1 and vn+1 in the same region from Eqs. (30) and (31) based on known wn+1. 4. Calculate the temperature Tn+1 in the region [h, 1  h]  [h, 1  h] using the Eq. (32) and the temperatures on the bottom and the top walls using Eqs. (33) and (34).

1135

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

Table 1 Comparison of the second-order finite difference solutions by Rudraiah et al.[8], the control volume solutions by Kandaswamy et al. [16], and the differential quadrature solutions by Lo [14] for different Grashof numbers (Gr = Ra/Pr) and Hartmann numbers at Pr = 0.733. Gr

Ha

Nu0 Ref. [8]

Ref. [16]

Ref. [14]

Present

41  41

51  51

31  31

33  33

65  65

129  129

2  104

0 10 50 100

2.5188 2.2234 1.0856 1.0110

2.6237 2.3234 1.0987 1.0245

2.5303 2.2381 1.0756 1.0048

2.5226 2.2319 1.0758 1.0061

2.5305 2.2381 1.0762 1.0062

2.5308 2.2383 1.0761 1.0062

2  105

0 10 50 100

4.9198 4.8053 2.8442 1.4317

5.1876 4.9825 2.9784 1.6318

5.0814 4.9752 2.9966 1.4644

4.9648 4.8588 2.9371 1.4510

5.0641 4.9556 2.9867 1.4617

5.0795 4.9712 2.9942 1.4629

1

(a)

1

(b)

0.8

0.6

0.6

y

y

0.8

0.4

0.4

0.2

0.2

0

0

0.2

0.4

x

0.6

0.8

1

(c)

0

1

0.2

0.4

0

0.2

0.4

x

0.6

0.8

1

0.6

0.8

1

0.6

0.8

1

1

(d)

0.8

0.8

0.6

0.6

y

y

0

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

0

1

x 1

(f)

1

0.8

0.8

0.6

0.6

y

y

(e)

x

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

x

0.8

1

0

0

0.2

0.4

x

Fig. 3. The streamfunction contours for MHD natural convection under horizontal magnetic field: (a) Ra = 104, Ha = 0, (b) Ra = 105,Ha = 0, (c) Ra = 104,Ha = 10, (d) Ra = 105, Ha = 30, (e) Ra = 104,Ha = 100 and (f) Ra = 105,Ha = 320.

1136

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

5. Determine the gradient of temperature T nþ1 in the region x [h, 1  h]  [h, 1  h] by solving the Eq. (35) and the values on the left and the right walls by solving Eqs. (36) and (37). 6. Repeat above steps from 2 to 5 by setting n = 0, 1, 2, . . . , until kun+1  unk1 < 1010, where u represents streamfunction w and temperature T. In order to examine correctness of the current codes, grid dependence studies are performed by a series of grid systems for simulations of the heat transfer characteristics in a square enclosure. The comparisons of the average Nusselt numbers along with the hot wall are shown in Table 1. It is clear that there is good agreement between the present computations and those of Rudraiah [8], Kandaswamy [16] and Lo [14]. Notice that our results are closer to those of Lo [14], who proposed a high-resolutions scheme. The maximum deviations of the calculated average Nusselt

(a)

numbers between current results and Lo [14] are found to be within 0.4%, meanwhile, the differences on 65  65 and 129  129 grids are no more than 0.3% as well. It should be mentioned that although the solutions of Lo [14] can use fewer grid points, it would take much more time to calculate the first and second derivatives and also should deal with the numerical boundary condi1 tions. In summary, our algorithm is trusted and a grid with h ¼ 64 would be used in coming investigation.

4. Results and discussion 4.1. Effects of Hartmann number The effect of Hartmann number has been extensively investigated by many researchers and it was found that the flow and heat

(b)

1

1

0.8

0.6

0.6

y

y

0.8

0.4

0.4

0.2

0.2

0

0

0.2

0.4

x

0.6

0.8

1

(c)

0

1

0.2

0.4

0

0.2

0.4

0.2

0.4

x

0.6

0.8

1

0.6

0.8

1

0.6

0.8

1

1

(d)

0.8

0.8

0.6

0.6

y

y

0

0.4

0.4

0.2

0.2

0

0.2

0.4

x

0.6

0.8

0

1

(f)

1

x

1

0.8

0.8

0.6

0.6

y

y

(e)

0

0.4

0.4

0.2

0.2

0

0

0.2

0.4

x

0.6

0.8

1

0

0

x

Fig. 4. The isothermal contours contours for MHD natural convection under horizontal magnetic field: (a) Ra = 104, Ha = 0, (b) Ra = 105,Ha = 0, (c) Ra = 104,Ha = 10, (d) Ra = 105, Ha = 30, (e) Ra = 104,Ha = 100 and (f) Ra = 105,Ha = 320.

1137

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

transfer would be suppressed with an increase of Hartmann numbers no matter what Pandtl numbers or Rayleigh numbers were. Here, we present our numerical results for Ra = 104 and 105 respectively in the presence of horizontal magnetic field (h = 0) in a square cavity. Figs. 3 and 4 exhibit the contours of streamfunction and isothermal for different Hartmann numbers. As mentioned by Ozoe and Maruo [6], the combination of the dimensionless parameter Ha2/Ra would be applicable at a critical state, we also find this parameter plays great role for the MHD natural convection and define this parameter as c. In Figs. 3 and 4, the left part stands for the results of different Ha at Ra = 104, while the right shows the solutions for the approximate c at Ra = 105. It can

Table 2 Summary of the results with different Ha at Ra = 104 and Ra = 105.

c

jwjmin

jujmax

jvjmax

E

Numax

Numin

Nu0

Ra = 10 0 0.000 10 0.010 50 0.250 100 1.000

4.759 3.446 0.478 0.125

15.656 11.762 2.489 0.961

15.534 11.883 1.905 0.483

110.234 60.550 1.672 0.156

0.6696 0.5812 0.7665 0.9179

3.0858 2.7394 1.3344 1.0926

2.0212 1.7862 1.0374 1.0032

Ra = 105 0 0.000 30 0.009 160 0.256 320 1.024

7.522 4.753 0.490 0.126

44.595 24.383 4.657 1.707

44.725 28.702 2.230 0.597

600.573 196.259 2.807 0.243

0.8618 0.5354 0.7198 0.9061

5.9335 4.7660 1.4677 1.1092

3.4704 2.7901 1.0506 1.0037

Ha

4

Ra=10 4 Ra=10 5 6 Ra=10 Fitting Curve

-4

log10(velocity/Ra)

(b)

-3

log10(u/Ra)

-5 log10(v/Ra)

-6

-7

-8

-3 Ra=10 4 Ra=10 5 6 Ra=10 Fitting Curve

-4

log10(velocity/Ra)

(a)

log10(v/Ra)

-5 log10(u/Ra)

-6

-7

1

1.5

2

2.5

3

3.5

-8

4

1

1.5

log10Ha

2

2.5

3

3.5

4

log10Ha

Fig. 5. The variation of the velocities with different magnetic fields at (a) h = 0° and (b) h = 90°.

(a) 1

(b)

-0.4 -0.8 -1.2 -1.6 -2 -2.4 -2.8

0.8

-0.4 -0.8 -1.2 -1.6

0.8

-2 -2.4

-2.8

-3.2

-3.2

0.6

y

y

0.6

1

0.4

0.4

0.2

0.2

0

0

0.2

0.4

(c) 1

x

-0.4 -0.8

0.8

0.6

0.8

0

1

(d) -1.2 -1.6 -2 -2.4

0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

1 -0.4 -0.8 -1.2 -1.6

0.8

-2 -2.4

-2.8

-2.8

-3.2

-3.2

y

0.6

y

0.6

x

0.4

0.4

0.2

0.2

0

0

0.2

0.4

x

0.6

0.8

1

0

0

0.2

0.4

x

Fig. 6. The streamfunction contours for Ra = 104 and Ha = 10 (c = 0.01): (a) h = 0°, (b) h = 45°, (c) h = 90° and (d) h = 135°.

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

be seen that the flow structures and temperature distributions for the approximate c at different Rayleigh numbers are very similar. For a weak magnetic field (c = 0.01), the flow structures are almost similar to those in a purely thermal case (Ha = 0). It is observed that there is a circulating flow in the enclosure with central streamline distorted slightly into an elliptic shape. When c increases, which means the magnetic field is strengthened, the stratification of the temperature field in the interior begins to diminish. Especially, when c is as large as one, it is noted that the isotherm stratification in the core diminishes and the isotherms are almost parallel to the vertical walls, indicating that most of the heat transfer is by heat conduction not by heat convection. In addition, two

(b)

.0 1 -0 25

0.8

0.6

.0 4

-0. 1

-0

0.8

1

.0 2

.1

-0 . 12

2 .0

-0.0 4 -0.06

-0

-0

1

-0.08

(a)

secondary vertices appear along the vertical centerline. In fact, comparing the items in the momentum equation, when Ra and Ha2 is large enough, the convection terms just play small effects. When c is large, it means the Lorentz force is dominant whereas the buoyancy is dominant. Quantitative comparison of the results for different Hartmann numbers is listed in Table 2. The quantities presented here are the maximum absolute value of streamfunction jwjmax, horizontal velocity jujmax and vertical velocity jvjmax, the total kinetic energy E, the average Nusselt number Nu0 on the hot wall, and the maximum and minimum values Numax and Numin of the local Nusselt number on the hot wall. For temperature, when Ha increases, it

-0

1138

0.6

-0

6 .0 8 14 1 6 8 .0 1 -0 -0. .12 -0. -0. -0.1 -0 2 -0 . 2 .2

y

y

-0

0.4

0.4

0.2

0

1 -0.

0

(c) 1

0.2

0.4

25

x

0.2

0.6

0.8

0

1

(d)

-0.02 -0.04

0

0.2

0.4

-0. -0.06 -0 . 0 8 1 -0. 0.12 1 -0 . 16 4 -0 .18 -0 .2

0.8

-0.08

0.6

0.6

-0.12

0.4

0.4

0.2

0.2

0.2

1

-0.0 -0.0 1 2

4

y

y

-0.125

0

0.8

-0.0

-0.1

0

0.6

1

-0.06

0.8

x

0.4

x

0.6

0.8

0

1

0

0.2

0.4

x

0.6

0.8

1

Fig. 7. The streamfunction contours for Ra = 104 and Ha = 100 (c = 1.0): (a) h = 0°, (b) h = 45°, (c) h = 90° and (d) h = 135°.

(a) 1.81

2.85 4

Ra=10 (LHS) Ra=10 5(RHS)

1.805 1.8

1.006 2.8

1.0055

2.75

1.0045

Nu

Nu

1.0065 4

Ra=10 Ra=10 5

1.005

1.795 1.79 1.785

1.004

2.7

1.0035

2.65

1.0025

2.6

1.0015

1.003

1.78 1.775 1.77

(b)

1.002 0

50

θ

100

150

0

50

θ

100

150

Fig. 8. Effect of inclination angle h on the average Nusselt number for: (a) c = 0.01 and (b) c = 1.0.

1139

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

is found the average Nusselt number approaches to one and the difference between the maximum and minimum of local Nusselt number becomes small. It also reveals that with the increase of Ha, the velocities and total kinetic energies decrease. Because the Lorentz force is perpendicular to the magnetic field, the vertical velocities are damped

(a)

1

(b)

θ=0 0 θ=45 0 0 θ=90 0 θ=135

0.8

more significantly when Ha is large enough. In [4], an analytical solution for the natural convection under the vertical uniform magnetic field with sufficient large Ha was proposed by Garandet et al. In their theory, the components of velocity u⁄ (=Pru) and v⁄ (=Prv) vary as Gr Ha2 and Gr Ha1.5. In our non-dimensional definition, the correlation can be rewritten as u / Ra Ha2 and

1 θ=0 0 θ=45 0 0 θ=90 0 θ=135

0.8 0.6 0.4 0.2 0

v

y

0.6

-0.2

0.4

-0.4 0.2

-0.6 -0.8

0 -1

-0.8 -0.6 -0.4 -0.2

0

0.2

u

0.4

0.6

0.8

1

-1

0

0.2

0.4

x

0.6

0.8

1

Fig. 9. Comparisons of (a) u-velocity along the vertical centerline and (b) v-velocity along the horizontal centerline for different inclination angles at Ra = 104 and Ha = 100 (c = 1.0).

-7 Ra=10 4 Ra=10 5 6 Ra=10 Fitting Curve

-8

2

log10(E/Ra )

-9

θ=0

-10

θ=45 o

o

-11 -12 -13 -14 -15

1

1.5

2

2.5

3

3.5

4

log10Ha Fig. 10. The variation of the total kinetic energies with different magnetic fields.

(a)

62.5

200

4

(b)

Ra=10 5(LHS) Ra=10 (RHS) 62

0.34

4

Ra=10 5 Ra=10

0.32 0.3

195

0.28 0.26

190

185

61

180

60.5

log10E

log10E

61.5

0.24 0.22 0.2 0.18 0.16 0.14

175

60 0

50

θ

100

150

0.12 0.1

0

50

100

θ

Fig. 11. Effect of inclination angle h on the total energy for: (a) c = 0.01 and (b) c = 1.0.

150

200

1140

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

v / Ra Ha1.5. To determine the different index a of Ra Haa, we plot log10(jujmax/Ra) and log10(jvjmax/Ra) with respect to log10Ha with h = 0° and h = 90° respectively in Fig. 5 and the slopes of the fitted curves can be regarded as the approximation of a. The functions of the fitting curves in Fig. 5 using least squares method are given by

log10

jujmax ¼ 1:49log10 Ha  1:05; Ra

ð38Þ

log10

jv jmax ¼ 1:92log10 Ha  0:45 Ra

ð39Þ

for h = 0°, and

log10

jujmax ¼ 1:90log10 Ha  0:49; Ra

ð40Þ

log10

jv jmax ¼ 1:49log10 Ha  1:03 Ra

ð41Þ

for h = 90°. According to Eqs. (40) and (41), it is clear that our result is close to the analytical one. However, it is noted, when the magnetic field is imposed in the horizontal direction, the relationship can be approximate to the contrary relationships as u / Ra Ha1.5 and v / Ra Ha2 following Eqs. (38) and (39). And thereby, we can conclude that the magnetic field can suppress the perpendicular

(a)

(b)

(c)

Fig. 12. Comparison of the streamfunction contours in the absence of magnetic field (left) and in the presence of magnetic field (right): (a) A = 0.5 (b) A = 0.75 and (c) A = 2.0.

1141

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

4.2. Effects of inclination angle It can be found from Eq. (25) that Fb(p + h) is equivalent to Fb(h), which indicates that the Lorentz force item Fb has a period of p. And then, in order to examine the effects of inclination angles, computations are carried out for the inclined angles varying from 0 to p, while other geometric parameters remain unchanged. As mentioned in Section 4.1, c = Ha2/Ra plays a great role in the magnetoconvection problem, and therefore, we implement numerical experiments with c = 0.01 and 1 at Ra = 104 and Ra = 105 in the square cavity. Figs. 6 and 7 contain the contours for the stream and temperature functions with different inclination angles for c = 0.01 and 1.0, respectively. It can be seen that the effect of inclination angles on the flow is considerable significant with an increase of magnetic fields. For c = 1.0, two vertical eddies occur when h = 0° and separate in the horizontal direction when h = 90°, while a primary vortex converges in the flow center when h = 45° and h = 135°. These phenomena can be explained by the combined effect of the magnetic force and the buoyancy force. When c is sufficient high, the flow is mainly driven by the Lorentz force, which is perpendicular to the direction of the magnetic field. For example, when h is zero, the Lorentz force only acts in the y-direction. Considering the symmetry, the u-velocity is dominant at the vertical centerline due to v = 0, which means the total speed is almost parallel to the magnetic field and therefore the magnetic force is close to zero. Because the magnetic force can reduce the effect of buoyancy, the total force at the vertical centerline is larger than that in the other areas. As the result, the conducting fluids are stretched more significant near the vertical centerline. Fig. 8 shows variation of the average Nusselt number with the different inclination angles. Apparently, for the low c, which the magnetic filed is relatively weak, compared with buoyancy force, the effect of the Lorentz force is not obvious and there exists only a single peak within a period. However, with the strengthening of the magnetic field, the influence of the Lorentz force also will be reinforced increasingly and the capabilities of the interaction of Lorentz force and buoyancy force become stronger so that the double peaks phenomenon appeared for Ra = 104 and Ra = 105 in the vicinity of h = 45° and h = 135°, respectively. And it also provides further related evidences that appropriate magnetic fields can weaken the ability of heat transfer of natural convection of conducting fluids. Fig. 9 exhibits the variations of the profiles of velocity component u along vertical centerline (x = 0.5) and velocity component v along horizontal centerline (y = 0.5) with different angles. From this figure, we can see that the effect of inclination angles is quite apparent for variations of velocity. For horizontal velocity, as the inclination angles increasing from 0° to 90°, the changing trend of the curve of horizontal velocity u becomes gentle relatively and then gradually reaches steep from 90° to 135°. Particularly, the velocity gradient ð@u j Þ along centerline in the cavity except @y x¼0:5 near the upper and lower wall tends to constant for h = 90°, which the Lorentz force is in the horizontal direction. We can see from Fig. 9(b) that vertical velocity v has the opposite trends, which this phenomenon coincides with the conclusion obtaining in 4.1. Meanwhile, there exists two thin Hartmann boundary layers near the walls, where the effect of viscosity should be considered [4].

According to the above analysis, the changing extent of velocity is not identical with various tilt angles for the different magnitudes of the external magnetic. The perpendicular and parallel components of velocity are damped with the magnetic field following different exponential rates. Here, in order to further investigate the flow characteristic as a whole, the total kinetic energy with various inclination angles is employed. The data of the numerical simulative total energy are plotted with the inclination angles at 0° and 45° in Fig. 10, and the correlations for the total energy for natural convection are presented by

log10

E Ra2

¼ 3:58log10 Ha  1:64

ð42Þ

for h = 0°, and

log10

E Ra2

¼ 3:90log10 Ha  0:81

ð43Þ

for h = 45°. It should also be mentioned that the fitting curve for h = 90° is similar with that of h = 0°. It is evident from (42) and (43) that log10(E/Ra2) and log10Ha are basically satisfied with linear relationship, and its variation also tends to more gentle as increasing of the inclination angles, and E is proportion to Ra2Hab, where b always maintains between 3 and 4, as expected, no matter how the inclination angles change. To elucidate further the effect of the magnetic fields and inclination angles, the changes of the total energy with the different inclination angles and magnetic fields are displayed in Fig. 11. And a similar changing tendency also can be observed from this figure, and more peaks can be found for Ra = 105 and c = 1.0, showing that the effect of inclination angles on flow characteristic is more notable for higher strength of the magnetic field. 4.3. Effects of aspect ratio In the previous researches, the investigations on effects of aspect ratio A are less than those of Hartman number. Fig. 12 compares the streamlines of steady flow motion in the absence of (left) and in the presence of (right) magnetic field with aspect ratio A from 0.5 to 2. Without the magnetic field, it can be seen that there is only one primary clockwise vortex in fluid field for A = 0.5 and A = 0.75, while several secondary vortices appear in

1.02

A=0.5 A=0.75 A=1 A=1.5 A=2

1.015

Nu

velocities u\ more significantly as the function u\ / Ra Ha2 than the parallel velocities uk which vary as Ra Ha1.5. In general, the parameter c is important to reflect the effects of Hartmann number in the MHD natural convection for different Rayleigh numbers. In the next parts, in order to examine systematically the characteristics of the flow and heat transfer by other factors, c rather than Ha will be treated as fixed to represent the strength of magnetic field.

1.01

1.005

1

0

30

60

90

θ

120

150

180

Fig. 13. Effect of inclination angle h and aspect ratio on the average Nusselt number at c = 1.

1142

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

(a)

(b) A=0.5 A=0.75 A=1 A=1.5 A=2

0.015

0.0004

A=0.5 A=0.75 A=1 A=1.5 A=2

0.00035

0.0003

0.00025 0.01

λ

λ

0.0002

0.00015 0.005

0.0001

5E-05

0

0

0

30

60

90

120

150

180

0

30

60

90

120

θ

θ

150

180

Fig. 14. Effect of inclination angle h and aspect ratio on the kinetic energy restrain ratio: (a) c = 1 and (b) c = 9.

Table 3 Summary of the results in the rectangular cavity with different aspect ratios in the absence of magnetic field at Ra = 104. A

jujmax

jvjmax

Eu

Ev

E

Eu/Ev

0.5 0.75 1 1.5 2

7.970 12.84 15.66 17.61 18.35

4.110 9.665 15.53 21.01 21.91

5.7995 26.827 55.629 88.445 91.890

1.7144 16.076 54.605 116.65 128.97

7.5139 42.902 110.23 200.09 220.86

3.383 1.669 1.019 0.715 0.713

the cavity for A = 2. However, inhibited by a magnetic field along the x-axis with c = 1, the convection structures of fluid field for different aspect ratios all become simple as shown in the right part of Fig. 12. It is observed that a pair of clockwise vortices exist in the middle of cavity and align perpendicularly to the direction of the magnetic field with the exception for A = 0.5. This flow structure can be also obtained for A = 0.5 by strengthen the magnetic field. Hence we conclude that the effects of a given magnetic field on the natural convection are similar over different aspect ratios. Next, the interaction effects of inclination angle and aspect ratio on heat transfer are considered. Fig. 13 has shown that the influence of inclination angle becomes much more significant on the average Nusselt number than that of A = 1. It is obvious that with A > 1,Nu0 reaches a peak around 70° < h < 110°, and the larger A is, the steeper the peak is. On the contrary, when A < 1,Nu0 reaches a valley around 60° < h < 120°. For aspect ratio closer to 1, the fluctuation in Nu0 is smaller. In general, Nu0 depends on the different aspect ratio A and the inclination angle h. To investigate the interaction effects of inclination angle and aspect ratio further, the attention is paid to the total kinetic energy.

(a)

Fig. 14 displays the restrain ratio k due to different aspect ratio A and inclination angle h for c = 1 and c = 9.

k¼

0.5

Eu ¼

Z

1

u2 dx dy;

Ev ¼

0

Z 0

A

Z

1

v 2 dx dy:

3 2

A=0.5 A=1 A=2 Fitting Curve

1

log10E

-1

-1.5

-2

-2

-3

-2.5

-4

1.6

1.8

2

2.2

log10Ha

2.4

2.6

ð45Þ

0

0

-1

-3

A

It indicates that Eu/Ev > 1 at A < 1 and Eu/Ev < 1 at A > 1. When A change from 0.5 to 2, the ratio of Eu to Ev decreases. In 4.1, it is noted that flow perpendicular to the magnetic field will be suppressed more than that along the magnetic field. For h = 0°, Ev is suppressed more than Eu, and therefore the total inhibition effect will be enhanced with aspect ratio A increases. On the contrary, the total inhibition effect will be reduced for h = 90° when A changes from 0.5 to 2. This trend has always been established for different Hartmann numbers, even for high Ha. For example, considering a specific Ra = 104, Fig. 15 exhibits the correlations between Ha and E and the detail fitting curves are given by

-0.5

log10E

Z 0

A=0.5 A=1 A=2 Fitting Curve

0

ð44Þ

where EB represents the total kinetic energy under the inhibition of magnetic field B while E0 stands for the total kinetic energy in the absence of magnetic field. Comparing to Fig. 13, we can see the pattern of restrain ratio is very similar to that of the average Nusselt number in Fig. 13, a peak around 70° < h < 110° for A > 1, and a valley around 60° < h < 120° for A < 1. It can also be observed that the curves are more smooth with smaller c, which represents a weaker magnetic field or a smaller Ha. This phenomenon is connected with variation of magnetic field’s inhibition with different inclination angles on convection. Table 3 lists the initial state for the different A without the magnetic field, here Eu and Ev are introduced as

(b)

1

EB ; E0

-5

1.6

1.8

2

2.2

log10Ha

Fig. 15. Effect of Ha on the total kinetic energy: (a) h = 0° and (b) h = 90°.

2.4

2.6

1143

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144 Table 4 Value of b for different A and h at Ra = 104. A

0°

15°

30°

45°

60°

75°

90°

105°

120°

135°

150°

165°

0.5 0.75 1 1.5 2

3.44 3.49 3.53 3.60 3.65

3.52 3.67 3.68 3.73 3.76

3.50 3.72 3.81 3.83 3.84

3.65 3.73 3.75 3.87 3.87

3.68 3.75 3.80 3.80 3.82

3.64 3.67 3.69 3.74 3.77

3.58 3.57 3.54 3.49 3.45

3.63 3.65 3.65 3.63 3.59

3.67 3.71 3.71 3.60 3.57

3.64 3.67 3.61 3.70 3.75

3.48 3.63 3.71 3.76 3.79

3.50 3.62 3.65 3.71 3.74

8 > < log10 E ¼ 3:44log10 Ha þ 5:94; A ¼ 0:5; log10 E ¼ 3:53log10 Ha þ 6:24; A ¼ 1; > : log10 E ¼ 3:65log10 Ha þ 6:66; A ¼ 2

ð46Þ

for h = 0°, and

8 > < log10 E ¼ 3:58log10 Ha þ 5:20; A ¼ 0:5; log10 E ¼ 3:54log10 Ha þ 6:26; A ¼ 1; > : log10 E ¼ 3:45log10 Ha þ 7:31; A ¼ 2

ð47Þ

for h = 90°. Obviously, jbj increases with the increase of A for h = 0°, decreases with the increase of A for h = 90°, as expected. Table 4 lists the value of b for changing A and h, and indicates that 4 < b < 3, which agrees well with the results in 4.2. Furthermore, it is notable from Fig. 14 that for different aspect ratios A, the minimums of the restrain ratios, which also means the maximal inhibitory effects, are close to each other. And hence, relative to the strength of the magnetic field, the inclination angle is only the secondary factor. However, by changing the angle, we can always get almost the same maximum inhibitory effect for different aspect ratios.

the magnetic field is dominant. For different aspect ratios, the inclination angle of the magnetic field may lead to the different inhibitory effects. Although the total kinetic energy is inhibited proportionally to Ra2Hab, where b is between 4 and 3, it is found that it is reduced more significantly if the magnetic field is imposed in the direction where the perpendicular velocity component is dominant.  When the magnetic field is given, as seen from Fig. 14, the maximum inhibitory effect, viz. the smallest value of k is basically the same no matter what aspect ratios, though it may happen at the different inclination angles. Finally, it is necessary to point out that the present numerical method still have some limitations. For the 2D problems, the extension of the present method to the complicated zone beyond rectangular and the unsteady cases are in progress. After that, the method can be chosen as an efficient tool to study more interesting phenomena in MHD convection for theoretical research or industrial applications. For the 3D cases, it is noted that some early work for the streamfunction (vector potential)-velocity formulation has been taken [24], and therefore the thorough research would also be considered in the future.

5. Conclusions Acknowledgments The present study considers steady laminar two-dimensional MHD natural convection within an electrically conducting fluid filled rectangular cavity with adiabatic horizontal walls and differentially heated vertical walls, in the presence of the inclined magnetic field. As the discovery of the previous study, it is verified that application of magnetic field can reduce the convective heat transfer rate in the cavity and suppress different components of velocity with different rates. When the parameter c = Ha2/Ra is close to one, the isotherms are almost parallel to the vertical walls and the average Nusselt number Nu0 tends to one. When Hartmann number is sufficient high, the velocity component perpendicular to the magnetic field is damped following the approximation function u\ / Ra Ha2, whereas the parallel component varies as Ra Ha1.5. The main contributions of the present paper are as follows:  We have extended the second-order accurate compact difference method based on the streamfunction-velocity formulation of Navier–Stokes equations to the MHD convection problem, which avoids the drawback having to develop numerical boundary conditions for the vorticity due to the lack of the physical boundary conditions in the streamfunction–vorticity formulation. To the best of our knowledge, this is originally report about the implementation of the streamfunction–velocity formulation-based compact difference to solve the MHD convection problems.  We found that the inclination angles of the magnetic field and the aspect ratios of the cavity also play great role for the MHD convection problems. For large Hartmann numbers, the flow structure is highly determined by the angles, while the temperature distribution is seen to be no significant difference due to

This work was supported in part by the National Natural Science Foundation of China under Grant 11372075, 91330112 and 10972058, Research Fund for the Doctoral Program of Higher Education of China (No. 20100071110017), China Postdoctoral Science Foundation, the Teaching and Research Award Program for Outstanding Young Teachers in Higher Education Institutions of MOE, PR China.

References [1] C. Vives, C. Perry, Effects of magnetically damped convection during the controlled solidification of metals and alloys, Int. J. Heat Mass Transfer 30 (1987) 479–496. [2] H.P. Utech, M.C. Flemmings, Elimination of solute banding in indium antimonide crystals by growth in a magneticfield, J. Appl. Phys. 37 (1996) 2021–2024. [3] K. Okada, H. Ozoe, Experimental heat transfer rates of natural convection of molten gallium suppressed under an external magnetic field in either the X, Y, or Z direction, J. Heat Transfer 114 (1992) 107–114. [4] J.P. Garandet, J.P. Alboussiere, T. Moreau, Buoyancy driven convection in a rectangular cavity with a transverse magnetic field, Int. J. Heat Mass Transfer 35 (1992) 741–748. [5] P. Vasseur, M. Hasnaoui, E. Bilgen, L. Robillard, Natural convection in an inclined fluid layer with a transverse magnetic field: analogy with a porous medium, Trans. ASME J. Heat Transfer 117 (1995) 121–129. [6] H. Ozoe, M. Maruo, Magnetic and gravitational natural convection of melted silicon-two dimensional numerical computations for the rate of heat transfer, JSME Int. J. 30 (1987) 774–784. [7] H. Ozoe, K. Okada, The effect of the direction of the external magnetic field on the three-dimensional natural convection in a cubical enclosure, Int. J. Heat Mass Transfer 32 (1989) 1939–1954. [8] N. Rudraiah, R.M. Barron, M. Venkatachalappa, C.K. Subbaraya, Effect of a magnetic field on free convection in a rectangular enclosure, Int. J. Eng. Sci. 33 (1995) 1075–1084.

1144

P.X. Yu et al. / International Journal of Heat and Mass Transfer 67 (2013) 1131–1144

[9] H. Ben Hadid, D. Henry, S. Kaddeche, Numerical study of convection in the horizontal Bridgman configuration under the action of a constant magnetic field. Part 1: Two-dimensional flow, J. Fluid Mech. 333 (1997) 23–56. [10] H. Ben Hadid, D. Henry, Numerical study of convection in the horizontal Bridgman configuration under the action of a constant magnetic field. Part 2: Three-dimensional flow, J. Fluid Mech. 333 (1997) 57–83. [11] N.M. Al-Najem, K.M. Khanafer, M.M. El-Refaee, Numerical study of laminar natural convection in tilted enclosure with transverse magnetic field, Int. J. Numer. Methods Fluids 8 (1998) 651–672. [12] M. Pirmohammadi, M. Ghassemi, G.A. Sheikhzadeh, Effect of a magnetic field on buoyancy-driven convection in differentially heated square cavity, IEEE Trans. Magn. 45 (2009) 407–411. [13] M. Pirmohammadi, M. Ghassemi, Effect of magnetic field on convection heat transfer inside a tilted square enclosure, Int. Commun. Mass 36 (2009) 776– 780. [14] D.C. Lo, High-resolution simulations of magnetohydrodynamic free convection in an enclosure with a transverse magnetic field using a velocity–vorticity formulation, Int. Commun. Mass 37 (2010) 514–523. [15] M. Venkatachalappa, C.K. Subbaraya, Natural convection in a rectangular enclosure in the presence of a magnetic field with uniform heat flux from the side wall, Acta Mech. 96 (1993) 13–26. [16] P. Kandaswamy, S. Malliga Sundari, N. Nithyadevi, Magnetoconvection in an enclosure with partially active vertical walls, Int. J. Heat Mass Transfer 51 (2008) 1946–1954.

[17] M. Sathiyamoorthy, A.J. Chamkha, Effect of magnetic field on natural convection flow in a liquid gallium filled square cavity for linearly heated side wall(s), Int. J. Therm. Sci. 49 (2010) 1856–1865. [18] R. Kupferman, A central-difference scheme for a pure stream function formulation of incompressible viscous flows, SIAM J. Sci. Comput. 23 (2001) 1–18. [19] M. Ben-Artzi, J.-P. Croisille, D. Fishelov, S. Trachtenberg, A pure-compact scheme for the streamfunction formulation of Navier–Stokes equations, J. Comput. Phys. 205 (2005) 640–664. [20] M.M. Gupta, J.C. Kalita, A new paradigm for solving Navier–Stokes equations: streamfunction-velocity formulation, J. Comput. Phys. 207 (2005) 52–68. [21] M. Ben-Artzi, J.P. Croisille, D. Fishelov, A high order compact scheme for the pure-streamfunction formulation of the Navier–Stokes equations, J. Sci. Comput. 42 (2010) 216–250. [22] Z.F. Tian, P.X. Yu, An efficient compact difference scheme for solving the streamfunction formulation of the incompressible Navier–Stokes equations, J. Comput. Phys. 230 (2011) 6404–6419. [23] P.X. Yu, Z.F. Tian, Compact computations based on a stream-function–velocity formulation of two-dimensional steady laminar natural convection in a square cavity, Phys. Rev. E 85 (2012) 036703. [24] J.T. Holdeman, A velocity-stream function method for three-dimensional incompressible fluid flow, Comput. Methods Appl. Mech. Eng. 209–212 (2012) 66–73.