- Email: [email protected]
Experimental studies and CFD calculations for buoyancy driven mixing phenomena
Nuclear Engineering and Design 240 (2010) 2185–2193
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
Experimental studies and CFD calculations for buoyancy driven mixing phenomena Marco Jose da Silva ∗ , Sebastian Thiele, Thomas Höhne, Roman Vaibar, Uwe Hampel Institute of Safety Research, Forschungszentrum Dresden-Rossendorf e.V., Dresden, Germany
a r t i c l e
i n f o
Article history: Received 20 January 2009 Received in revised form 3 August 2009 Accepted 7 August 2009
a b s t r a c t In nuclear reactor safety the mixing of borated and deborated water is a critical issue that needs investigation, assessment and prediction. Such mixing is buoyancy driven and numerical codes must correctly model momentum transfer between fluids of different density. To assess and develop CFD models for buoyancy driven mixing we set up a simple vertical mixing test facility (VeMix) and equipped it with a newly developed planar electrical imaging sensor. This imaging sensor acquires conductivity images of the liquid at the rear channel wall with a speed of 2,500 frames/s. By adding NaCl tracer to the denser fluid we were able to visualize the mixing process in high spatial and temporal detail. Furthermore, an image processing algorithm based on the optical flow concept was implemented and tested which allows the measurement of flow pattern velocities. Selected experiments at different Richardson numbers were run with two components of different density (pure water and glucose–water mixture) simulating borated and deborated water in a light water reactor scenario. These experiments were compared to CFD calculations using standard turbulence models. Good agreement between experimental data and CFD simulations was found. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Buoyancy driven flow is found in many engineering applications such as mixing of fluids with different densities. In nuclear reactors, fluid mixing induced by the buoyancy effect is, amongst others, relevant for • boron dilution issues, when highly borated water with different temperature mixes with the ambient water in the reactor pressure vessel, and • pressurized thermal shock scenarios, when cold emergency cooling water is injected into the reactor and is contacted with the vessel wall. Especially the problems related to the mixing of borated and unborated water have recently been in the scope of investigations (Umminger et al., 2001; Rohde et al., 2005; Höhne et al., 2006; Kliem et al., 2008a,b). In certain situations such as when boron concentration in the reactor drops below a critical level, additional borated water is added to the core via an emergency injection system. The additional borated water usually has a lower temperature than the
∗ Corresponding author. E-mail address: [email protected] (M.J. da Silva). 0029-5493/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2009.11.023
ambient coolant. Both boron concentration and temperature affect the density difference of the fluids and with that the intensity of the mixing processes. Since the mixing of weakly and highly borated coolant is a critical issue with respect to the reactivity insertion into the reactor core, its understanding and modeling is essential for the design, operation and safety analysis of nuclear reactors but also for other processes in industry. In this article we present a test facility which is used for the study of the mixing of fluids with different densities. The VeMix facility is a simple vertical flow channel equipped with a new imaging device—a planar conductivity imaging sensor, which allows us to study mixing with high spatial and temporal detail. The generic geometry thereby eases model development and simulation but also flow measurement and therefore helps to gain an improved understanding of the fundamentals of mixing governed by buoyancy. We further implemented the optical flow method for the estimation of velocity fields in the conductivity images obtained with the planar sensor. Accuracy of velocity measurements was evaluated by a simple experiment beforehand. We further present experimental and simulation results for selected mixing experiments. For a first flow simulation with standard turbulence models we used the commercial CFD code CFX-11. Qualitatively there is good agreement between experiment and simulation. Eventually, it is believed that in the future CFD with improved turbulence models can tackle complex geometries, such as a nuclear reactor vessel.
2186
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
Fig. 1. Schematic diagram (a) and photograph (b) of the VeMix test facility used in the buoyancy driven flow experiment.
2. Experimental setup and CFD model
2.2. CFD model
2.1. VeMix test facility
As a CFD code for simulating the mixing studies we used CFX-11 (CFX11, 2007). CFX-11 is an element-based finite-volume method with second-order discretization schemes in space and time. It uses a coupled algebraic multigrid algorithm to solve the linear systems arising from discretization. The discretization schemes and the multigrid solver are scalably parallelized. CFX-11 works with unstructured hybrid grids consisting of tetrahedral, hexahedral, prism and pyramid elements. An extensive grid study was done and as a result of it, a uniform hexahedral node distribution was selected for the simulations, whereby the resulting grid has 420,000 nodes and 346,000 elements. The following parameters were used in the simulations. The SST buoyancy turbulence model was used. The SST model (Menter, 1993) works by solving a turbulence/frequency-based model (k–ω) at the wall and standard k–ω in the bulk flow. A blending function ensures a smooth transition between the two models. All transient calculations were done with the higher order scheme “high resolution” in space and “fully implicit 2nd order backward Euler” in time. For both discretization schemes the target variable does not change significantly for convergence criteria below 1 × 10−4 . Therefore this convergence criterion is used for all calculations. Double precision was used. The inlet boundary conditions were set depending on the chosen flow rate, i.e.:
Fig. 1 shows the design of the vertical mixing test facility (VeMix) which is essentially a simple upright channel of rectangular crosssection. It consists of five identical rectangular segments stacked together on top of each other. Each segment is made of acrylic glass and is 500 mm wide, 625 mm high, and 100 mm deep. Altogether, the test section is 3.32 m tall. A flow inlet is located in between segments 1 and 2, while an outlet is positioned in between segments 4 and 5. Segment number 2, right after the inlet, is instrumented with the planar sensor (see Section 3). This region is of major interest with respect to mixing. The facility ensures that the gravitational force acts along a significant length scale. For the experiments we used pure deionized water as the lighter component and a 6% water–glucose mixture as the heavier component. The latter has a density of 1020 kg/m3 , i.e. 2% density difference compared to pure water. The execution of one single experiment consists of the following steps:
1. Water is filled into the test section up to the bottom of segment 2. 2. Water–glucose solution is very slowly inserted from the bottom taking care not to premix both components. In this way, a layered structure is obtained with pure water in segments 1 to 2 and the denser water–glucose solution in segments 3 to 5. The formation of a well-defined boundary between both components is important for later comparison with CFD simulations, where this known initial distribution serves as an initial condition. 3. Finally, the mixing process is initiated by injecting water–glucose solution through the inlet with a constant flow rate.
• velocityIn = Q/Ain • velocityOut = Q/Aout where Q is the flow rate (constant in time) and Ain = Aout = 0.005 m2 the cross-section area of inlet and outlet sections. Furthermore, for all simulations and experiments water with w = 997 kg/m3 (T = 25 ◦ C) and water–glucose solution with a density of
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
2187
Fig. 2. Flow regime during the flow condition = 2%, L = 0.83 m and flow rate is indicated in each figure. The water–glucose solution was colored with a dye tracer: (a) falling down regime: flow rate = 0.4 l/s, Ri = 2.55 and (b) horizontal jet regime: flow rate = 0.7 l/s, Ri = 0.83.
wgs = 1020 kg/m3 were used, i.e. with a density difference of = 2%. The simulation domain was filled up to 1.8 m, i.e. up to bottom of segment 2, with heavy liquid in order to meet the experimental conditions as described in Section 2.1. Transient experimental values of velocity were used. No specific velocity profile is given. As an initial guess of the turbulent kinetic energy and the dissipation rate the code standard is used (turbulent intensity = 5%, eddy viscosity ratio = 10). The outlet boundary condition was an opening using outflow velocities. The multicomponent buoyant flow option using different fluid densities were used to describe the mixing processes (Höhne et al., 2006). Transient calculations with 300 s simulation time with time steps: • t = 0.1 s, for 0 ≤ t ≤ 30, and • t = 0.5 s for 30 < t ≤ 300 were performed on the FZD LINUX cluster. The cluster runs at Linux Scientific 64 bit operational system and comprises 32 AMD OPTERON Computer Nodes. Each computer node is formed by 2 × AMD OPTERON 285 (2.6 GHz, dual-core) with 16 GB memory. Simulations took approximately 3 days each to complete. Four processors were used for the simulations presented here in a parallel mode with the message passing protocol PVM. 2.3. Flow characterization First validations of the simulated flow were performed by visual inspection of the mixing behavior using a video camera. For this purpose, water–glucose solution was colored with a blue dye tracer. Fig. 2 shows snap shots of both simulation results and images taken with the camera for the indicated conditions which show different flow regimes with typical vertical and horizontal flow patterns. The observed flow behavior can be characterized by an internal Richardson number, a non-dimensional number which is given by
the ratio between the buoyancy and the momentum forces: Ri =
gL w v2
(1)
Here, g = 9.81 m/s2 is the acceleration due to gravity, L = Hinlet − Hfill is the characteristic length scale, which is the height of the pure water layer above the water–glucose mixture (Hinlet = 2.63 m and Hfill = 1.8 m for the VeMix facility), is the density difference between the two components, w is the density of pure water and v the characteristic velocity of the VeMix test facility, which describes the movement of the interface between the liquids and is determined by
v=
Q A
(2)
Q denotes the flow rate and A = 0.05 m2 is the cross-section area of the VeMix test section. Other characteristic length and velocity scales could be applied in defining the Richardson number. For instance, the hydraulic inlet area in (2) for the definition of the characteristic velocity could also be used. However, the described scales were chosen because they enhance the sensitivity of the obtained Richardson numbers regarding to the flow pattern transition of regimes and thus better describe the mixing pattern developed in the facility. The Froude number has been widely used to characterize buoyant mixing which is inversely proportional to the Richardson number (Ri = Fr−2 ). However, we have chosen the Ri number in this study because it has a greater sensitivity to describe the transitions of characteristic flow patterns occurring at VeMix facility. Basically two different flow regimes can be identified from the simulation and experiments and are properly characterized by the Richardson number: • For Ri > 2 a regime with a falling jet pattern is observed (Fig. 2a). The gravity force acts along the vertical axis and is significantly stronger than momentum forces of the inlet flow. The jet drops near the wall below the inlet nozzle.
2188
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
Fig. 3. Planar sensor with 64 × 64 interdigital sensing structures fitted to a segment of VeMix test facility. The size of each sensing structure is 7.8 mm × 6.6 mm and the size of the whole sensor is 620 mm by 500 mm.
• For Ri < 1 a regime with a horizontal jet pattern is observed (Fig. 2b). The momentum force from the inlet flow is stronger than the gravitational force, which results in the formation of a horizontal jet. The horizontal jet influences the transient mixing behavior, as oscillating waves are formed on interface between the heavier and the lighter component. For intermediate values between 1 and 2 a transition regime occurs. The gravity and momentum forces are in near equilibrium and the resulting flow pattern strongly depends on the force ratio. In this regime it is not easy to predict the flow pattern. 3. Planar array sensor and velocity field estimation 3.1. Planar array sensor The most interesting location in the VeMix vessel was found to be just below the inlet nozzle. Therefore, the planar sensor was placed here in order to study the mixing behavior of the two liquids. The planar array sensor used to investigate the mixing behavior can be understood as a fast conductivity camera, realized as a pixellated
planar sensor surface aligned with the rear wall of VeMix. The sensor was manufactured using standard printed-circuit board (PCB) fabrication technology. It was built in a rectangular geometry with 4096 interdigital sensing structures evenly distributed over an area of 620 mm × 500 mm (Fig. 3). The sensor board fits exactly into a segment of the VeMix flow channel (see Section 2). Each sensing structure of the sensor consists of a small transmitter electrode and a small receiver electrode, opposing each other in an interlocking way (Fig. 3). The transmitter electrodes are connected row-wise whereas the receiver electrodes are connected column-wise. Electrical conductivity at the interdigital sensing structures is consecutively, repeatedly and rapidly interrogated by a special electronics. The transmitter electrodes are therefore sequentially excited by a bipolar voltage signal while the non-active electrodes are grounded. The electrical currents flowing from the activated transmitter electrode to the receiver electrodes are measured synchronously while the transmitter electrodes are being activated consecutively. After all transmitter electrodes have been activated once a complete frame of the conductivity distribution at all sensor structures has been acquired. Fig. 4 depicts a schematic diagram of the electronics and the timing diagram of the electrical signals. The maximal frame rate achievable for this sensor is
Fig. 4. (a) Simplified block diagram of the measuring electronics. (b) Timing diagram for the multiplexed excitation and measuring scheme for the sensor. Voltages are sampled at the rising edge of the sample-and-hold signal. One frame is complete when the 64 transmitter electrodes have been activated. The signal SP denotes the polarity switch and S1–S64 the analog switches.
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
2189
2,500 images/s (Da Silva et al., 2007). It was demonstrated that the penetration depth of each sensing structure composing the whole sensor is of about 1 mm (Da Silva et al., 2007). That means, the sensor is able to visualize the flow only close to the wall it is attached. The electrical current flowing from the transmitter to the receiver electrode of a sensing structure depends on the electrical conductivity of the liquid surrounding the sensing structure. The current is converted into a voltage by amplifier electronics and the voltage value V therefore has general dependency on electrical conductivity in the form (Da Silva et al., 2007): V (i, j, k) = a(i, j) + b(i, j)
(3)
where a and b are slope and offsets describing electronics’ characteristics, i and j are the grid coordinates and k is the temporal index. To obtain an electrical contrast between pure water and water–glucose mixture one of the fluids needs to be marked with a conductivity tracer. This was done by adding a small amount of NaCl to the water–glucose mixture which increases electrical conductivity but only marginally changes density. It has been shown previously that the conductivity measurement can be linked to density gradients by (Höhne et al., 2006): (i, j, k) =
(i, j, k) − L ∼ (i, j, k) − L V (i, j, k) − VL (i, j) = = H − L H − L VH (i, j) − VL (i, j)
(4)
where is called mixing scalar with values between 0 and 1. The subscripts denote the heavy (H) and light (L) component. Thus, in order to obtain mixing scalar distributions (i,j,k) from the planar sensor, calibration measurements must be taken prior to experiments. This calibration involves the acquisition of voltage readings for the sensor fully covered with pure water VL (i,j) and for the sensor fully covered with NaCl–water–glucose solution VH (i,j).
As a first step, we implemented a differential optical flow algorithm to estimate the flow vector f(i,j) for each pixel i and j at time t of the 2D-image from the spatio-temporal derivatives of the gray values in the image sequence. For the present sensor, i and j vary from 1 to 64, since the sensor has 64 × 64 sensing structures. The algorithm is based on a concept similar to the continuity equation in fluid dynamics, where mass is conserved. Here, the gray value distribution g(i,j) in the image is assumed to be nearly constant for a short period, which yields the optical flow constraint equation (Jähne, 2002, pp. 387):
∇ g f + gt = 0
⎛
With increasing computer performance computation of velocity fields from image sequences have gained increasing attention. Optical flow is a prominent method (amongst other, such as cross-correlation algorithms) and is being widely used for motion tracking. Details of this method are found in Beauchemin and Barron (1995), Jähne (2002), as well as in Corpetti et al. (2004). The optical flow method has been developed especially for processing of images from optical acquisition systems. There, the projected motion field from the real 3D scene onto the image scene is extracted from the spatio-temporal change of the gray values in the corresponding image sequence which is defined as optical flow (Jähne, 2002). The optical flow method evaluates the change of brightness (or intensity) in a sequence of images caused, for example, by the movement of an object, and associates these intensity gradients with a velocity field. The planar array sensor is evidently no optical acquisition system and the gray values in the images represent an electrical parameter, i.e. the measured conductivity value at each electrode pair, which is represented by a gray value of the corresponding pixel in the image. However, the measured 2D conductivity value distribution represents the projected flow scene on the sensor surface in a given experiment and thus carries information of the near-wall velocity field. That is, the velocity field can be calculated from the spatio-temporal gradients of the gray values. The velocity vector v(x, y) can be inferred from the 2D optical flow vector f(i,j) according to
kg f (i, j) = fi /fj t
3.3. Optical flow for planar sensor data
(6)
where g is the local gradient of gray value, given by
3.2. Optical flow
v(x, y) =
Fig. 5. Illustration of the aperture problem in optical flow. The displacement vector of the edge in (a) cannot be determined uniquely, while for the movement of the corner and (b) it is possible. Adapted from Jähne (2002).
(5)
where kg is an scalar value which describes the geometry factor of the system and t is the time between two consecutive images.
∇g = ⎝
∂g ∂x ∂g ∂y
⎞
⎛
⎠≈⎝
g x g y
⎞
⎠=
gi gj
(7)
Further, gt is the discrete temporal derivative of the local gray value: ∂g g ≈ = gt t ∂t
(8)
However, the spatio-temporal derivatives are local operators, which act only in small regions of the image. This can be illustrated as a mask or aperture applied to the image (see Fig. 5). The solid line in Fig. 5a corresponds to an edge that is moving to the position of the dashed line in the second image. Due to the limited view of the local operators, the displacement vector of the optical flow may be connected from one point on the solid line to any point on the dashed line, i.e. there is no unambiguous correspondence. Thus, only the normal component f⊥ of the displacement vector can be correctly determined, whereas the parallel f component remains unknown. This phenomenon is well known as aperture problem (Jähne, 2002). Thus, the optical flow constraint equation (6) is nonunique. From the local operators, the optical flow vector can only be exactly determined when a corner of the object falls within the mask of the local operators as depicted in Fig. 5b. Due to the aperture problem, further conditions have to be taken into account to determine the vector of the optical flow. We have applied a global smoothness constraint first introduced by Horn and Schunck (1981) for regularization of the optical flow constraint equation. The global smoothness constraint assumes that optical flow vectors of neighboring pixels are nearly identical, i.e. the difference between the local average value of the optical flow f¯ and the optical flow of the neighboring pixels is assumed to be minimal.
2190
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
Fig. 6. Experimental setup for the validation of the optical flow algorithm.
This hypothesis yields the quadratic difference: ı2f = (f − f )
2
process is repeated until the difference between two iterations is below a given value in the form: (9)
Combining equation (9) with the uncertainty in the constraint equation resulting from signal noise: ın = ∇ g f + gt
(10)
yields the overall uncertainty constraint: ı2 = ˛2 ı2f + ı2n
(11)
which has to be minimized to achieve an estimation of the optical flow vector field. Thereby, ˛ is a fixed weighting factor being estimated empirically to adjust the impact of both terms of uncertainty. In our calculations the value of 0.1 was used. Minimization of ı2 is achieved by demanding ∂ı2 ! =0 ∂fi
(12)
n+1 n
f −f ≤ε
(16)
which is the norm of the difference between two consecutive iterations. In the calculations we empirically chose ε as 10−3 . In this way, the progress of the algorithm is adapted to the individual composition of the local spatio-temporal derivatives. In fact, the calculation of the local spatio-temporal derivatives in the standard algorithm by Horn and Schunck (1981) achieves sufficient accuracy in case of linear gray value distributions and oversampled image data only (Beauchemin and Barron, 1995). To achieve improved performance of the algorithm we implemented a method that evaluates five consecutive images for estimation of spatio-temporal derivatives for each step in time, which are then applied in the iterative calculation of the current optical flow vector field according to Eqs. (14) and (15) (Simoncell, 1994; Thiele, 2007). Once the optical flow is calculated the velocity can be easily achieved by implying the spatial dimensions of the sensing elements and the data acquisition rate as described in Eq. (5).
and 4. Validation experiment for velocity measurement
∂ı2 ! =0 ∂fj
(13)
from which one obtains two equations to be solved for each individual component of f. Since analytical solution of these both equations for all pixels of the image is highly computational intensive, an iterative set of two equations is derived from (12) and (13) solved with the Gauss–Seidel scheme (Horn and Schunck, 1981): n
f n+1 = f¯ i − gi i
n
n
gi f¯ i + gj f¯ j + gt ˛2 + gi + gj
(14)
and n
f n+1 = f¯ j − gj j
n n gi f¯ i + gj f¯ j + gt
˛2 + gi + gj
(15)
Eqs. (14) and (15) are solved iteratively for every step in time t and the optical flow in the iteration (n + 1) is estimated by the spatiotemporal derivatives and the previous locally averaged optical flow n f¯ (Horn and Schunck, 1981). The iterative algorithm for each step in time t is briefly described below. Initially, the spatio-temporal derivatives are calculated for every pixel and the averaged values f¯i0 and f¯j0 are set to zero. The first iteration is carried out yielding the optical flow components fi1 and fj1 . The neighboring pixels are locally averaged with a 2D averaging filter for every pixel of the image to obtain f¯i1 and f¯j1 which are used in the next iteration. This
In order to verify the feasibility of the implemented algorithm for velocity measurement from planar sensor data, a simple validation experiment was performed. A smaller prototype array sensor was used in this experiment. The electrodes dimensions are accordingly smaller than in the VeMix sensor, but the results obtained in this experiment are representative for the 64 × 64 sensor. Fig. 6 depicts the experimental setup. It consists of the sensor mounted in a horizontal position and a brush which was mounted on a 200 mm long fixation arm and mechanically moved across the sensor by means of a stepper motor. The motor was driven at constant angular speed of ω = 5.934 rad/s giving a tangential velocity of vT = 1.116 m/s for the moving brush. The sensor and a shallow vessel filled with tap water were placed horizontally within the circular path of the brush. In this way, the wet brush periodically passed across the sensor surface at constant velocity in a circular path. The brush movement along the planar sensor surface was captured with the electronics at 10,000 fps. This velocity is possible due to the fact that the prototype sensor used has only 32 × 32 electrodes (Da Silva et al., 2008). The motion of the wet brush in the image sequence was analyzed with the optical flow algorithm described in Section 3.3. From calculated optical flow vectors, the velocity vectors were determined using Eq. (5) with the spatial electrode dimensions (as shown in Fig. 6) and the time step between two consecutive images t = 0.1 ms. Fig. 7 shows four snap shoots of the liquid distribution
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
2191
Table 1 Results of the validation experiment for velocity measurement. Calculation type
Measured velocity (m/s)
Deviation from theoretical vT (1.116 m/s)
Maximum value over whole sensor Mean value of velocity vectors at the moving front
1.138 1.101
−2.0% 1.3%
on the sensor surface where the brush motion can clearly be identified. Furthermore, the figure also displays images of the temporal gradient values (Eq. (8)) which clearly show the front of the moving brush. For accuracy assessment an image sequence of 450 single images representing one passage of the wet brush over the whole sensor surface was analyzed. In order to compare the theoretical value of vT = 1.116 m/s and a measured velocity value, two approaches have been chosen. First, the maximal value of the calculated velocity vectors over the whole sensor surface was taken to represent the brush velocity for each time step. Second, the analyzed velocity vectors were confined to the moving front as indicated in Fig. 7 and the mean value within this region was chosen to represent the brush velocity for each time step. From the two data series containing 450 data points, the average value was calculated to determine the brush velocity, called vmax and vfront to indicate
the prior and latter calculation method. Table 1 shows the obtained values. Both calculated values show good agreement with theoretical velocity value vT with deviations below 2%. The value vmax is slightly higher than vfront . This can be explained by the fact that, when calculating maximal values, some image noise is picked up by the algorithm. 5. Vemix experiment and simulation In this section, the experimental and simulation results for an exemplary flow condition are presented. The essential experimental parameters are Q = 0.5 l/s (flow rate), Hfill = 2.12 m (i.e. L = 0.51 m), = 2%. This gives a Richardson number of 1.004 from Eq. (1) and, with a horizontal jet pattern flow regime (Section 2.3). Simulations were performed with the CFD model explained in Section 2.2. The planar array sensor was set to acquire 200 images per second during a 100 s time interval. In the experiment, pure water with w = 997 kg/m3 and water–glucose solution with wgs = 1017 kg/m3 and electrical conductivity of 60 and 200 S/cm, respectively, were used. The velocity fields were computed by the optical flow method as presented in Section 3.3. As previously discussed, the measured data with the planar sensor represents the flow close to the wall, while from CFD simulations a threedimensional data set of the flow is obtained. Thus, only the CFD data of an axial plane located at 1 mm distance from the wall were extracted and compared with the planar sensor’s experimental results. 5.1. Fluid distribution In Fig. 8 images of the distribution of mixing scalar (Eq. (4)) from the measured and simulated results are depicted. As the water–glucose solution flows into the channel a circular jet forms both in the experiment and the simulation. It seems that in the experiment mixing happens more rapidly and therefore the jet is more pronounced in the simulated distributions. In the experiment, a transition layer at the bottom of the images can be observed. It is caused by the pre-mixing of water and water–glucose solution prior to the start of the experiment. That means, the initial condition in the experiment is slightly different to that used in the simulation. Nevertheless, the agreement between both methods is very good. The use of the Richardson number could also correctly predict the flow regime for this experiment. 5.2. Velocity fields
Fig. 7. Selected images of (a) the liquid distribution on the sensor surface and (b) the temporal gradient for the validation experiment of velocity measurement. In (b) the location of the moving brush can be clearly recognized.
Further analysis of the measured data was carried out by computation of the velocity field. Fig. 9 illustrates four selected images of measured and the simulated velocity fields. The magnitude of the velocity is indicated by the background color and velocity vectors indicating only the direction are shown in the foreground. The simulation images show a dominating vortex at the bottom of the jet which bends the initial straight jet and creates the circular liquid distribution as observed in Fig. 8. In the velocity field computed from the measured data the vortex is not as perfect but still visible. Comparison of velocity values for experiment and simulation turns out to be difficult. It was found that reliable velocity vectors from optical flow are only found in regions of high conductivity contrast, i.e. high spatial gradient. Due to micro-mixing
2192
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
Fig. 8. Distributions of mixing scalar at the rear wall of the VeMix test facility for the mixing experiment at Ri = 1.004: (a) planar array sensor measurement and (b) simulation with ANSYS CFX.
Fig. 9. Velocity fields at the rear wall of the VeMix test facility for the mixing experiment for Ri = 1.004: (a) velocity information calculated by the optical flow method from experimental sensor data and (b) velocity information extracted from CFD data.
the boundary between high and low density liquid is smeared out which decreases local contrast and thus leads to noise pickup by the optical flow algorithm. A good example is the image at t = 18 s where grayscale gradients are smaller and the algorithm accuracy tends to be lower. The jet almost disappears in this image and a velocity spot appears which is not physical. The spot occurs because the gray value gradient in that region is high. In spite of these limitations, the optical flow method is still able to evaluate the velocity field assuming the gradients in the grayscale are present. For the VeMix experiment described here, this occurs predominantly in the jet region and at the beginning of the experiment. Another finding is that the experimental results show much finer turbulence details than predicted by the applied SST model in the simulation. This indicates that modification of the turbulence model is necessary to improve agreement. In a qualitative sense the agreement in the values of the velocity magnitude is quite good, indicating that the
general features of the flow are already well captured by the CFD simulation. 6. Conclusions The VeMix test facility was constructed to investigate buoyancy driven flow with an emphasis on CFD code development and validation. In the facility, the mixing of weakly and highly borated coolant can be simulated. The use of a new developed planar array sensor allows the measurement of the mixing scalar with high spatial and temporal resolution which can be directly and quantitatively compared with density distributions predicted from CFD simulations. In this paper, we also introduced a new method to obtain velocity fields by image processing of the planar sensor images with the optical flow method. A simple experiment with a moving brush was carried out to validate this velocity measure-
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
ment method. A mixing experiment was performed and analyzed. Good agreement between experimental data and first CFD simulations with a simple SST turbulence model was found. The velocity measurement by optical flow method proved to be applicable in general but exposed some problems in regions with progressed micro-mixing and accordingly small conductivity gradients. Experimental data was found to disclose much more turbulence detail than the simulations which indicates that improvement in turbulence modeling is necessary. In the future, more experiments will be carried out in the VeMix test. Further work with respect to the new measurement modality will be directed towards validation of optical flow method by simultaneous measurements with PIV. Acknowledgement This work was in part supported by the German Federal Ministry of Economy and Labor, grant 150 1287. References Beauchemin, S.S., Barron, J.L., 1995. The computation of optical flow. ACM Computing Surveys 27, 433–467. CFX11, 2007. User Manual. ANSYS-CFX. Corpetti, T., Heitz, D., Arroyol, G., Memin, E., Santa-Cruz, A., 2004. Estimation of motion using a PIV correlation-based method and an ‘optical-flow’ one for two experimental flows: quantitative and qualitative comparison. In: Proceedings of 12th International Symposium on Applications of Laser Techniques to Fluid Mechanics, Lisbon, Portugal.
2193
Da Silva, M.J., Suhnel, T., Schleicher, E., Lucas, D., Vaibar, R., Hampel, U., 2007. Planar array sensor for high-speed component distribution imaging in fluid flow applications. Sensors 7, 2430–2445. Da Silva, M.J., Lu, Y., Suhnel, T., Schleicher, E., Thiele, S., Kernchen, R., Diele, K.-H., Hampel, U., 2008. Autonomous planar conductivity array sensor for fast liquid distribution imaging in a fluid coupling. Sensors and Actuators A 147, 508–515. Höhne, T., Kliem, S., Bieder, U., 2006. Modeling of a buoyancy-driven flow experiment at the ROCOM test facility using the CFD-codes CFX-5 and TRIO U. Nuclear Engineering and Design 236, 1309–1325. Horn, B.K.P., Schunck, B.G., 1981. Determining optical flow. Artificial Intelligence 17, 185–203. Jähne, B., 2002. Digital Image Processing, 5th revised and extended Edition. Springer. Kliem, S., Suhnel, T., Rohde, U., Hohne, T., Prasser, H.-M., Weiss, F.-P., 2008a. Experiments at the mixing test facility ROCOM for benchmarking of CFD codes. Nuclear Engineering and Design 238, 566–576. Kliem, S., Prasser, H.-M., Suhnel, T., Weiss, F.-P., Hansen, A., 2008b. Experimental determination of the boron concentration distribution in the primary circuit of a PWR after a postulated cold leg small break loss-of-coolant-accident with cold leg safety injection. Nuclear Engineering and Design 238, 1788–1801. Menter, F.R., 1993. Zonal two equation k–ω turbulence models for aerodynamic flows. In: Proc. of 24th AIAA Fluid Dynamics Conference. Orlando, FL, AIAA paper 93-2906, July 1993. Rohde, U., Kliem, S., Höhne, T., Karlsson, R., Hemström, B., Lillington, J., Toppila, T., Elter, J., Bezrukov, Y., 2005. Fluid mixing and flow distribution in the reactor circuit, measurement data base. Nuclear Engineering and Design 235, 421–443. Simoncell, E.P., 1994. Design of multi-dimensional derivative filters. In: Proc. of 1st IEEE International Conference on Image Processing, Austin, TX, USA, pp. 790–793. Thiele, S., 2007. Entwicklung und Aufbau eines kapazitiven Oberflächensensors für die Phasenverteilung- und Geschwindigkeitsmessung zweiphasiger Strömungen. Diploma Thesis, Technische Universität Dresden, Dresden, Germany (in German). Umminger, K., Kastner, W., Liebert, J., Mull, T., 2001. Thermal hydraulics of PWRS with respect to boron dilution phenomena. Experimental results from the test facilities PKL and UPTF. Nuclear Engineering and Design 204, 191–203.
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
Experimental studies and CFD calculations for buoyancy driven mixing phenomena Marco Jose da Silva ∗ , Sebastian Thiele, Thomas Höhne, Roman Vaibar, Uwe Hampel Institute of Safety Research, Forschungszentrum Dresden-Rossendorf e.V., Dresden, Germany
a r t i c l e
i n f o
Article history: Received 20 January 2009 Received in revised form 3 August 2009 Accepted 7 August 2009
a b s t r a c t In nuclear reactor safety the mixing of borated and deborated water is a critical issue that needs investigation, assessment and prediction. Such mixing is buoyancy driven and numerical codes must correctly model momentum transfer between fluids of different density. To assess and develop CFD models for buoyancy driven mixing we set up a simple vertical mixing test facility (VeMix) and equipped it with a newly developed planar electrical imaging sensor. This imaging sensor acquires conductivity images of the liquid at the rear channel wall with a speed of 2,500 frames/s. By adding NaCl tracer to the denser fluid we were able to visualize the mixing process in high spatial and temporal detail. Furthermore, an image processing algorithm based on the optical flow concept was implemented and tested which allows the measurement of flow pattern velocities. Selected experiments at different Richardson numbers were run with two components of different density (pure water and glucose–water mixture) simulating borated and deborated water in a light water reactor scenario. These experiments were compared to CFD calculations using standard turbulence models. Good agreement between experimental data and CFD simulations was found. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Buoyancy driven flow is found in many engineering applications such as mixing of fluids with different densities. In nuclear reactors, fluid mixing induced by the buoyancy effect is, amongst others, relevant for • boron dilution issues, when highly borated water with different temperature mixes with the ambient water in the reactor pressure vessel, and • pressurized thermal shock scenarios, when cold emergency cooling water is injected into the reactor and is contacted with the vessel wall. Especially the problems related to the mixing of borated and unborated water have recently been in the scope of investigations (Umminger et al., 2001; Rohde et al., 2005; Höhne et al., 2006; Kliem et al., 2008a,b). In certain situations such as when boron concentration in the reactor drops below a critical level, additional borated water is added to the core via an emergency injection system. The additional borated water usually has a lower temperature than the
∗ Corresponding author. E-mail address: [email protected] (M.J. da Silva). 0029-5493/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2009.11.023
ambient coolant. Both boron concentration and temperature affect the density difference of the fluids and with that the intensity of the mixing processes. Since the mixing of weakly and highly borated coolant is a critical issue with respect to the reactivity insertion into the reactor core, its understanding and modeling is essential for the design, operation and safety analysis of nuclear reactors but also for other processes in industry. In this article we present a test facility which is used for the study of the mixing of fluids with different densities. The VeMix facility is a simple vertical flow channel equipped with a new imaging device—a planar conductivity imaging sensor, which allows us to study mixing with high spatial and temporal detail. The generic geometry thereby eases model development and simulation but also flow measurement and therefore helps to gain an improved understanding of the fundamentals of mixing governed by buoyancy. We further implemented the optical flow method for the estimation of velocity fields in the conductivity images obtained with the planar sensor. Accuracy of velocity measurements was evaluated by a simple experiment beforehand. We further present experimental and simulation results for selected mixing experiments. For a first flow simulation with standard turbulence models we used the commercial CFD code CFX-11. Qualitatively there is good agreement between experiment and simulation. Eventually, it is believed that in the future CFD with improved turbulence models can tackle complex geometries, such as a nuclear reactor vessel.
2186
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
Fig. 1. Schematic diagram (a) and photograph (b) of the VeMix test facility used in the buoyancy driven flow experiment.
2. Experimental setup and CFD model
2.2. CFD model
2.1. VeMix test facility
As a CFD code for simulating the mixing studies we used CFX-11 (CFX11, 2007). CFX-11 is an element-based finite-volume method with second-order discretization schemes in space and time. It uses a coupled algebraic multigrid algorithm to solve the linear systems arising from discretization. The discretization schemes and the multigrid solver are scalably parallelized. CFX-11 works with unstructured hybrid grids consisting of tetrahedral, hexahedral, prism and pyramid elements. An extensive grid study was done and as a result of it, a uniform hexahedral node distribution was selected for the simulations, whereby the resulting grid has 420,000 nodes and 346,000 elements. The following parameters were used in the simulations. The SST buoyancy turbulence model was used. The SST model (Menter, 1993) works by solving a turbulence/frequency-based model (k–ω) at the wall and standard k–ω in the bulk flow. A blending function ensures a smooth transition between the two models. All transient calculations were done with the higher order scheme “high resolution” in space and “fully implicit 2nd order backward Euler” in time. For both discretization schemes the target variable does not change significantly for convergence criteria below 1 × 10−4 . Therefore this convergence criterion is used for all calculations. Double precision was used. The inlet boundary conditions were set depending on the chosen flow rate, i.e.:
Fig. 1 shows the design of the vertical mixing test facility (VeMix) which is essentially a simple upright channel of rectangular crosssection. It consists of five identical rectangular segments stacked together on top of each other. Each segment is made of acrylic glass and is 500 mm wide, 625 mm high, and 100 mm deep. Altogether, the test section is 3.32 m tall. A flow inlet is located in between segments 1 and 2, while an outlet is positioned in between segments 4 and 5. Segment number 2, right after the inlet, is instrumented with the planar sensor (see Section 3). This region is of major interest with respect to mixing. The facility ensures that the gravitational force acts along a significant length scale. For the experiments we used pure deionized water as the lighter component and a 6% water–glucose mixture as the heavier component. The latter has a density of 1020 kg/m3 , i.e. 2% density difference compared to pure water. The execution of one single experiment consists of the following steps:
1. Water is filled into the test section up to the bottom of segment 2. 2. Water–glucose solution is very slowly inserted from the bottom taking care not to premix both components. In this way, a layered structure is obtained with pure water in segments 1 to 2 and the denser water–glucose solution in segments 3 to 5. The formation of a well-defined boundary between both components is important for later comparison with CFD simulations, where this known initial distribution serves as an initial condition. 3. Finally, the mixing process is initiated by injecting water–glucose solution through the inlet with a constant flow rate.
• velocityIn = Q/Ain • velocityOut = Q/Aout where Q is the flow rate (constant in time) and Ain = Aout = 0.005 m2 the cross-section area of inlet and outlet sections. Furthermore, for all simulations and experiments water with w = 997 kg/m3 (T = 25 ◦ C) and water–glucose solution with a density of
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
2187
Fig. 2. Flow regime during the flow condition = 2%, L = 0.83 m and flow rate is indicated in each figure. The water–glucose solution was colored with a dye tracer: (a) falling down regime: flow rate = 0.4 l/s, Ri = 2.55 and (b) horizontal jet regime: flow rate = 0.7 l/s, Ri = 0.83.
wgs = 1020 kg/m3 were used, i.e. with a density difference of = 2%. The simulation domain was filled up to 1.8 m, i.e. up to bottom of segment 2, with heavy liquid in order to meet the experimental conditions as described in Section 2.1. Transient experimental values of velocity were used. No specific velocity profile is given. As an initial guess of the turbulent kinetic energy and the dissipation rate the code standard is used (turbulent intensity = 5%, eddy viscosity ratio = 10). The outlet boundary condition was an opening using outflow velocities. The multicomponent buoyant flow option using different fluid densities were used to describe the mixing processes (Höhne et al., 2006). Transient calculations with 300 s simulation time with time steps: • t = 0.1 s, for 0 ≤ t ≤ 30, and • t = 0.5 s for 30 < t ≤ 300 were performed on the FZD LINUX cluster. The cluster runs at Linux Scientific 64 bit operational system and comprises 32 AMD OPTERON Computer Nodes. Each computer node is formed by 2 × AMD OPTERON 285 (2.6 GHz, dual-core) with 16 GB memory. Simulations took approximately 3 days each to complete. Four processors were used for the simulations presented here in a parallel mode with the message passing protocol PVM. 2.3. Flow characterization First validations of the simulated flow were performed by visual inspection of the mixing behavior using a video camera. For this purpose, water–glucose solution was colored with a blue dye tracer. Fig. 2 shows snap shots of both simulation results and images taken with the camera for the indicated conditions which show different flow regimes with typical vertical and horizontal flow patterns. The observed flow behavior can be characterized by an internal Richardson number, a non-dimensional number which is given by
the ratio between the buoyancy and the momentum forces: Ri =
gL w v2
(1)
Here, g = 9.81 m/s2 is the acceleration due to gravity, L = Hinlet − Hfill is the characteristic length scale, which is the height of the pure water layer above the water–glucose mixture (Hinlet = 2.63 m and Hfill = 1.8 m for the VeMix facility), is the density difference between the two components, w is the density of pure water and v the characteristic velocity of the VeMix test facility, which describes the movement of the interface between the liquids and is determined by
v=
Q A
(2)
Q denotes the flow rate and A = 0.05 m2 is the cross-section area of the VeMix test section. Other characteristic length and velocity scales could be applied in defining the Richardson number. For instance, the hydraulic inlet area in (2) for the definition of the characteristic velocity could also be used. However, the described scales were chosen because they enhance the sensitivity of the obtained Richardson numbers regarding to the flow pattern transition of regimes and thus better describe the mixing pattern developed in the facility. The Froude number has been widely used to characterize buoyant mixing which is inversely proportional to the Richardson number (Ri = Fr−2 ). However, we have chosen the Ri number in this study because it has a greater sensitivity to describe the transitions of characteristic flow patterns occurring at VeMix facility. Basically two different flow regimes can be identified from the simulation and experiments and are properly characterized by the Richardson number: • For Ri > 2 a regime with a falling jet pattern is observed (Fig. 2a). The gravity force acts along the vertical axis and is significantly stronger than momentum forces of the inlet flow. The jet drops near the wall below the inlet nozzle.
2188
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
Fig. 3. Planar sensor with 64 × 64 interdigital sensing structures fitted to a segment of VeMix test facility. The size of each sensing structure is 7.8 mm × 6.6 mm and the size of the whole sensor is 620 mm by 500 mm.
• For Ri < 1 a regime with a horizontal jet pattern is observed (Fig. 2b). The momentum force from the inlet flow is stronger than the gravitational force, which results in the formation of a horizontal jet. The horizontal jet influences the transient mixing behavior, as oscillating waves are formed on interface between the heavier and the lighter component. For intermediate values between 1 and 2 a transition regime occurs. The gravity and momentum forces are in near equilibrium and the resulting flow pattern strongly depends on the force ratio. In this regime it is not easy to predict the flow pattern. 3. Planar array sensor and velocity field estimation 3.1. Planar array sensor The most interesting location in the VeMix vessel was found to be just below the inlet nozzle. Therefore, the planar sensor was placed here in order to study the mixing behavior of the two liquids. The planar array sensor used to investigate the mixing behavior can be understood as a fast conductivity camera, realized as a pixellated
planar sensor surface aligned with the rear wall of VeMix. The sensor was manufactured using standard printed-circuit board (PCB) fabrication technology. It was built in a rectangular geometry with 4096 interdigital sensing structures evenly distributed over an area of 620 mm × 500 mm (Fig. 3). The sensor board fits exactly into a segment of the VeMix flow channel (see Section 2). Each sensing structure of the sensor consists of a small transmitter electrode and a small receiver electrode, opposing each other in an interlocking way (Fig. 3). The transmitter electrodes are connected row-wise whereas the receiver electrodes are connected column-wise. Electrical conductivity at the interdigital sensing structures is consecutively, repeatedly and rapidly interrogated by a special electronics. The transmitter electrodes are therefore sequentially excited by a bipolar voltage signal while the non-active electrodes are grounded. The electrical currents flowing from the activated transmitter electrode to the receiver electrodes are measured synchronously while the transmitter electrodes are being activated consecutively. After all transmitter electrodes have been activated once a complete frame of the conductivity distribution at all sensor structures has been acquired. Fig. 4 depicts a schematic diagram of the electronics and the timing diagram of the electrical signals. The maximal frame rate achievable for this sensor is
Fig. 4. (a) Simplified block diagram of the measuring electronics. (b) Timing diagram for the multiplexed excitation and measuring scheme for the sensor. Voltages are sampled at the rising edge of the sample-and-hold signal. One frame is complete when the 64 transmitter electrodes have been activated. The signal SP denotes the polarity switch and S1–S64 the analog switches.
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
2189
2,500 images/s (Da Silva et al., 2007). It was demonstrated that the penetration depth of each sensing structure composing the whole sensor is of about 1 mm (Da Silva et al., 2007). That means, the sensor is able to visualize the flow only close to the wall it is attached. The electrical current flowing from the transmitter to the receiver electrode of a sensing structure depends on the electrical conductivity of the liquid surrounding the sensing structure. The current is converted into a voltage by amplifier electronics and the voltage value V therefore has general dependency on electrical conductivity in the form (Da Silva et al., 2007): V (i, j, k) = a(i, j) + b(i, j)
(3)
where a and b are slope and offsets describing electronics’ characteristics, i and j are the grid coordinates and k is the temporal index. To obtain an electrical contrast between pure water and water–glucose mixture one of the fluids needs to be marked with a conductivity tracer. This was done by adding a small amount of NaCl to the water–glucose mixture which increases electrical conductivity but only marginally changes density. It has been shown previously that the conductivity measurement can be linked to density gradients by (Höhne et al., 2006): (i, j, k) =
(i, j, k) − L ∼ (i, j, k) − L V (i, j, k) − VL (i, j) = = H − L H − L VH (i, j) − VL (i, j)
(4)
where is called mixing scalar with values between 0 and 1. The subscripts denote the heavy (H) and light (L) component. Thus, in order to obtain mixing scalar distributions (i,j,k) from the planar sensor, calibration measurements must be taken prior to experiments. This calibration involves the acquisition of voltage readings for the sensor fully covered with pure water VL (i,j) and for the sensor fully covered with NaCl–water–glucose solution VH (i,j).
As a first step, we implemented a differential optical flow algorithm to estimate the flow vector f(i,j) for each pixel i and j at time t of the 2D-image from the spatio-temporal derivatives of the gray values in the image sequence. For the present sensor, i and j vary from 1 to 64, since the sensor has 64 × 64 sensing structures. The algorithm is based on a concept similar to the continuity equation in fluid dynamics, where mass is conserved. Here, the gray value distribution g(i,j) in the image is assumed to be nearly constant for a short period, which yields the optical flow constraint equation (Jähne, 2002, pp. 387):
∇ g f + gt = 0
⎛
With increasing computer performance computation of velocity fields from image sequences have gained increasing attention. Optical flow is a prominent method (amongst other, such as cross-correlation algorithms) and is being widely used for motion tracking. Details of this method are found in Beauchemin and Barron (1995), Jähne (2002), as well as in Corpetti et al. (2004). The optical flow method has been developed especially for processing of images from optical acquisition systems. There, the projected motion field from the real 3D scene onto the image scene is extracted from the spatio-temporal change of the gray values in the corresponding image sequence which is defined as optical flow (Jähne, 2002). The optical flow method evaluates the change of brightness (or intensity) in a sequence of images caused, for example, by the movement of an object, and associates these intensity gradients with a velocity field. The planar array sensor is evidently no optical acquisition system and the gray values in the images represent an electrical parameter, i.e. the measured conductivity value at each electrode pair, which is represented by a gray value of the corresponding pixel in the image. However, the measured 2D conductivity value distribution represents the projected flow scene on the sensor surface in a given experiment and thus carries information of the near-wall velocity field. That is, the velocity field can be calculated from the spatio-temporal gradients of the gray values. The velocity vector v(x, y) can be inferred from the 2D optical flow vector f(i,j) according to
kg f (i, j) = fi /fj t
3.3. Optical flow for planar sensor data
(6)
where g is the local gradient of gray value, given by
3.2. Optical flow
v(x, y) =
Fig. 5. Illustration of the aperture problem in optical flow. The displacement vector of the edge in (a) cannot be determined uniquely, while for the movement of the corner and (b) it is possible. Adapted from Jähne (2002).
(5)
where kg is an scalar value which describes the geometry factor of the system and t is the time between two consecutive images.
∇g = ⎝
∂g ∂x ∂g ∂y
⎞
⎛
⎠≈⎝
g x g y
⎞
⎠=
gi gj
(7)
Further, gt is the discrete temporal derivative of the local gray value: ∂g g ≈ = gt t ∂t
(8)
However, the spatio-temporal derivatives are local operators, which act only in small regions of the image. This can be illustrated as a mask or aperture applied to the image (see Fig. 5). The solid line in Fig. 5a corresponds to an edge that is moving to the position of the dashed line in the second image. Due to the limited view of the local operators, the displacement vector of the optical flow may be connected from one point on the solid line to any point on the dashed line, i.e. there is no unambiguous correspondence. Thus, only the normal component f⊥ of the displacement vector can be correctly determined, whereas the parallel f component remains unknown. This phenomenon is well known as aperture problem (Jähne, 2002). Thus, the optical flow constraint equation (6) is nonunique. From the local operators, the optical flow vector can only be exactly determined when a corner of the object falls within the mask of the local operators as depicted in Fig. 5b. Due to the aperture problem, further conditions have to be taken into account to determine the vector of the optical flow. We have applied a global smoothness constraint first introduced by Horn and Schunck (1981) for regularization of the optical flow constraint equation. The global smoothness constraint assumes that optical flow vectors of neighboring pixels are nearly identical, i.e. the difference between the local average value of the optical flow f¯ and the optical flow of the neighboring pixels is assumed to be minimal.
2190
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
Fig. 6. Experimental setup for the validation of the optical flow algorithm.
This hypothesis yields the quadratic difference: ı2f = (f − f )
2
process is repeated until the difference between two iterations is below a given value in the form: (9)
Combining equation (9) with the uncertainty in the constraint equation resulting from signal noise: ın = ∇ g f + gt
(10)
yields the overall uncertainty constraint: ı2 = ˛2 ı2f + ı2n
(11)
which has to be minimized to achieve an estimation of the optical flow vector field. Thereby, ˛ is a fixed weighting factor being estimated empirically to adjust the impact of both terms of uncertainty. In our calculations the value of 0.1 was used. Minimization of ı2 is achieved by demanding ∂ı2 ! =0 ∂fi
(12)
n+1 n
f −f ≤ε
(16)
which is the norm of the difference between two consecutive iterations. In the calculations we empirically chose ε as 10−3 . In this way, the progress of the algorithm is adapted to the individual composition of the local spatio-temporal derivatives. In fact, the calculation of the local spatio-temporal derivatives in the standard algorithm by Horn and Schunck (1981) achieves sufficient accuracy in case of linear gray value distributions and oversampled image data only (Beauchemin and Barron, 1995). To achieve improved performance of the algorithm we implemented a method that evaluates five consecutive images for estimation of spatio-temporal derivatives for each step in time, which are then applied in the iterative calculation of the current optical flow vector field according to Eqs. (14) and (15) (Simoncell, 1994; Thiele, 2007). Once the optical flow is calculated the velocity can be easily achieved by implying the spatial dimensions of the sensing elements and the data acquisition rate as described in Eq. (5).
and 4. Validation experiment for velocity measurement
∂ı2 ! =0 ∂fj
(13)
from which one obtains two equations to be solved for each individual component of f. Since analytical solution of these both equations for all pixels of the image is highly computational intensive, an iterative set of two equations is derived from (12) and (13) solved with the Gauss–Seidel scheme (Horn and Schunck, 1981): n
f n+1 = f¯ i − gi i
n
n
gi f¯ i + gj f¯ j + gt ˛2 + gi + gj
(14)
and n
f n+1 = f¯ j − gj j
n n gi f¯ i + gj f¯ j + gt
˛2 + gi + gj
(15)
Eqs. (14) and (15) are solved iteratively for every step in time t and the optical flow in the iteration (n + 1) is estimated by the spatiotemporal derivatives and the previous locally averaged optical flow n f¯ (Horn and Schunck, 1981). The iterative algorithm for each step in time t is briefly described below. Initially, the spatio-temporal derivatives are calculated for every pixel and the averaged values f¯i0 and f¯j0 are set to zero. The first iteration is carried out yielding the optical flow components fi1 and fj1 . The neighboring pixels are locally averaged with a 2D averaging filter for every pixel of the image to obtain f¯i1 and f¯j1 which are used in the next iteration. This
In order to verify the feasibility of the implemented algorithm for velocity measurement from planar sensor data, a simple validation experiment was performed. A smaller prototype array sensor was used in this experiment. The electrodes dimensions are accordingly smaller than in the VeMix sensor, but the results obtained in this experiment are representative for the 64 × 64 sensor. Fig. 6 depicts the experimental setup. It consists of the sensor mounted in a horizontal position and a brush which was mounted on a 200 mm long fixation arm and mechanically moved across the sensor by means of a stepper motor. The motor was driven at constant angular speed of ω = 5.934 rad/s giving a tangential velocity of vT = 1.116 m/s for the moving brush. The sensor and a shallow vessel filled with tap water were placed horizontally within the circular path of the brush. In this way, the wet brush periodically passed across the sensor surface at constant velocity in a circular path. The brush movement along the planar sensor surface was captured with the electronics at 10,000 fps. This velocity is possible due to the fact that the prototype sensor used has only 32 × 32 electrodes (Da Silva et al., 2008). The motion of the wet brush in the image sequence was analyzed with the optical flow algorithm described in Section 3.3. From calculated optical flow vectors, the velocity vectors were determined using Eq. (5) with the spatial electrode dimensions (as shown in Fig. 6) and the time step between two consecutive images t = 0.1 ms. Fig. 7 shows four snap shoots of the liquid distribution
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
2191
Table 1 Results of the validation experiment for velocity measurement. Calculation type
Measured velocity (m/s)
Deviation from theoretical vT (1.116 m/s)
Maximum value over whole sensor Mean value of velocity vectors at the moving front
1.138 1.101
−2.0% 1.3%
on the sensor surface where the brush motion can clearly be identified. Furthermore, the figure also displays images of the temporal gradient values (Eq. (8)) which clearly show the front of the moving brush. For accuracy assessment an image sequence of 450 single images representing one passage of the wet brush over the whole sensor surface was analyzed. In order to compare the theoretical value of vT = 1.116 m/s and a measured velocity value, two approaches have been chosen. First, the maximal value of the calculated velocity vectors over the whole sensor surface was taken to represent the brush velocity for each time step. Second, the analyzed velocity vectors were confined to the moving front as indicated in Fig. 7 and the mean value within this region was chosen to represent the brush velocity for each time step. From the two data series containing 450 data points, the average value was calculated to determine the brush velocity, called vmax and vfront to indicate
the prior and latter calculation method. Table 1 shows the obtained values. Both calculated values show good agreement with theoretical velocity value vT with deviations below 2%. The value vmax is slightly higher than vfront . This can be explained by the fact that, when calculating maximal values, some image noise is picked up by the algorithm. 5. Vemix experiment and simulation In this section, the experimental and simulation results for an exemplary flow condition are presented. The essential experimental parameters are Q = 0.5 l/s (flow rate), Hfill = 2.12 m (i.e. L = 0.51 m), = 2%. This gives a Richardson number of 1.004 from Eq. (1) and, with a horizontal jet pattern flow regime (Section 2.3). Simulations were performed with the CFD model explained in Section 2.2. The planar array sensor was set to acquire 200 images per second during a 100 s time interval. In the experiment, pure water with w = 997 kg/m3 and water–glucose solution with wgs = 1017 kg/m3 and electrical conductivity of 60 and 200 S/cm, respectively, were used. The velocity fields were computed by the optical flow method as presented in Section 3.3. As previously discussed, the measured data with the planar sensor represents the flow close to the wall, while from CFD simulations a threedimensional data set of the flow is obtained. Thus, only the CFD data of an axial plane located at 1 mm distance from the wall were extracted and compared with the planar sensor’s experimental results. 5.1. Fluid distribution In Fig. 8 images of the distribution of mixing scalar (Eq. (4)) from the measured and simulated results are depicted. As the water–glucose solution flows into the channel a circular jet forms both in the experiment and the simulation. It seems that in the experiment mixing happens more rapidly and therefore the jet is more pronounced in the simulated distributions. In the experiment, a transition layer at the bottom of the images can be observed. It is caused by the pre-mixing of water and water–glucose solution prior to the start of the experiment. That means, the initial condition in the experiment is slightly different to that used in the simulation. Nevertheless, the agreement between both methods is very good. The use of the Richardson number could also correctly predict the flow regime for this experiment. 5.2. Velocity fields
Fig. 7. Selected images of (a) the liquid distribution on the sensor surface and (b) the temporal gradient for the validation experiment of velocity measurement. In (b) the location of the moving brush can be clearly recognized.
Further analysis of the measured data was carried out by computation of the velocity field. Fig. 9 illustrates four selected images of measured and the simulated velocity fields. The magnitude of the velocity is indicated by the background color and velocity vectors indicating only the direction are shown in the foreground. The simulation images show a dominating vortex at the bottom of the jet which bends the initial straight jet and creates the circular liquid distribution as observed in Fig. 8. In the velocity field computed from the measured data the vortex is not as perfect but still visible. Comparison of velocity values for experiment and simulation turns out to be difficult. It was found that reliable velocity vectors from optical flow are only found in regions of high conductivity contrast, i.e. high spatial gradient. Due to micro-mixing
2192
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
Fig. 8. Distributions of mixing scalar at the rear wall of the VeMix test facility for the mixing experiment at Ri = 1.004: (a) planar array sensor measurement and (b) simulation with ANSYS CFX.
Fig. 9. Velocity fields at the rear wall of the VeMix test facility for the mixing experiment for Ri = 1.004: (a) velocity information calculated by the optical flow method from experimental sensor data and (b) velocity information extracted from CFD data.
the boundary between high and low density liquid is smeared out which decreases local contrast and thus leads to noise pickup by the optical flow algorithm. A good example is the image at t = 18 s where grayscale gradients are smaller and the algorithm accuracy tends to be lower. The jet almost disappears in this image and a velocity spot appears which is not physical. The spot occurs because the gray value gradient in that region is high. In spite of these limitations, the optical flow method is still able to evaluate the velocity field assuming the gradients in the grayscale are present. For the VeMix experiment described here, this occurs predominantly in the jet region and at the beginning of the experiment. Another finding is that the experimental results show much finer turbulence details than predicted by the applied SST model in the simulation. This indicates that modification of the turbulence model is necessary to improve agreement. In a qualitative sense the agreement in the values of the velocity magnitude is quite good, indicating that the
general features of the flow are already well captured by the CFD simulation. 6. Conclusions The VeMix test facility was constructed to investigate buoyancy driven flow with an emphasis on CFD code development and validation. In the facility, the mixing of weakly and highly borated coolant can be simulated. The use of a new developed planar array sensor allows the measurement of the mixing scalar with high spatial and temporal resolution which can be directly and quantitatively compared with density distributions predicted from CFD simulations. In this paper, we also introduced a new method to obtain velocity fields by image processing of the planar sensor images with the optical flow method. A simple experiment with a moving brush was carried out to validate this velocity measure-
M.J. da Silva et al. / Nuclear Engineering and Design 240 (2010) 2185–2193
ment method. A mixing experiment was performed and analyzed. Good agreement between experimental data and first CFD simulations with a simple SST turbulence model was found. The velocity measurement by optical flow method proved to be applicable in general but exposed some problems in regions with progressed micro-mixing and accordingly small conductivity gradients. Experimental data was found to disclose much more turbulence detail than the simulations which indicates that improvement in turbulence modeling is necessary. In the future, more experiments will be carried out in the VeMix test. Further work with respect to the new measurement modality will be directed towards validation of optical flow method by simultaneous measurements with PIV. Acknowledgement This work was in part supported by the German Federal Ministry of Economy and Labor, grant 150 1287. References Beauchemin, S.S., Barron, J.L., 1995. The computation of optical flow. ACM Computing Surveys 27, 433–467. CFX11, 2007. User Manual. ANSYS-CFX. Corpetti, T., Heitz, D., Arroyol, G., Memin, E., Santa-Cruz, A., 2004. Estimation of motion using a PIV correlation-based method and an ‘optical-flow’ one for two experimental flows: quantitative and qualitative comparison. In: Proceedings of 12th International Symposium on Applications of Laser Techniques to Fluid Mechanics, Lisbon, Portugal.
2193
Da Silva, M.J., Suhnel, T., Schleicher, E., Lucas, D., Vaibar, R., Hampel, U., 2007. Planar array sensor for high-speed component distribution imaging in fluid flow applications. Sensors 7, 2430–2445. Da Silva, M.J., Lu, Y., Suhnel, T., Schleicher, E., Thiele, S., Kernchen, R., Diele, K.-H., Hampel, U., 2008. Autonomous planar conductivity array sensor for fast liquid distribution imaging in a fluid coupling. Sensors and Actuators A 147, 508–515. Höhne, T., Kliem, S., Bieder, U., 2006. Modeling of a buoyancy-driven flow experiment at the ROCOM test facility using the CFD-codes CFX-5 and TRIO U. Nuclear Engineering and Design 236, 1309–1325. Horn, B.K.P., Schunck, B.G., 1981. Determining optical flow. Artificial Intelligence 17, 185–203. Jähne, B., 2002. Digital Image Processing, 5th revised and extended Edition. Springer. Kliem, S., Suhnel, T., Rohde, U., Hohne, T., Prasser, H.-M., Weiss, F.-P., 2008a. Experiments at the mixing test facility ROCOM for benchmarking of CFD codes. Nuclear Engineering and Design 238, 566–576. Kliem, S., Prasser, H.-M., Suhnel, T., Weiss, F.-P., Hansen, A., 2008b. Experimental determination of the boron concentration distribution in the primary circuit of a PWR after a postulated cold leg small break loss-of-coolant-accident with cold leg safety injection. Nuclear Engineering and Design 238, 1788–1801. Menter, F.R., 1993. Zonal two equation k–ω turbulence models for aerodynamic flows. In: Proc. of 24th AIAA Fluid Dynamics Conference. Orlando, FL, AIAA paper 93-2906, July 1993. Rohde, U., Kliem, S., Höhne, T., Karlsson, R., Hemström, B., Lillington, J., Toppila, T., Elter, J., Bezrukov, Y., 2005. Fluid mixing and flow distribution in the reactor circuit, measurement data base. Nuclear Engineering and Design 235, 421–443. Simoncell, E.P., 1994. Design of multi-dimensional derivative filters. In: Proc. of 1st IEEE International Conference on Image Processing, Austin, TX, USA, pp. 790–793. Thiele, S., 2007. Entwicklung und Aufbau eines kapazitiven Oberflächensensors für die Phasenverteilung- und Geschwindigkeitsmessung zweiphasiger Strömungen. Diploma Thesis, Technische Universität Dresden, Dresden, Germany (in German). Umminger, K., Kastner, W., Liebert, J., Mull, T., 2001. Thermal hydraulics of PWRS with respect to boron dilution phenomena. Experimental results from the test facilities PKL and UPTF. Nuclear Engineering and Design 204, 191–203.