
Full Citation 
Material Information 

Title: 
14.4.5  Two Phase Cross Jet in a Fuel Rod Assembly Using DNS/LevelSet Method Computational Techniques for Multiphase Flows 

Series Title: 
7th International Conference on Multiphase Flow  ICMF 2010 Proceedings 

Physical Description: 
Conference Papers 

Creator: 
Behafarid, F. Bolotnov, I.A. Podowski, M.Z. Jansen, K.E. 

Publisher: 
International Conference on Multiphase Flow (ICMF) 

Publication Date: 
June 4, 2010 
Subjects 

Subject: 
twophase flow stabilized finite element level set method DNS 
Notes 

Abstract: 
One possible accident scenario in a hexagonal fuel assembly of a Gen. IV sodiumcooled fast reactor (SFR) is associated with
local heat degradation around the hottest fuel element due to the coolant flow blockage by foreign objects and/or mechanical
failure of the hottest fuel rods. As a result, high pressure fission gases would be injected into the liquid sodium coolant
channels as a high speed twophase crossflow jet.
A model of gasjet/liquidsodium flow evolution has been modeled using the Parallel Hierarchic Adaptive Stabilized
Transient Analysis (PHASTA) software. This new model of turbulent, unsteady, threedimensional twophase flow uses
unstructured grids and multiple (more than 8000) processors. The PHASTA code uses a direct numerical simulation (DNS)
method and is capable of capturing the shape of gas/liquid interfaces using the level set method.
An intrinsic flaw was found in the standard approach to the modeling of gas/liquid/solid interfaces using the level set method,
and appropriate steps to correct it are discussed in the paper. It has also been demonstrated that the pressure outflow boundary
conditions (BCs) commonly used in singlephase flow simulations are not applicable to twophase flows with resolved
gas/liquid interfaces, since they cannot handle the surfacetensioninduced sharp pressure gradients across the interface. A
new approach to treating the outflow BC in twophase flow problems has been proposed. 

General Note: 
The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: BioFluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and NanoScale Multiphase Flows; Microgravity in TwoPhase Flow; Multiphase Flows with Heat and Mass Transfer; NonNewtonian Multiphase Flows; ParticleLaden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows 
Record Information 

Bibliographic ID: 
UF00102023 

Volume ID: 
VID00354 

Source Institution: 
University of Florida 

Holding Location: 
University of Florida 

Rights Management: 
All rights reserved by the source institution and holding location. 

Resource Identifier: 
1445BehafaridICMF2010.pdf 

Full Text 
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Two Phase Cross Jet in a Fuel Rod Assembly Using DNS/LevelSet Method
F. Behafarid, I.A. Bolotnov, M.Z. Podowski and K.E. Jansen
Center for Multiphase Research,
Rensselaer Polytechnic Institute,
110 8th St., Troy, NY 12180, USA
behaff@rpi.edu
Keywords: Two phase flow, stabilized finite element, level set, DNS
Abstract
One possible accident scenario in a hexagonal fuel assembly of a Gen. IV sodiumcooled fast reactor (SFR) is associated with
local heat degradation around the hottest fuel element due to the coolant flow blockage by foreign objects and/or mechanical
failure of the hottest fuel rods. As a result, high pressure fission gases would be injected into the liquid sodium coolant
channels as a high speed twophase crossflow jet.
A model of gasjet/1iquidsodium flow evolution has been modeled using the Parallel Hierarchic Adaptive Stabilized
Transient Analysis (PHASTA) software. This new model of turbulent, unsteady, threedimensional twophase flow uses
unstructured grids and multiple (more than 8000) processors. The PHASTA code uses a direct numerical simulation (DNS)
method and is capable of capturing the shape of gas/1iquid interfaces using the level set method.
An intrinsic flaw was found in the standard approach to the modeling of gas/1iquid/solid interfaces using the level set method,
and appropriate steps to correct it are discussed in the paper. It has also been demonstrated that the pressure outflow boundary
conditions (BCs) commonly used in singlephase flow simulations are not applicable to twophase flows with resolved
gas/1iquid interfaces, since they cannot handle the surfacetensioninduced sharp pressure gradients across the interface. A
new approach to treating the outflow BC in twophase flow problems has been proposed.
V, Volume of each phase within a cell (m )
w Pseudo velocity in redistancing
Greek letters
e Half of the thickness of the interface (m)
pu Dynamic viscosity (kg/m s)
p Density (kg/m )
r Pseudo time in redistancing
E.Stress tensor (N/m )
Subscripts
1 first phase: Liquid
2 Second phase: Gas
i or j Coordinates (Defines spatial components)
,i or j Spatial derivative in that direction
,t Time derivative
PHASTA Description
PHASTA is a parallel, hierarchic (higher order accurate,
from the 2nd to 5th order accuracy, depending on function
choice), unstructured/adaptive, stabilized finite element
transient analysis flow solver for both incompressible and
Introduction
Development of next generation nuclear reactors, in
particular Gen. IV reactors, will require a new generation
of advanced modeling tools. Thus, it is important that
advanced computational capabilities be developed to
capture the multiscale physics phenomena governing
reactor transients and accidents.
The purpose of this paper is to discuss a novel
mechanistic approach to 3D simulations of the injection
of gaseous fission products into a partially blocked SFR
coolant channel following locahized cladding overheat and
failure. The proposed methodology has been developed
for, and implemented in, the PHASTA DNS code,
coupled with the Level Set Method (LSM) solver.
The computational mesh used in the simulation is highly
refined around the injection zone, to let PHASTA capture
the effect of turbulence and accurately predict the
evolution of interface shape.
Nomenclature
d Corrected distance field (m)
SBody forces (N)
H, Heaviside kernel function
p Pressure (N/m )
S, Strain rate tensor (s')
ui Velocity (m/s)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
incompressible flow of a Newtonian fluid, the stress
tensor is related to the strain rate tensor, SU, as:
.i. = 2pSij.= p(ut.j+ ujp)
PHASTA
domain
Figure 2. Top view of the PHASTA computational
domain.
Using the Continuum Surface Tension (CST) model of
Brackbill et al. (1992), the surface tension force is
computed as a local interfacial force density.
The finite element formulation is based on the so called
weak torm of equations which are generated by dotting
continuity and momentum equations on the left by the
weight functions and integrating over the domain.
Following the stabilized formulation of Taylor et al
(1998) the residuals of the continuity and momentum
equations are added together.
The stabilized finite element formulation is described in
detail in Rodriguez (2009) and Jansen et al (1999). The
stabilization and its parameters were described by Franca
and Frey (1992) and by Whiting and Jansen (2001).
Level Set Method
The level set method of Sussman (Sussman et al, 1998,
Sussman and Fatemi, 1999 and Sussman et al., 1999) and
Sethian (1999) involves modeling the interface as the zero
level set of a smooth function, 9, which represents the
signed distance from the interface. Hence, the interface is
defined by p = 0. The scalar 9, the socalled first scalar, is
convected within a moving fluid according to,
D =fv +.V= 0
Phase1, typically the liquid phase, is indicated by a
positive level set, p > 0, and phase2 by a negative level
set, p < 0. Evaluating the jump in physical properties
using a step change across the interface leads to poor
computational results. Therefore, the properties near an
interface are defined using a smoothed Heaviside kernel
function, H,, given by:
HE ( ) = 1+ + i
The fluid properties are thus defined as:
compressible flows (Jansen et al., 2000, Whiting and
Jansen, 2001)
PHASTA was the first unstructured grid LES code
(Jansen, 1993) and has been applied to turbulent flows
ranging from validation benchmarks (channel flow, decay
of isotropic turbulence) to complex flows (airfoils at
maximum lift, flow over a cavity, near lip jet engine flows
and fintube heat exchangers).
PHASTA uses advanced anisotropic adaptive algorithms
(Sahni et al., 2006) and the most advanced LES/DES
models (TejadaMartinez and Jansen, 2005). Note that
DES, LES, and DNS are computationally intensive even
for single phase flows. The twophase version of
PHASTA utilizes the Level Set method to define the
interface between the gas and liquid phases (Sethian,
1999). Here, it has been used to simulate gas jet injection
into the coolant region.
Simulation Domain
PHASTA simulations were performed for a small region
in front of the crack. Figures 1 and 2 show the
computational domain sharing in the multiscale approach
used in the simulations.
PHASTA outlet twophase
mixture outflow
Inlet
velocity
profile
Figure 1. Computational domain (side view).
PHASTA Formulation
Governing Equations
PHASTA is a parallel, hierarchic, unstructured adaptive,
stabilized (finite element) transient analysis flow solver
(both incompressible and compressible). The spatial and
temporal discretization of the Incompressible NS
equations within PHASTA has been described in Whiting
(1999) and Nagrath !1lI H~ This section follows their
discussion.
The strong form of the INS equations is given by:
Continuity: ui j 0
momentum: pui~t +pu ui, j ,i +9r, j i
where ui is ith component of velocity and fi represents
body forces along the ith coordinate. For the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where H' is aHE (d) Bd and is given by:
HE()=11 1 a
Hence, to minimize the deviation in volume each new di~
is constrained to satisfy the following:
J] H'/(do)(dk Odo x= 0
This is accomplished by projecting the value of the
current level set, dy onto the new level set field, dio ,
where di is computed from:
di = di + f fk () )HE (d )
The scaling function /1U is assumed constant within each
cell and is computed from:
HE (d ) ~k ()d
k O
p(#) = PlHE ~ P2(1 HE ())
pU( ) = pl1HE(~ 2 (1 HE ())
Although the solution may be reasonably good in the
immediate vicinity of the interface, the distance field may
not be correct throughout the domain since the varying
fluid velocities throughout the flow field distorts the level
set contours. Thus the level set is corrected with a re
distancing operation by solving the following PDE:
ad
= S( ) 1 Vd
where d is a scalar that represents the corrected distance
field and r is the pseudo time over which the PDE is
solved to steadystate. This may be more clearly
expressed as the following transport equation:
+ E:~Vd = S(4)
The second scalar, d, is originally assigned the level set
field, p, and is convected with a pseudo velocity Ei,
where
Vd
G = S( ) Vd
and S(0) is defined as:
S() )= + ]1si ,<
Note that the zeroth level set or interface, p = 0, does not
move since its convecting velocity Ei is zero. Solving the
second scalar to steadystate restores the distance field
to Vd = f1but does not alter the location of the interface.
First scalar, 9, is updated using the steady solution of d
Volume Constraint on Level Set Redistance Step
Sussman proposed an additional constraint to be applied
during the redistance step to help ensure the interface
(p=0) does not move (Sussman and Fatemi, 1999 and
Susman et al., 1999). They found imposing this constraint
improves the convergence of the redistance step while
maintaining the original zero level set. The essence of the
constraint is to preserve the original volume of each phase
during the redistance step. The volume of each phase
within a cell, V , is defined as:
where dk is the value of the distance field, d, at iteration k.
Small changes in the zero level set from pseudo time step
O to Zk may be approximated by:
k () k () He(d )
k a
= J~ H:(d ) d d )dv
e e
"j HE(d ~) dx
Note when the distance field has converged,
dk = ,) and
/1U= 0. In incompressible PHASTA, the volume constraint
correction occurs at the global level after each re
distancing iteration. The calculation of /1U occurs at the
element level and is projected back to the global level.
The Level set method was implemented into PHASTA by
Nagrath !II II 1, and Rodriguez !II II **,s For the simulations
presented herein the authors used the version of the code
provided by Dr. Rodriguez with several modifications
described below.
LSM Shadowing Problem
Level set method models the exact shape of the interface
as the zero level set of a smooth function which represents
the signed distance from the interface. The convection of
the distance field by the liquid flow moves the interface
properly but distorts the distance field. To recover it, a
few redistancing (or reinitiating) iterations should be
applied after each time step. Redistancing has an intrinsic
flaw which can cause errors and introduce a false phase in
some parts of the domain.
This problem is related to the lack of boundary conditions
for distance field and has not been discussed in the
literature before. In redistancing, the information flows in
the direction normal and away from both sides of the
interface. If the interface is not attached to the boundary
or is attached in the normal direction to it (see Figure 3),
the level set method does not cause errors.
However, when the interface is attached to the wall
(creating a triple point) at an arbitrary angle, there are
some regions that cannot get the correct information about
the distance field due to the location of those regions in
relation to the interface and to the wall. We call these
regions "shadowed". One example of a shadowed region
is shown in Figure 4.
If angle is 900 or less, there is
O ~no shadowing problem
Figure 3. Schematic of conditions when LSM has no
problem even if the interface touches the boundary
In the shadowed regions, the information should come
from the boundary which has no boundary condition for
the distance field; therefore, meaningless information may
propagate into, and pollute, the domain. In this case, more
redistancing iterations make the distance field worse
Figure 4. Schematic of two typical conditions when
interface: touches the boundary (Solid) and LSM
encounters the shadowing problem.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
\h,,,1. as. J region
Figure 6. Distance filed after redistancing. The distance
field in shadowed region is getting worse.
In most cases it is difficult to set a boundary condition for
the distance field in order to prevent the shadowing
problem. The current approach is to minimize the impact
of this error by using the limited number of redistancing
operations.
Future work will include a level set method enhancement
which will allow the enforcement of a physically based
value of a contact angle at the triple point.
New BC for Two Phase Flow Outlet
Three types of BCs commonly used for outflow for
incompressible fluids are "Natural Pressure BC",
"Convective Pressure BC" and "No Pressure BC". None
of the three can properly handle the lateral pressure
gradient due to surface tension force in two phase flow
simulations.
The common way to avoid problems in outflow is to use
periodic boundary conditions. The issue can also be
avoided in cases where the pressure jump through
interface is low; either due to a large interface curvature
or a small value of the surface tension. In such cases, the
"Natural Pressure Boundary Condition" can handle this
small pressure difference with little or no effect on the
interface far from the outlet. In other similar cases, the
main stream velocity is high enough to convect the
unphysical data out of the domain. In the case of interest
to this work, none of the above mentioned conditions are
satisfied. Since no solutions are available in the open
literature a new method has been developed to solve the
problem.
A simple explanation of the problem is as follows. The
surface tension force acts on, and is balanced by the
pressure gradient across, the interface. At the outlet, the
"Natural Pressure BC" tries to make the pressure uniform
throughout the outlet and eliminates the pressure
difference. Therefore, the surface tension cannot be
balanced by pressure difference and moves the interface
in a nonphysical manner. Figures 7 9 show a coflow
case with the "natural pressure outlet BC In order to
solve this problem, we added an artificial body force at
each nodal point in a very thin layer adjacent to the outlet,
to maintain the surface tension/pressure gradient balance.
The simulation results mn this region can be considered as
a numerical extension to the real domain and have no
physical meaning. Smnce this region is very thmn (of the
order of 0.1 mm), the additional computational expenses
are negligible. To properly apply the body forces, this thin
layer should have at least three to five elements across.
Figure 5. Distance field after convection by flow and
before redistancing. The arrows show the direction of
information, normal to the interface used to rebuild the
distance field.
Figures 5 and 6 show a shadowing problem observed in
the current PHASTA simulations. In Figure 5, the
distance field is shown after it has been convected by the
fluid flow solution. At this stage the interface has moved
and the distance field is distorted. A series of re
distancing iterations is applied to reconstruct the distance
field.
Figure 6 shows the distance field after redistancing
iterations have been performed. As discussed above, the
distance field degrades after redistancing in the
shadowed regions only, where the distance field
information flows from the boundary to the interior of the
domain. In this region, no matter how many redistancing
iterations are used, the distance field cannot be recovered.
I tI I I I
Figure 10. Schematic of the artificial body forces applied
to modify the outlet BC when theof plume of gas exits the
computational domain.
Outlet: Natural Pressure
Inlet: Velocity poi
0s~i 0200ir
00150"
Figure 11. Test case geometry and velocity profile.
velocity Magnitude
000 0.0500 0.100 0.150 0.200
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
SPressume gradient
Surface tension
SInterface
SBody forces inside of the interface
SBody forces in the gas region
Artificial body
forces are added
after this line
SOutlet
Figure 7. Geometry colored by velocity magnitude. The
small opening in the middle reflects gas injection region.
Outlet
Natural Pressure
Liquid bIlet 7
Gas bIlet
Figure 8. Pressure field for the coflow case before
gas/1iquid interface reaches domain outlet.
Figure 9. Effect of surface tension on pressure imbalance
when the interface reaches domain outlet.
Figure 10 shows a schematic of the plume of gas reaching
the outlet and shows the size and direction of the artificial
body forces. The violet arrows show the body forces
applied inside the interface to balance the difference
between surface tension and pressure gradient. On the
other hand, in the gas region the pressure goes down in
the direction normal to the outlet. To avoid gas
acceleration, another body force has applied which is
shown as black arrows.
Figure 11 shows the geometry, velocity profile and the
boundary conditions of a simple case used to test the new
outlet BC.
Figures 12 and 13 present the results of simulations using
the new modified boundary condition for the case of co
flow with a cylindrical interface.
Since the initial interface is defined as a perfect steady
cylinder, even very small changes in the interface shape at
the outlet can be readily noticed. As shown in Figure 13,
pressure inside the gas plume goes down to zero but the
artificial body forces balance the surface tension.
Figure 12. Cross section of the velocity field in a coflow
test case. The initial interface is a cylindrical plume of
gas (white lines) all the way to the outlet.
In test calculations, the interface maintained its original
shape after 2000 time steps. Similar tests were performed
for other interfaces including bubbles and more complex
shapes, and all of them shown desired results.
t t i t 1 i t
c
c
c
c
c
c
1 lit f f 1
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Different gas injection velocities have been used to test
their effect on the interface shape. A parabolic velocity
profile was used at the crack inlet, corresponding to three
different values of maximum velocity: 5 m/s, 15 m/s and
120 m/s.
070s 10
Figure 13. Pressure field after an artificial body force has
been applied to the right of the white line. A uniform
pressure is observed at the outlet, whereas the interface
shape stays unchanged.
Results
Gas injection
The first geometry under consideration, shown in Figure
14, included two fuel rods. Figure 15 shows the 3D
geometry colored by distance field. The crack has a
rectangular 0.1 mm by 0.1 mm shape, and it is initially
filled with gas (as shown by a negative sign of the
distance field variable).
/ conywasutroal
r
/ *
Figure 14. Computational domain (top view).
The following fluid properties have been used in the
simulations presented below:
kg ; 1431nkg
p, = 754.9 ,:p .8x0 
p, = 4.481 k ,=4051 g
n2 n2 s
where the liquid sodium properties are taken at 800oC and
the fission gas is assumed to have the properties of
nitrogen. The surface tension coefficient used in the
simulations was, a = 0.19 N/2.
sea cano e
33 01 06
scrl
MO.OO89
0.00430
0.00270
0.00110
0 000501
Figure 15. Computational domain including a small
crack colored by distance field (having negative values in
the gas region. Initially, the interface was at the
intersection of the crack and the main domain.
Figures 16 18 show the velocity and pressure fields, as
well as the 3D interface shape for the maximum gas
injection velocity of 5 m/s after three different numbers of
time steps. The left part of each figure shows the velocity
magnitude distribution in a 2D cut of the computational
domain. The cut plane was chosen to include the
directions of mean jet velocity and mean liquid velocity.
The solid white line shows the gas liquid interface. The
center part of each figure presents the pressure
distribution in the same plane cut. The interface is shown
with a black solid line. Higher pressure values can be
observed in the gas region due to the surface tension
effect. Finally, the right part of each figure presents a 3D
view of the gas/1iquid interface. The color of the interface
represents the local velocity magnitude.
As discussed before, in these simulations a parabolic
velocity profile has been used at the gas injection
boundary condition. It should be noticed that the
maximum liquid velocity before the gas injection is
almost 10 times lower than the maximum gas injection
velocity, but the effect of gravity is strong enough to shed
bubbles from the inlet gas stream. In this simulation and
the following ones, the gas density was almost 168 times
less than the liquid sodium.
In figures 19 21, the maximum gas injection velocity
increased from 5 m/s to 15 m/s. In this case enough gas is
provided to avoid bubble shedding. Also, the gas flow is
less stable and exhibits somewhat turbulent behavior
compared to the previous case.
Figures 22 and 23 show the results for a maximum gas
injection velocity of 120 n2/s. The gas flow is very
turbulent and starts to expand the interface close to the
wall. For this special geometry and gas injection velocity,
the interface is expected to create a lot of small bubbles.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
essure
1 500+03
938.
375
188
750
11.2
7.50
3 7
2 50
11.2
Figure 16. Simulation of twophase cross jet injection.
The gas and liquid velocity fields are colored by pressure,
via; = 5 m/s, of time step number: 2,000.
Figure 19. Simulation of twophase cross jet injection.
The gas and liquid velocity fields are colored by pressure,
vi4j= 15 m/s, of time step number: 5,000.
11.2 938. 11.2
7.50 1 = 376./ 7.50
Figure 20. Simulation of twophase cross jet injection.
The gas and liquid velocity fields are colored by pressure,
vi4j= 15 m/s, time step number: 8,500.
11.2 938. 11.2
7.50 376. 7.50
H .75 188 37
Figure 21. Simulation of twophase cross jet injection.
The gas and liquid velocity fields are colored by pressure,
vi,j = 15 m/s, time step number: 11,500.
Figure 17. Simulation of twophase cross jet injection.
The gas and liquid velocity fields are colored by pressure,
v~i = 5 m/s, time step number: 4,000.
3mi
2 50
25:
Figure 18. Simulation of twophase cross jet injection.
The gas and liquid velocity fields are colored by pressure,
vin; = 5 m/s, time step number: 5,500.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
coarse region is eight times bigger than in the refined
region, the computational cost of the newly added domain
is about 500 times (8 ) smaller per unit volume.
Figure 24. Top view of the velocity field for the
computational domain including 2 fuel rods, vmax=0.462
Tegsadlqiveoiyfedarcooebypressure,
v;,,y= 120m/s, time step onume:100
a/ I
Figure 23. Simulation of twophase cross jet injection.
The gas and liquid velocity fields are colored by pressure,
vi,j = 120 m/s, time step number: 2,000.
.
mesh resoltio which cain be tpaffred inos jea pareticula
Inthe prevn ious aeapaaoi a velocity fed r ord proie was ,
s et y asa oundar condition capt the ga tint To obtain a
motrephyical vhnoeloit profle winside the crack, the
comutatisonald doain wtas late extmenddto inclusdiae gas
(oresult 00) Tu, the rc oehsbeent etene areosse the
Iths allso bubeenosre in the simulations thateen uo the
syme etry faces which are cloe tfo the interace patcanaffct
dmai ie h odified geometry ishonnFgue
Iet shoul be menti oned here tat the compultationa ostai of
thre adde reion i veriy smafll since the mehcrank be
loars in the regionnls farfrom the interfacdte whchpertin .
Ito hsinglehs ben fowFr exape if the smeltos h szei the
Figure 25. Top view of the velocity field for the
computational domain including several fuel rods
surrounding one intact fuel element, vmax= 0.008 m/s.
Azimuthal flow instability
The inclusion of a larger domain also yielded an
interesting finding during the stage of simulation where
the fully developed liquid flow was occurred prior to the
start of gas injection. Since the geometry is quite
complex, no analytical methods could be used to
formulate the initial conditions for the liquid coolant
inflow velocity profile. Instead, a stagnant flow was used
at the starting point, combined with periodic boundary
conditions in the streamwise coolant direction and a
proper body force to drive the flow and to create the fully
developed velocity profile. The resultant velocity profile
was then used as the initial and coolant inflow conditions.
One of the major differences between the geometries
shown in Figure 24 and Figure 25 is that the latter
geometry does not impose the symmetry boundary
condition which was an inherent part of the former,
As shown in Figures 25 and 26, after the symmetry lines
have been removed and a complete fuel rod was includes
inside the domain, an azimuthal instability have occurred,
0.00592~
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
quite a stable behavior. Due to the relatively slow liquid
flow, no bubble shedding is expected for this case. It can
be noticed from the pressure plots in Figures 28 and 29,
that a higher gas pressure pushes the jet out of the crack
(the left part of the central pressure contour plot). Also,
as the size of the gas region in the channel grows, the
incoming gas jet starts to show a more turbulent behavior,
as it is no more constrained by the relatively heavy fluid.
Figure 30 shows the velocity distribution and the interface
evolution at various downstream distances from the
injection zone. It is interesting to notice two nearly
symmetrical recirculation regions in the middle of the gas
bubble in Figure 30 (b, c).
enhancing turbulence development at a much lower mass
flow rate than before. This instability first starts rotating
the fluid in one direction, and then in the opposite
direction, so that the flow eventually becomes turbulent,
10 .0 2
Figure 26. Unstable velocity field around the failed rod,
vmax = 0.07 m/s. A counterclockwise rotation can be
noticed.
The maximum velocity to create a laminar velocity profile
in Figure 25 is almost 57 times smaller than that of the
Figure 24. This can be explained by the fact that the
effective characteristic diameter of the flow shown in
Figure 25 is much larger due to the enclosed liquid loop
around the fuel rod. This fact results in reaching the
critical Reynolds number at much lower liquid velocities
compared to the case shown in Figure 24. Figure 27
shows the mesh corresponding to the case of a larger
computational domain with the ring. This mesh has
around 15 million elements and more than 8000
processors were used in the simulations.
Figure 28. Sideview cuts of the pressure and velocity
fields and a top view of the velocity field: vinj, max= 8.4
m/s: crack size: 1 by 2mm, time step number: 12,000:
transient duration (physical time): 0.021 s.
Figure 29: Sideview cuts of the pressure and velocity
fields and a top view of the velocity field: vinj, max = 8.4
m/s: crack size: 1 by 2 mm, time step number: 24,000:
transient duration (physical time): 0.043 s.
Summary and Conclusions
A DNS approach with Level Set Method for interface
tracking has been used to numerically simulate a
hypothetical nuclear reactor accident. It was shown that
the standard Level Set Method has an intrinsic problem on
the boundaries near the gas/1iquid/solid interfaces during
the redistancing stage, a ways to overcome this problem
were discussed. It was also was shown that no common
outlet boundary condition can handle the sharp pressure
jump across the interface due to surface tension force. A
Figure 27. Crosssection of the mesh. The mesh is refined
inside and around the crack, and is coarse far from it.
Figures 28 and 29 show the evolution of gas flow for the
new case, based on the modified geometry and the new
velocity profile at the crack. The gas/1iquid interface is
shown as a white solid line in the left and center part of
the figures. The right parts of Figures 28 and 29 show the
domain slices across the coolant flow direction, with the
interface shown again as a white line.
It has been observed that with a larger size crack and a
relatively small jet velocity, the gas/1iquid interface is
well controlled by the surface tension force and it exhibits
velocity Magnitude
velocity Magn tude
4.2~113
velocity Magiitude
.169
.056
U
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Jansen, K. E., Collis, S. S., Whiting, C., Shakib, F., "A
better consistency for loworder stabilized finite element
methods" Computer methods in applied mechanics and
engineering, Vol. 174, 153170 (1999).
Jansen, K. E., Whiting, C.H. and Hulbert, G.M., "A
generalizedalpha method for integrating the filtered
NavierStokes equations with a stabilized finite element
method", Computer Medrods in Applied Mechcaries card
Engineering, Vol. 190, 34, 305319 (2000).
Nagrath, S. "Adaptive stabilized finite element analysis of
multiphase flows using level set approach", Ph.D. Thesis,
MANE, Rensselaer Polytechnic Institute (2004).
Nagrath, S., Jansen, K.E. and Lahey, Jr., R.T.,
"Computation of incompressible bubble dynamics with a
stabilized finite element level set method", Computer
Medrods in Applied Mechcaries and Engineering, Vol.
194(4244) 45654587, (2005).
Rodriguez J. M., "Numerical simulation of twophase
annular flow", Ph.D. Thesis, MANE, Troy, NY,
Rensselaer Polytechnic Institute (2009).
Sahni, O., Mueller, J., Jansen, K.E., Shephard, M.S. and
Taylor, C.A., "Efficient Anisotropic Adaptive
Discretization of the Cardiovascular System", Computer
Medrods in Applied Mechcaries and Engineering, Vol.
195, 4143, 56345655 (2006).
Sethian, J. A., "Level Set Methods and Fast Marching
Methods", Cambridge University Press. (1999)
Sussman, M., E. Fatemi, P. Smereka and S. Osher., "An
improved level set method for incompressible twophase
flows", Journal of Computers & Fluids, Vol. 27, Issue 5
6, 663680 (1998).
Sussman, M. and E. Fatemi., "An efficient, interface
preserving level set redistancing algorithm and its
application to interfacial incompressible fluid flow."
Siam Journal on Scientific Computing, Vol. 20, Issue 4,
Pages 11651191 (1999).
Sussman, M., A. S. Almgren, J. B. Bell, P. Colella, L. H.
Howell and M. L. Welcome., "An adaptive level set
approach for incompressible twophase flows.", Journal of
Computational Physics, Vol. 148, Issue 1, 81124 (1999).
Taylor, C. A., T. J. R. Hughes and C. K. Zarins, "Finite
element modeling of blood flow in arteries", Computer
Methods in Applied Mechanics and Engineering, Vol.
158, Issues 12, 155196 (1998).
TejadaMartinez, A.E. and Jansen, K.E., "A parameter
free dynamic subgridscale model for largeeddy
simulation", Computer Medrods in Applied Mechanics
and Engineering, Vol. 194, No. 9, 12251248 (2005).
Whiting, C. H., "Stabilized finite element methods for
fluid dynamics using a hierarchical basis", Ph.D. Thesis,
MEAE, Troy, NY, Rensselaer Polytechnic Institute
(1999).
Whiting, C.H. and Jansen, K.E., "A Stabilized Finite
Element Formulation For The Incompressible Navier
Stokes Equations Using A Hierarchical Basis,"
International Journal of Numerical Medrods in Fluids,
Vol. 35, 93116 (2001).
new method was presented to provide a proper outlet BC
for twophase flow simulations.
The results of the gas injection into the liquid sodium for
different crack sizes and gas velocities were shown and
discussed. The selection of the main domain geometry
and its effect on the simulation were investigated.
Future work will meclude the simulations for extended
domains with timedependent gas mnflow rates and a time
varying crack size.
;eocty Magntuc~
PO 22 43 65 8.6
Figure 30: Top views of velocity field: (a) at the center of
the crack, (b) 2 mm, (c) 4 mm, and (d) 6 mm downstream
from the center of the crack. Time step number: 13,000.
Acknowledgments
The authors would like to acknowledge the financial
support provided to this study by the U.S. Department of
Energy and the computation resources provided by the
Computational Center for Nanotechnology Innovations
(CCNI) at Rensselaer Polytechnic Institute (RPI).
References
Brackbill, J. U., D. B. Kothe and C. Zemach, "A
continuum method for modeling surface tension", Journal
of Computational Physics, Vol. 100, Issue 2, Pages 335
354 (1992).
Franca, L. P. and S. L. Frey. "Stabilized FiniteElement
Methods. II. The Incompressible NavierStokes
Equations", Computer Methods in Applied Mechanics
and Engineering, Vol. 99, Issues 23, Pages 209233
(1992).
Jansen, K.E., "Unstructured Grid Large Eddy Simulations
of Wall Bounded Flows", Annual Research Briefs, Center
for Turbulence Research, NASA Anzes/Stcozford
University, 151 (1993).
U U

