- Email: [email protected]
Development of a simulation method for evaluating Marangoni convection with free surface for tungsten divertor
Fusion Engineering and Design xxx (xxxx) xxx–xxx
Contents lists available at ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Development of a simulation method for evaluating Marangoni convection with free surface for tungsten divertor ⁎
Kohei Hamaguchi , Yu Teramoto, Eiji Hoashi, Takafumi Okita, Kenzo Ibano, Yoshio Ueda Graduate School of Engineering, Osaka University, Osaka, Japan
A R T I C L E I N F O
A B S T R A C T
Keywords: Tungsten divertor Molten tungsten Marangoni convection Free surface Heat flux
A tungsten (W) divertor has been developed for tokamak-type nuclear fusion reactors, such as ITER. Because the divertor will be repetitively subjected to high heat flux when plasma disruption or edge localized mode occurs, W at the surface of the divertor will be melted and solidified during reactor operation. The surface shape of the W divertor after or during such a heat load strongly influences its lifetime. Thus, evaluating the behavior of W in such conditions is important. An important driving force in molten metal is Marangoni convection due to a difference in surface tension resulting from a temperature gradient at the surface. To clarify the mechanism of Marangoni convection in molten W, a simulation method that can simultaneously treat surface deformation and convection in molten W, including Marangoni convection, is required. We have been developing such a simulation code based on the constrained interpolation profile combined unified procedure method. In this paper, we simulated molten W under a given temperature gradient or heat load and validated the model by simulating the flow in a molten W pool with a temperature gradient. The development process of Marangoni convection in molten W with a free surface was evaluated.
1. Introduction Tungsten (W) has been studied for the candidate material of plasma facing materials, such as divertor, in tokamak-type nuclear fusion reactor. W has the highest melting point among all metals, a high thermal conductivity, and a low sputtering rate, thus making it a promising divertor surface material. In tokamak-type nuclear fusion reactors, such as ITER, because the divertor will be subjected to a high heat flux of at least 100 MW/m2 within several milliseconds or less when plasma disruption or type-I edge localized mode (ELM) occurs, the W in the divertor will undergo repeated melting and solidification [1]. Thus, the surface shape of the W divertor may change after solidification. If the surface shape would be smooth, it is expected the W will stand up to subsequent repetitive heat load. Because surface roughness can decrease the thermal resistance, the surface status after heat load influences on the lifetime of the W divertor. Because controlling plasma disruption and ELM is very difficult, knowledge of the surface status after solidification is required to ensure the lifetime. To acquire such knowledge, it is necessary to investigate and understand the details of melting and solidification process of W subjected to high heat flux. Previous research has indicated that molten state flow is one of the main factors that determine the surface shape [1]. One primary factor that controls flow in a molten metal is the
⁎
Marangoni effect, which is caused by a difference in the surface tension resulting from a temperature gradient at the surface. Convection due to the Marangoni effect is called Marangoni convection and has been previously investigated [2–4]. A liquid layer with a uniform surface tension maintains a stationary state. However, when the temperature or concentration at the liquid surface is not uniform, movement of the liquid surface occurs and leads to Marangoni convection. In the analysis of Marangoni convection, a dimensionless number, which is called the Marangoni number (Ma), represents the strength and instability of the convection; it is expressed as follows:
Ma = −
(∂σ / ∂T )ΔT d, μα
(1)
where T, ΔΤ, μ, α, and d are temperature, characteristic temperature difference, viscosity coefficient, thermal diffusion coefficient, and the characteristic length, respectively. σ represents the surface tension coefficient depending on the temperature. The characteristic length d corresponds to the depth of the liquid layer. Nield has reported that Marangoni convection becomes more dominant than buoyancy in a thin liquid layer, such as a layer that is a few millimeters thick [5]. Extensive research on Marangoni convection for fluids such as water, which have relatively high Prandtl numbers (Pr) than molten W, has been reported [6]. For example, the Pr of water and molten W are 7
Corresponding author. E-mail address: [email protected] (K. Hamaguchi).
https://doi.org/10.1016/j.fusengdes.2018.02.007 Received 29 September 2017; Received in revised form 2 February 2018; Accepted 2 February 2018 0920-3796/ © 2018 Elsevier B.V. All rights reserved.
Please cite this article as: Hamaguchi, K., Fusion Engineering and Design (2018), https://doi.org/10.1016/j.fusengdes.2018.02.007
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
component Fsa(t) is expressed as follows:
and 0.015, respectively. However, the literature contains few reports on the process of Marangoni convection in low-Pr fluids. Li et al. evaluated the relevance of Marangoni convection and buoyancy in a circular pool using molten silicon (Pr = 0.011) as a low-Pr fluid [7]. Smith and Davis indicated the instabilities of flows for high- and low-Pr fluids [8] in twodimensional (2D) simulations. In addition, the influence of impurities on Marangoni convection flow patterns in a molten layer was reported by Tsotridis et al. [1]. However, research on Marangoni convection using the molten W has not been reported. Ueda and Inoue et al. have conducted experiments on W irradiated by a laser with a duration time from microseconds to milliseconds [9,10]. According to their observation of the W surface after the laser irradiation, the center of the irradiation area was found to have swollen like a cone. However, the formation mechanism of the cone by heat flux was not fully understood. W has a high melting point and must be maintained in an ultra high purity state to prevent surface oxidation; it is thus, are difficult to perform experiments with the molten W. In addition, measuring inner flow is difficult because the molten W is metal with an opaque surface. Thus, we have been aiming to develop a three-dimensional (3D) unified simulation method that can evaluate the melting and solidification process of W, including phase change. The surface of the molten W divertor becomes a free surface; however, free surfaces have not been considered in the research conducted by researches such as Smith and Davis [8], Tsotridis et al. [1]. Thus, our simulation code newly considers free surface. In our simulation, we use the properties of molten W reported by Ishikawa et al., who conducted measurements using an electrostatic levitator [11]. In this study, we develop the simulation code that can reproduce the Marangoni effect in surface tension model for future evaluation of melting and solidification process of W, and then we clarify the development process of Marangoni convection of W with the free surface using this code.
Fsa(t) = ∇σ.
(3)
Temperature dependence of surface tension coefficient σ is defined as following empirical equation obtained from experimental data [11]:
σ (T ) = 2.48 − 0.31 × (T − Tm ) × 10−3
(4)
where Tm is the melting point. 2.2. Governing equations As the governing equations, the equation of continuity (5), momentum equation (6), equation of energy (7), and EOS (8) are used in our simulations as follows:
Dρ ∂u = −ρ i , Dt ∂x i
(5)
Dui 1 ∂P 1 ∂Sij + gi , =− + Dt ρ ∂x i ρ ∂x i
(6)
Sij ∂ui DT 1 ∂ ⎛ ∂T ⎞ Q P ∂u = − TH i + + κ + , and Dt ρCV ∂x i ρCV ∂x i ρCV ∂x i ⎝ ∂x i ⎠ ρCV
(7)
P = P(ρ, T)
(8)
⎜
⎟
where ρ, P, S, Cv, κ, and Q indicate the density, pressure, stress or viscous stress, specific heat at constant volume, thermal conductivity, and the amount of volumetric heat generated, respectively. Parameter PTH is calculated according to the expression PTH = T(∂P/∂T). The GRAY-EOS developed at Lawrence Livermore National Laboratory [17] is used as EOS.
2. Material and methods
3. Results and discussion
2.1. Simulation method
3.1. Comparison of our results with those of previous studies
The constrained interpolation profile (CIP) method, which enables the accurate solving of the advective term, was proposed by Yabe et al. [12]. In addition, the CIP combined unified procedure (C-CUP) method based on the CIP method has been developed. The C-CUP method is suitable for managing both compressible and incompressible fluids [13]. However, temperature is not treated in the C-CUP method. Hoashi et al. developed a method that can treat temperature with the C-CUP method by determining it using the equation of energy (EOS) [14]. In the analysis of the flow in which multiple phases coexist, Yabe et al. used a density function to capture the interface [12]. We considered a method that can solve the density function with tangent conversion to decrease numerical diffusion to capture the interface accurately [15]. Surface tension should be considered because it is one of the primary factors, such as Marangoni convection, that affect inner flow. Thus, we constructed the surface tension model based on the continuous surface force (CSF) model suggested by Brackbill et al. [16]. The total surface forces Fsa is expressed as follows:
Because the heat flux that W undergoes at the surface of the divertor will vary for different locations on the divertor, a temperature gradient parallel to the surface will occur. However, the mechanism by which Marangoni convection develops in molten W with a free surface has not been clarified. Therefore, we considered the 3D simulation model shown in Fig. 1. As an initial condition, molten W was placed in the lower region and air was placed in the upper region. Air was assumed to behave as an ideal gas, and the initial velocity in all regions was set to zero. The plane at x = 0 was treated as the cold wall with a temperature of 4000 K, and the plane opposite the cold wall in direction x was treated as the hot wall with a temperature of 4100 K. Under these
Fsa = Fsa(n) + Fsa(t) Fsa(n)
(2) (t)
is the normal component to the interface and Fsa is the where tangential component. Fsa(n) is obtained by a curvature and surface tension coefficient according to CSF model. A curvature is obtained by a density function. Calculating the tangential component Fsa(t) is required so that Marangoni convection occurs due to the contribution of Fsa(t). Thus, it is necessary to specify the direction in which Fsa(t) works at the interface. We used the polar coordinate transformation for specifying the direction. Then, the cell, where the interface was adjacent, was determined. The difference in the surface tension between those cells becomes the tangential component Fsa(t). Therefore, the tangential
Fig. 1. Schematic of the simulation model.
2
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
in the x direction for liquid phase normalized by the velocity at the surface. The results of our simulation are consistent with those of Smith and Davis. Therefore, our simulation is sufficiently accurate to reproduce the inner flow created by surface tension due to a temperature gradient. 3.2. Simulation of Marangoni convection with free surface Our simulation was conducted using the model shown in Fig. 1. The cross-sectional streamlines at the center in the y direction and the surface of the molten W every 0.1 s for as long as 0.8 s are shown in Fig. 3. The temperature distribution and velocity vectors at the same cross section but in molten W at 0.2 s and 0.6 s are shown in Fig. 4. In the initial stage of this simulation, the inner flow is generated in molten W by a surface tension gradient due to the temperature gradient, as shown in Fig. 3(a). Because molten W has a negative surface tension coefficient, the surface tension increases with a decrease in temperature and becomes strongest near the cold wall. Therefore, the direction of the inner flow is from the hot wall to the cold wall near the surface, as shown in Fig. 4(a). As shown in Fig. 3(b)–(d), a small vortex is also generated at the hot wall side of a large vortex in the inner flow. This small vortex is a separation vortex due to the inner flow and large vortex at the cold wall side. Thus, the directions of both flows’ vortices are counterclockwise. As the large vortex at the cold wall side develops, the separation vortex also becomes larger. However, from ∼0.5 s, the separation vortex disappears via the development of the large vortex at the cold wall side, and the flow at the hot wall side is weakened, as shown in Fig. 3(e). Then, the vortex is generated at the hot wall side again; however, the flow direction of the new vortex at the hot wall side
Fig. 2. Comparison of the velocity distribution of an inner liquid.
conditions, (∂σ ⁄ ∂T), μ, α, and ΔΤ are 3.1 × 10−4 N/m/K, 5.0 × 10−3 Pa s, 2.2 × 10−5 m2/s [11,18], and 1000 K respectively, and the Ma becomes 400. A slip condition was applied to the xz planes at y = 0 and to the end of the xz planes in the y direction; other planes were nonslip wall conditions. For verification, the results of our simulation were compared with those of the numerical simulation reported by Smith and Davis [8]. Surface deformation was not considered in their model, in which the velocity perpendicular to the surface was set to zero and the only velocity parallel to it was considered at the surface. They evaluated the inner flow in a liquid layer in a case in which the temperature gradient parallel to the surface was given, as with our model. It is considered that the comparison of the velocities in the case of Ma = 400 is shown in Fig. 2. The vertical axis indicates the position in the z direction normalized by surface position; thus, zero means the bottom position of the liquid phase. The horizontal axis indicates the velocity distribution
Fig. 3. Streamlines of the inner molten W and the surface shape of the molten W at each time.
3
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
Fig. 4. Temperature distribution and velocity vectors for the inner molten W.
Fig. 5. Temperature gradient and velocity vector at each time, and streamlines and surface shape.
is clockwise. The flow at the center in the x direction becomes vertical upward flow, and both vortices develop, as shown in Figs. Fig. 33(f) and Fig. 44(b). The vertical upward flow causes surface deformation, and the two vortices inner molten W are disturbed because of the flow instability, as shown in Fig. 3(g) and (h). The inner flow pattern then becomes complex, and such disturbed vortices lead to subsequent surface deformation. If the temperature difference varies, the flow in molten W is considered to change. Here, another simulation corresponding to a temperature difference of 10 K was performed in the same simulation model as that shown in Fig. 1; the Ma in this case became 40. The crosssectional streamline, temperature distribution, and the velocity vector at the center in the y direction are shown as Fig. 5. In addition, Fig. 6 shows the profiles of the z component of the velocity plotted along the dotted line at the half depth of molten W in Fig. 5(a).
Fig. 6. Profiles of the z component of velocity plotted along the dotted line in Fig. 4(a) at each time.
Fig. 7. Schematic of the depth and flow patterns of a molten metal pool loaded by high heat flux [1].
4
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
lowest in the center region of the surface. Therefore, surface flow is driven toward the outside, as shown Fig. 7(a). However, in the case of the surface tension coefficient, the driving force at the surface is in the opposite direction, as shown in Fig. 7(b). Because the surface geometry was fixed as a boundary condition in their report, Tsotridis et al. [1] could not consider surface deformation. However, actually, the surface of materials loaded by high heat flux can be deformed with molten-phase flow; thus, the model of Tsotridis et al. cannot reproduce the actual behaviors of such materials. Therefore, we constructed a simulation code that can treat two-phase flow including the interface tracking for molten W and air. We conducted the simulation of laser irradiation at molten W, as shown in Fig. 8. In our simulation, the diameter of the laser was 0.6 mm, the pulse duration was 1 ms, and the energy density on the surface was 10 MW/m2. The results of the velocity vector and temperature distribution in molten W at the center in the y direction are shown in Fig. 9(a). Those corresponding to the cross-sectional streamline and the surface shape of molten W at 0.8 s are shown in Fig. 9(b). Because of the heat flux, a temperature distribution arose at the surface of molten W. ΔΤ meant the temperature difference between the highest point and the lowest point at the surface; the Ma then became 280. Marangoni convection occurred, and the surface flow was driven toward the outside, as shown in Fig. 9(a). This result is in qualitatively good agreement with the findings of Tsotridis et al. [1], as shown in Fig. 7(a). The direction of Marangoni convection is similar to the case in Section 3.2. Surface deformation due to Marangoni convection could also be observed in this case. Fig. 10 shows profiles of the z component of the velocity at the surface of molten W from 0.6 s to 0.8 s after every 0.05 s. The vertical and horizontal axes respectively indicate the z component of the velocity and the position in the x direction. The z component of the velocity in the periphery with a large surface tension changes substantially. From 0.7 s, the z component of the velocity increases remarkably at the center of the surface because of the vertical upward flow, as shown in Fig. 9(a). Thus, surface deformation occurs at the center of the surface, as shown in Fig. 9(b). The warmed molten W is transported to the bottom via the periphery of the surface while cooling by releasing heat. Therefore, the heat transferred to the interior by heat conduction could melt the solid layer under the liquid layer, but the heat transported by Marangoni convection in molten W could not promote the melting of the solid layer. However, the molten region could spread in the radial direction through Marangoni convection.
Fig. 8. Simulation model of laser irradiation.
In Fig. 6, the main vertical axis represents the z component of the velocity at 0.5, 1.0, and 1.2 s. The second vertical axis is for the case at 1.5 s. The horizontal axis is the position in the x direction. The inner flow with separated flow is also generated in this case. Two vortices with counterclockwise direction are then observed as shown in Fig. 5(a) and (b). As the inner flow develops further, the vortex at the cold wall side becomes large and the separation vortex disappears, as shown in Fig. 5(c). Finally, the vortex at the hot wall side with the clockwise direction is induced, as shown in Fig. 5(d). In Fig. 6, the z component of the velocity is shown to change rapidly from 1.2 to 1.5 s. This rapid change causes flow instability, but surface deformation was not observed at Ma = 40. 3.3. Marangoni convection in a molten W pool irradiated by high heat flux Tsotridis et al. reported that surface tension coefficients of materials influence the flow patterns in molten pools of the materials, as shown in Fig. 7 [1]. In their model, high heat flux, assumed to occur via plasma disruption, was simulated by laser irradiation at the solid. The heat flux melted the surface, and a temperature distribution arose in the molten pool. Because of the temperature distribution, the surface tension decreased or increased radially from the center to the edge; thus, Marangoni convection occurred. In the case of a negative surface tension coefficient, as observed for many pure metals, the surface tension decreases with increasing temperature. The surface tension becomes highest in the colder region of the surface, i.e., at the periphery of the molten metal, and becomes
Fig. 9. Temperature distribution, velocity vector, streamlines, and surface shape at 0.8 s.
5
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
conducted for the W divertor. References [1] G. Tsotridis, H. Rother, E.D. Hondros, On modelling of Marangoni convection flows in simulated plasma disruptions, Fusion Eng. Des. 15 (1991) 155–162. [2] J. Thomson, On certain curious motions observable at the surfaces of wine and other alcoholic liquors, Philos. Mag. Ser. 4 (1855) 330. [3] C. Marangoni, Ann. Phys. Chem. (Poggendorff) 143 (1871) 337. [4] C. Marangoni, J. Plateau, N. Cimento Giornale di Fisica, C. e storia Naturale, Ser.2, V-VI, (1871), p. 239. [5] D.A. Nield, Surface tension and buoyancy effects in cellular convection, J. Fluid Mech. 19 (1964) 341–352. [6] N. Imanishi, Fundamental of the marangoni convection, Int. J. Microgravity Sci. (Supplement) (2014). [7] Y.R. Li, L. Peng, S.Y. Wu, N. Imaishi, D.L. Zeng, Thermocapillary-buoyancy flow of silicon melt in a shallow annular pool, Cryst. Res. Technol. 39 (2004) 1055–1062. [8] M.K. Smith, S.H. Davis, Instabilities of dynamic thermocapillary liquid layers, J. Fluid Mech. 132 (1983) 114–144. [9] Y. Ueda, et al., Effects of repetitive ELM-like heat pulses on tungsten surface morphology, Fusion Eng. Des. (2007) 1904–1910. [10] D. Inoue, K. Ibano, S. Yoshikawa, T. Maeji, Y. Ueda, Molten layer characteristics of W materials and film coated W by pulsed laser irradiation, Fusion Eng. Des. 124 (2017) 316–320. [11] T. Ishikawa, Thermophysical properties of molten tungsten measured with an electrostatic levitator, Netsu Bussei 19 (2005) 61–66 (in Japanese). [12] T. Yabe, T. Ishikawa, P.Y. Wang, T. Aoki, Y. Kadota, F. Ikeda, Comp. Phys. Commun. 66 (1991) 223. [13] T. Yabe, T. Aoki, A universal solver for hyperbolic equations by cubic-polynomial interpolation I. One-dimensional solver, Comp. Phys. Commun. 66 (1991) 219. [14] E. Hoashi, T. Yokomine, A. Shimizu, T. Kunugi, Numerical analysis of wave-type heat transfer propagating in a thin foil irradiated by short-pulsed laser, Int. J. Heat Mass Transf. 46 (2003) 4083–4095. [15] T. Yabe, F. Xiao, Description of complex and sharp interface with fixed grids in incompressible and compressible fluid, Comput. Math. Appl. 29 (1995) 15–25. [16] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comp. Phys. 100 (1992) 335. [17] E.B. Royce, Gray, A Three-Phase Equation of State for Metals, (1971). [18] M.W. Chase, NIST-JANAF thermochemical tables for oxygen fluorides, J. Phys. Chem. Ref. Data 25 (1996) 551–603.
Fig. 10. Profiles of the z component of the velocity at the surface of molten W at each time.
4. Conclusions Numerical simulations that could evaluate Marangoni convection in molten W with free surface were conducted to clarify the process of Marangoni convection. When a temperature gradient parallel to the surface of molten W was applied, Marangoni convection occurred and two vortices were generated. According to the Ma, surface deformation was observed as Marangoni convection developed. However, the convection after the two vortices is disturbed owing to the flow instability has not been clarified; thus, it should therefore be evaluated in the future. In addition, the evaluation of surface deformation due to Marangoni convection in molten W would be required subsequently for determining surface shape after solidification. In the situation where heat flux was applied to the surface of the molten W, the vortices were generated and surface deformation was caused by Marangoni convection. Hence, a detailed evaluation of the influence of Ma on Marangoni convection and on the surface deformation for molten W should be
6
Contents lists available at ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Development of a simulation method for evaluating Marangoni convection with free surface for tungsten divertor ⁎
Kohei Hamaguchi , Yu Teramoto, Eiji Hoashi, Takafumi Okita, Kenzo Ibano, Yoshio Ueda Graduate School of Engineering, Osaka University, Osaka, Japan
A R T I C L E I N F O
A B S T R A C T
Keywords: Tungsten divertor Molten tungsten Marangoni convection Free surface Heat flux
A tungsten (W) divertor has been developed for tokamak-type nuclear fusion reactors, such as ITER. Because the divertor will be repetitively subjected to high heat flux when plasma disruption or edge localized mode occurs, W at the surface of the divertor will be melted and solidified during reactor operation. The surface shape of the W divertor after or during such a heat load strongly influences its lifetime. Thus, evaluating the behavior of W in such conditions is important. An important driving force in molten metal is Marangoni convection due to a difference in surface tension resulting from a temperature gradient at the surface. To clarify the mechanism of Marangoni convection in molten W, a simulation method that can simultaneously treat surface deformation and convection in molten W, including Marangoni convection, is required. We have been developing such a simulation code based on the constrained interpolation profile combined unified procedure method. In this paper, we simulated molten W under a given temperature gradient or heat load and validated the model by simulating the flow in a molten W pool with a temperature gradient. The development process of Marangoni convection in molten W with a free surface was evaluated.
1. Introduction Tungsten (W) has been studied for the candidate material of plasma facing materials, such as divertor, in tokamak-type nuclear fusion reactor. W has the highest melting point among all metals, a high thermal conductivity, and a low sputtering rate, thus making it a promising divertor surface material. In tokamak-type nuclear fusion reactors, such as ITER, because the divertor will be subjected to a high heat flux of at least 100 MW/m2 within several milliseconds or less when plasma disruption or type-I edge localized mode (ELM) occurs, the W in the divertor will undergo repeated melting and solidification [1]. Thus, the surface shape of the W divertor may change after solidification. If the surface shape would be smooth, it is expected the W will stand up to subsequent repetitive heat load. Because surface roughness can decrease the thermal resistance, the surface status after heat load influences on the lifetime of the W divertor. Because controlling plasma disruption and ELM is very difficult, knowledge of the surface status after solidification is required to ensure the lifetime. To acquire such knowledge, it is necessary to investigate and understand the details of melting and solidification process of W subjected to high heat flux. Previous research has indicated that molten state flow is one of the main factors that determine the surface shape [1]. One primary factor that controls flow in a molten metal is the
⁎
Marangoni effect, which is caused by a difference in the surface tension resulting from a temperature gradient at the surface. Convection due to the Marangoni effect is called Marangoni convection and has been previously investigated [2–4]. A liquid layer with a uniform surface tension maintains a stationary state. However, when the temperature or concentration at the liquid surface is not uniform, movement of the liquid surface occurs and leads to Marangoni convection. In the analysis of Marangoni convection, a dimensionless number, which is called the Marangoni number (Ma), represents the strength and instability of the convection; it is expressed as follows:
Ma = −
(∂σ / ∂T )ΔT d, μα
(1)
where T, ΔΤ, μ, α, and d are temperature, characteristic temperature difference, viscosity coefficient, thermal diffusion coefficient, and the characteristic length, respectively. σ represents the surface tension coefficient depending on the temperature. The characteristic length d corresponds to the depth of the liquid layer. Nield has reported that Marangoni convection becomes more dominant than buoyancy in a thin liquid layer, such as a layer that is a few millimeters thick [5]. Extensive research on Marangoni convection for fluids such as water, which have relatively high Prandtl numbers (Pr) than molten W, has been reported [6]. For example, the Pr of water and molten W are 7
Corresponding author. E-mail address: [email protected] (K. Hamaguchi).
https://doi.org/10.1016/j.fusengdes.2018.02.007 Received 29 September 2017; Received in revised form 2 February 2018; Accepted 2 February 2018 0920-3796/ © 2018 Elsevier B.V. All rights reserved.
Please cite this article as: Hamaguchi, K., Fusion Engineering and Design (2018), https://doi.org/10.1016/j.fusengdes.2018.02.007
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
component Fsa(t) is expressed as follows:
and 0.015, respectively. However, the literature contains few reports on the process of Marangoni convection in low-Pr fluids. Li et al. evaluated the relevance of Marangoni convection and buoyancy in a circular pool using molten silicon (Pr = 0.011) as a low-Pr fluid [7]. Smith and Davis indicated the instabilities of flows for high- and low-Pr fluids [8] in twodimensional (2D) simulations. In addition, the influence of impurities on Marangoni convection flow patterns in a molten layer was reported by Tsotridis et al. [1]. However, research on Marangoni convection using the molten W has not been reported. Ueda and Inoue et al. have conducted experiments on W irradiated by a laser with a duration time from microseconds to milliseconds [9,10]. According to their observation of the W surface after the laser irradiation, the center of the irradiation area was found to have swollen like a cone. However, the formation mechanism of the cone by heat flux was not fully understood. W has a high melting point and must be maintained in an ultra high purity state to prevent surface oxidation; it is thus, are difficult to perform experiments with the molten W. In addition, measuring inner flow is difficult because the molten W is metal with an opaque surface. Thus, we have been aiming to develop a three-dimensional (3D) unified simulation method that can evaluate the melting and solidification process of W, including phase change. The surface of the molten W divertor becomes a free surface; however, free surfaces have not been considered in the research conducted by researches such as Smith and Davis [8], Tsotridis et al. [1]. Thus, our simulation code newly considers free surface. In our simulation, we use the properties of molten W reported by Ishikawa et al., who conducted measurements using an electrostatic levitator [11]. In this study, we develop the simulation code that can reproduce the Marangoni effect in surface tension model for future evaluation of melting and solidification process of W, and then we clarify the development process of Marangoni convection of W with the free surface using this code.
Fsa(t) = ∇σ.
(3)
Temperature dependence of surface tension coefficient σ is defined as following empirical equation obtained from experimental data [11]:
σ (T ) = 2.48 − 0.31 × (T − Tm ) × 10−3
(4)
where Tm is the melting point. 2.2. Governing equations As the governing equations, the equation of continuity (5), momentum equation (6), equation of energy (7), and EOS (8) are used in our simulations as follows:
Dρ ∂u = −ρ i , Dt ∂x i
(5)
Dui 1 ∂P 1 ∂Sij + gi , =− + Dt ρ ∂x i ρ ∂x i
(6)
Sij ∂ui DT 1 ∂ ⎛ ∂T ⎞ Q P ∂u = − TH i + + κ + , and Dt ρCV ∂x i ρCV ∂x i ρCV ∂x i ⎝ ∂x i ⎠ ρCV
(7)
P = P(ρ, T)
(8)
⎜
⎟
where ρ, P, S, Cv, κ, and Q indicate the density, pressure, stress or viscous stress, specific heat at constant volume, thermal conductivity, and the amount of volumetric heat generated, respectively. Parameter PTH is calculated according to the expression PTH = T(∂P/∂T). The GRAY-EOS developed at Lawrence Livermore National Laboratory [17] is used as EOS.
2. Material and methods
3. Results and discussion
2.1. Simulation method
3.1. Comparison of our results with those of previous studies
The constrained interpolation profile (CIP) method, which enables the accurate solving of the advective term, was proposed by Yabe et al. [12]. In addition, the CIP combined unified procedure (C-CUP) method based on the CIP method has been developed. The C-CUP method is suitable for managing both compressible and incompressible fluids [13]. However, temperature is not treated in the C-CUP method. Hoashi et al. developed a method that can treat temperature with the C-CUP method by determining it using the equation of energy (EOS) [14]. In the analysis of the flow in which multiple phases coexist, Yabe et al. used a density function to capture the interface [12]. We considered a method that can solve the density function with tangent conversion to decrease numerical diffusion to capture the interface accurately [15]. Surface tension should be considered because it is one of the primary factors, such as Marangoni convection, that affect inner flow. Thus, we constructed the surface tension model based on the continuous surface force (CSF) model suggested by Brackbill et al. [16]. The total surface forces Fsa is expressed as follows:
Because the heat flux that W undergoes at the surface of the divertor will vary for different locations on the divertor, a temperature gradient parallel to the surface will occur. However, the mechanism by which Marangoni convection develops in molten W with a free surface has not been clarified. Therefore, we considered the 3D simulation model shown in Fig. 1. As an initial condition, molten W was placed in the lower region and air was placed in the upper region. Air was assumed to behave as an ideal gas, and the initial velocity in all regions was set to zero. The plane at x = 0 was treated as the cold wall with a temperature of 4000 K, and the plane opposite the cold wall in direction x was treated as the hot wall with a temperature of 4100 K. Under these
Fsa = Fsa(n) + Fsa(t) Fsa(n)
(2) (t)
is the normal component to the interface and Fsa is the where tangential component. Fsa(n) is obtained by a curvature and surface tension coefficient according to CSF model. A curvature is obtained by a density function. Calculating the tangential component Fsa(t) is required so that Marangoni convection occurs due to the contribution of Fsa(t). Thus, it is necessary to specify the direction in which Fsa(t) works at the interface. We used the polar coordinate transformation for specifying the direction. Then, the cell, where the interface was adjacent, was determined. The difference in the surface tension between those cells becomes the tangential component Fsa(t). Therefore, the tangential
Fig. 1. Schematic of the simulation model.
2
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
in the x direction for liquid phase normalized by the velocity at the surface. The results of our simulation are consistent with those of Smith and Davis. Therefore, our simulation is sufficiently accurate to reproduce the inner flow created by surface tension due to a temperature gradient. 3.2. Simulation of Marangoni convection with free surface Our simulation was conducted using the model shown in Fig. 1. The cross-sectional streamlines at the center in the y direction and the surface of the molten W every 0.1 s for as long as 0.8 s are shown in Fig. 3. The temperature distribution and velocity vectors at the same cross section but in molten W at 0.2 s and 0.6 s are shown in Fig. 4. In the initial stage of this simulation, the inner flow is generated in molten W by a surface tension gradient due to the temperature gradient, as shown in Fig. 3(a). Because molten W has a negative surface tension coefficient, the surface tension increases with a decrease in temperature and becomes strongest near the cold wall. Therefore, the direction of the inner flow is from the hot wall to the cold wall near the surface, as shown in Fig. 4(a). As shown in Fig. 3(b)–(d), a small vortex is also generated at the hot wall side of a large vortex in the inner flow. This small vortex is a separation vortex due to the inner flow and large vortex at the cold wall side. Thus, the directions of both flows’ vortices are counterclockwise. As the large vortex at the cold wall side develops, the separation vortex also becomes larger. However, from ∼0.5 s, the separation vortex disappears via the development of the large vortex at the cold wall side, and the flow at the hot wall side is weakened, as shown in Fig. 3(e). Then, the vortex is generated at the hot wall side again; however, the flow direction of the new vortex at the hot wall side
Fig. 2. Comparison of the velocity distribution of an inner liquid.
conditions, (∂σ ⁄ ∂T), μ, α, and ΔΤ are 3.1 × 10−4 N/m/K, 5.0 × 10−3 Pa s, 2.2 × 10−5 m2/s [11,18], and 1000 K respectively, and the Ma becomes 400. A slip condition was applied to the xz planes at y = 0 and to the end of the xz planes in the y direction; other planes were nonslip wall conditions. For verification, the results of our simulation were compared with those of the numerical simulation reported by Smith and Davis [8]. Surface deformation was not considered in their model, in which the velocity perpendicular to the surface was set to zero and the only velocity parallel to it was considered at the surface. They evaluated the inner flow in a liquid layer in a case in which the temperature gradient parallel to the surface was given, as with our model. It is considered that the comparison of the velocities in the case of Ma = 400 is shown in Fig. 2. The vertical axis indicates the position in the z direction normalized by surface position; thus, zero means the bottom position of the liquid phase. The horizontal axis indicates the velocity distribution
Fig. 3. Streamlines of the inner molten W and the surface shape of the molten W at each time.
3
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
Fig. 4. Temperature distribution and velocity vectors for the inner molten W.
Fig. 5. Temperature gradient and velocity vector at each time, and streamlines and surface shape.
is clockwise. The flow at the center in the x direction becomes vertical upward flow, and both vortices develop, as shown in Figs. Fig. 33(f) and Fig. 44(b). The vertical upward flow causes surface deformation, and the two vortices inner molten W are disturbed because of the flow instability, as shown in Fig. 3(g) and (h). The inner flow pattern then becomes complex, and such disturbed vortices lead to subsequent surface deformation. If the temperature difference varies, the flow in molten W is considered to change. Here, another simulation corresponding to a temperature difference of 10 K was performed in the same simulation model as that shown in Fig. 1; the Ma in this case became 40. The crosssectional streamline, temperature distribution, and the velocity vector at the center in the y direction are shown as Fig. 5. In addition, Fig. 6 shows the profiles of the z component of the velocity plotted along the dotted line at the half depth of molten W in Fig. 5(a).
Fig. 6. Profiles of the z component of velocity plotted along the dotted line in Fig. 4(a) at each time.
Fig. 7. Schematic of the depth and flow patterns of a molten metal pool loaded by high heat flux [1].
4
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
lowest in the center region of the surface. Therefore, surface flow is driven toward the outside, as shown Fig. 7(a). However, in the case of the surface tension coefficient, the driving force at the surface is in the opposite direction, as shown in Fig. 7(b). Because the surface geometry was fixed as a boundary condition in their report, Tsotridis et al. [1] could not consider surface deformation. However, actually, the surface of materials loaded by high heat flux can be deformed with molten-phase flow; thus, the model of Tsotridis et al. cannot reproduce the actual behaviors of such materials. Therefore, we constructed a simulation code that can treat two-phase flow including the interface tracking for molten W and air. We conducted the simulation of laser irradiation at molten W, as shown in Fig. 8. In our simulation, the diameter of the laser was 0.6 mm, the pulse duration was 1 ms, and the energy density on the surface was 10 MW/m2. The results of the velocity vector and temperature distribution in molten W at the center in the y direction are shown in Fig. 9(a). Those corresponding to the cross-sectional streamline and the surface shape of molten W at 0.8 s are shown in Fig. 9(b). Because of the heat flux, a temperature distribution arose at the surface of molten W. ΔΤ meant the temperature difference between the highest point and the lowest point at the surface; the Ma then became 280. Marangoni convection occurred, and the surface flow was driven toward the outside, as shown in Fig. 9(a). This result is in qualitatively good agreement with the findings of Tsotridis et al. [1], as shown in Fig. 7(a). The direction of Marangoni convection is similar to the case in Section 3.2. Surface deformation due to Marangoni convection could also be observed in this case. Fig. 10 shows profiles of the z component of the velocity at the surface of molten W from 0.6 s to 0.8 s after every 0.05 s. The vertical and horizontal axes respectively indicate the z component of the velocity and the position in the x direction. The z component of the velocity in the periphery with a large surface tension changes substantially. From 0.7 s, the z component of the velocity increases remarkably at the center of the surface because of the vertical upward flow, as shown in Fig. 9(a). Thus, surface deformation occurs at the center of the surface, as shown in Fig. 9(b). The warmed molten W is transported to the bottom via the periphery of the surface while cooling by releasing heat. Therefore, the heat transferred to the interior by heat conduction could melt the solid layer under the liquid layer, but the heat transported by Marangoni convection in molten W could not promote the melting of the solid layer. However, the molten region could spread in the radial direction through Marangoni convection.
Fig. 8. Simulation model of laser irradiation.
In Fig. 6, the main vertical axis represents the z component of the velocity at 0.5, 1.0, and 1.2 s. The second vertical axis is for the case at 1.5 s. The horizontal axis is the position in the x direction. The inner flow with separated flow is also generated in this case. Two vortices with counterclockwise direction are then observed as shown in Fig. 5(a) and (b). As the inner flow develops further, the vortex at the cold wall side becomes large and the separation vortex disappears, as shown in Fig. 5(c). Finally, the vortex at the hot wall side with the clockwise direction is induced, as shown in Fig. 5(d). In Fig. 6, the z component of the velocity is shown to change rapidly from 1.2 to 1.5 s. This rapid change causes flow instability, but surface deformation was not observed at Ma = 40. 3.3. Marangoni convection in a molten W pool irradiated by high heat flux Tsotridis et al. reported that surface tension coefficients of materials influence the flow patterns in molten pools of the materials, as shown in Fig. 7 [1]. In their model, high heat flux, assumed to occur via plasma disruption, was simulated by laser irradiation at the solid. The heat flux melted the surface, and a temperature distribution arose in the molten pool. Because of the temperature distribution, the surface tension decreased or increased radially from the center to the edge; thus, Marangoni convection occurred. In the case of a negative surface tension coefficient, as observed for many pure metals, the surface tension decreases with increasing temperature. The surface tension becomes highest in the colder region of the surface, i.e., at the periphery of the molten metal, and becomes
Fig. 9. Temperature distribution, velocity vector, streamlines, and surface shape at 0.8 s.
5
Fusion Engineering and Design xxx (xxxx) xxx–xxx
K. Hamaguchi et al.
conducted for the W divertor. References [1] G. Tsotridis, H. Rother, E.D. Hondros, On modelling of Marangoni convection flows in simulated plasma disruptions, Fusion Eng. Des. 15 (1991) 155–162. [2] J. Thomson, On certain curious motions observable at the surfaces of wine and other alcoholic liquors, Philos. Mag. Ser. 4 (1855) 330. [3] C. Marangoni, Ann. Phys. Chem. (Poggendorff) 143 (1871) 337. [4] C. Marangoni, J. Plateau, N. Cimento Giornale di Fisica, C. e storia Naturale, Ser.2, V-VI, (1871), p. 239. [5] D.A. Nield, Surface tension and buoyancy effects in cellular convection, J. Fluid Mech. 19 (1964) 341–352. [6] N. Imanishi, Fundamental of the marangoni convection, Int. J. Microgravity Sci. (Supplement) (2014). [7] Y.R. Li, L. Peng, S.Y. Wu, N. Imaishi, D.L. Zeng, Thermocapillary-buoyancy flow of silicon melt in a shallow annular pool, Cryst. Res. Technol. 39 (2004) 1055–1062. [8] M.K. Smith, S.H. Davis, Instabilities of dynamic thermocapillary liquid layers, J. Fluid Mech. 132 (1983) 114–144. [9] Y. Ueda, et al., Effects of repetitive ELM-like heat pulses on tungsten surface morphology, Fusion Eng. Des. (2007) 1904–1910. [10] D. Inoue, K. Ibano, S. Yoshikawa, T. Maeji, Y. Ueda, Molten layer characteristics of W materials and film coated W by pulsed laser irradiation, Fusion Eng. Des. 124 (2017) 316–320. [11] T. Ishikawa, Thermophysical properties of molten tungsten measured with an electrostatic levitator, Netsu Bussei 19 (2005) 61–66 (in Japanese). [12] T. Yabe, T. Ishikawa, P.Y. Wang, T. Aoki, Y. Kadota, F. Ikeda, Comp. Phys. Commun. 66 (1991) 223. [13] T. Yabe, T. Aoki, A universal solver for hyperbolic equations by cubic-polynomial interpolation I. One-dimensional solver, Comp. Phys. Commun. 66 (1991) 219. [14] E. Hoashi, T. Yokomine, A. Shimizu, T. Kunugi, Numerical analysis of wave-type heat transfer propagating in a thin foil irradiated by short-pulsed laser, Int. J. Heat Mass Transf. 46 (2003) 4083–4095. [15] T. Yabe, F. Xiao, Description of complex and sharp interface with fixed grids in incompressible and compressible fluid, Comput. Math. Appl. 29 (1995) 15–25. [16] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comp. Phys. 100 (1992) 335. [17] E.B. Royce, Gray, A Three-Phase Equation of State for Metals, (1971). [18] M.W. Chase, NIST-JANAF thermochemical tables for oxygen fluorides, J. Phys. Chem. Ref. Data 25 (1996) 551–603.
Fig. 10. Profiles of the z component of the velocity at the surface of molten W at each time.
4. Conclusions Numerical simulations that could evaluate Marangoni convection in molten W with free surface were conducted to clarify the process of Marangoni convection. When a temperature gradient parallel to the surface of molten W was applied, Marangoni convection occurred and two vortices were generated. According to the Ma, surface deformation was observed as Marangoni convection developed. However, the convection after the two vortices is disturbed owing to the flow instability has not been clarified; thus, it should therefore be evaluated in the future. In addition, the evaluation of surface deformation due to Marangoni convection in molten W would be required subsequently for determining surface shape after solidification. In the situation where heat flux was applied to the surface of the molten W, the vortices were generated and surface deformation was caused by Marangoni convection. Hence, a detailed evaluation of the influence of Ma on Marangoni convection and on the surface deformation for molten W should be
6