- Email: [email protected]
Code development and validation for analyzing liquid metal MHD flow in rectangular ducts
Fusion Engineering and Design 85 (2010) 1736–1741
Contents lists available at ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Code development and validation for analyzing liquid metal MHD flow in rectangular ducts Tao Zhou a,c,∗ , Zhiyi Yang a , Mingjiu Ni c , Hongli Chen a,b a b c
Institute of Plasma Physics, Chinese Academy of Sciences, Hefei, Anhui 230031, China School of Nuclear Science and Technology, University of Science and Technology of China, Hefei, Anhui 230027, China Graduate University of Chinese Academy of Sciences, Beijing 100049, China
a r t i c l e
i n f o
Keywords: Numerical simulations Liquid metal magnetohydrodynamics Current density conservative scheme
a b s t r a c t A code named MTC-H 1.0 which can simulate 3D magnetohydrodynamics (MHD) flow in rectangular ducts has been developed by FDS Team. In this code, a conservative scheme of the current density was employed for calculation of the induced current and the Lorentz force, and the consistent projection method was employed for solving the incompressible Navier–Stokes equations with the Lorentz force included as a source term. The code was developed on a structured collocated grid, on which velocity, pressure, and electrical potential were located in the cell center, while current fluxes were located on the cell faces. The code MTC-H 1.0 was test by simulating MHD rectangular ducts flows, and the Shercliff’s insulated walls and Hunt’s conductive walls were used as benchmark models. The validation cases were conducted with Hartmann number from 102 to 104 on structured collocated grids. The results matched well with Hunt’s and Shercliff’s analytical solutions and it showed good accuracy of the code. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Liquid metal breeder blanket concept has been a topic of great interest in fusion reactor blanket design [1–5] because of many advantages. A series of fusion reactors (named FDS series) have been designed and assessed by the FDS Team in China [6–10]. And a series of liquid metal blanket concepts [11–13] have been developed correspondingly. In these blankets, lead–lithium alloy is considered as tritium breeder only or as breeder and coolant. But the motion of liquid metal in fusion reactor strong magnetic field cause serious magnetohydrodynamic (MHD) effects, which have dramatic impacts on velocity distribution, heat transfer characteristics, pressure drop and the required pumping power for the cooling system. Therefore, the liquid metal MHD effect under fusion relevant condition is one of the key issues in making an optimal design for LM blanket [14–18]. Up to now, liquid metal MHD duct flow has been extensively studied by theory [19–22] and numerical simulation [23–28]. Due to the complex geometry of flow channel in blankets, it is difficult to study the detail of MHD flow characteristics by theoretical analysis. As the development of computer science, the numerical
∗ Corresponding author at: Institute of Plasma Physics, Chinese Academy of Sciences, Hefei, Anhui 230031, China. Tel.: +86 551 5593296; fax: +86 551 5591397. E-mail address: [email protected] (T. Zhou). 0920-3796/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.fusengdes.2010.05.034
simulation has become an effective tool for analyzing MHD effects and optimizing blanket design. For MHD duct flow, the Hartmann layer perpendicular to the magnetic field scales with Ha−1 , and the side layer parallel to the magnetic field scales with Ha−1/2 . In fusion blankets, Hartmann number can be as high as 104 –105 , which leads to some special numerical difficulty. Therefore, simulating and studying liquid metal MHD flow at high Hartmann number has become a hot research topic for fusion blanket design. A number of numerical codes have been developed for analyzing liquid metal MHD effects [23–28]. Smolentsev developed a 2D code based on structured grid and induced magnetic equation [25]. A 3D code named HIMAG [26] and another 3D code based on the commercial package CFX 5.6 [27] were developed based on unstructured grid and electrical potential equation. By solving induced magnetic equations, it is difficult to construct a proper boundary condition for induced magnetic field. For fusion reactor conditions, the magnetic Reynolds numbers is small, the electrical potential formula can be employed for MHD with good accuracy [29–31]. Ni developed a consistent and conservative scheme [32] for calculate the current density and the Lorentz force by solving electrical potential equations, this scheme can solve LM MHD problem at high Hartmann number. In order to research and evaluate MHD effect in liquid metal blanket as well as make an optimal design for FDS series blankets, FDS team has been developing the MHD simulation code named MTC (Magnetic Thermo-hydraulics Coupling Code) for years. MTC-F 1.0 and MTC-F 2.0 are developed
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741
based on the commercial code FLUENT, which can solve LM MHD problem only at low Hartman number (Ha < 500). In this paper, a numerical code named MTC-H 1.0 for simulating MHD flow in rectangular duct is developed based on the structure collocated grid, in which consistent and conservative scheme are used to calculate the current density and the Lorentz force. We validated the numerical code by simulating 2D and 3D MHD flows with analytical solutions existed. The validation cases are conducted with Hartmann number from 102 to 104 .
1737
Table 1 Computed flow rate. Hartmann number
Calculated flow rate
Error
Shercliff’s case
500 1000
3.979 3.964
0.5% 0.9%
Hunt’s case
500 1000
4.008 4.016
0.2% 0.4%
Ohm’s law: J = −∇ ϕ + u × B
2. Governing equations and numerical methods
Conservation of current density:
2.1. Governing equations LM is a kind of incompressible fluid. Under fusion blanket conditions, the magnetic Reynolds numbers is small, this means the induced magnetic field can be negligible compared to the apply magnetic field, the governing equations for LM MHD flow in dimensionless form can be expressed as follows: Conservation of momentum: ∂u 1 2 + u · ∇ u = −∇ p + ∇ u + NJ × B Re ∂t
(1)
Conservation of mass (Navier–Stokes equation):
∇ ·u=0
(3)
(2)
∇ ·J =0
(4)
From Eqs. (3) and (4),we can get the Poisson’s equation for electrical potential:
∇ · (∇ ϕ) = ∇ · (u × B)
(5)
Here t, u, p, J, B, ϕ are dimensionless time, velocity, pressure, current density, applied magnetic field, and electrical potential. They are scaled with L/u0 , u0 , 0 u20 , u0 B0 , B0 and Lu0 B0 , respectively, where L is a characteristic length, u0 is characteristic velocity, and 0 are the conductivity and density of fluid. The Reynolds number is Re = u0 L/, the Hartmann number is Ha = B0 L /0 , and the interaction parameter is N = Ha2 /Re.
Fig. 1. Calculating results for Shercliff’s fully developed flow at Ha = 1000: (a) convergence of flow rate, (b) current in the duct, (c) velocity distribution along the middle line normal to the Hartmann walls, and (d) velocity distribution along the middle line normal to the side walls.
1738
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741
Fig. 2. Calculating results for Hunt’s fully developed flow at Ha = 1000: (a) convergence of flow rate, (b) current in the duct, (c) velocity distribution along the middle line normal to the Hartmann walls, and (d) velocity distribution along the middle line normal to the side walls.
2.2. Numerical scheme and boundary conditions The code was developed based on the structured grid. Since the Hartmann layer and the side layer are very thin, and require at least 4–5 grid points to resolve these layers, a non-uniform grid system is needed to improve the computational efficiency. The control volume method on the structured collocated grid system was used to ensure the conservation of the discrete equations. All the unknown variables were located at the cell center, and the unknown fluxes were located at cell faces. For the current and Lorentz force calculation, the consistent scheme [32] was employed to calculate the current flux at cell face, and the current fluxes obtained were divergence free, which means the obtained current fluxes on cell faces are conservative in a control volume. Then, the Lorentz force at cell center can be calculated base on the divergence free current flux using a conservative interpolation of current densities to the cell centers. For the incompressible Navier–Stokes equations, fourstep projection method [33] was used to calculate the velocity and pressure. This kind of projection method is a general second order projection method, and the SIMPLE-type methods, such as the standard SIMPLE method in and the SIMPLEC method in [34,35], can be acquired from this projection method [36,37]. Both the convection and diffusion term of Navier–Stokes equation were discretized with the central difference scheme in spatial discretization, and for temporal updating, explicit and Crank–Nicholson scheme were utilized respectively.
For the velocity at walls, no-slip boundary conditions were set. We have assumed a fully developed flow at outlet, and then, we can set Dirichlet boundary condition u = c(y,z)and Neumann boundary condition ∂u/∂n = 0, respectively at inlet and outlet of the duct. For electrical potential, the Neumann boundary condition ∂ϕ/∂n = 0 was set at all boundaries.
2.3. Code and solving procedure The code was written in Fortran 90. An input file is provided to define the input data, such as geometry size of calculating domain, physical properties of fluid and duct, the number of grid points and time step size. To begin with, the code reads the input data and initializes the field of velocity, pressure, current, Lorentz force and electrical potential. In a time level, the four-step projection method is conducted to calculate the flow field. The Navier–Stokes system is solved by using Approximation factorization (AF) technique, and alternative direction implicit (ADI) method is employed to solve the pressure Poisson’s equation. Then, we solve the electrical Poisson’s equation and use the consistent and conservative scheme to calculate current density and Lorentz force. The obtained Lorentz force is considered as source term for Navier–Stokes system in next time level. In order to avoid the oscillations as well as to accelerate the convergence in the solution, a relaxation technique [38] is used for solving the pressure and electrical potential Poisson’s equations.
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741
1739
Fig. 3. Results for Hunt’s fully developed flow at Ha = 15,000: (a) velocity profile, (b) velocity distribution near the side walls, and (c) current distribution.
3. Validating cases Shercliff’s case and Hunt’s case are the exact solutions for fully developed incompressible laminar flows in rectangular ducts with transverse magnetic fields. These two cases were usually used to test numerical code for MHD duct flows. In Shercliff’s problem, all of the duct walls are insulating. In Hunt’s problem, the Hartmann walls perpendicular to the field are conducting and the side walls parallel to the field are insulating. We also used these two cases to validate the MTC-H 1.0 code. All the calculations were performed using a PC (Intel core2 Quad, 2.66 GHz, 4 GB RAM) with double precision. 3.1. 2D fully developed flow 2D fully developed flow of Shercliff’s and Hunt’s case were simulated to validate the numerical scheme of the code. The calculations were conducted on a non-uniform grid system at Ha = 500, Re = 10,
and Ha = 1000, Re = 10 with 65 × 65 nodes for resolving fluid region and 9 nodes for resolving duct walls. The dimension of fluid domain is 2 × 2, and the wall thickness is 0.1. The wall conductance ratio is 0.0 for Shercliff’s case and 0.05 for Hunt’s case. For Shercliff’s case, the pressure gradient −52.08 and −102.87 were given at mass flow rate 4 obtained form Shercliff’s analysis results at Ha = 500 and Ha = 1000, respectively, and for Hunt’s case, the pressure gradient were −953.48 and −3379.82, respectively. The results were compared with the exact solutions. Table 1 shows the computed flow rate for the case mentioned above, the error of the calculated results are less than 1%. Figs. 1(a) and 2(a) shows the convergence of mass flow rate. Induced current in the quarter duct at Ha = 1000 are shown in Figs. 1(b) and 2(b), The current density in the duct traces a closed streamline, which means the current density is conservative. Figs. 1(c) and (d) and 2(c) and (d) show the computed velocity distribution compared with the exact solution at Ha = 1000 for Shercliff’s and Hunt’s cases, respectively, we can see the computed results match well with exact solutions.
Table 2 The results of flow rate and pressure gradient. Shercliff’s case Hartmann number Calculated flow rate at outlet The given flow rate at inlet Error of flow rate Calculated pressure gradient Analytical pressure gradient Error of pressure gradient
1000 4.0069 4 0.17% 2.0582 2.0578 0.01%
Hunt’s case 2000 4.0135 4 0.3% 4.0956 4.0800 0.38%
1000 4.0087 4 0.22% 67.6525 67.5424 0.16%
2000 4.0082 4 0.21% 236.395 236.089 0.13%
1740
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741
3.3. 3D MHD flow To investigate the code in 3D simulation, the calculations of Shercliff’s and Hunt’s case were conducted at Ha = 1000, Re = 500 and Ha = 2000, Re = 500. We choose here the same cross section as above and extend it in the x-direction to 50 units, and the length in x direction is 5. A non-uniform mesh containing 51 × 55 × 55 nodes and 145,800 cells was used in this simulation. The wall conductance ratio is 0.05. Table 2 illustrates the results of flow rate and pressure gradient, and the relative errors are also given in this table. The results match well with the analytical solution. 3.4. 3D MHD flow at high Hartmann number A 3D MHD flow case was tested at high Hartmann number Ha = 10,000. The non-uniform mesh containing 62 × 65 × 65 nodes and 249,865 cells was used in this simulation, and wall conductance ratio is 0.05. The dimension of fluid domain is 5 × 2 × 2, and the wall thickness is 0.1. A parabolic velocity profile was given as velocity boundary condition at inlet, and the flow rate was 4. Fig. 4(a) shows the velocity profile in the duct at x = 0.05, x = 2.0, x = 3.95. Fig. 4(b) shows the comparison of velocity distribution between computed result and exact solution at the fully developed region near the side wall. We can see the numerical result match well with Hunt’s analytical solution. Fig. 4(c) shows the current distribution in the duct at x = 0.05, x = 2.0, x = 3.95. Since the flow is not fully developed, a 3D current occurs at x = 0.1. The calculated pressure gradient at the fully developed region is 19450.62, this results is very close to the Hunt’s analytical result of 19688.50. And the flow rate at outlet is 4.0092. 4. Conclusions
Fig. 4. Results for Hunt’s 3D case at Ha = 10,000: (a) velocity profiles at x = 0.05, x = 2.0, x = 3.95, (b) velocity distribution near the side walls, and (c) current distribution at x = 0.05, x = 2.0, x = 3.95.
In order to research and evaluate MHD effect in liquid metal blanket as well as make an optimal design for FDS series blankets, a code named MTC-H 1.0 which can simulate 3D MHD flow in rectangular ducts has been developed by FDS Team. The code was developed based on the structured grid. The current density conservative scheme was employed for calculation of the induced current and the Lorentz force, and the consistent projection method was employed for solving the incompressible Navier–Stokes equations with the Lorentz force included as a source term. We validated MTC-H 1.0 code by simulating MHD rectangular ducts flows, and the Shercliff’s insulated walls and Hunt’s conductive walls were used as benchmark models. The results showed that accuracy results can be obtained for duct flow at Ha = 104 . In the next step, in order to improve the capability for dealing with the complex geometry problem, the unstructured grid system will be needed, and the ability of parallel computing will be developed to promote efficiency of calculation. Acknowledgements
3.2. A 2D case at fusion relevant Hartmann number Since the Hartmann number at fusion reactor conditions is very large. We validated the code by Hunt’s case at a fusion relevant Hartmann number Ha = 15,000, Re = 10. A non-uniform mesh with 75 × 75 nodes was used for this calculation. The dimension of fluid is 2 × 2 and wall thickness is 0.1. The wall conductance ratio is 0.05. The pressure gradient −390809.68 were given at mass flow rate 4. Fig. 3(a) shows the velocity distribution, a great jet occurs near the side walls, and Fig. 3(b) illustrate that the computed velocity matches well with the analytical solution. Fig. 3(c) shows the conservative current distribution in the duct.
This work was supported partly by the National Natural Science Foundation of China with the grants No. 10675123, 10875145, 50936006, 50706015, 50676108. References [1] T. Ihli, T.K. Basu, L.M. Giancarli, S. Konishi, S. Malang, F. Najmabadi, et al., Review of blanket designs for advanced fusion reactors, Fusion Engineering and Design 83 (7–9) (2008) 912–919. [2] N.B. Morley, Y. Katoh, S. Malang, B.A. Pint, A.R. Raffray, S. Sharafat, et al., Recent research and development for the dual-coolant blanket concept in the US, Fusion Engineering and Design 83 (7–9) (2008) 920–927. [3] I.R. Kirillov, R.D. Team, Lithium cooled blanket of RF DEMO reactor, Fusion Engineering and Design 49 (2000) 457–465.
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741 [4] P. Sardain, D. Maisonnier, L. Di Pace, L. Giancarli, A.L. Puma, P. Norajitra, et al., The European power plant conceptual study: helium-cooled lithium-lead reactor concept, Fusion Engineering and Design 81 (23/24) (2006) 2673–2678. [5] E.R. Kumar, C. Danani, I. Sandeep, C. Chakrapani, N.R. Pragash, V. Chaudhari, et al., Preliminary design of Indian Test Blanket Module for ITER, Fusion Engineering and Design 83 (7–9) (2008) 1169–1172. [6] Y. Wu, FDS Team, Conceptual design activities of FDS series fusion power plants in China, Fusion Engineering and Design 81 (23/24) (2006) 2713–2718. [7] Y.C. Wu, J.P. Qian, J.N. Yu, The fusion-driven hybrid system and its material selection, Journal of Nuclear Materials 307–311 (Part 2) (2002) 1629–1636. [8] Y. Wu, FDS Team, Conceptual design of the China fusion power plant FDS-II, Fusion Engineering and Design 83 (10–12) (2008) 1683–1689. [9] Y. Wu, FDS Team, Fusion-based hydrogen production reactor and its material selection, Journal of Nuclear Materials 386–388 (2009) 122–126. [10] Y. Wu, L. Qiu, Y. Chen, Conceptual study on liquid metal center conductor post in spherical tokamak reactors, Fusion Engineering and Design 51–52 (2000) 395–399. [11] Y. Wu, A fusion neutron source driven sub-critical nuclear energy system: a way for early application of fusion technology, Plasma Science and Technology 6 (2001) 1085. [12] Y. Wu, FDS Team, Design status and development strategy of China liquid lithium–lead blankets and related material technology, Journal of Nuclear Materials 367–370 (Part 2) (2007) 1410–1415. [13] H. Chen, Y. Wu, S. Konishi, J. Hayward, A high temperature blanket concept for hydrogen production, Fusion Engineering and Design 83 (7–9) (2008) 903–911. [14] S. Malang, P. Leroy, G.P. Casini, R.F. Mattas, Y. Strebkov, Crucial issues on liquidmetal blanket design, Fusion Engineering and Design 16 (1991) 95–109. [15] L. Barleon, V. Casal, L. Lenhart, Mhd flow in liquid-metal-cooled blankets, Fusion Engineering and Design 14 (3/4) (1991) 401–412. [16] I.R. Kirillov, C.B. Reed, L. Barleon, K. Miyazaki, Present understanding of Mhd and heat-transfer phenomena for liquid-metal blankets, Fusion Engineering and Design 27 (1995) 553–569. [17] N.B. Morley, S. Smolentsev, L. Barleon, I.R. Kirillov, M. Takahashi, Liquid magnetohydrodynamics—recent progress and future directions for fusion, Fusion Engineering and Design 51-2 (2000) 701–713. [18] N.B. Morley, S. Malang, I. Kirillov, Thermofluid magnetohydrodynamic issues for liquid breeders, Fusion Science and Technology 47 (3) (2005) 488–501. [19] J.A. Shercliff, Steady motion of conducting fluids in pipes under transverse magnetic fields, Proceedings of the Cambridge Philosophical Society 49 (1) (1953) 136–144. [20] J.A. Shercliff, The flow of conducting fluids in circular pipes under transverse magnetic fields, Journal of Fluid Mechanics 1 (6) (1956) 644–666. [21] C.C. Chang, T.S. Lundgren, Duct flow in magnetohydrodynamics, Zeitschrift Fur Angewandte Mathematik Und Physik 12 (2) (1961) 100.
1741
[22] J.C.R. Hunt, Magnetohydrodynamic flow in rectangular ducts, Journal of Fluid Mechanics 21 (1965) 577. [23] A. Sterl, Numerical-simulation of liquid-metal Mhd flows in rectangular ducts, Journal of Fluid Mechanics 216 (1990) 161–191. [24] N. Umeda, M. Takahashi, Numerical analysis for heat transfer enhancement of a lithium flow under a transverse magnetic field, Fusion Engineering and Design 51-2 (2000) 899–907. [25] 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 (1) (2005) 83–93. [26] M.J. Ni, R. Munipalli, N.B. Morley, M.A. Abdou, Validation strategies of HIMAG in interfacial flow computation for fusion applications, Fusion Engineering and Design 81 (8–14) (2006) 1535–1541. [27] C. Mistrangelo, L. Buhler, Numerical investigation of liquid metal flows in rectangular sudden expansions, Fusion Engineering and Design 82 (15–24) (2007) 2176–2182. [28] M.J. Pattison, K.N. Premnath, N.B. Morley, M.A. Abdou, Progress in lattice Boltzmann methods for magnetohydrodynamic flows relevant to fusion applications, Fusion Engineering and Design 83 (4) (2008) 557–572. [29] R.J. Moreau, Magnetohydrodynamics, Kluwer Academic Publishers, 1990. [30] U. Müller, L. Bühler, Magnetofluiddynamics in Channels and Containers, Springer, 2001. [31] P.A. Davidson, An Introduction to Magnetohydrodynamics, Cambridge University Press, 2001. [32] M.-J. Ni, R. Munipalli, N.B. Morley, P. Huang, M.A. Abdou, A current density conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. Part I: on a rectangular collocated grid system, Journal of Computational Physics 227 (1) (2007) 174–204. [33] M.J. Ni, S. Komori, N. Morley, Projection methods for the calculation of incompressible unsteady flows, Numerical Heat Transfer Part B-Fundamentals 44 (6) (2003) 533–551. [34] S.V. Patankar, D.B. Spalding, Calculation procedure for heat, mass and momentum-transfer in 3-dimensional parabolic flows, International Journal of Heat and Mass Transfer 15 (10) (1972) 1787. [35] J.P. Vandoormaal, G.D. Raithby, Enhancements of the simple method for predicting incompressible fluid-flows, Numerical Heat Transfer 7 (2) (1984) 147–163. [36] M.J. Ni, Projection and simple methods for single-fluid and multiple-fluid incompressible unsteady flows—general formula, Numerical Heat Transfer Part B-Fundamentals 50 (5) (2006) 395–408. [37] M.J. Ni, M.A. Abdou, A bridge between projection methods and SIMPLE type methods for incompressible Navier–Stokes equations, International Journal for Numerical Methods in Engineering 72 (12) (2007) 1490–1512. [38] W.Q. Tao, Numerical Heat Transfer, Xi’an Jiaotong Univercity Press, 2001.
Contents lists available at ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Code development and validation for analyzing liquid metal MHD flow in rectangular ducts Tao Zhou a,c,∗ , Zhiyi Yang a , Mingjiu Ni c , Hongli Chen a,b a b c
Institute of Plasma Physics, Chinese Academy of Sciences, Hefei, Anhui 230031, China School of Nuclear Science and Technology, University of Science and Technology of China, Hefei, Anhui 230027, China Graduate University of Chinese Academy of Sciences, Beijing 100049, China
a r t i c l e
i n f o
Keywords: Numerical simulations Liquid metal magnetohydrodynamics Current density conservative scheme
a b s t r a c t A code named MTC-H 1.0 which can simulate 3D magnetohydrodynamics (MHD) flow in rectangular ducts has been developed by FDS Team. In this code, a conservative scheme of the current density was employed for calculation of the induced current and the Lorentz force, and the consistent projection method was employed for solving the incompressible Navier–Stokes equations with the Lorentz force included as a source term. The code was developed on a structured collocated grid, on which velocity, pressure, and electrical potential were located in the cell center, while current fluxes were located on the cell faces. The code MTC-H 1.0 was test by simulating MHD rectangular ducts flows, and the Shercliff’s insulated walls and Hunt’s conductive walls were used as benchmark models. The validation cases were conducted with Hartmann number from 102 to 104 on structured collocated grids. The results matched well with Hunt’s and Shercliff’s analytical solutions and it showed good accuracy of the code. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Liquid metal breeder blanket concept has been a topic of great interest in fusion reactor blanket design [1–5] because of many advantages. A series of fusion reactors (named FDS series) have been designed and assessed by the FDS Team in China [6–10]. And a series of liquid metal blanket concepts [11–13] have been developed correspondingly. In these blankets, lead–lithium alloy is considered as tritium breeder only or as breeder and coolant. But the motion of liquid metal in fusion reactor strong magnetic field cause serious magnetohydrodynamic (MHD) effects, which have dramatic impacts on velocity distribution, heat transfer characteristics, pressure drop and the required pumping power for the cooling system. Therefore, the liquid metal MHD effect under fusion relevant condition is one of the key issues in making an optimal design for LM blanket [14–18]. Up to now, liquid metal MHD duct flow has been extensively studied by theory [19–22] and numerical simulation [23–28]. Due to the complex geometry of flow channel in blankets, it is difficult to study the detail of MHD flow characteristics by theoretical analysis. As the development of computer science, the numerical
∗ Corresponding author at: Institute of Plasma Physics, Chinese Academy of Sciences, Hefei, Anhui 230031, China. Tel.: +86 551 5593296; fax: +86 551 5591397. E-mail address: [email protected] (T. Zhou). 0920-3796/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.fusengdes.2010.05.034
simulation has become an effective tool for analyzing MHD effects and optimizing blanket design. For MHD duct flow, the Hartmann layer perpendicular to the magnetic field scales with Ha−1 , and the side layer parallel to the magnetic field scales with Ha−1/2 . In fusion blankets, Hartmann number can be as high as 104 –105 , which leads to some special numerical difficulty. Therefore, simulating and studying liquid metal MHD flow at high Hartmann number has become a hot research topic for fusion blanket design. A number of numerical codes have been developed for analyzing liquid metal MHD effects [23–28]. Smolentsev developed a 2D code based on structured grid and induced magnetic equation [25]. A 3D code named HIMAG [26] and another 3D code based on the commercial package CFX 5.6 [27] were developed based on unstructured grid and electrical potential equation. By solving induced magnetic equations, it is difficult to construct a proper boundary condition for induced magnetic field. For fusion reactor conditions, the magnetic Reynolds numbers is small, the electrical potential formula can be employed for MHD with good accuracy [29–31]. Ni developed a consistent and conservative scheme [32] for calculate the current density and the Lorentz force by solving electrical potential equations, this scheme can solve LM MHD problem at high Hartmann number. In order to research and evaluate MHD effect in liquid metal blanket as well as make an optimal design for FDS series blankets, FDS team has been developing the MHD simulation code named MTC (Magnetic Thermo-hydraulics Coupling Code) for years. MTC-F 1.0 and MTC-F 2.0 are developed
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741
based on the commercial code FLUENT, which can solve LM MHD problem only at low Hartman number (Ha < 500). In this paper, a numerical code named MTC-H 1.0 for simulating MHD flow in rectangular duct is developed based on the structure collocated grid, in which consistent and conservative scheme are used to calculate the current density and the Lorentz force. We validated the numerical code by simulating 2D and 3D MHD flows with analytical solutions existed. The validation cases are conducted with Hartmann number from 102 to 104 .
1737
Table 1 Computed flow rate. Hartmann number
Calculated flow rate
Error
Shercliff’s case
500 1000
3.979 3.964
0.5% 0.9%
Hunt’s case
500 1000
4.008 4.016
0.2% 0.4%
Ohm’s law: J = −∇ ϕ + u × B
2. Governing equations and numerical methods
Conservation of current density:
2.1. Governing equations LM is a kind of incompressible fluid. Under fusion blanket conditions, the magnetic Reynolds numbers is small, this means the induced magnetic field can be negligible compared to the apply magnetic field, the governing equations for LM MHD flow in dimensionless form can be expressed as follows: Conservation of momentum: ∂u 1 2 + u · ∇ u = −∇ p + ∇ u + NJ × B Re ∂t
(1)
Conservation of mass (Navier–Stokes equation):
∇ ·u=0
(3)
(2)
∇ ·J =0
(4)
From Eqs. (3) and (4),we can get the Poisson’s equation for electrical potential:
∇ · (∇ ϕ) = ∇ · (u × B)
(5)
Here t, u, p, J, B, ϕ are dimensionless time, velocity, pressure, current density, applied magnetic field, and electrical potential. They are scaled with L/u0 , u0 , 0 u20 , u0 B0 , B0 and Lu0 B0 , respectively, where L is a characteristic length, u0 is characteristic velocity, and 0 are the conductivity and density of fluid. The Reynolds number is Re = u0 L/, the Hartmann number is Ha = B0 L /0 , and the interaction parameter is N = Ha2 /Re.
Fig. 1. Calculating results for Shercliff’s fully developed flow at Ha = 1000: (a) convergence of flow rate, (b) current in the duct, (c) velocity distribution along the middle line normal to the Hartmann walls, and (d) velocity distribution along the middle line normal to the side walls.
1738
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741
Fig. 2. Calculating results for Hunt’s fully developed flow at Ha = 1000: (a) convergence of flow rate, (b) current in the duct, (c) velocity distribution along the middle line normal to the Hartmann walls, and (d) velocity distribution along the middle line normal to the side walls.
2.2. Numerical scheme and boundary conditions The code was developed based on the structured grid. Since the Hartmann layer and the side layer are very thin, and require at least 4–5 grid points to resolve these layers, a non-uniform grid system is needed to improve the computational efficiency. The control volume method on the structured collocated grid system was used to ensure the conservation of the discrete equations. All the unknown variables were located at the cell center, and the unknown fluxes were located at cell faces. For the current and Lorentz force calculation, the consistent scheme [32] was employed to calculate the current flux at cell face, and the current fluxes obtained were divergence free, which means the obtained current fluxes on cell faces are conservative in a control volume. Then, the Lorentz force at cell center can be calculated base on the divergence free current flux using a conservative interpolation of current densities to the cell centers. For the incompressible Navier–Stokes equations, fourstep projection method [33] was used to calculate the velocity and pressure. This kind of projection method is a general second order projection method, and the SIMPLE-type methods, such as the standard SIMPLE method in and the SIMPLEC method in [34,35], can be acquired from this projection method [36,37]. Both the convection and diffusion term of Navier–Stokes equation were discretized with the central difference scheme in spatial discretization, and for temporal updating, explicit and Crank–Nicholson scheme were utilized respectively.
For the velocity at walls, no-slip boundary conditions were set. We have assumed a fully developed flow at outlet, and then, we can set Dirichlet boundary condition u = c(y,z)and Neumann boundary condition ∂u/∂n = 0, respectively at inlet and outlet of the duct. For electrical potential, the Neumann boundary condition ∂ϕ/∂n = 0 was set at all boundaries.
2.3. Code and solving procedure The code was written in Fortran 90. An input file is provided to define the input data, such as geometry size of calculating domain, physical properties of fluid and duct, the number of grid points and time step size. To begin with, the code reads the input data and initializes the field of velocity, pressure, current, Lorentz force and electrical potential. In a time level, the four-step projection method is conducted to calculate the flow field. The Navier–Stokes system is solved by using Approximation factorization (AF) technique, and alternative direction implicit (ADI) method is employed to solve the pressure Poisson’s equation. Then, we solve the electrical Poisson’s equation and use the consistent and conservative scheme to calculate current density and Lorentz force. The obtained Lorentz force is considered as source term for Navier–Stokes system in next time level. In order to avoid the oscillations as well as to accelerate the convergence in the solution, a relaxation technique [38] is used for solving the pressure and electrical potential Poisson’s equations.
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741
1739
Fig. 3. Results for Hunt’s fully developed flow at Ha = 15,000: (a) velocity profile, (b) velocity distribution near the side walls, and (c) current distribution.
3. Validating cases Shercliff’s case and Hunt’s case are the exact solutions for fully developed incompressible laminar flows in rectangular ducts with transverse magnetic fields. These two cases were usually used to test numerical code for MHD duct flows. In Shercliff’s problem, all of the duct walls are insulating. In Hunt’s problem, the Hartmann walls perpendicular to the field are conducting and the side walls parallel to the field are insulating. We also used these two cases to validate the MTC-H 1.0 code. All the calculations were performed using a PC (Intel core2 Quad, 2.66 GHz, 4 GB RAM) with double precision. 3.1. 2D fully developed flow 2D fully developed flow of Shercliff’s and Hunt’s case were simulated to validate the numerical scheme of the code. The calculations were conducted on a non-uniform grid system at Ha = 500, Re = 10,
and Ha = 1000, Re = 10 with 65 × 65 nodes for resolving fluid region and 9 nodes for resolving duct walls. The dimension of fluid domain is 2 × 2, and the wall thickness is 0.1. The wall conductance ratio is 0.0 for Shercliff’s case and 0.05 for Hunt’s case. For Shercliff’s case, the pressure gradient −52.08 and −102.87 were given at mass flow rate 4 obtained form Shercliff’s analysis results at Ha = 500 and Ha = 1000, respectively, and for Hunt’s case, the pressure gradient were −953.48 and −3379.82, respectively. The results were compared with the exact solutions. Table 1 shows the computed flow rate for the case mentioned above, the error of the calculated results are less than 1%. Figs. 1(a) and 2(a) shows the convergence of mass flow rate. Induced current in the quarter duct at Ha = 1000 are shown in Figs. 1(b) and 2(b), The current density in the duct traces a closed streamline, which means the current density is conservative. Figs. 1(c) and (d) and 2(c) and (d) show the computed velocity distribution compared with the exact solution at Ha = 1000 for Shercliff’s and Hunt’s cases, respectively, we can see the computed results match well with exact solutions.
Table 2 The results of flow rate and pressure gradient. Shercliff’s case Hartmann number Calculated flow rate at outlet The given flow rate at inlet Error of flow rate Calculated pressure gradient Analytical pressure gradient Error of pressure gradient
1000 4.0069 4 0.17% 2.0582 2.0578 0.01%
Hunt’s case 2000 4.0135 4 0.3% 4.0956 4.0800 0.38%
1000 4.0087 4 0.22% 67.6525 67.5424 0.16%
2000 4.0082 4 0.21% 236.395 236.089 0.13%
1740
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741
3.3. 3D MHD flow To investigate the code in 3D simulation, the calculations of Shercliff’s and Hunt’s case were conducted at Ha = 1000, Re = 500 and Ha = 2000, Re = 500. We choose here the same cross section as above and extend it in the x-direction to 50 units, and the length in x direction is 5. A non-uniform mesh containing 51 × 55 × 55 nodes and 145,800 cells was used in this simulation. The wall conductance ratio is 0.05. Table 2 illustrates the results of flow rate and pressure gradient, and the relative errors are also given in this table. The results match well with the analytical solution. 3.4. 3D MHD flow at high Hartmann number A 3D MHD flow case was tested at high Hartmann number Ha = 10,000. The non-uniform mesh containing 62 × 65 × 65 nodes and 249,865 cells was used in this simulation, and wall conductance ratio is 0.05. The dimension of fluid domain is 5 × 2 × 2, and the wall thickness is 0.1. A parabolic velocity profile was given as velocity boundary condition at inlet, and the flow rate was 4. Fig. 4(a) shows the velocity profile in the duct at x = 0.05, x = 2.0, x = 3.95. Fig. 4(b) shows the comparison of velocity distribution between computed result and exact solution at the fully developed region near the side wall. We can see the numerical result match well with Hunt’s analytical solution. Fig. 4(c) shows the current distribution in the duct at x = 0.05, x = 2.0, x = 3.95. Since the flow is not fully developed, a 3D current occurs at x = 0.1. The calculated pressure gradient at the fully developed region is 19450.62, this results is very close to the Hunt’s analytical result of 19688.50. And the flow rate at outlet is 4.0092. 4. Conclusions
Fig. 4. Results for Hunt’s 3D case at Ha = 10,000: (a) velocity profiles at x = 0.05, x = 2.0, x = 3.95, (b) velocity distribution near the side walls, and (c) current distribution at x = 0.05, x = 2.0, x = 3.95.
In order to research and evaluate MHD effect in liquid metal blanket as well as make an optimal design for FDS series blankets, a code named MTC-H 1.0 which can simulate 3D MHD flow in rectangular ducts has been developed by FDS Team. The code was developed based on the structured grid. The current density conservative scheme was employed for calculation of the induced current and the Lorentz force, and the consistent projection method was employed for solving the incompressible Navier–Stokes equations with the Lorentz force included as a source term. We validated MTC-H 1.0 code by simulating MHD rectangular ducts flows, and the Shercliff’s insulated walls and Hunt’s conductive walls were used as benchmark models. The results showed that accuracy results can be obtained for duct flow at Ha = 104 . In the next step, in order to improve the capability for dealing with the complex geometry problem, the unstructured grid system will be needed, and the ability of parallel computing will be developed to promote efficiency of calculation. Acknowledgements
3.2. A 2D case at fusion relevant Hartmann number Since the Hartmann number at fusion reactor conditions is very large. We validated the code by Hunt’s case at a fusion relevant Hartmann number Ha = 15,000, Re = 10. A non-uniform mesh with 75 × 75 nodes was used for this calculation. The dimension of fluid is 2 × 2 and wall thickness is 0.1. The wall conductance ratio is 0.05. The pressure gradient −390809.68 were given at mass flow rate 4. Fig. 3(a) shows the velocity distribution, a great jet occurs near the side walls, and Fig. 3(b) illustrate that the computed velocity matches well with the analytical solution. Fig. 3(c) shows the conservative current distribution in the duct.
This work was supported partly by the National Natural Science Foundation of China with the grants No. 10675123, 10875145, 50936006, 50706015, 50676108. References [1] T. Ihli, T.K. Basu, L.M. Giancarli, S. Konishi, S. Malang, F. Najmabadi, et al., Review of blanket designs for advanced fusion reactors, Fusion Engineering and Design 83 (7–9) (2008) 912–919. [2] N.B. Morley, Y. Katoh, S. Malang, B.A. Pint, A.R. Raffray, S. Sharafat, et al., Recent research and development for the dual-coolant blanket concept in the US, Fusion Engineering and Design 83 (7–9) (2008) 920–927. [3] I.R. Kirillov, R.D. Team, Lithium cooled blanket of RF DEMO reactor, Fusion Engineering and Design 49 (2000) 457–465.
T. Zhou et al. / Fusion Engineering and Design 85 (2010) 1736–1741 [4] P. Sardain, D. Maisonnier, L. Di Pace, L. Giancarli, A.L. Puma, P. Norajitra, et al., The European power plant conceptual study: helium-cooled lithium-lead reactor concept, Fusion Engineering and Design 81 (23/24) (2006) 2673–2678. [5] E.R. Kumar, C. Danani, I. Sandeep, C. Chakrapani, N.R. Pragash, V. Chaudhari, et al., Preliminary design of Indian Test Blanket Module for ITER, Fusion Engineering and Design 83 (7–9) (2008) 1169–1172. [6] Y. Wu, FDS Team, Conceptual design activities of FDS series fusion power plants in China, Fusion Engineering and Design 81 (23/24) (2006) 2713–2718. [7] Y.C. Wu, J.P. Qian, J.N. Yu, The fusion-driven hybrid system and its material selection, Journal of Nuclear Materials 307–311 (Part 2) (2002) 1629–1636. [8] Y. Wu, FDS Team, Conceptual design of the China fusion power plant FDS-II, Fusion Engineering and Design 83 (10–12) (2008) 1683–1689. [9] Y. Wu, FDS Team, Fusion-based hydrogen production reactor and its material selection, Journal of Nuclear Materials 386–388 (2009) 122–126. [10] Y. Wu, L. Qiu, Y. Chen, Conceptual study on liquid metal center conductor post in spherical tokamak reactors, Fusion Engineering and Design 51–52 (2000) 395–399. [11] Y. Wu, A fusion neutron source driven sub-critical nuclear energy system: a way for early application of fusion technology, Plasma Science and Technology 6 (2001) 1085. [12] Y. Wu, FDS Team, Design status and development strategy of China liquid lithium–lead blankets and related material technology, Journal of Nuclear Materials 367–370 (Part 2) (2007) 1410–1415. [13] H. Chen, Y. Wu, S. Konishi, J. Hayward, A high temperature blanket concept for hydrogen production, Fusion Engineering and Design 83 (7–9) (2008) 903–911. [14] S. Malang, P. Leroy, G.P. Casini, R.F. Mattas, Y. Strebkov, Crucial issues on liquidmetal blanket design, Fusion Engineering and Design 16 (1991) 95–109. [15] L. Barleon, V. Casal, L. Lenhart, Mhd flow in liquid-metal-cooled blankets, Fusion Engineering and Design 14 (3/4) (1991) 401–412. [16] I.R. Kirillov, C.B. Reed, L. Barleon, K. Miyazaki, Present understanding of Mhd and heat-transfer phenomena for liquid-metal blankets, Fusion Engineering and Design 27 (1995) 553–569. [17] N.B. Morley, S. Smolentsev, L. Barleon, I.R. Kirillov, M. Takahashi, Liquid magnetohydrodynamics—recent progress and future directions for fusion, Fusion Engineering and Design 51-2 (2000) 701–713. [18] N.B. Morley, S. Malang, I. Kirillov, Thermofluid magnetohydrodynamic issues for liquid breeders, Fusion Science and Technology 47 (3) (2005) 488–501. [19] J.A. Shercliff, Steady motion of conducting fluids in pipes under transverse magnetic fields, Proceedings of the Cambridge Philosophical Society 49 (1) (1953) 136–144. [20] J.A. Shercliff, The flow of conducting fluids in circular pipes under transverse magnetic fields, Journal of Fluid Mechanics 1 (6) (1956) 644–666. [21] C.C. Chang, T.S. Lundgren, Duct flow in magnetohydrodynamics, Zeitschrift Fur Angewandte Mathematik Und Physik 12 (2) (1961) 100.
1741
[22] J.C.R. Hunt, Magnetohydrodynamic flow in rectangular ducts, Journal of Fluid Mechanics 21 (1965) 577. [23] A. Sterl, Numerical-simulation of liquid-metal Mhd flows in rectangular ducts, Journal of Fluid Mechanics 216 (1990) 161–191. [24] N. Umeda, M. Takahashi, Numerical analysis for heat transfer enhancement of a lithium flow under a transverse magnetic field, Fusion Engineering and Design 51-2 (2000) 899–907. [25] 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 (1) (2005) 83–93. [26] M.J. Ni, R. Munipalli, N.B. Morley, M.A. Abdou, Validation strategies of HIMAG in interfacial flow computation for fusion applications, Fusion Engineering and Design 81 (8–14) (2006) 1535–1541. [27] C. Mistrangelo, L. Buhler, Numerical investigation of liquid metal flows in rectangular sudden expansions, Fusion Engineering and Design 82 (15–24) (2007) 2176–2182. [28] M.J. Pattison, K.N. Premnath, N.B. Morley, M.A. Abdou, Progress in lattice Boltzmann methods for magnetohydrodynamic flows relevant to fusion applications, Fusion Engineering and Design 83 (4) (2008) 557–572. [29] R.J. Moreau, Magnetohydrodynamics, Kluwer Academic Publishers, 1990. [30] U. Müller, L. Bühler, Magnetofluiddynamics in Channels and Containers, Springer, 2001. [31] P.A. Davidson, An Introduction to Magnetohydrodynamics, Cambridge University Press, 2001. [32] M.-J. Ni, R. Munipalli, N.B. Morley, P. Huang, M.A. Abdou, A current density conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. Part I: on a rectangular collocated grid system, Journal of Computational Physics 227 (1) (2007) 174–204. [33] M.J. Ni, S. Komori, N. Morley, Projection methods for the calculation of incompressible unsteady flows, Numerical Heat Transfer Part B-Fundamentals 44 (6) (2003) 533–551. [34] S.V. Patankar, D.B. Spalding, Calculation procedure for heat, mass and momentum-transfer in 3-dimensional parabolic flows, International Journal of Heat and Mass Transfer 15 (10) (1972) 1787. [35] J.P. Vandoormaal, G.D. Raithby, Enhancements of the simple method for predicting incompressible fluid-flows, Numerical Heat Transfer 7 (2) (1984) 147–163. [36] M.J. Ni, Projection and simple methods for single-fluid and multiple-fluid incompressible unsteady flows—general formula, Numerical Heat Transfer Part B-Fundamentals 50 (5) (2006) 395–408. [37] M.J. Ni, M.A. Abdou, A bridge between projection methods and SIMPLE type methods for incompressible Navier–Stokes equations, International Journal for Numerical Methods in Engineering 72 (12) (2007) 1490–1512. [38] W.Q. Tao, Numerical Heat Transfer, Xi’an Jiaotong Univercity Press, 2001.