- Email: [email protected]
Numerical analysis of liquid metal MHD flows through circular pipes based on a fully developed modeling
Fusion Engineering and Design 88 (2013) 226–232
Contents lists available at SciVerse ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Numerical analysis of liquid metal MHD flows through circular pipes based on a fully developed modeling Xiujie Zhang ∗ , Chuanjie Pan, Zengyu Xu Southwestern Institute of Physics, Chengdu, Sichuan, China
h i g h l i g h t s
2D MHD code based on a fully developed modeling is developed and validated by Samad analytical results. The results of MHD effect of liquid metal through circular pipes at high Hartmann numbers are given. M type velocity profile is observed for MHD circular pipe flow at high wall conductance ratio condition. Non-uniform wall electrical conductivity leads to high jet velocity in Robert layers.
a r t i c l e
i n f o
Article history: Received 28 November 2011 Received in revised form 20 July 2012 Accepted 8 February 2013 Available online 20 March 2013 Keywords: MHD flow Circular pipe Numerical simulation Velocity profile
a b s t r a c t Magnetohydrodynamics (MHD) laminar flows through circular pipes are studied in this paper by numerical simulation under the conditions of Hartmann numbers from 18 to 10000. The code is developed based on a fully developed modeling and validated by Samad’s analytical solution and Chang’s asymptotic results. After the code validation, numerical simulation is extended to high Hartmann number for MHD circular pipe flows with conducting walls, and numerical results such as velocity distribution and MHD pressure gradient are obtained. Typical M-type velocity is observed but there is not such a big velocity jet as that of MHD rectangular duct flows even under the conditions of high Hartmann numbers and big wall conductance ratio. The over speed region in Robert layers becomes smaller when Hartmann numbers increase. When Hartmann number is fixed and wall conductance ratios change, the dimensionless velocity is through one point which is in agreement with Samad’s results, the locus of maximum value of velocity jet is same and effects of wall conductance ratio only on the maximum value of velocity jet. In case of Robert walls are treated as insulating and Hartmann walls as conducting for circular pipe MHD flows, there is big velocity jet like as MHD rectangular duct flows of Hunt’s case 2. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Liquid metal blanket have many attractive features such as low operating pressure, design simplicity and a convenient tritium breeding cycle, but MHD effect is a remaining key issue to be resolved. The mainly purpose of MHD effect study is to reduce the high MHD pressure drop and understand its velocity distribution. It is necessary to model MHD flows through Maxwell’s equations and the Navier–Stokes equation, in order to better understand the MHD effects of liquid metal under strong magnetic field. One of the most basic cases is the laminar MHD flow through circular pipe under a uniform strong magnetic field. Although there are some theoretical results available, such as analytical solutions for MHD pipe flows with insulating or conducting walls [1–3], but all those solutions are under the form of infinite series expansions
∗ Corresponding author. Tel.: +86 28 82850419. E-mail address: [email protected] (X. Zhang). 0920-3796/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.fusengdes.2013.02.032
involving modified Bessel functions, which limit those solutions to small Hartmann number (Ha). The square of Hartmann number is the ratio of Lorentz force comparing to viscous force, while the Hartmann number normally ranges from 103 to 105 in magnetic fusion reactor. There are some approximate solutions for high Hartmann number based on the asymptotic method [4–6], but this method is not full solution, there has some difference with analytical solution at small Ha condition which can be seen from Fig. 1. In the analytical solution Samad [3] considered the case of circular pipe with finite electric conductivity and finite wall thickness, which obtained the M-type velocity profile under the conditions of small Ha and big wall conductance ratio. However, in approximate approaches, Chang [4] did not observe the M-type velocity profile. Since the analytical solution of Samad is limited to small Ha which is about 30, does there exist big velocity jet like as MHD rectangular duct flows [9] when the Hartmann number becomes bigger to 104 ? Therefore it is important to extend study to high Hartmann numbers through numerical simulations.
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
227
where a is the inner radius of the circular pipe and v0 is the mean velocity and P = (ˇa2 )/v0 , Ha = aB0 /, and ˇ = −(∂p/∂z), where the non-dimensional quantity Ha is the Hartmann number and P is the Poiseuille number. For numerical simulations, we convert Eqs. (3) and (4) into cylinder coordinate where z-axis is along the axis of the circular cylinder: ∂2 V 1 ∂V 1 ∂2 V + + 2 + Ha 2 R ∂R ∂R R ∂ 2 ∂ ∂R
∂B ¯ ∂R
∂B 1 ∂ + + R ∂R R ∂
+ Ha
∂B sin ∂B cos − ∂R ∂ R ∂B R ∂
+1=0
(5)
∂V ∂V sin cos − ∂R ∂ R
=0
(6)
where is the ratio of the electrical conductivity of liquid to its local value and it is assumed that electrical conductivities of liquid and solid are isotropy. The dimensionless electrical current is calculated through the following equations: Fig. 1. The difference of velocity profiles between analytical solutions and approximate results.
In this paper, high resolution numerical simulation of MHD flows through circular pipes based on a fully developed modeling under the cylindrical coordinate is done to obtain the results of velocity distribution and MHD pressure drop under the conditions of Hartmann numbers from 18 to 10000. Numerical results with small Hartmann number are validated by Samad analytical results, and then give the results of velocity and induced electrical current distribution, MHD pressure drop under high Hartmann number conditions. 2. Mathematical modeling The unidirectional incompressible laminar flows are considered in this modeling as shown in Fig. 2a, where the magnetic field is along X direction and the flow is driven by a uniform pressure gradient in the Z direction. It can be assumed that there is only one component vz of velocity and only one component of induced magnetic field Bz due to the fluid is only flow along Z direction. The equations [3,5,13] governing this electrical conducting flow in the Cartesian coordinates are shown as below:
∂ 2 vz ∂2 vz + ∂x2 ∂y2
1 ∂ ∂x
1 ∂Bz x ∂x
−
∂p B0 ∂BZ + =0 ∂x ∂z
1 ∂ + ∂y
1 ∂Bz y ∂y
+ B0
(1) ∂ vZ =0 ∂x
(2)
where is the constant viscosity, is the magnetic permeability of vacuum, x and y are the electrical conductivities in x and y directions respectively. By introducing the dimensional variables: X=
x y z , Y= , Z= a a a
R=
r vz Bz , V= , B= √ a v0 P v0 P
Eqs. (1) and (2) can be rewritten in a dimensionless form as follows: ∂2 V ∂2 V ∂B + + Ha +1=0 ∂X ∂X 2 ∂Y 2
(3)
∂2 B ∂2 B ∂V =0 + + Ha ∂X ∂X 2 ∂Y 2
(4)
jx =
1 ∂Bz − → x Ha ∂y
(7)
jy =
1 ∂Bz − → y Ha ∂x
(8)
3. Numerical methods 3.1. Numerical scheme Similar to Reference [13], a control volume technique based on non-uniform collocated meshes is used to get the finite differential equations. The velocity and induced magnetic field are defined at the center of the control volume cell, and is taken at the sides of the cell. Eq. (5) is only solved in liquid area while Eq. (6) is solved in all area which contains liquid and solid. The code is written in Fortran90, uses Alternative Direction Implicit (ADI) method to solve the finite differential equations and contains an effective convergence acceleration technique same as Reference [13]. 3.2. Mesh and treatment of cylinder coordinate singularities For liquid metal MHD pipe flows the boundary layer is very thin with about Ha−1 non-dimensionless thickness. Under high Hartmann conditions it cannot use the uniform meshes for very large computation, so the non-uniform meshes related to Ha is employed in numerical simulations, which can insure that there are at least 7 mesh points within the boundary layer and 342 × 342 mesh points in the radial and azimuthal directions. The illustration of the computational mesh is shown in Fig. 2b. Due to the symmetry of the circular pipe it is reasonable to solve only one quarter of the pipe area, i.e. theta () from /2 to . Because of the singularity at the origin point where need special treatment, we choose the grids which are half-integer in radial direction and integer in azimuthal direction, i.e. Ri = (i − 3/2)R and j = (j − 1), where i = 2,3,. . .,N + 1 and j = 1,2,3,. . .,M. When R1 = 0 and Ri is computed from 2 to N in radial direction and from References [10,11], it can be seen that this method can avoid the singularity in radial direction. The integered mesh in radial direction ranges only from 0 to 0.8 since the non-uniform mesh is required to solve the thin boundary layer at high Hartmann numbers. If we compute from /2 to in azimuthal direction, the symmetry boundary conditions are V(/2 + )= V(/2) and V( − ) = V(), which is not correct from physical consideration and cannot get reasonable numerical results from computational test. Therefore, the computational area
228
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
Fig. 2. Sketch of the computational area and mesh: (a) the cross-section of the computational area, where the red area represents liquid metal, (b) non-uniform computational mesh used in numerical simulations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
is changed from /2 − to + and thereafter the symmetry boundary conditions are:
4. Validation of the code 4.1. Samad’s case
V (/2 − ) = V (/2 + ),
V ( − ) = V ( + ),
B(/2 − ) = −B(/2 + )
B( − ) = B( + )
This treatment is reasonable from physical understanding and validated to be correct for high Hartman numbers even Ha = 104 . There is no slip condition for the velocity at the interface of liquid metal and solid. At the interface of solid wall and air, the induced magnetic field is set to be zero.
To validate the code, the numerical results are compared with Samad analytical solutions under small ha conditions. Due to the appearance of small differences between large numbers, Samad limit his computation to Ha = 18. With the recent development of computing hardware, we can extend the computation to Ha = 30. The numerical results are compared with those of Samad analytical solution for Ha = 18 and Ha = 30 respectively, as seen in Fig. 3. As mentioned above, a is the inner radius of the circular pipe, alpha*a is defined as the outer radius of the pipe, and C takes the electrical
Fig. 3. Comparison of the velocity profiles between numerical results and those of Samad analytical solutions under the conditions of small Ha: (a) Ha = 18 and (b) Ha = 30.
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
229
Fig. 4. Comparison of the velocity profiles between numerical results and those of Chang’s asymptotic solution at Ha = 1000 with all insulating walls: (a) compare the velocity profiles at = /2, (b) velocity distribution in the cross section (3D display).
conductivity of the solid to liquid metal as same as that in Samad paper. From the comparison in Fig. 3, it can be seen that numerical results match well with those of Samad analytical solutions and the M-type velocity is observed from both results.
Even the Ha changes the velocity profile in this direction is always parabola-shape.
4.2. MHD pipe flow with all insulating walls
Because of the computational difficulty of Bessel function in Samad’s analytical solution at Hartmann number greater than 30, the numerical method is validated by Samad analytical solution at small Ha and then extended for simulation of MHD flows in a circular duct with velocity and pressure drop obtained at high Hartmann number.
In the case of MHD flow through circular pipe with all insulating walls, there is no difference between Samad analytical and Chang asymptotic results. Then the code is validated against Chang’s results at high Ha [12]. It is seen from the comparison in Fig. 4a that numerical results match well with those of asymptotic solution. Fig. 4b shows the non-dimensional velocity distribution in the cross section of circular pipe. There is an interesting result that in the line of x = 0 ( = /2) the velocity profile is parabola distribution, which is not similar to that of MHD rectangular duct flow with all insulating walls but that of normal fluid duct flows.
5. Results and discussion
5.1. Electrical current and velocity distribution in the cross section In Figs. 5 and 6, Cw is the wall conductance ratio defined as follows:
Fig. 5. Induced electrical current distribution in the cross-section of the MHD pipe flows.
230
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
Fig. 6. Velocity distribution in the cross section of the MHD pipe flows at different Ha (3D display).
Fig. 7. Velocity profiles in the line of = /2 changing with Hartmann numbers.
Cw = (w tw )/(f a), where w , f , tw are electrical conductivity of solid wall and liquid metal, solid wall thickness, respectively. It can be seen from electrical current distribution results in Fig. 5 that there are two circuits if we convert the current distribution to the whole circular pipe domain by symmetry, which is reasonable in physical understanding. Small changes of the velocity over-speed regions in the Robert layers [7,8,12] are seen when Ha increases, as shown in Figs. 6 and 7. There has no big velocity jet in MHD circular pipe flows even when Ha = 104 and Cw = 1.0, which is different from that of MHD rectangular duct flows, because the difference of electrical current distributions in the cross section between rectangular and circular pipes, there is no obvious area where the electrical current is parallel to the imposed magnetic field in MHD circular pipe flows. The most interesting result is that the velocity distribution changes with different Cw at Ha = 1000, as shown in Fig. 8. All dimensionless velocities go through one point, which is in agreement with Samad results at Ha = 18. Another interesting result is that with the fixed Ha number, the locus of the maximum non-dimensional velocity does not change with the wall conductance ratio (Cw ), which only affects the value of maximum jet velocity.
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
Fig. 8. Velocity distribution in the line of = /2 vs. wall conductance ratio (Cw ) at Ha = 1000.
231
Fig. 9. Numerical results of dimensionless pressure gradient compared with theoretical results.
compare the pressure gradient with numerical results, the following formula [14] is used for the theoretical results: 5.2. Compare MHD pressure drop with theory − In numerical simulation the pressure gradient can be computed using:
−
dp v0 = 4Q dz
derived from Q =
(9)
1 /2
0
V (R, )dRd, where Q is the dimension-
less flow rate related to the dimensional pressure gradient. To
Cw dp (dimensionless) = 1 + Cw dz theory
(10)
Formula (9) is then transformed to the following dimensionless format: −
dp (dimensionless) = dz numerical 4Qf B02
(11)
where f and B0 are electrical conductivity of liquid metal and imposed magnetic field, respectively. It can be seen from comparison shown in Fig. 9 that the dimensionless pressure gradient
Fig. 10. MHD flows through circular pipe with non-uniform wall electrical conductivity: (a) velocity distribution in the cross-section (3D display), (b) electrical current distribution in the cross-section, Ha = 1000, Cw = 0.016.
232
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
change smaller with Hartmann number becoming bigger, when Hartmann number is equal to 100 the dimensionless pressure gradient is much close to that of theory values, while under Ha = 18 and Ha = 30 conditions there are big differences between numerical and theoretical results. The reason is because the theoretical formula (9) is suitable when Ha 1. In addition, it is obvious that the theoretical formula only considers the effects of wall electrical conductivity and the pipe shape on the pressure gradient, not considers effect of the changing Lorentz force on it, i.e. the dimensionless pressure gradient should also be the function of Hartmann number. 6. Effect of non-uniform wall electrical conductivity In Hunt’s case 2 [9] of rectangular duct flows, Hartmann walls were treated as conducting and side walls as insulating. There also have Hartmann and Robert walls for MHD circular pipe flows, the Robert walls in circular pipe are like as the side walls in rectangular duct. When we treated Robert walls as insulating and Hartmann walls as conducting in MHD circular pipe flows, does there have big velocity jet like as Hunt case 2? In numerical simulation the wall from = /2 to = 1.175/2 (i.e. 1/5.7 of whole computational wall area, approximately in Robert wall) is treated as electrical insulating and other walls as conducting. It can be seen from numerical results in Fig. 10 that in Robert layer there is big velocity jet like as Hunt’s case 2, which is because there has obvious electric current distribution which is parallel to the imposed magnetic field in Robert layer where the Lorentz force is very small, so in this area form big velocity jet like as MHD rectangular duct flows with conducting walls. 7. Conclusions In the case of MHD flows through circular pipes at high values of Hartmann number and wall conductance ratio (Cw ), there is a typical M-type velocity profile as observed by Hunt in Reference [9]. Compared with the case of MHD rectangular duct flows, the maximum of jet velocity in MHD circular pipe flows is smaller than that in rectangular duct flows. Another difference lies in the effect of wall conductance ratio on velocity profile. In MHD circular pipe flows, if the Hartmann number is fixed, for different wall conductance ratios the dimensionless velocity profiles all through one point and the locus of maximum velocity jet is same. Cw only has the effects on the maximum value of the velocity jet in Robert layers.
From the comparison of MHD pressure gradients, it can be seen that numerical results have some differences from theoretical results since the theory does not consider the effect of magnetic field on the pressure gradient, which should be the function of wall conductance ratio and Hartmann numbers. When the electrical conductivity approximately in Robert walls are treated as insulating and other walls as conducting, there has big velocity jet like as MHD rectangular duct flows with conducting walls. Acknowledgments Part of this work is supported by China National Nature Science Fund Grant No. 11105044; the authors would like to express gratitude to professor Zongze Mu and Dr Jianzhou Zhu for their helpful suggestions. References [1] R.R. Gold, Magnetohydrodynamic pipe flow. Part 1, Journal of Fluid Mechanics 13 (1962) 505–512. [2] S. Ihara, T. Kiyohiro, A. Matsushima, The flow of conducting fluids in circular pipes with finite conductivity under uniform transverse magnetic fields, Journal of Applied Mechanics 34 (1) (1967) 29–36. [3] S. Samad, The flow of conducting fluids through circular pipes having finite conductivity and finite thickness under uniform transverse magnetic fields, International Journal of Engineering Science 19 (1981) 1221–1232. [4] C. Chang, S. Lundgren, Duct flow in magnetohydrodynamics, Zeitschrift für angewandte Mathematik und Physik XII 10 (1961) 100–114. [5] A. Shercliff, The flow of conducting fluids in circular pipes under transverse magnetic fields, Journal of Fluid Mechanics 1 (1956) 644–666. [6] A. Shercliff, Magnetohydrodynamic pipe flow. Part 2. High Hartmann number, Journal of Fluid Mechanics 13 (1962) 513–518. [7] P.H. Roberts, An Introduction to Magnetohydrodynamics, American Elsevier Publishing Company, Inc., New York, 1967. [8] P.H. Roberts, Singularities of Hartmann layers, Proceedings of the Royal Society A 300 (1967) 94–107. [9] J.C.R. Hunt, Magnetohydrodynamic flow in rectangular ducts, Journal of Fluid Mechanics 21 (4) (1965) 577–590. [10] K. Mohseni, T. Colonius, Numerical treatment of polar coordinate singularities, Journal of Computational Physics 157 (2000) 787–795. [11] M.-C. Lai, W.-C. Wang, Fast direct solvers for Poisson equation on 2D polar and spherical geometries, Numerical Methods for Partial Differential Equations 18 (2002) 56–68. [12] S. Vantieghem, X. Albets-Chico, B. Knaepen, The velocity profile of laminar MHD flows in circular conducting pipes, Theoretical and Computational Fluid Dynamics 23 (2009) 525–533. [13] S. Smolentsev, N. Morley, M. Abdou, Code development for analysis of MHD pressure drop reduction in a liquid metal blanket using insulation technique based on a fully developed flow model, Fusion Engineering and Design 73 (2005) 83–93. [14] I.R. Kirillov, C.B. Reed, et al., Present understanding of MHD and heat transfer phenomena for liquid metal blankets, Fusion Engineering and Design 27 (1995) 553–569.
Contents lists available at SciVerse ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Numerical analysis of liquid metal MHD flows through circular pipes based on a fully developed modeling Xiujie Zhang ∗ , Chuanjie Pan, Zengyu Xu Southwestern Institute of Physics, Chengdu, Sichuan, China
h i g h l i g h t s
2D MHD code based on a fully developed modeling is developed and validated by Samad analytical results. The results of MHD effect of liquid metal through circular pipes at high Hartmann numbers are given. M type velocity profile is observed for MHD circular pipe flow at high wall conductance ratio condition. Non-uniform wall electrical conductivity leads to high jet velocity in Robert layers.
a r t i c l e
i n f o
Article history: Received 28 November 2011 Received in revised form 20 July 2012 Accepted 8 February 2013 Available online 20 March 2013 Keywords: MHD flow Circular pipe Numerical simulation Velocity profile
a b s t r a c t Magnetohydrodynamics (MHD) laminar flows through circular pipes are studied in this paper by numerical simulation under the conditions of Hartmann numbers from 18 to 10000. The code is developed based on a fully developed modeling and validated by Samad’s analytical solution and Chang’s asymptotic results. After the code validation, numerical simulation is extended to high Hartmann number for MHD circular pipe flows with conducting walls, and numerical results such as velocity distribution and MHD pressure gradient are obtained. Typical M-type velocity is observed but there is not such a big velocity jet as that of MHD rectangular duct flows even under the conditions of high Hartmann numbers and big wall conductance ratio. The over speed region in Robert layers becomes smaller when Hartmann numbers increase. When Hartmann number is fixed and wall conductance ratios change, the dimensionless velocity is through one point which is in agreement with Samad’s results, the locus of maximum value of velocity jet is same and effects of wall conductance ratio only on the maximum value of velocity jet. In case of Robert walls are treated as insulating and Hartmann walls as conducting for circular pipe MHD flows, there is big velocity jet like as MHD rectangular duct flows of Hunt’s case 2. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Liquid metal blanket have many attractive features such as low operating pressure, design simplicity and a convenient tritium breeding cycle, but MHD effect is a remaining key issue to be resolved. The mainly purpose of MHD effect study is to reduce the high MHD pressure drop and understand its velocity distribution. It is necessary to model MHD flows through Maxwell’s equations and the Navier–Stokes equation, in order to better understand the MHD effects of liquid metal under strong magnetic field. One of the most basic cases is the laminar MHD flow through circular pipe under a uniform strong magnetic field. Although there are some theoretical results available, such as analytical solutions for MHD pipe flows with insulating or conducting walls [1–3], but all those solutions are under the form of infinite series expansions
∗ Corresponding author. Tel.: +86 28 82850419. E-mail address: [email protected] (X. Zhang). 0920-3796/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.fusengdes.2013.02.032
involving modified Bessel functions, which limit those solutions to small Hartmann number (Ha). The square of Hartmann number is the ratio of Lorentz force comparing to viscous force, while the Hartmann number normally ranges from 103 to 105 in magnetic fusion reactor. There are some approximate solutions for high Hartmann number based on the asymptotic method [4–6], but this method is not full solution, there has some difference with analytical solution at small Ha condition which can be seen from Fig. 1. In the analytical solution Samad [3] considered the case of circular pipe with finite electric conductivity and finite wall thickness, which obtained the M-type velocity profile under the conditions of small Ha and big wall conductance ratio. However, in approximate approaches, Chang [4] did not observe the M-type velocity profile. Since the analytical solution of Samad is limited to small Ha which is about 30, does there exist big velocity jet like as MHD rectangular duct flows [9] when the Hartmann number becomes bigger to 104 ? Therefore it is important to extend study to high Hartmann numbers through numerical simulations.
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
227
where a is the inner radius of the circular pipe and v0 is the mean velocity and P = (ˇa2 )/v0 , Ha = aB0 /, and ˇ = −(∂p/∂z), where the non-dimensional quantity Ha is the Hartmann number and P is the Poiseuille number. For numerical simulations, we convert Eqs. (3) and (4) into cylinder coordinate where z-axis is along the axis of the circular cylinder: ∂2 V 1 ∂V 1 ∂2 V + + 2 + Ha 2 R ∂R ∂R R ∂ 2 ∂ ∂R
∂B ¯ ∂R
∂B 1 ∂ + + R ∂R R ∂
+ Ha
∂B sin ∂B cos − ∂R ∂ R ∂B R ∂
+1=0
(5)
∂V ∂V sin cos − ∂R ∂ R
=0
(6)
where is the ratio of the electrical conductivity of liquid to its local value and it is assumed that electrical conductivities of liquid and solid are isotropy. The dimensionless electrical current is calculated through the following equations: Fig. 1. The difference of velocity profiles between analytical solutions and approximate results.
In this paper, high resolution numerical simulation of MHD flows through circular pipes based on a fully developed modeling under the cylindrical coordinate is done to obtain the results of velocity distribution and MHD pressure drop under the conditions of Hartmann numbers from 18 to 10000. Numerical results with small Hartmann number are validated by Samad analytical results, and then give the results of velocity and induced electrical current distribution, MHD pressure drop under high Hartmann number conditions. 2. Mathematical modeling The unidirectional incompressible laminar flows are considered in this modeling as shown in Fig. 2a, where the magnetic field is along X direction and the flow is driven by a uniform pressure gradient in the Z direction. It can be assumed that there is only one component vz of velocity and only one component of induced magnetic field Bz due to the fluid is only flow along Z direction. The equations [3,5,13] governing this electrical conducting flow in the Cartesian coordinates are shown as below:
∂ 2 vz ∂2 vz + ∂x2 ∂y2
1 ∂ ∂x
1 ∂Bz x ∂x
−
∂p B0 ∂BZ + =0 ∂x ∂z
1 ∂ + ∂y
1 ∂Bz y ∂y
+ B0
(1) ∂ vZ =0 ∂x
(2)
where is the constant viscosity, is the magnetic permeability of vacuum, x and y are the electrical conductivities in x and y directions respectively. By introducing the dimensional variables: X=
x y z , Y= , Z= a a a
R=
r vz Bz , V= , B= √ a v0 P v0 P
Eqs. (1) and (2) can be rewritten in a dimensionless form as follows: ∂2 V ∂2 V ∂B + + Ha +1=0 ∂X ∂X 2 ∂Y 2
(3)
∂2 B ∂2 B ∂V =0 + + Ha ∂X ∂X 2 ∂Y 2
(4)
jx =
1 ∂Bz − → x Ha ∂y
(7)
jy =
1 ∂Bz − → y Ha ∂x
(8)
3. Numerical methods 3.1. Numerical scheme Similar to Reference [13], a control volume technique based on non-uniform collocated meshes is used to get the finite differential equations. The velocity and induced magnetic field are defined at the center of the control volume cell, and is taken at the sides of the cell. Eq. (5) is only solved in liquid area while Eq. (6) is solved in all area which contains liquid and solid. The code is written in Fortran90, uses Alternative Direction Implicit (ADI) method to solve the finite differential equations and contains an effective convergence acceleration technique same as Reference [13]. 3.2. Mesh and treatment of cylinder coordinate singularities For liquid metal MHD pipe flows the boundary layer is very thin with about Ha−1 non-dimensionless thickness. Under high Hartmann conditions it cannot use the uniform meshes for very large computation, so the non-uniform meshes related to Ha is employed in numerical simulations, which can insure that there are at least 7 mesh points within the boundary layer and 342 × 342 mesh points in the radial and azimuthal directions. The illustration of the computational mesh is shown in Fig. 2b. Due to the symmetry of the circular pipe it is reasonable to solve only one quarter of the pipe area, i.e. theta () from /2 to . Because of the singularity at the origin point where need special treatment, we choose the grids which are half-integer in radial direction and integer in azimuthal direction, i.e. Ri = (i − 3/2)R and j = (j − 1), where i = 2,3,. . .,N + 1 and j = 1,2,3,. . .,M. When R1 = 0 and Ri is computed from 2 to N in radial direction and from References [10,11], it can be seen that this method can avoid the singularity in radial direction. The integered mesh in radial direction ranges only from 0 to 0.8 since the non-uniform mesh is required to solve the thin boundary layer at high Hartmann numbers. If we compute from /2 to in azimuthal direction, the symmetry boundary conditions are V(/2 + )= V(/2) and V( − ) = V(), which is not correct from physical consideration and cannot get reasonable numerical results from computational test. Therefore, the computational area
228
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
Fig. 2. Sketch of the computational area and mesh: (a) the cross-section of the computational area, where the red area represents liquid metal, (b) non-uniform computational mesh used in numerical simulations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
is changed from /2 − to + and thereafter the symmetry boundary conditions are:
4. Validation of the code 4.1. Samad’s case
V (/2 − ) = V (/2 + ),
V ( − ) = V ( + ),
B(/2 − ) = −B(/2 + )
B( − ) = B( + )
This treatment is reasonable from physical understanding and validated to be correct for high Hartman numbers even Ha = 104 . There is no slip condition for the velocity at the interface of liquid metal and solid. At the interface of solid wall and air, the induced magnetic field is set to be zero.
To validate the code, the numerical results are compared with Samad analytical solutions under small ha conditions. Due to the appearance of small differences between large numbers, Samad limit his computation to Ha = 18. With the recent development of computing hardware, we can extend the computation to Ha = 30. The numerical results are compared with those of Samad analytical solution for Ha = 18 and Ha = 30 respectively, as seen in Fig. 3. As mentioned above, a is the inner radius of the circular pipe, alpha*a is defined as the outer radius of the pipe, and C takes the electrical
Fig. 3. Comparison of the velocity profiles between numerical results and those of Samad analytical solutions under the conditions of small Ha: (a) Ha = 18 and (b) Ha = 30.
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
229
Fig. 4. Comparison of the velocity profiles between numerical results and those of Chang’s asymptotic solution at Ha = 1000 with all insulating walls: (a) compare the velocity profiles at = /2, (b) velocity distribution in the cross section (3D display).
conductivity of the solid to liquid metal as same as that in Samad paper. From the comparison in Fig. 3, it can be seen that numerical results match well with those of Samad analytical solutions and the M-type velocity is observed from both results.
Even the Ha changes the velocity profile in this direction is always parabola-shape.
4.2. MHD pipe flow with all insulating walls
Because of the computational difficulty of Bessel function in Samad’s analytical solution at Hartmann number greater than 30, the numerical method is validated by Samad analytical solution at small Ha and then extended for simulation of MHD flows in a circular duct with velocity and pressure drop obtained at high Hartmann number.
In the case of MHD flow through circular pipe with all insulating walls, there is no difference between Samad analytical and Chang asymptotic results. Then the code is validated against Chang’s results at high Ha [12]. It is seen from the comparison in Fig. 4a that numerical results match well with those of asymptotic solution. Fig. 4b shows the non-dimensional velocity distribution in the cross section of circular pipe. There is an interesting result that in the line of x = 0 ( = /2) the velocity profile is parabola distribution, which is not similar to that of MHD rectangular duct flow with all insulating walls but that of normal fluid duct flows.
5. Results and discussion
5.1. Electrical current and velocity distribution in the cross section In Figs. 5 and 6, Cw is the wall conductance ratio defined as follows:
Fig. 5. Induced electrical current distribution in the cross-section of the MHD pipe flows.
230
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
Fig. 6. Velocity distribution in the cross section of the MHD pipe flows at different Ha (3D display).
Fig. 7. Velocity profiles in the line of = /2 changing with Hartmann numbers.
Cw = (w tw )/(f a), where w , f , tw are electrical conductivity of solid wall and liquid metal, solid wall thickness, respectively. It can be seen from electrical current distribution results in Fig. 5 that there are two circuits if we convert the current distribution to the whole circular pipe domain by symmetry, which is reasonable in physical understanding. Small changes of the velocity over-speed regions in the Robert layers [7,8,12] are seen when Ha increases, as shown in Figs. 6 and 7. There has no big velocity jet in MHD circular pipe flows even when Ha = 104 and Cw = 1.0, which is different from that of MHD rectangular duct flows, because the difference of electrical current distributions in the cross section between rectangular and circular pipes, there is no obvious area where the electrical current is parallel to the imposed magnetic field in MHD circular pipe flows. The most interesting result is that the velocity distribution changes with different Cw at Ha = 1000, as shown in Fig. 8. All dimensionless velocities go through one point, which is in agreement with Samad results at Ha = 18. Another interesting result is that with the fixed Ha number, the locus of the maximum non-dimensional velocity does not change with the wall conductance ratio (Cw ), which only affects the value of maximum jet velocity.
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
Fig. 8. Velocity distribution in the line of = /2 vs. wall conductance ratio (Cw ) at Ha = 1000.
231
Fig. 9. Numerical results of dimensionless pressure gradient compared with theoretical results.
compare the pressure gradient with numerical results, the following formula [14] is used for the theoretical results: 5.2. Compare MHD pressure drop with theory − In numerical simulation the pressure gradient can be computed using:
−
dp v0 = 4Q dz
derived from Q =
(9)
1 /2
0
V (R, )dRd, where Q is the dimension-
less flow rate related to the dimensional pressure gradient. To
Cw dp (dimensionless) = 1 + Cw dz theory
(10)
Formula (9) is then transformed to the following dimensionless format: −
dp (dimensionless) = dz numerical 4Qf B02
(11)
where f and B0 are electrical conductivity of liquid metal and imposed magnetic field, respectively. It can be seen from comparison shown in Fig. 9 that the dimensionless pressure gradient
Fig. 10. MHD flows through circular pipe with non-uniform wall electrical conductivity: (a) velocity distribution in the cross-section (3D display), (b) electrical current distribution in the cross-section, Ha = 1000, Cw = 0.016.
232
X. Zhang et al. / Fusion Engineering and Design 88 (2013) 226–232
change smaller with Hartmann number becoming bigger, when Hartmann number is equal to 100 the dimensionless pressure gradient is much close to that of theory values, while under Ha = 18 and Ha = 30 conditions there are big differences between numerical and theoretical results. The reason is because the theoretical formula (9) is suitable when Ha 1. In addition, it is obvious that the theoretical formula only considers the effects of wall electrical conductivity and the pipe shape on the pressure gradient, not considers effect of the changing Lorentz force on it, i.e. the dimensionless pressure gradient should also be the function of Hartmann number. 6. Effect of non-uniform wall electrical conductivity In Hunt’s case 2 [9] of rectangular duct flows, Hartmann walls were treated as conducting and side walls as insulating. There also have Hartmann and Robert walls for MHD circular pipe flows, the Robert walls in circular pipe are like as the side walls in rectangular duct. When we treated Robert walls as insulating and Hartmann walls as conducting in MHD circular pipe flows, does there have big velocity jet like as Hunt case 2? In numerical simulation the wall from = /2 to = 1.175/2 (i.e. 1/5.7 of whole computational wall area, approximately in Robert wall) is treated as electrical insulating and other walls as conducting. It can be seen from numerical results in Fig. 10 that in Robert layer there is big velocity jet like as Hunt’s case 2, which is because there has obvious electric current distribution which is parallel to the imposed magnetic field in Robert layer where the Lorentz force is very small, so in this area form big velocity jet like as MHD rectangular duct flows with conducting walls. 7. Conclusions In the case of MHD flows through circular pipes at high values of Hartmann number and wall conductance ratio (Cw ), there is a typical M-type velocity profile as observed by Hunt in Reference [9]. Compared with the case of MHD rectangular duct flows, the maximum of jet velocity in MHD circular pipe flows is smaller than that in rectangular duct flows. Another difference lies in the effect of wall conductance ratio on velocity profile. In MHD circular pipe flows, if the Hartmann number is fixed, for different wall conductance ratios the dimensionless velocity profiles all through one point and the locus of maximum velocity jet is same. Cw only has the effects on the maximum value of the velocity jet in Robert layers.
From the comparison of MHD pressure gradients, it can be seen that numerical results have some differences from theoretical results since the theory does not consider the effect of magnetic field on the pressure gradient, which should be the function of wall conductance ratio and Hartmann numbers. When the electrical conductivity approximately in Robert walls are treated as insulating and other walls as conducting, there has big velocity jet like as MHD rectangular duct flows with conducting walls. Acknowledgments Part of this work is supported by China National Nature Science Fund Grant No. 11105044; the authors would like to express gratitude to professor Zongze Mu and Dr Jianzhou Zhu for their helpful suggestions. References [1] R.R. Gold, Magnetohydrodynamic pipe flow. Part 1, Journal of Fluid Mechanics 13 (1962) 505–512. [2] S. Ihara, T. Kiyohiro, A. Matsushima, The flow of conducting fluids in circular pipes with finite conductivity under uniform transverse magnetic fields, Journal of Applied Mechanics 34 (1) (1967) 29–36. [3] S. Samad, The flow of conducting fluids through circular pipes having finite conductivity and finite thickness under uniform transverse magnetic fields, International Journal of Engineering Science 19 (1981) 1221–1232. [4] C. Chang, S. Lundgren, Duct flow in magnetohydrodynamics, Zeitschrift für angewandte Mathematik und Physik XII 10 (1961) 100–114. [5] A. Shercliff, The flow of conducting fluids in circular pipes under transverse magnetic fields, Journal of Fluid Mechanics 1 (1956) 644–666. [6] A. Shercliff, Magnetohydrodynamic pipe flow. Part 2. High Hartmann number, Journal of Fluid Mechanics 13 (1962) 513–518. [7] P.H. Roberts, An Introduction to Magnetohydrodynamics, American Elsevier Publishing Company, Inc., New York, 1967. [8] P.H. Roberts, Singularities of Hartmann layers, Proceedings of the Royal Society A 300 (1967) 94–107. [9] J.C.R. Hunt, Magnetohydrodynamic flow in rectangular ducts, Journal of Fluid Mechanics 21 (4) (1965) 577–590. [10] K. Mohseni, T. Colonius, Numerical treatment of polar coordinate singularities, Journal of Computational Physics 157 (2000) 787–795. [11] M.-C. Lai, W.-C. Wang, Fast direct solvers for Poisson equation on 2D polar and spherical geometries, Numerical Methods for Partial Differential Equations 18 (2002) 56–68. [12] S. Vantieghem, X. Albets-Chico, B. Knaepen, The velocity profile of laminar MHD flows in circular conducting pipes, Theoretical and Computational Fluid Dynamics 23 (2009) 525–533. [13] S. Smolentsev, N. Morley, M. Abdou, Code development for analysis of MHD pressure drop reduction in a liquid metal blanket using insulation technique based on a fully developed flow model, Fusion Engineering and Design 73 (2005) 83–93. [14] I.R. Kirillov, C.B. Reed, et al., Present understanding of MHD and heat transfer phenomena for liquid metal blankets, Fusion Engineering and Design 27 (1995) 553–569.