A Radioactive High Level Waste (HLW) is produced as a byproduct of nuclear power plant operations and it is important to develop methods to reprocess HLW. In reprocessing, the HLW is dissolved and becomes a High Level Liquid Waste (HLLW), and the HLLW is poured into molten borosilicate glass in a glass melter to make a stable mixture of HLLW and glass for geological disposal. In Japan, a Liquid-Fed Ceramic Melter (LFCM) type glass melter is being developed for the reprocessing operation. The LFCM glass melter at Tokai Vitrification Facility (TVF) is composed of a cuboid upper part and a lower pyramid-shape part as shown in Figure 1 . The glass melter applies Joule-heating to melt a glass material and HLLW is poured into the molten glass. The melter can mix HLLW and molten glass by convective flow mainly induced by Joule-heating and radiation from the upper surface of the melter. The volumetric heating in the lower part of the melter and cooling on the top surface cause continuous chaotic flow behavior; this is known as a “chaotic steady state” .
Many factors affect the flow behavior such as electrode cooling , a cold cap on the top surface , the presence of platinum group particles , the existence of foaming reagent , etc. Hence, the flow in the glass melter involves a multiphase flow of the molten glass containing the HLLW, a foaming reagent (gas),
Figure 1. LFCM Joule-heating glass melter .
the solid phase cold cap and unmelted platinum group particles. Probably due to the complex multiphase flow, sudden temperature changes or other events sometimes occur in the melter. Since an understanding of the chaotic flow behavior is inadequate, the melter has to be stopped to avoid severe accidents. Therefore, to operate the glass melter effectively, flow behavior in the melter should be better understood. However, details of this multiphase flow are difficult to obtain because many phenomena occur in the melter and they affect the flow behavior. Thus, the melter geometry was simplified to better understand the chaotic flow behavior and a cubic cavity was used in a simplified case .
Since the glass melter has a pair of electrodes on the cuboid upper part , the chaotic flow induced by the Joule-heating will occur mainly in the upper part. Changes in the heating area may affect the behavior of Joule-heating induced flow. The effects of cavity geometry on the flow behavior would also depend on the shape of the cavity, so the model cavity has a similar 2-D shape as that of the glass melter and is called a “sloping bottom cavity”.
In measuring the flow behavior of molten glass, the high temperature of the molten glass makes it difficult to observe the flow. The Ultrasound Velocity Profiler (UVP) method is one way to observe this high temperature flow by applying a buffer rod . This method is non-invasive, simple to use and applicable to opaque fluids such as molten glass and time-series velocity measurements. However, a difficulty in the UVP method is the temperature distribution which changes the sound speed in the medium. Therefore, the UVP method is applied with another optical method, particle image velocimetry (PIV) in the case of Joule-heating induced flow.
Besides measuring the molten glass flow directly, CFD is a useful tool applied to predict the Joule-heating flow behavior . However, the previous numerical studies focused on the flow and temperature fields and the magnetic field was ignored. In the glass melter, the flow behavior needs to be considered based on the flow, temperature and magnetic fields. To develop flow analysis and fluid measurement techniques for a molten glass, a CFD code is required to account for these three fields. A GSMAC-FEM (Generalized Simplified Marker and Cell-Finite Element Method) code has a feature that can cross-couple an electromagnetic field and a flow field. Since the glass melter uses a high electrical current, the fully coupled analysis of the three fields is valuable when the electromagnetic field effect becomes large. Therefore, GSMAC-FEM is used in the present numerical simulation.
In this paper, the velocity profile on the vertical center line is measured by the UVP method, and the PIV method is used to observe the two-dimensional flow behavior, while the GSMAC method is applied to calculate the 3-D flow in the cavity. The flow fluctuations measured by the UVP method and predicted by the GSMAC method are compared. The fluctuation frequency is analyzed by using the Fast Fourier Transform (FFT) method. Numerical predictions of the flow behaviors are compared with the experimental data.
2. Experimental Apparatus and Method
2.1. Experimental Apparatus
A schematic drawing of the sloping bottom cavity is shown in Figure 2. The cavity has dimensions of 126 × 152 × 154 (mm) in the x, y, z-direction, respectively, and is a 1/5 scale model of a TVF melter. The representative length of the cavity is 126 mm and the cavity has a cuboid upper part and a lower trapezoidal part, both of which are of the same shape as in the glass melter. A pair of electrode plates made of graphite is set on two opposite side walls of the cavity. Below each electrode, the cavity has a sloping surface. The ratio of the heating area to the non-heating area is designed to be similar to that of the glass melter. The top surface of the cavity is cooled which simulates the cold cap in the glass melter. Heat sinks made of copper are placed on the top surface of the electrodes to cool them. The heat sinks can be changed to an adiabatic acrylic board to change the thermal boundary condition for the electrodes from isothermal to adiabatic. The working fluid is a 80 wt% glycerin-water solution with LiCl as the electrolyte. However, the viscosity of molten glass is much higher than the glycerol solution (3.81 × 10−2 m2/s against 4.96 × 10−5 m2/s). Thus, the internal Rayleigh number should be of the similar order between the experiment and actual glass melter. The internal Rayleigh number is defined as follows:
The experimental apparatus is shown in Figure 3. AC power is applied to generate Joule-heating. The temperature profile is measured in the center of the cavity using a sheathed K-type thermocouple probe. To reduce the effect of the thermocouple on the flow, the thermocouple probe of 0.5 mm diameter is utilized. The temperature and current data are sent to a PC using a data logger. A chiller is connected to a heat sink on the top surface of the test section separately to maintain the coolant temperature. The side electrodes were adiabatic in the experiment. A UVP transducer is set underneath the cavity while the PIV’s laser sheet is set slightly away from the UVP beam. Nylon particles (d = 80 μm; ρ = 1.020 g/cm3) are scattered in the working fluid as tracer particles for both UVP and PIV methods.
In the sloping bottom cavity, because the laser light is reflected by the sloping wall, the laser is set to the bottom left side of the cavity, so the flow behavior near the wall is observed by the PIV method. The UVP transducer measures the vertical flow behavior on the measurement line in the test section. Results of the UVP method are the focus of the present experiments.
2.2. Measurement Methods for Joule-Heating Induced Flow
A schematic of the UVP measurement method is depicted in Figure 4. In the UVP method developed by Takeda , an ultrasound transducer emits a pulse and receives the echo signal reflected from the particles suspended in the liquid on the measurement line X where the ultrasound beam path in the liquid. Based on an echo signal analysis of ultrasonic pulses reflected by particles suspended in the fluid at each position along the measurement line X, instantaneous velocities Vx of the suspended particles on the measurement line can be obtained. The particle position in each channel can be extracted from the time delay τPRF as follows:
Figure 2. Sloping bottom cavity.
Figure 3. Experimental apparatus.
Figure 4. Schematic of measurement by the UVP.
where c is the sound speed in the medium. At the same time, an instantaneous local velocity, vuvp(x), is derived from the Doppler shift frequency (fD) as:
where f0 is the basic ultrasound frequency. The UVP method measures the flow velocity Vx on the measurement line X using the Doppler shift in the ultrasound frequency. Therefore, the UVP method can measure velocities in opaque fluids. However, the signal analysis in the UVP method requires the sound speed which depends on the temperature of the medium. Therefore, the temperature measurement in the medium is necessary. In the present work, the maximum temperature variation was less than 20˚C, so the error in the velocity measurement by UVP due to the fluid temperature variation was small .
In the experiment, a commercial UVP system (UVP-DUO) was used. The measurement conditions for the UVP method are shown in Table 1. The basic frequency of the ultrasound was set at 4 MHz, which can pass through a glycerin-water solution easily. Since the Joule-heating induced flow is very slow, the velocity resolution was chosen at the maximum value. The sound speed was calculated using the median fluid temperature in the cavity. Because the flow in the cavity was chaotic, the measurement channel width and sampling interval were set at the minimum values to obtain accurate measurements.
The PIV method was used in the 2-D flow behavior measurement. A laser light sheet was set to illuminate the tracer particles dispersed in the fluid. In the experiment, particle images were recorded by a CCD camera at 30 frames per second (fps). Tracer particle motions in the same mesh in each frame were calculated by PIV lab using MATLAB to obtain 2-D velocity distributions in the plane of the laser light sheet. Two-dimensional velocity was visualized and calculated by the PIV software in the middle part of the cavity between the two electrodes. The measurement conditions for the PIV method are summarized in Table 2. Because the laser light was reflected by the sloping shape, the flow near the cavity wall could be observed by the PIV method. The width of the test window was 50 mm and the width of the mesh size was 2.5 mm to ensure the accuracy of the calculated velocity.
2.3. Experimental Methods
In this experiment, the room temperature was kept constant at around 20˚C. The initial temperature of the fluid was also 20˚C. The cooling temperature of the heat sink was 20˚C by using a copper heat sink at the top surface and a water circulator. The experiment was started by applying a voltage of 95 V between the electrodes, and it led to Joule-heating induced flow of the working fluid (a 80 wt% glycerol-water solution) in the sloping bottom cavity. The experimental conditions are shown in Table 3. After a sufficient time from the start of heating, the temperature tends to reach a steady state. In the cavity, the temperature became stable at about 6000 s after Joule-heating was started. Therefore, flow measurements by UVP and PIV methods were started 7200 s after the experiment was begun.
Table 1. Measurement conditions for the UVP method.
Table 2. Measurement conditions for the PIV method.
Table 3. Experimental conditions.
3. Numerical Simulations Using GSMAC
3.1. Governing Equations
Governing equations in the original GSMAC are shown below, derived from the formulas for the velocity field, electromagnetic field and temperature field. The velocity field formula consists of the continuity and Naiver-Stokes equations. The electromagnetic field formula is derived from the conservation laws of current and Ohm’s law. The temperature field formula is used in the energy equation. GSMAC executes fully-coupled analyses of the flow field, thermal field, and electromagnetic field. More details of the calculation method in GSMAC are given in , and .
Equation of continuity:
Thermal energy and electric field equations:
The revision of Ampere’s rule is derived from the introduction of the A − φ method. The magnetic vector potential A is introduced with the definition as:
From the Gauss’s law and Faraday’s law, the fundamental equations for the electromagnetic field are given by:
And then using Equation (15),
Equation (16) can be written as the gradient of some scalar function. Thus, a scalar potential can be calculated as:
The electric current density is then derived as:
Using the Coulomb gauge , the electric current density is written as follows:
Using Equation (20) and Equation (12), GSMAC-FEM solves for the electric current distribution.
3.2. Computational Method
As a numerical method for incompressible viscous flow, the Generalized-Simplified Marker and Cell Finite Element Method (GSMAC-FEM) are applied. The energy equation and Ampere-Maxwell equation are computed by the node based FEM and edge based FEM, respectively. Solvers of three fields are coupled through a staggered scheme. The coupling method enables to analyze natural convection under Joule-heating induced by an electromagnetic field.
The governing equations for the velocity field are solved by the GSMAC-FEM. The GSMAC-FEM is extended from Highly Simplified Marker and Cell (HSMAC). The velocity and pressure are calculated at a predictor step and a simultaneous relaxation step. The Naiver-Stokes equations are discretized semi-implicitly, where the velocity and the pressure are discretized explicitly and implicitly, respectively, as follows:
where a forward difference is used for the time derivative, n denotes the n-th time step and t is the time increment. Using the fractional step method, Equation (22) is divided into two parts:
where is the predictor, and is the modified pressure. Equation (14) is solved at the predictor step. When Equation (4) is satisfied at the next time step, the following Poisson equation is solved with Equation (24) at the simultaneous relaxation step:
The thermal energy equation is explicitly discretized as follows:
The above-mentioned governing equations have been simplified by the following assumptions:
1) The fluid is Newtonian;
2) The Boussinesq approximation may be applied to the buoyancy force term in Equation (22);
3) The flow is laminar and incompressible;
4) The displacement current is neglected.
3.3. Numerical Model
The analytical model of Joule-heated cavity is shown in Figure 5. Three-dimensional numerical analyses of GSMAC were executed using this model. The blue area contains the analysis mesh for the thermoelectrically conducting field and the white area is the analysis mesh for the electromagnetic field. The dimensions of this fluid domain are 126 × 152 × 154 mm. This volume is modeled by hexahedral elements with 27,000 meshes. The top surface is modeled as an isothermal surface and the rest of the cavity surfaces are treated as adiabatic surfaces. At the initial stage, the liquid is at rest with a constant temperature of 20˚C. The electrodes are placed on the upper side of opposite side walls (Y = 76 mm - 152 mm) of the fluid domain and set under an adiabatic condition. The voltage between the electrodes was 95 V, which was the same as in the experiments. The electromagnetic field is modeled to cover the fluid domain which is three times larger than the analysis mesh for the thermoelectrically conducting field. The dimensions of the electromagnetic domain are 378 × 456 × 462 mm and modeled by 48,000 hexahedral elements. The outer boundary of the MHD region satisfies electro-magnetic condition as follows:
where n is the outward unit vector. Equation (27) represents the boundary condition for the magnetic field calculation.
4. Results and Discussion
4.1. Experimental Results
The electrodes were set at the upper cuboid part of the cavity and induced Joule-heating. Figure 6 shows the 2-D flow behavior at 7200 s after the experiment started between the electrode surface and the center of the sloping bottom cavity. In the upper part of the cavity, a vortex was formed by Joule-heating. There was only a small flow in the lower sloping part as the fluid in the bottom part was not heated by Joule-heating. The temperature of the fluid in this sloping part of the cavity was lower than that of the fluid in the upper part. Therefore, a chaotic flow occurred just in the upper part, where the fluid was heated by Joule-heating.
Figure 7 depicts the color graphs of the velocity obtained by the UVP method. In these profiles, 0 - 76 mm was the lower part of the cavity and 76 - 152 mm was the upper part. The figure shows that the downflow only existed in the upper part of the cavity and did not reach the lower sloping part. There was only a small flow in the lower sloping part of the cavity.
Figure 8 shows a comparison of the 30 s average velocity profiles between the UVP and PIV methods. Because the bottom of the cavity is a non-flow area and the Joule-heating induced flow only occurred in the upper part of the cavity, a comparison between the UVP and PIV results could be made only in the upper part of the cavity, where Y = 80 mm - 152 mm. In addition, for comparisons of the PIV and UVP results, the laser light sheet set at the bottom of the cavity is shown in Figure 9. Therefore, the interrogation area of the PIV method decreased to 10 × 75 mm due to the reflection of the light by sloping wall. Measurement results from UVP and PIV were in good agreement. However, between Y = 130 mm and Y = 140 mm, the velocities measured by the two methods were different by about 0.5 mm/s. This difference could be due to the difference in the measurement positions. Both the laser sheet and the ultrasound sensor were set at the bottom of the cavity. However, if the laser sheet and ultrasound sensor were set at the same position, the ultrasound sensor would hide the laser sheet, so the laser sheet and the ultrasound sensor were offset by a little distance (2.5 mm). Accordingly, it is revealed that the UVP method could measure the Joule-heating induced flow.
4.2. Numerical Results
Figure 10 depicts the snapshot of the flow behavior in the sloping bottom cavity calculated by the GSMAC code. The flow was unstable in the upper part of the cavity where Joule-heating occurred and the downflow did not reach the lower sloping part. The fluid in the bottom part was not heated by Joule-heating. The flow behavior calculated by the GSMAC code was similar to that observed in the experiment.
To investigate the chaotic flow behavior in the cavity, the long-term flow behavior is required because the flow changes at every position in the cavity. In order to observe the long-term flow behavior using a numerical analysis, the vertical component of the velocity profile was calculated at the same position of the experimental configuration (Y = 78 mm) and plotted in color. Figure 11 shows a color plot of the vertical velocity magnitude as in Figure 7. The chaotic flow could be observed in the upper part of the cavity, which was also observed in Figure10 at t = 2000 s. As a result, GSMAC could reproduce the chaotic flow behavior in the cavity with Joule-heating. However, the Joule-heating area calculated by the GSMAC was a little larger than that measured by the UVP method. In the experiment, the flow occurred around Y = 80 - 152 mm, but under the simulation, the flow can be observed from Y = 60 mm to 152 mm. There was some discrepancy between the simulation and experimental results.
4.3. Comparison of Numerical Results with Experimental Data
Figure 12 shows the time series change in the temperature at the center of the cavity in the simulation and experimental results. Fluctuations in temperature were induced by the downward flows from the top of the cavity which was cooled and maintained at a constant temperature of 20˚C. However, the frequency and magnitude of the fluctuations were not coincident with each other between the simulation and the experiment. The temperature shown after 3000 s in Figure 12 was 38˚C in the simulation result, whereas it was 41˚C in the experimental result. In addition, the amplitude of temperature fluctuations in the experiment was larger than that in the simulation. It could be considered to be due to the effect of the electromagnetic field in the cavity. Because AC power was applied in the cavity to generate Joule-heating, an electromagnetic field was generated in the cavity. As the temperature was measured by a K-type thermocouple in the experiment, the magnetic field could have affected the response of the thermocouple resulting in the measured temperature being a little higher than the actual temperature. Another reason could be due to the wall of the cavity in the experiment not absolutely adiabatic. Therefore, to keep the balance of the temperature field, the center temperature in the cavity would be higher than that calculated by the numerical method.
A frequency analysis could be done with the velocity data obtained from the experiments and numerical simulations for validation purposes. The frequency spectra or Power Spectrum Density (PSD) was calculated using the Fast Fourier Transform (FFT) method on the time-series data. Figure 13 shows the PSD of
Figure 5. Calculation model and mesh layout for numerical analysis.
Figure 6. Velocity vectors in the cavity.
Figure 7. Color graphs showing vertical velocity in the cavity obtained by the UVP method.
Figure 8. Comparison of velocity profile.
Figure 9. Laser light sheet setting for comparison of the UVP and PIV results.
Figure 10. Simulation result of velocity field. (a) t = 500 s; (b) t = 1000 s; (c) t = 2000 s; (d) t = 5000 s.
Figure 11. Vertical component of velocity profile using GSMAC model.
Figure 12. History of temperature from experimental data and simulation result measured at center of the cavity using GSMAC-FEM.
Figure 13. PSD analysis of velocity from the experimental and numerical results.
the velocity data at the same position of z = 125 mm from the experimental results and the numerical simulations. They show the same slope in both PSDs of experiment and simulation results. The slopes range from around 10−4 Hz to 10−2 Hz and have almost the same inclination. The PSD of velocity describes the flow regime. The flow is turbulent when the PSD shows no clear peaks in the oscillation frequency . Thus, the simulation results are in qualitative agreement with the experimental results.
The chaotic Joule-heating induced flow behavior was examined by using the GSMAC-FEM method. The computational results have elucidated the flow field in the sloping bottom cavity, and the predicted flow regime in the cavity was similar to that observed in the experiments. The PSD results calculated from the velocity in the vertical direction show a good agreement between the experiments and simulations.
Laboratory-scale experiments and numerical simulations were conducted to study the flow behavior in a glass melter used in the HLLW vitrification process where multiphase flow occurs involving many phenomena. To clarify the flow behavior, a simplified model cavity was used with a glycerol-water solution as the working fluid.
To understand the flow behavior in the glass melter more clearly, a model cavity was designed to be closer to the design of an actual glass melter than that used in our previous study. The new cavity was named a sloping bottom cavity and had a different Joule-heating area and bottom shape from those of the previously studied cubic cavity.
In the model cavity, a pair of electrodes was placed on opposite side walls of the cavity in just the upper part where the Joule-heating area was, and two sloping surfaces were set in a non-Joule-heating area. The flow behavior was observed by the UVP and PIV methods in the experiment and calculated by the GSMAC method in the simulation. Conclusions are given as follows.
1) In the sloping bottom cavity, a chaotic flow occurred just in the upper part and only a small flow was in the sloping bottom part. Another result revealed that UVP technique could be applicable to the observation of the chaotic Joule-heating induced flow. The applicability of the UVP measurement was verified by means of the PIV method.
2) The GSMAC calculations were conducted for the same condition of the experiment. Accordingly, the Joule-heating induced flow as unstable flow appearing in the upper part of the cavity was evaluated. Qualitative agreement was obtained between the simulation results and experimental data.
3) Comparing the experimental data and the numerical results, the Joule-heating induced flow area and temperature in the center showed some discrepancies because the cavity wall in the experiment may not have been perfectly adiabatic. However, the PSD’s calculated from the velocity data showed the same slopes in both the experiment and simulation. Simulation results qualitatively agreed with the experimental results.
A: Magnetic Vector potential [T m]
B: Magnetic flux density [T]
Cv: Specific heat at constant volume [J/kg K]
d: Diameter [m]
D: Deformation velocity tensor [1/s]
Da: Damkohler number [-]
E: Electric field [V/m]
g: Gravitational acceleration [m/s2]
Gr: Grashof number [-]
H: Magnetic field intensity [A/m]
I: Identity tensor [-]
Jc: Induced current density [A/m2]
Jf: Coil current density [A/m2]
L: Characteristic length [m]
n: Unity vector [-]for magnetic field
p: Pressure [Pa]
Pr: Prandlt number [-]
Q: Volumetric heat density [W/m3]
q: Heat flux [W/m2]
Ra: Rayleigh number [-]
Rai: Internal Rayleigh number [-]
Tem: Magnetic stress tensor [1/s]
T: Temperature [K]
T0: Reference temperature [K]
v: Velocity [m/s]
v: Velocity vector [m/s]
x: Channel of UVP measurement [-]
X: Distance to the UVP transducer [m]
Y: Height of the sloping bottom cavity [mm]
α: Thermal diffusivity [m2/s]
β: Thermal expansion coefficient [1/K]
k: Thermal conductivity [W/m K]
μ: Coefficient of viscosity [Pa s]
μm: Magnetic permeability [H/m]
ν: Kinematic viscosity [m2/s]
ρ: Density [kg/m3]
σe: Electrical conductivity[S/m]
φ: Electrical potential [V]
 Tsuzuki, N., Kikura, H., Kawai, H., Kawaguchi, T., Inagaki, A. and Ochi, E. (2011) The Numerical Analysis on Unsteady Flow in the Cavity with Joule-Heating. 48th Heat Transfer Symposium of Japan, Vol. 3, Okayama, June 2011, 727-728.
 Duong, T.T., Tanaka, H., Tsuzuki, N., Kawai, H. and Kikura, H. (2015) Effect of Cooling Temperature of Electrodes on Joule-Heating Flow in Cubic Cavity. Progress in Nuclear Energy, 82, 165-175.
 Duong, T.T., Tanaka, H., Tsuzuki, N., Kawai, H. and Kikura, H. (2014) Experimental and Numerical Investigation of Joule-Heating Flow in a Square Cavity Effect of Cold Cap Condition. International Topical Meeting on Advances in Thermal Hydraulics 2014 (ATH14), Reno, NV, June 2014, 1-6.
 Igarashi, H., Kawamura, K. and Takahashi, T. (1991) Effect of Noble Metal Elements on Viscosity and Electrical Resistivity of Simulated Vitrified Products for High-Level Liquid Waste. MRS Online Proceedings Library, Vol. 257, 169.
 Ihara, T., Tsuzuki, N. and Kikura, H. (2016) Application of Ultrasonic Doppler Velocimetry to Molten Glass by Using Broadband Phase Difference Method. Flow Measurement and Instrumentation, 48, 90-96.
 Matsuno, S., Iso, Y., Uchida, H., Oono, I., Fukui, T. and Ooba, T. (2008) CFD Modeling Coupled with Electric Field Analysis for Joule-Heated Glass Melters. Journal of Power and Energy Systems, 2, 447-455.
 Tanaka, H., Duong, T.T., Tsuzuki, N. and Kikura, H. (2013) Flow Measurement in an Electrode Cooling Joule-Heating Cavity by UVP Method. 41st Visualization Symposium of Japan, Tokyo, July 2013, 415-418.
 Tanahashi, T., Okanaga, H. and Saito, T. (1990) GSMAC Finite Element Method for Unsteady Incompressible Naiver-Stokes Equation at High Reynolds Number. International Journal for Numerical Method, 11, 479-499.
 Kohno, H. and Tanahashi, T. (2005) An Application of GSMAC-FEM to Coupled Natural and Marangoni Convection in a Square Cavity. International Journal of Computational Fluid Dynamics, 19, 329-335.
 Bathchelor, G.K. (1959) Small-Scale Variation of Convected Quantities like Temperature in Turbulent Fluid Part 1. General Discussion and the Case of Small Conductivity. Journal of Fluid Mechanics, 5, 113-133.