Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0003720/00001
## Material Information- Title:
- A Computational Thermal-Fluid Modeling Approach to Pulsed and Steady Gas Core Reactors
- Creator:
- CHARMEAU, ANNE (
*Author, Primary*) - Copyright Date:
- 2008
## Subjects- Subjects / Keywords:
- Charge flow devices ( jstor )
Gas density ( jstor ) Inlets ( jstor ) Mach number ( jstor ) Modeling ( jstor ) Nozzles ( jstor ) Pressure pulses ( jstor ) Shock tubes ( jstor ) Shock waves ( jstor ) Uranium ( jstor )
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Anne Charmeau. Permission granted to University of Florida to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 4/30/2004
- Resource Identifier:
- 55893158 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

PAGE 1 A COMPUTATIONAL THERMAL-FLUID MODELING APPROACH TO PULSED AND STEADY GAS CORE REACTORS By ANNE CHARMEAU A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2004 PAGE 2 Copyright 2004 by Anne Charmeau PAGE 3 To my great-grandmother Mathilde PAGE 4 ACKNOWLEDGMENTS I would like to thank Dr. Samim Anghaie, Dr. Edward Dugan and Dr. Bruce Carroll for being members of my supervisory committee. I would like to particularly express my sincere gratitude to Dr. Samim Anghaie, who provided excellent tutoring as my advisor and showed constant patience and support. I would like to acknowledge Dr. Travis Knight for his assistance on the coupling between FLUENT and MCNP and Dr. Blair Smith for his help and his amazing sketches. Innovative Nuclear Space Power and Propulsion Institute (INSPI) at the University of Florida provided support for this research. I would like to thank my family for their love, support and understanding, and my friends for always being here for me. iv PAGE 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES............................................................................................................vii LIST OF FIGURES.........................................................................................................viii NOMENCLATURE..........................................................................................................xi CHAPTER 1 INTRODUCTION........................................................................................................1 Introduction...................................................................................................................1 History of the Gas Core Nuclear Rocket Why Choosing a Propulsion based on Nuclear Energy........................................................................................................1 Shock Wave Driven Gas Core Reactor........................................................................3 Thesis Objectives..........................................................................................................5 Thesis Organization......................................................................................................5 2 COMPUTATIONAL FLUID DYNAMICS.................................................................7 What is Computational Fluid Dynamics?.....................................................................7 Governing Equations....................................................................................................8 Solver Formulation.....................................................................................................10 Turbulence Modeling..................................................................................................10 FLUENT 6.1...............................................................................................................13 Creating Grids and Meshing.......................................................................................14 User-Defined Functions (UDFs).................................................................................14 Solution Controls........................................................................................................15 Convergence Check-Residuals...................................................................................15 3 MODELING AND PROPERTY DATA....................................................................17 General Assumptions..................................................................................................17 Assumptions on the Fluids Used in the Study............................................................18 Nature of Pressure Pulses Generated..........................................................................24 v PAGE 6 Summary of Approximations and Assumptions.........................................................25 4 STUDY OF PRESSURE PULSES EVOLUTION IN SHOCK TUBE.....................26 Design of the Core......................................................................................................26 Fluid Used...................................................................................................................27 Creation of the Grids..................................................................................................27 Preliminary Studies.....................................................................................................28 Set-up of the Cases.....................................................................................................35 Results and discussion................................................................................................38 Conclusion on the Shock Wave Interaction Analysis.................................................48 5 COUPLING OF CFD CODE WITH MCNP..............................................................49 Input............................................................................................................................49 Processing of FLUENTs Output...............................................................................51 Results from the Neutronic Calculations....................................................................54 Conclusion..................................................................................................................57 6 DESIGN ANALYSIS OF EXPANSION NOZZLE...................................................58 General Description of the Nozzle.............................................................................58 Meshing of the Nozzle................................................................................................59 Specificity of the Nozzle Case Inlet Boundary Condition and Initialization Step...61 Cases Considered and their Set-up.............................................................................61 Results and Discussion...............................................................................................62 Non-Adiabatic Case....................................................................................................65 Conclusion..................................................................................................................66 7 SUMMARY AND CONCLUSIONS.........................................................................71 Summary of Results....................................................................................................71 Future Work................................................................................................................72 APPENDIX A EVOLUTION OF DENSITY IN THE SHOCK TUBE.............................................74 B POWER DENSITY EVOLUTION IN THE SHOCK TUBE....................................77 LIST OF REFERENCES...................................................................................................80 BIOGRAPHICAL SKETCH.............................................................................................82 vi PAGE 7 LIST OF TABLES Table page 3.1: Properties of Natural Uranium gas at 10 bars and 5,000 K........................................24 4.1: Pressure ratios depending on the distance traveled in the gas when the pulse is generated with a pressure ratio of 4.........................................................................31 4.2: Extensive description of the settings for the shock waves interaction problem in the case of Uranium Tetrafluoride at an operating pressure of 10 bars.........................37 4.3: Description of the initialization step...........................................................................38 4.4: Relevant cases run during the shock tube study.........................................................39 4.5: Maximum ratios obtained when varying the height of the pressure pulses in Uranium Tetrafluoride.............................................................................................................42 5.1: Maximum power generation and its location at every time step................................56 6.1: Maximum Mach number reached in the nozzle for different back pressures.............65 vii PAGE 8 LIST OF FIGURES Figure page 1.1: Principle of Nuclear Thermal Propulsion.....................................................................3 1.2: Conceptual design for a high-pressure ionized shock generator pulsed magnetic fission reactor with magnetocumulative flux compression power generator.............4 2.1: Plot of the six residuals for 5 time steps. Each time step requires about 75 iterations.16 3.1: Curve of liquefaction of Uranium Tetrafluoride depending on temperature and pressure.....................................................................................................................19 3.2: Heat Capacity of UF 4 for the real formulation and the polynomial approximation...21 3.3: Viscosity of UF 4 for the real formulation and the polynomial approximation...........21 3.4: Thermal Conductivity of UF 4 for the real formulation and the polynomial approximation...........................................................................................................22 3.5: Comparison of density given by real formulation and by ideal gas formulation a) as a function of Pressure at T=3,000 K b) as a function of Temperature at P=50 bars.23 3.6: Shape of the pulse of pressure generated by a pulsed magnetic field as a function of time. Here the pressure ratio is 5..............................................................................24 4.1: Sketch of the two designs considered for shock duct. a) One side is closed b) Opened on each end.................................................................................................27 4.2: Grids generated with GAMBIT 2.0 and read by FLUENT 6.1 for the two geometries studied......................................................................................................................28 4.3A): Pressure wave generated on the left side of the tube propagates in the pipe..........29 4.3B): Pressure wave is reflected on the right closed side of the tube...............................30 4.4: Absolute pressure profile on a section of the pipe as a function of normalized position for different diameters................................................................................33 4.5: Density profile on a section of the pipe as a function of normalized position for different diameters....................................................................................................33 viii PAGE 9 4.6: Absolute Pressure profile near the wall when the pipe is 2 and 20 cm diameter.......34 4.7: Relative difference in pressure along the axis for 2 and 20 centimeter pipes............34 4.8: User Define Function used to set time dependent pressure pulses at inlet of the pipes.36 4.9: Variation of the Specific Heat Ratio on the axis of the shock tube............................40 4.10: Profile of normalized density, pressure and temperature of Uranium Tetrafluoride when the density reaches a maximum......................................................................41 4.11: Density profile on axis of the shock tube for different initial pressure ratios..........43 4.12: Evolution of density with time at the center of the shock tube.................................44 4.13: Density profile on axis of the shock tube for initial pressure pulses of different widths.......................................................................................................................45 4.14: Density normalized to the operating density on axis of the tube for different operating pressures...................................................................................................46 4.15: At time t = 15 micros, maximum overlapping, profile of properties on the axis of the shock tube a) Absolute pressure b) Density c) Temperature..............................47 4.16: At time t = 12.5 micros, maximum overlapping, profile of properties on the axis of the shock tube a) Absolute pressure b) Density.......................................................48 5.1: Example input for MCNP based on FLUENTs results.............................................51 5.2: Shock tube models used in the combined simulation codes. (a) Illustration of the shocktube model output by FLUENT with many cells. b) Illustration of the shocktube model after cell collapse by MHD2MCNP where many of the cells with nearly the same density/pressure are collapsed to simplify the MCNP model........52 5.3: Flow diagram of operating script linking FLUENT and MCNP................................53 5.4: Example output from MCNP......................................................................................55 5.5: Profile of density as it is input in MCNP and profile of power density as calculated by MCNP at t=1.659 ms...........................................................................................56 5.6: Evolution of total neutron flux in the tube with time.................................................56 5.7: Evolution of k eff with time..........................................................................................57 6.1: Design and boundary conditions for the expansion nozzle........................................59 6.2: Meshing of the expansion nozzle................................................................................60 ix PAGE 10 6.3: Refinement of the size of the grid in the throat of the converging-diverging nozzle.60 6.4: Absolute pressure on the exit plane of the nozzle when the boundary condition sets it to 10 bars..................................................................................................................63 6.5: Profile of Mach Numbers along vertical axis of nozzle for different exit pressures..64 6.6: Difference in Mach number profiles on the axis of the C-D nozzle for adiabatic walls and cooled walls.......................................................................................................66 6.7: Contours of absolute pressure in the nozzle (only the right part is represented) when the exit pressure is set to 10 bars..............................................................................67 6.8: Contours of Mach number when the exit pressure is set to 10 bars...........................68 6.9: Contours of temperature when the exit pressure is set to 10 bars...............................69 6.10: In the case of cooled walls (T w =1,800 K), contours of temperature when the exit pressure is set to 10 bars...........................................................................................70 x PAGE 11 NOMENCLATURE The following is a list of definitions of the main symbols used in this thesis. SI units are considered in the study. c p c v Heat capacity at constant pressure, volume J/(kg.K) h Enthalpy J/kg k Thermal conductivity W/(m.K) keff Effective conductivity composed of thermal conductivity of the material plus the turbulent thermal conductivity W/(m.K) P Pressure Pa q Volumetric heat source J/m 3 t Time s T Temperature K u, v, w Components of velocity m/s Specific heat ratio Dynamic viscosity Pa.s t Turbulent viscosity Pa.s Density kg/m 3 Shear stress Pa F External body forces xi PAGE 12 Tv Effects of volume dilatation I Unit tensor J Diffusion flux of the gas. xii PAGE 13 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science A COMPUTATIONAL THERMAL-FLUID MODELING APPROACH TO PULSED AND STEADY GAS CORE REACTORS By Anne Charmeau May 2004 Chair: Samim Anghaie Major Department: Nuclear and Radiological Engineering Steady and pulsed Gas Core Reactors (GCR) are considered for terrestrial and space power generation. The main characterizing feature of this kind of reactor is its fuel that is in gaseous or vapor of uranium fluorides. Steady Gas Core Reactors are among concepts considered for Generation IV nuclear power systems. The Shock Wave Driven Gas Core Reactors are proposed for very high power density Nuclear Electrical Propulsion applications that are intended to power spacecrafts for manned or robotic exploration of the solar system. In Shock Wave Driven Gas Core Reactors, criticality is achieved by constructive interference of two opposing shock waves. The main goal of this research is to develop a Computational Fluid Dynamics (CFD) model for steady and pulsed Gas Core Reactors. The GCR CFD model is used to analyze the interaction of the two counter approaching pressure pulses in a shock tube that simulates the Gas Core Reactor. The main reason for the interference of the shock wave is to achieve fissile gas densities in the reactor core that is several times larger than xiii PAGE 14 the ambient gas density. The significant increase in the fissile gas density is expected to ignite prompt fission criticality, which will be followed by a large pulse of energy release. The CFD model also is used to simulate fissile gas flow in the nozzle region of steady GCR. The FLUENT CFD code is used to develop the GCR model. Real fissile gas properties such as uranium tetrafluoride (UF 4 ) are incorporated in the FLUENT code to simulate the evolution of the flow, density, pressure and temperature in both steady and pulsed GCRs. Results of CFD analysis in pulsed GCR indicate that the attenuation of the approaching shock waves are larger than expected. The second part of the analysis involves fissile gas flow through a two-stage nozzle of the steady GCR. The two-stage nozzle is used to convert the enthalpy to very high kinetic energy that is need for power conversion in magnetohydrodynamic (MHD) power generator. Results of the CFD analysis show that the fissile gas approaches the chocking condition (Mach number equal one) at the throat of the second stage nozzle. This confirms the validity of the nozzle design for the steady GCR. xiv PAGE 15 CHAPTER 1 INTRODUCTION Introduction New concepts of nuclear power generation under investigation involve high temperature and high-pressure cores, as well as non-traditional fuel and coolant. One example is the Gas Core Reactor (GCR). The fuel in a GCR remains gaseous through the whole cycle of the core. Typically, the fuel is enriched uranium gas, Uranium Tetrafluoride (UF 4 ) for example. New tools are utilized to predict the behavior of fluids in these types of reactors. One of those tools is the Computational Fluid Dynamics (CFD). CFD codes are used to solve numerically the fundamental governing equations of thermodynamics and heat transfer. For the research presented in this thesis CFD has been applied to two concepts of nuclear power generation: Gas Core Reactors and Shock Wave Driven Gas Core Reactors (SWD-GCR). Gas Core Reactors are among concepts developed for generation IV nuclear power systems. Shock Wave Driven Gas Core Reactors are an application of the GCR technology to nuclear space power. This introductory chapter discusses briefly the context of nuclear propulsion for space rockets as well as the principle of SWD-GCR. History of the Gas Core Nuclear Rocket Why Choosing a Propulsion based on Nuclear Energy Chemical rockets use the heat produced by burning fuel and oxidizer in a trust chamber, which requires a huge amount of fuel and represents a lot of space and weight. In the mid 60s, American government and industry looked for a cheaper, reliable 1 PAGE 16 2 alternative to chemical rocket engines. It was in the era of President Kennedy when the Cold War and the race to the moon were at their peak. Project Rover, the first project of this kind, aimed at building a nuclear reactor to power a rocket in space. The concept of nuclear thermal propulsion-NTP-is simple: the exiting gas of a solid-core hydrogen cooled reactor expands through a rocket nozzle to discharge in space, producing massive thrust, as sketched on Figure 1.1. The motivation for development of this kind of propulsion is that it would create twice as much specific impulse 1 as for the best chemical rocket, with a reduction by a factor of 5 in the ratio of take-off mass to final mass. The project terminated in January 1973, but saw the creation of several advanced designed engines and was appraised as a technical success. Some new studies on nuclear thermal propulsion were lead during the past several years since the next giant leap for mankind will be human exploration of Mars. Missions to Mars will be longer than usual ones and thus astronauts will be exposed to greater amounts of radiation in space. To face up to these issues, the travel time needs to be shortened, and more shielding added to the rocket. In the case of a trip to Mars, the nuclear rocket would provide a 200-day trip instead of 600 days with a traditional chemical rocket. One of the possibilities studied is that the nuclear rocket engine would accelerate for most of the trip and then decelerate until it reaches Mars: turning the rocket would provide the deceleration as thrust would be used as a break. The next generation of space propulsion is the Nuclear Electrical Propulsion NEP. In the case of this futuristic propulsion electricity is produced by a Shock Wave 1 Specific impulse: thrust per unit flow rate of propellant PAGE 17 3 Driven Gas Core Reactor. The principle of a gas heated by fission processes is abandoned. Figure 1.1: Principle of Nuclear Thermal Propulsion Shock Wave Driven Gas Core Reactor The idea behind Shock Wave Driven Gas Core Reactor (SWD-GCR) is that a set of strong ionizing twin shock waves would be generated in a shock tube filled with fissile gaseous fuel, Uranium Tetrafluoride for example [1]. These ionizing waves would overlap in the vicinity of a neutron moderator material. The collision would result in a brief explosive burst of energy ionizing the fissile gas. Solenoidal coils would surround the chamber where the collision occurs, and the ionizing wave would be strong enough to create magnetoculmulative flux-compression power, leading to extraction of pulsed electrical power. The entire system involves three coupled stages: PAGE 18 4 shock generation shock interaction and fission release magnetoculmulative flux-compression generator power extraction Figure 1.2 shows a conceptual design of the core of the Shock Wave Driven Gas Core Reactor. The central chamber, surrounded by the BeO reflector is where the shock waves collide, the energy is released and the plasma is created. The present investigation considers this chamber to be 100 cm long and 80 cm radius. Shock interaction is the main focus of the research presented in this thesis. However the scope is targeted on the aeroacoustic of the gas, and the fission release process is not studied. This study will not discuss power extraction either. Reactors of this type have never been prototyped or even tested on a small scale owing to the exceeding stringent materials requirements, analysis can only be done by modeling the system and using appropriate codes. Figure 1.2: Conceptual design for a high-pressure ionized shock generator pulsed magnetic fission reactor with magnetocumulative flux compression power generator Shock Wave Driven Gas Core Reactor is a futuristic concept brought into serious considerations for future space power applications. The main advantage of this core design is that it does not use destructive explosives and is intended to be a fully PAGE 19 5 functioning power cycle. Also, it appears that a great deal of effort is required to lead the core to criticality; this system is therefore quite safe from a neutronic point of view. Viability of Shock Wave Driven Gas Core Reactor concept depends on the ability and efficiency of creating strong pressure pulses on both sides of a duct with such large dimensions. The idea under investigation is to generate these pressure pulses as shock waves, created by a pulsed magnetic field. Pulse generation is still one of the biggest challenges of the SWD-GCR conception. Thesis Objectives The goal of this analysis is to simulate the flow in several parts of cores involving gaseous fuel. To model the different flows, the commercial CFD code FLUENT 6.1 has been run. One of the aims of the thesis is to present the approximations, assumptions and models used to lead the different studies. The thesis mainly focuses on presenting the results of two studies that are described below. Interaction of shock waves : In Shock Wave Driven Gas Core Reactors pressure pulses are generated on both sides of the shock tube. Behavior of UF 4 is unknown in the conditions for which the core is run and when the two shock waves collide. Thanks to the CFD simulation we get a better idea of the properties of the gas when the pressure pulses overlap, and then predict the conditions in which the core becomes critical. This study is challenging since it is a time-dependant problem presenting a huge pressure discontinuity. Expansion nozzle : In the design of Gas Core Reactors currently under investigation the GCR is coupled to a high efficiency MHD-generator. Coupling of the two components is done through an expansion nozzle. Design of the nozzle is performed modeling the flow with the CFD tool, and making sure that the gas does not choke anywhere in the system. Thesis Organization A general overview of Computational Fluid Dynamics is given in Chapter 2. Chapter 3 presents the assumptions and approximations used during the research on the modeling and property data. Simulation of the interaction of pressure waves in the shock PAGE 20 6 tube core is described in chapter 4. The results obtained in Chapter 4 are then coupled to a neutronic code (MCNP) to get an evaluation of the energy released by the shock wave collision. The first results of that coupled analysis are compiled in Chapter 5. Design analyses of the expansion nozzle linking gas core reactors to MHD-generator are included in Chapter 5. Finally, Chapter 6 presents a summary of the results as well as some recommendations for future work. PAGE 21 CHAPTER 2 COMPUTATIONAL FLUID DYNAMICS What is Computational Fluid Dynamics? Generally, a flow can be described as solving the Navier-Stoke equations derived from the four conservation equations: Conservation of mass Conservation of momentum Conservation of energy Equation of state When geometries start being complex and the flow cannot be considered one-dimensional, Navier-Stokes equations become highly non-linear. They are a real challenge since there is no analytical solution and computational methods are required to solve thermodynamics and heat transfer problems. Codes that solve the conservation equations for flows are called Computational Fluid Dynamics codes, referred as CFD codes. With advances in computer technology, CFD has become a common tool for equipment design, R&D and process optimization. As a numerical tool, CDF codes solve the conservation equations for microscopic regions of the flow defined by a grid, using finite differences methods. A given CFD model contains a grid or geometry, physical models (turbulence, multiphase, etc.) and boundary conditions. The domain is divided into discrete control volumes using a computational grid. The governing equations are integrated on the individual control volume to construct algebraic equations for the unknowns. The governing equations are then linearized to update the values of the 7 PAGE 22 8 dependant variables or unknowns. This chapter presents the governing equations and differences between solvers and models used by CFD codes. Governing Equations To describe the flow properly, the following four intensive variables have to be known at every fluid location: Velocity v Density Pressure P Temperature T In this thesis, the Cartesian two-dimensional domain is considered. Flows are described in the (x,z) plane. The core is supposed to be run in conditions such as fluid always remains in the gaseous phase. Hence, we can consider only single-phase flows and their associated equations. Mass Conservation Equation Continuity Equation The general form of the continuity equation is given by 0.v t It can be reduced to 0)()( z w x u t (Eq. 2.1) Momentum Conservation Equation Newtons Second Law The general form of Newtons second law is given by FPvvtv .. As discussed in Chapter 3, gravity is not considered in the cases studied: there are no external body forces 0 F and the general equation can be reduced to the xand z-components: PAGE 23 9 zxzPxwwxwutwzxxPzuwxuutuzzxzzxxx (Eq. 2.2) The stress tensor is given by IvvvT.32 (Eq. 2.21) Energy Conservation Equation First Law of Thermodynamics vJhTkPEvEteffeff ..))(.()( (Eq. 2.3) Components of the deviatoric stress tensor, energy E and enthalpy h are defined respectively by ijiieffjiijeffijxuxuxu32 (Eq. 2.31) 22vPhEdTcdhp (Eq. 2.32). (Eq. 2.33) When equation (2.33) is integrated, the temperature of reference considered is usually 298.15 K. Equation of State In equations 2.1, 2.2 and 2.3, there are five unknowns, pressure, temperature, density, x-velocity and z-velocity but only a set of four equations. The fifth equation is not a conservation equation, but the equation of state that links density to intensive parameters of the domain. Typically, the equation of state can be a constant if the flow is considered incompressible ( = constant ), an ideal-gas law if the PAGE 24 10 behavior of the fluid can be considered perfect ( = P/RT ), or any other empirical relation defined between density, pressure and temperature. Solver Formulation Solvers in CFD codes have thre different formulations: Segregated Coupled implicit Coupled explicit The three solvers give accurate results for a broad range of flows but one formulation may perform better for a given case. Segregated vs. Coupled Segregated approach solves for one single variable field (e.g., pressure, x-velocity) by considering all cells at the same time whereas coupled approach solves for all variables at a time. Implicit vs. Explicit If the coupled solver is selected, equations can be linearized two ways: using an implicit or an explicit scheme. With an implicit scheme, for a given variable, the unknown variable in each cell is computed using both existing and unknown values from neighboring cells whereas with an explicit scheme, the unknown is computed using only existing values. Hence, with the implicit method, the unknown appears in more than one equation and all the equations are solved simultaneously, but with the explicit method, the unknown quantity appears in only one equation and the equations can be solved one at a time to give all the unknowns. Turbulence Modeling Flows considered all through this work have a high Reynolds number, meaning that all the problems presented in the thesis consider fully developed turbulent flows. One of PAGE 25 11 the most difficult hurdles to the proper use of CFD codes is the turbulence modeling. In a turbulent flow, particles of fluid move unsteadily in an unpredictable path. Unfortunately, there is not a single universal turbulence formulation to model all turbulent flows. In most CFD codes different turbulence models are proposed. Actually, a distinct turbulence theory may need to be developed for each geometry, as it depends on the fluid type and properties, the geometry and other parameters. Turbulent flows are characterized by fluctuating velocity fields. These fluctuations can be of small scale and high frequency and then become too computational time consuming when simulated in practical engineering calculations. For the research presented here, the Spalart-Allmaras model and the standard kmodel have been used. They are described below. Other viscous models are inviscid and laminar (1 equation) for laminar flows and k-omega (2 equations) and Reynolds stress (4 equations) for turbulent flows. Spalart-Allmaras Model Spalart-Allmaras model is a relatively simple one-equation model [2] that solves a simplified transport equation for the kinematic eddy viscosity. In the Spalart-Allmaras model, the equation is generated from scratch using empiricism and arguments of dimensional analysis, Galilean invariance and selective dependence on the molecular viscosity. This model allows saving computational time for problems where the local shear layer thickness is small compared to the dimensions of the system. The Spalart-Allmaras model was designed specifically for aerospace applications involving wall-bounded flows [3]. On a cautionary note, however, the Spalart-Allmaras model is still relatively new, and no claim is made regarding its suitability to all types of complex engineering flows. PAGE 26 12 Standard kModel Two-equation models are the most widely used in engineering simulations [4-5]. The k-epsilon model is used as a better alternative to the one-equation Spalart-Allmaras model. K-epsilon model is a very common turbulence model that has been proven to give good results without requiring too much computational time. This model, first proposed by Jones and Launder [6], solves the transport equations of turbulent kinetic energy (k) and its dissipation rate () based on the classical balance equation: Rate of change + Convection = Diffusion + {Generation Destruction} As it is explained case by case in the chapters that follow, Mach numbers during this research were relatively low, smaller than 3, and then it can be considered that turbulence has the same characteristics as incompressible turbulence [7]. Hence, only the k-epsilon model for incompressible flow will be discussed here, adding terms to the initial equations would treat incompressibility. K-epsilon model starts from the Reynolds hypothesis, which postulates a linear relation with the symmetric part of the velocity gradient. The transport equations for the standard k-epsilon model are kCGCGkCxxDtDGGxkxDtDkbkitibkikti2231 (Eq. 2.4) ijjikuuuuG'' is the production of turbulent kinetic energy and ittibxTgGPr is the generation term due to buoyancy. 1C Cand Care constants that can be determined experimentally. 2 3 PAGE 27 13 The main difficulty when using a k-epsilon turbulence modeling is to choose boundary conditions and the initial values for k and This issue is detailed in Chapter 4, when presenting the shock tube case. FLUENT 6.1 For several years INSPI has been developing expertise in the domain of CFD using the commercial code FLUENT developed by Fluent, Inc Lebanon, NH. Today, FLUENT is widely used and its use is only increasing in the nuclear industry and National Laboratories (General Atomics, Westinghouse, INEEL). It has proved to give good results in past in-house studies as well as in industry, on problems like modeling complex gas reactor coolant systems. For the previous researches FLUENT 4 was used but this time finer study is required, then the use of FLUENT 6, and its latest release 6.1 (June 2004) is necessary. FLUENT is a state-of-the-art computer program for modeling fluid flow and heat transfer in complex geometries. It solves problems for unstructured meshes that can be generated using the software GAMBIT. This code is detailed later in this chapter. All functions to set up a case or extract the results of calculations are accessible through a user interface. Grids created with gambit are read in FLUENT. The grid can be checked by FLUENT to make sure that all the volumes created are positive. Then the user defines the models for solver, energy, viscosity or turbulence, radiation. Operating conditions and boundary conditions need to be specified, as well as the properties of the different fluids present in the system. For a better control of solution different algorithms can be selected and the level of convergence desired by the user is selected. Before starting to iterate, properties of the flow field such as ambient pressure, velocity temperature and viscosity parameters are initialized from a specified boundary. Once the solution has converged PAGE 28 14 FLUENT allows plotting contours and profiles of all properties. Profiles can also be saved in a text file where the value of the property is given for each spatial node. Creating Grids and Meshing Although discretization of the partial differential equations is powerful, it has important limitations concerning stability. The finer the mesh is the less stability problems occur. However, a larger grid requires a higher Computational Time (CPU). Attention will be first focused on the building of the grid, finding a compromise between good stability and fast calculations. Since the meshes in the study are slightly complex, an external code, GAMBIT 2.0 is used to generate grids. Creating a grid with GAMBIT is done in three steps: Defining the geometry by creating volumes, surfaces, edges and points. Creating a structured grid. The choice of the number of meshes, type and sizing functions is up to the user. Meshes can be quadrilateral or triangular. Specifying boundaries and continuum zone types depending on the solver chosen (in our case, FLUENT 5/6). Once the mesh is created it is saved as a binary file and can be read by FLUENT. User-Defined Functions (UDFs) A user-defined function, or UDF, is a programmed function that can be dynamically loaded with the FLUENT solver to enhance the standard features of the code. UDFs are written in the C or C++ programming language. They are defined using macros that are supplied by Fluent Inc. It is for example used to set a time-dependant boundary condition or a material property as a function of temperature or pressure for example. PAGE 29 15 Solution Controls To solve the conservation equations, FLUENT uses a control-volume-based technique. By default, FLUENT stores discrete values of the scalar phi at the cell centers. However, face values phif are required for the convection terms and must be interpolated from the cell center values. This is accomplished using an upwind scheme. Upwinding means that the face value phif is derived from quantities in the cell upstream, or upwind, relative to the direction of the normal velocity. FLUENT allows choosing from several upwind schemes: first-order upwind, second-order upwind, power law, and QUICK. Before running a case, it has to be determined which scheme is the better adapted to it. Convergence Check-Residuals For every iteration, the residual for continuity, components of velocity, energy and turbulence parameters is calculated. The residual for each variable is the difference between the value calculated at iteration N-1 and iteration N divided by its physical time step. FLUENT plots the evolution of the residuals vs iteration number as showed on Figure 2.1. As long as the residuals are greater than the convergence criterions the code keeps iterating. The convergence criterions can be decreased to obtain more accurate results but it also increases dramatically the CPU time. Convergence criterions are typically of 10 -3 but can be lowered to 10 -6 in the case of the analysis of the expansion nozzle. PAGE 30 16 Figure 2.1: Plot of the six residuals for 5 time steps. Each time step requires about 75 iterations. PAGE 31 CHAPTER 3 MODELING AND PROPERTY DATA In the analysis described all dimensions are expressed in SI units. General Assumptions Two-Dimensional Representation of the Tube and Nozzle Both tube and nozzle are rotationally symmetric. To save computational time, the study considers a two dimensional section of the systems. The 2D version of FLUENT 6.1 is run. Gravitational force In the whole study, the influence of gravitation is neglected. This assumption is pretty accurate since the core is run in a non-gravitational environment. In the Operating Conditions panel of FLUENT the Gravity option is turned off. Adiabatic Walls This study is performed in order to get a global idea of the phenomenon that can happen when the experiment is run. In a first approach of the problem the walls are considered adiabatic, so heat transfer between walls and fluid do not interfere. Since the inlet gas is supposed to be at 3,000 K, the first idea was to set the wall at constant temperature, i.e. 3,000 K. The losses of energy, measured, as total heat transfer through the walls are not negligible. In the case of the expansion nozzle, described in more details in Chapter 5, the total heat transfer rate reaches 133 kW. 17 PAGE 32 18 A heat flux of zero on all the walls is then considered. The heat transfer rate is reduced down to 194 W. This loss of heat can be considered negligible: the total length of the walls is more than 10 meters and the heat transfer at inlet is more than 2 GW. In the study, adiabatic walls are modeled when the heat flux through the walls is set to zero. Assumptions on the Fluids Used in the Study The core of the GCR uses two main fluids: Helium gas and Uranium Tetrafluoride (UF 4 ) gas. Basic functions of FLUENT allow the user to set properties of any material either constant, using referenced values, or as a function of one extensive property. Another function of FLUENT is the User-Defined Real Gas Model (UDRGM) that is incorporated in the algorithm to solve the Navier-Stokes equations described Chapter 2. The UDRGM enables to describe properties of a material as a function of several extensive properties and should be written in c++ like any other User Defined Function. However, depending on the gas considered the default properties given by FLUENTs library can be very accurate. The assumption regarding the modeling of each gas is given below. Helium Helium does not play an important role in the criticality process. The aim is just to get a global idea of its behavior in systems where there is uranium tetrafluoride. The ideal gas law is selected to model density, otherwise, properties are set using values of FLUENTs library [8]. Uranium Tetrafluoride Properties of UF 4 were measured and compiled in Russia in the late 70s [9]. Correlations to describe the evolution of density, specific heat ratio, viscosity, thermal PAGE 33 19 conductivity and heat capacity of uranium tetrafluoride were derived at University of Florida in 1993 [10]. For UF 4 liquid phase boiling the following equations are adopted: )()(44gasUFliquidUF for T < 1600K TTatmPln0.7/37977217.74)(ln for T > 1600K TTatmPln05.7/3845388.74)(ln Figure 3.1 gives the limit between gaseous and liquid phase of UF 4 for temperature higher than 1,600 K. The fluid at the exit of a Gas Core Reactor has been evaluated to be at a temperature of 3,000 K. In a gas core reactor, pulsed or steady, the fissile gas should remain under a gaseous state. If the gas is considered to remain at a temperature near 3,000K the core should not be run for pressures higher than a hundred bars. This is in agreement with the preliminary designs of the expansion nozzle and the shock tube. 0.1110100100016001800200022002400260028003000Temperature (K)Pressure (atm) LIQUID VAPOR Figure 3.1: Curve of liquefaction of Uranium Tetrafluoride depending on temperature and pressure. PAGE 34 20 The different properties are given for monolithic UF 4 gaseous phase. It is considered that UF 4 is composed of one atom of natural Uranium and four atoms of fluorine. The molecular weight considered in Fluent is then M = 314.0225 g/mol. Correlations used in all the studies led at INSPI are the following: ).(*393.10))./((*2.3))./((10*04.267.350010*10.780.3500*10*357.3).(/10*06.3*10*24.25.121))./(()(77.13)(10*83.3)/(55639323sPamolKJcKmWkbaKTandbaKTwithbTaTsPaTTmolKJcTbarsPTPaPmkgpp (Eq. 3.1) Note: in FLUENT, heat capacity has to be expressed in J/(kg.K). FLUENT allows the user to set temperature dependant properties of a new material. From equations (3.1) it can be seen that heat capacity, viscosity and thermal conductivity are only temperature dependant. The cases are run at a temperature of 3,000 K and the variations of temperature in the system are small, since we consider adiabatic walls. To accelerate the calculations, we want to find out polynomial relations between the different properties and temperature for temperatures ranging from 2,500 K to 3,500 K. After running each case, we need to make sure that the temperatures obtained are in this range. Figures 3.2, 3.3 and 3.4 present the plots of the real evolution and the polynomial regression of respectively heat capacity, viscosity and thermal conductivity. PAGE 35 21 40440540640740840941041141241325002600270028002900300031003200330034003500Temperature (K)Heat Capacity (J/kg.K) Correlation Approximation Figure 3.2: Heat Capacity of UF 4 for the real formulation and the polynomial approximation 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-044.0E-0425002600270028002900300031003200330034003500Temperature (K)Viscosity (Pa.s) Correlation Approximation Figure 3.3: Viscosity of UF 4 for the real formulation and the polynomial approximation PAGE 36 22 00.020.040.060.080.10.120.140.160.1825002600270028002900300031003200330034003500Temperature (K)Thermal Conductivity (W/m.K) Correlation Approximation Figure 3.4: Thermal Conductivity of UF 4 for the real formulation and the polynomial approximation The polynomial approximations used are: (Eq. 3.2) Error on the values obtained with the approximations is small, no more than 1.7%, in the range of temperature considered. A linear variation of properties with temperature gives results accurate enough. TKmWkTsPaTKkgJcp585310*937.3020824.0))./((*10*910*4).(10*13.791.386))./(( Density is a function of both temperature and pressure. The cases were first run with the real formulation for density. In doing so, FLUENT ends up with a specific heat ratio in the system of 1.27. Now, the study led at INSPI [10] showed that the gamma of UF 4 should be about 1.08 in the range of temperature and pressure that we consider. This is the main feature of interest of Uranium Tetrafluoride. The real formulation for density in FLUENT raises some problems. The same study was then led using an ideal-gas PAGE 37 23 formulation for density. Once the change has been made, the specific heat ratio of the gas in the system matches the experimental values. The difference between the two formulations is shown Figure 3.5. The error when using the ideal gas density is 1.39%, which can be considered negligible. The ideal gas formulation for density is given by )(*)./()()/(3KTKkgJRPaPmkg (Eq. 3.3) where R is the gas constant for UF 4 )./(4768.26KkgJR a) 020406080010203040Pressure (bars)Density (kg/m3) 50 real formulation ideal gas b) 6062.56567.5702800290030003100Temperature (K)Density (kg/m3) real formulation ideal gas Figure 3.5: Comparison of density given by real formulation and by ideal gas formulation a) as a function of Pressure at T=3,000 K b) as a function of Temperature at P=50 bars. PAGE 38 24 Uranium gas In some preliminary studies, Natural Uranium gas is used in a first approximation. Its properties are kept constant in the range of temperature and pressure of the study and equal to the properties of Natural Uranium at 10 bars and 5,000 K [11]. They are given in table 3.1. Table 3.1: Properties of Natural Uranium gas at 10 bars and 5,000 K c p (kJ/(kg.K)) 0.2316 k (W/(m.K)) 0.5456 (mP) 1.6971 M (g/mol) 238.029 Nature of Pressure Pulses Generated Generation of pressure pulses has been described in the introduction, Chapter 1. The assumption was made that a step function, pressure vs. time, can represent the pulses generated at inlets of the shock tube. Pulses are expected to be 5 mm width, which correspond to a pulse time of approximately 5 s. Shape of pressure pulses is given Figure 3.6. Uniform profile of pressure on the inlet section is also assumed. 01234560.0E+005.0E-061.0E-051.5E-052.0E-052.5E-053.0E-05time (s)Pressure Figure 3.6: Shape of the pulse of pressure generated by a pulsed magnetic field as a function of time. Here the pressure ratio is 5. PAGE 39 25 In the study, it will often be referred to the pressure ratio. It is the ratio of the absolute pressure at a given point of the tube over the operating pressure. Summary of Approximations and Assumptions Two-dimensional mock up of the systems. Adiabatic walls. Gravitation negligible. Gases properties: Helium: Ideal gas law, other properties from Fluents library. Uranium Tetrafluoride Gas: Ideal gas law, other properties linearly dependent on temperature. Natural Uranium Gas: Constant properties. Pressure pulses described by a step function. PAGE 40 CHAPTER 4 STUDY OF PRESSURE PULSES EVOLUTION IN SHOCK TUBE Feasibility of manufacturing Shock Wave Driven Gas Core Reactors will depend on the ability to get the fissile gas to be critical. The condition when gas will ignite prompt fission supercriticality is reached if there is a sudden increase by order of magnitude of the density of the gas located near the BeO reflectors. This chapter studies the evolution of properties of the gas in the shock tube when subjected to pressure pulses. The real challenge of this part of the analysis is to capture pressure properly in the tube since there are huge pressure discontinuities upstream and downstream from the shock wave. This problem also requires a time-dependent treatment. The scope of this part of the study is to prove that when pressure pulses collide, the density of the fissile gas, uranium tetrafluoride in the case considered, reaches high values. Design of the Core The analysis considers two types of designs sketched, Figure 4.1.a and 4.1.b. In the first case, the tube is closed on one side, and the pressure pulses are only generated on the opened side. The first pulse is reflected on the opposite wall and the overlapping of the reflected pressure wave with an incident wave is studied. In the second case, a tube is opened on both ends, and identical pressure pulses are generated from each side of the tube. The phenomenon of interest is the behavior of the gas when the two pressure pulses overlap. 26 PAGE 41 27 Figure 4.1: Sketch of the two designs considered for shock duct. a) One side is closed b) Opened on each end Fluid Used The shock tube in the Shock Wave Driven Gas Core Reactor can be considered as the core of the reactor. It is filled with Uranium Tetrafluoride. The interesting property of UF 4 is that its specific heat constant is particularly low 1.08 compared to the gamma of 1.33 for a triatomic gas. Behaviors of natural uranium gas not enriched and helium have also been investigated. Creation of the Grids Grids are created for the two geometries using the software GAMBIT 2.0. The mesh is refined in the area where the phenomenons of interest are expected to occur: when the shocks are generated from each side of the tube, the area of interest is in the middle of the tube. If the tube is closed at one end, the reflection of the wave is the interesting property, meshing next to the wall is refined. The final grids are presented Figure 4.2. PAGE 42 28 Figure 4.2: Grids generated with GAMBIT 2.0 and read by FLUENT 6.1 for the two geometries studied. Preliminary Studies Quick Validation of the CFD Code FLUENT in Capturing Pressure. The very first step of this study is to model a pressure tube closed on one end, pressure waves are generated continuously on only one side of the tube referred previously as geometry a). The tube modeled is 10 cm long and 2 cm diameter. Contours of pressure in the tube are plotted and its evolution with time is recorded. With this study we want to check that FLUENT 6.1 is adapted to capture pressure when shock wave reflects on a solid wall or when shock waves interact. A similar study was led at Fluent, PAGE 43 29 Inc before the release of Fluent 6.1. 1 Operating pressure is set to 10 bars. A single continuous shock wave defined by a pressure of 20 bars is implemented at inlet, on the left part of section of the tube. The time-step is set to 0.1 s, and we look at the evolution of absolute pressure of the gas along the axis of the tube given Figure 4.3. The pressure wave gets attenuated when it travels in the gas, but FLUENT can handle well the propagation of the wave and its reflection on the wall. Results are similar from what was expected. From this study, we can consider that FLEUNT is an acceptable tool to capture pressure waves and shocks in the conditions of the analysis. A Figure 4.3A): Pressure wave generated on the left side of the tube propagates in the pipe. 1 Dr Samir Rida, Fluent Inc., private conversation. PAGE 44 30 Figure 4.3B): Pressure wave is reflected on the right closed side of the tube. Length of the Pipe The original design sets the length of the pipe to about 1 m. When the case is run for the first time with Uranium Tetrafluoride or Helium, the first noticeable fact is the large attenuation of the pressure when the pulse propagates in the tube. This can be seen for example Figure 4.3 A) when a continuous pressure raise is implemented on the inlet. The attenuation of the intensity of the generated pulse is even more sensitive if the pressure discontinuity is a single pulse of a few microseconds width. This dissipation is so important that when using a 10 cm long tube, after traveling 4 cm in the gas pressure pulses with a 5 s width are reduced to a tenth of their initial value. Table 4.1 lists the attenuation value recorded when pulses are generated on both ends of a 10 cm long 2 cm diameter pipe. The process of attenuation is more sensitive at the inlet of the shock tube. It can be interfered from Table 4.1 that it is not required to model the whole length of the pipe PAGE 45 31 since the pressure pulse is attenuated in the first millimeters after its generation but not so much during its propagation. To save computational time, only the central part of the tube is meshed. Table 4.1: Pressure ratios depending on the distance traveled in the gas when the pulse is generated with a pressure ratio of 4. Distance traveled (mm) Uranium Tetrafluoride Helium 0 2 2.9 10 1.4 2.3 45 1.3 2.1 Diameter of the Pipe Radius of the pipe should be about 80 cm. However, meshing a 1.6 m diameter pipe requires either a lot of cells thus a long computational time. One can choose to mesh cells of big size but this reduces the precision on the solutions obtained. This part studies the influence of the diameter of the pipe on the pressure and density distribution at a given section of the pipe. First, a tube of 5 mm long is considered and its diameter is changed using the scaling function of FLUENT. Diameters of 2, 5, 10 and 25 cm are studied. When the scaling function is utilized, the size of the geometry is changed by stretching or shrinking the cells. The number of cells remains the same for all the cases run. The absolute pressure profile at the center of the tube is given Figure 4.4 for the different radiuses considered. The absolute pressure is given as a function of normalized position, which is obtained by dividing the diametrical position by the diameter of the pipe. We consider the pressure profile at the time when there is a maximum overlapping of the pressure pulses. Pressure profile on a section of the shock tube is uniform, with a slight pressure drop near the wall. The extent of the pressure drop is comparatively the same for the different sizes of the pipes. The absolute value of pressure decreases when the size of the PAGE 46 32 pipe increases. However, the decrease is not really significant since there is less than 1% difference between the 2 cm and 25 cm diameter pipes. The exact same tendency is observed for any other property including density as can be seen on Figure 4.5. In the cases compared, the number of meshes is the same, only their size changes to get a larger or smaller pipe. In a next step, pressure profiles are compared for two different diameters, but with meshes of the same sizes in the two pipes. The pipes considered have diameters of 2 and 20 cm and the grids have 2,400 cells and 242,501 cells respectively. Profiles of absolute pressure are computed on the axis of the tubes and near the walls. Profiles near walls are given Figure 4.6. Radial size of the pipe does not affect the pressure distribution on the axis of the shock tube. The difference in absolute pressure obtained for the two pipes is showed Figure 4.7. Meshing the real size shock tube requires a lot of cells and thus a longer computational time. The results presented previously show that the diameter of the pipe does not influence the pressure profiles and distributions in the area of interest. Hence, instead of meshing a 1 m long and 80 cm radius core, only a portion of 2 cm diameter pipe with a length varying from one to several centimeters is considered in the analysis presented in the chapter. PAGE 47 33 4.66E+064.67E+064.68E+064.69E+064.70E+064.71E+064.72E+064.73E+064.74E+064.75E+0600.10.20.30.40.50.60.70.80.91Normalized positionPressure (Pa) d=2cm d=5 cm d=10 cm d=25 cmt Figure 4.4: Absolute pressure profile on a section of the pipe as a function of normalized position for different diameters. 53.25 Density (kg/m3) 53 2 cm 5 cm 10 cm 52.75 0.1 0.2 0.8 0.9 0 0.3 0.4 0.5 0.6 0.7 1 Normalized position Figure 4.5: Density profile on a section of the pipe as a function of normalized position for different diameters. PAGE 48 34 1.0E+062.0E+063.0E+064.0E+0600.10.20.30.40.50.60.70.80.91Position (mm)Pressure (Pa) 20 cm 2 cm Figure 4.6: Absolute Pressure profile near the wall when the pipe is 2 and 20 cm diameter. 1.E-041.E-031.E-021.E-011.E+0000.511.522.5Position (mm)Relative difference in pressure Figure 4.7: Relative difference in pressure along the axis for 2 and 20 centimeter pipes. Time Step Pressure shock waves are generated during 5 s. The time step to describe time dependant behavior in the pipe should be at least one order of magnitude smaller. When running the problem with an implicit time dependant formulation, the typical time step PAGE 49 35 chosen by the code is 5*10 -8 sec. In the study presented here, we choose an explicit time formulation so every time step is the same, but the explicit formulation gives an idea of the typical time step required to solve the cases considered. Typically, time steps are chosen to be 0.1 s when initializing the problem or when the wave gets close to an interaction and 0.5 s if the wave only propagates. Conclusion of the Preliminary Study In this preliminary study, we showed that Fluent seems to be appropriate to get a global description of the gas properties when shock waves propagate and interact in Shock Wave Driven Gas Core Reactors. We also showed that it is not necessary to mock-up the shock tube in its whole, but only a small portion of it is enough to describe properly the density evolution of the fissile gas. Last, we figured that time steps of tenth of a microsecond are well appropriate in areas of interest to treat the shock wave interaction problem. Set-up of the Cases Set-up of Inlet Boundary Conditions: Pressure Pulses As mentioned in Chapter 3, the pressure pulses at inlet of the tube can be described by a uniform time dependent step function. To set up this boundary condition, a User Defined Function -UDFis interpreted. The UDF is a small program written in C++ as shown Figure 4.8. In the analysis, the time-dependant profile of pressure is called press_pulse. In the menu Define/Boundary_Condition/inlet/pressure-inlet, Gauge Total Pressure and Supersonic/Initial Gauge Pressure are set to the UDF press_pulse. Inlet temperature is chosen to be equal to 3,000 K. PAGE 50 36 #include "udf.h" #define P0 3.e6 /* relative height of pressure pulse */ #define step 5.e-6 /* spatial width of the pressure pulse */ DEFINE_PROFILE(press_pulse, t inlet ) { double time, pressure; face_t f; begin_f_loop(f, t) { time=RP_Get_Real("flow-time"); if((time>=0)&&(time PAGE 51 37 Table 4.2: Extensive description of the settings for the shock waves interaction problem in the case of Uranium Tetrafluoride at an operating pressure of 10 bars. MODEL Solver Solver Coupled Formulation Explicit Time Unsteady Unsteadyformulation 2nd-orderImplicit Space 2D Energy EnergyequationON Viscous Model k-epsilon k-epsilonmodel Standard Near-Wall Treatment Standard Wall Functions Radiation OFF OPERATING CONDITION Pressure Operating Pressure 1000000 Reference Pressure Location x=0, y=0 Gravity OFF MATERIALS Density ideal-gas Heat Capacity polynomial (2 coefficients) 386.91 0.00713 Thermal Conductivity polynomial (2 coefficients) 2.08E-02 3.94E-05 Viscosity polynomial (2coefficients) 4.00E-05 9.00E-08 Molecular Mass constant 314.0225 BOUNDARY CONDITIONS Walls Thermal conditions Heat Flux all parameters to ZERO Inlet pressure-inlet Gauge TotalPressure udf press_pulse Supersonic/Initial Gauge Pressure udf press_pulse TotalTemperature 3000K Direction Specification Method Normal to Boundary Fluid Fluid Uranium Tetrafluoride PAGE 52 Turbulence Model The standard kmodel presented widely in Chapter 2 has proven in previous studies led at INSPI [4, 5] that this model is most of the time sufficient for the type of problems considered. It is also the turbulence model usually used when dealing with shock waves [12]. In the preliminary study to validate the use of Fluent for our problem, the k-epsilon model was used and it has proved to give good results. Settings Description of the models and tools used to set up the case of propagation and interaction of shock waves in the pipe is given Table 4.2. Properties and models that are kept to default values are not mentioned. The initialization step values are described Table 4.3. Table 4.3: Description of the initialization step INITIALIZE Compute from Inlet Reference frame Absolute Initial values Gauge Pressure 0 x-velocity 0 y-velocity 0 Turbulence kinetic energy 1 Turbulence dissipation rate 1 Temperature 3000.001 Results and discussion Cases Run The different cases run for this study are presented Table 4.4. This analysis deals with the influence of: height of the initial pressure wave initial temperature width of the shock wave length of the pipe, PAGE 53 39 and the differences in properties when the gas is changed from Uranium Tetrafluoride to Helium gas. Table 4.4: Relevant cases run during the shock tube study Case Ambient pressure Material Pressure ratio Width pulse (s) Length pipe Temperature 1 10 bars UF 4 3 5 5mm 3000 K 2 10 bars UF 4 2 5 5mm 3000 K 3 10 bars UF 4 4 5 5mm 3000 K 4 10 bars He 3 5 5mm 2000 K 5 1 bar UF 4 3 5 5mm 3000 K 6 5 bars UF 4 3 5 5mm 3000 K 7 10 bars UF 4 3 15 5mm 3000 K 8 10 bars UF 4 3 1 5mm 3000 K 9 10 bars UF 4 3 continuous 5mm 3000 K 10 10 bars UF 4 4 continuous 5mm 3000 K 11 10 bars UF 4 3 5 5mm 2900 K 12 10 bars UF 4 3 5 5mm 2800 K 13 10 bars UF 4 2.5 5 5mm 3000 K 14 10 bars UF 4 1.5 5 5mm 3000 K 15 10 bars UF 4 3 10 5mm 3000 K The reference case is the one called Case 1, where the gas is Uranium Tetrafluoride with an operating pressure of 10 bars, an ambient temperature of 3,000K and an initial pressure ratio of 3, i.e. the shock wave has an intensity of 30 bars. The reference density in these conditions is 12.5896 kg/m3. Properties are computed every microsecond on a 5mm long portion of the shock tube. The other cases are compared to that reference case and data is plotted every 5s. With the ideal-gas formulation used to model density, the specific heat ratio or gamma calculated by FLUENT remains almost constant at any point of the fluid and at any time as shown Figure 4.9. Values of gamma vary around an average value of 1.068 and the variation is only of a few percent. The error on the computed value compared to the experimental value given by Watanabe and al. [10] is of 1.1%. PAGE 54 40 1.08 Specific Heat Ratio 1.07 1.06 -2.5 -2 -1.5 1.5 2 -1 -0.5 0.5 2.5 0 1 Position (mm) Figure 4.9: Variation of the Specific Heat Ratio on the axis of the shock tube. With the adiabatic model used, the heat loss through the walls is negligible: there is a net heat transfer of 0.3 mW through the walls to be compared to 11 MW through the two ends of the pipe. When the shock wave propagates in the gas, the global shape of the pulse changes. Not only the intensity of the pulse is attenuated, but also the sharp edges are rounded. Yet this study can only be a qualitative study: we are not exactly sure at that point what the characteristics of the pressure pulses will be. It also seems difficult to generate identical pulses with the same width and the same height in two different locations at the same time. We are interested in the profile of properties in the shock tube when the overlapping of the pressure shock wave is maximum, meaning when the pressure reaches a maximum in the shock tube. Properties can be considered constant on the radius PAGE 55 41 therefore only properties on the axis of the shock pipe are computed. Figure 4.10 presents the profile of normalized density, pressure and temperature on the axis of the sock tube. Dividing the values of the properties by 10 bars, 12.5896 kg/m3 and 3,000 K for pressure, density and temperature respectively normalizes the values. Pressure and density reach their maximum at the same time. Fluctuations of temperature are small since the model has been set up for the walls to be adiabatic and then minimize the heat exchanges. No heat generation form fission is modeled then it makes sense that the variations in temperature remain of small proportions. 6 5 density pressure 4 Normalized properties temperature 3 2 1 0 -2.5 -2 -1.5 1.5 2 -1 -0.5 0.5 2.5 0 1 Position (mm) Figure 4.10: Profile of normalized density, pressure and temperature of Uranium Tetrafluoride when the density reaches a maximum. We are interested in the maximum ratio of pressure and density. They are given Table 4.5 for different initial pressure ratios. Since the specific heat ratio of Uranium Tetrafluoride is high, we expect the ratio to be greater than twice the initial ratio from a PAGE 56 42 pressure perspective. As mentioned previously, there is a lot of attenuation. The pressure jump is not as high as expected. Density is proportional to pressure and inversely proportional to temperature. The increase in density is less relevant than in pressure, because there is a small variation of temperature that makes the gas less dense. The process is not as efficient as one would have expected. For an initial pressure ratio of 3, it would have been of great interest to get a pressure ratio of at least 6 when in reality it barely excesses 5. However, we need to consider the dynamics of the process. A neutronic calculation needs to be run to see if the process leads to criticality. Not only are we interested in the absolute value of density of the gas, but we also need to consider how quickly the density raise occurs, how long the gas remains at a high density and how large is the area of high pressure. Table 4.5: Maximum ratios obtained when varying the height of the pressure pulses in Uranium Tetrafluoride. Initial pressure ratio Time of maximum density (sec) Maximum density ratio Maximum pressure ratio 1.5 8.5 1.929 2.019 2 8 2.863 3.085 2.5 7.5 3.752 4.127 3 7 4.550 5.084 When the initial pressure ratio increases, the maximum pressure is higher, but the area where the density is increased becomes narrower as seen Figure 4.11: The area where there is potentially fission and ionization in the gas is more localized. The starting idea was that, for example, if the intensity of the pulse was of 30 bars for an operating pressure of 10 bars, the maximum static pressure reached in the tube would be at least 40 bars. In the computed model the absolute pressure is attenuated to 22 bars before the to pulses overlap. We would expect the pressure to reach at least 34 bars. The maximum pressure is 51.3 bars. Even if the phenomenon does not have the extent expected we PAGE 57 43 would have think that the pressure might be multiplied by 5 or 6 the pressure and the density are increased in an important manner using the interaction of two shock waves. 60 50 ratio of 1.5 ratio of 2 ratio of 2.5 ratio of 3 40 Density (kg/m3) 30 20 10 0 0 0.5 1 1.5 2 2.5 -2.5 -2 -1.5 -1 -0.5 Position (mm) Figure 4.11: Density profile on axis of the shock tube for different initial pressure ratios. Figure 4.12 presents the typical evolution of density with time at the center of the pipe. When the two waves start overlapping, there is a quick increase in density: it is multiplied by 4.5 within 2 s. Density remains at its maximum for a period of 3 s, and then it decreases back to its equilibrium value following an exponential trend. The equilibrium state is reached 15 s after the beginning of the pulses overlapping. If the gas is relaxed, there is even a drop in density after 25 s which is due to a huge drop in pressure. This has no effect on the power generation. PAGE 58 44 60 50 40 Density (kg/m3) 30 20 10 0 0 5 25 10 15 20 30 Time (microsecond) Figure 4.12: Evolution of density with time at the center of the shock tube. From this Computational Fluid Dynamics analysis, a design optimization of the core of the SWD-GCR can be lead. The aim is to find the best operating conditions that would lead to a great increase in the density of Uranium Tetrafluoride, hoping to drive the gas, and the whole core to criticality. The parameters we analyze are the influence of the operating pressure, the inlet temperature of the gas, the pressure ratio in the shock wave and the period length during which the pulse is created. The influence of the initial pressure ratio has already been discussed previously and can be seen Figure 4.11 and Table 4.4. As the initial ratio increases, the difference between density before and during overlapping increases too. A higher pressure ratio leads to a greater and faster peak in density. The area where there is a discontinuity in density is narrowed, which is a plus for the generation of electricity since energy would be released in a very located space and then the losses would be minimized. PAGE 59 45 The shape of the pulses is then analyzed. The width of the pressure pulses is increased or decreased and its influence on the maximum density is the plot on Figure 4.13. When the pulse is too small, a few microseconds, there is a lot of attenuation of pressure as previously discussed. The resulting density remains small. When the width of the pulse reaches 10 s, we get a density maximum that cannot be increased by increasing the width of the pulse. The initial design plans on a pulse of 5 s. From Figure 4.11 it seems to be a good choice since density is high in a big portion of the pipe, even if it does not reach the maximum possible. 010203040506070-2.5-2-1.5-1-0.500.511.522.5Position (mm)Density (kg/m3) width of 1 microsec width of 5 microsec width of 10 microsec width of 15 microsec continuous pulse Figure 4.13: Density profile on axis of the shock tube for initial pressure pulses of different widths. The last analysis is done by varying the operating pressure of the core. The normalized density on the axis when overlapping of the pressure waves is at its maximum is plotted on Figure 4.14. The profiles are equivalent for operating pressures ranging from 1 to 10 bars. PAGE 60 46 5 4.5 4 10 bars 5 bars 3.5 1 bar Normalized density 3 2.5 2 1.5 1 0.5 0 -2.5 -2 -1.5 1.5 2 -1 -0.5 0.5 2.5 0 1 Position (mm) Figure 4.14: Density normalized to the operating density on axis of the tube for different operating pressures. Description of the Results with Helium Gas The shock wave in Helium gas propagates faster than Uranium Tetrafluoride gas. This can be expected since the molecule is lighter. The case cannot be processed at 3,000K. The solver diverges, due to a too high viscous turbulence. Problems with Helium can be run in longer tubes, since the shock wave propagates faster, here the pipe is 5 cm long. The inlet temperature of the gas is set to 2,000 K and the first noticeable particularity of the shock wave interaction in Helium is the high gradient of temperature obtained during the process. On Figure 4.15.c the temperature reaches 4,000 K when the pressure waves are at their maximum overlapping. The inlet temperature is maintained at 2,000 K, and the walls are adiabatic. Curiously, density does not peak at the center of the tube where the gradient of pressure is the biggest (Figure 4.15), and also, the highest PAGE 61 47 values of density are reached a moment before total overlapping of the shock waves (Figure 4.16). Since the specific heat ratio of Helium is higher than the gamma of UF 4 1.667 and 1.08 respectively the ratios of maximum density over operating density and maximum pressure over operating pressure are smaller: The density ratio is equal to 3.8 and the pressure ratio to 5.4 in Helium when these values for a comparable case were equal to 4.5 and 5.1 in the case of Uranium Tetrafluoride. The pressure jump in Helium is of greater importance but because of the variations in temperature the density ratio is not as important as for Uranium Tetrafluoride. a)0.E+002.E+064.E+066.E+06-25-15-551525Position (mm)Pressure (Pa) b)00.250.50.751-25-15-551525Position (mm)Density (kg/m3) c)010002000300040005000-25-15-551525Position (mm)Temperature (K) Figure 4.15: At time t = 15 micros, maximum overlapping, profile of properties on the axis of the shock tube a) Absolute pressure b) Density c) Temperature. PAGE 62 48 a)0.E+002.E+064.E+066.E+06-25-15-551525Position (mm)Pressure (Pa) b)00.250.50.751-25-15-551525Position (mm)Density (kg/m3) Figure 4.16: At time t = 12.5 micros, maximum overlapping, profile of properties on the axis of the shock tube a) Absolute pressure b) Density. Conclusion on the Shock Wave Interaction Analysis The previous analysis was performed on a purely thermodynamic aspect. It highlighted the high attenuation of the traveling shock waves. When the waves approach density increases and reaches a maximum at overlapping of the pressure pulses. However the highest densities obtained do not seem sufficient to ignite prompt fission supercriticality in the fissile gas. A neutronic study needs to be led to determine if the jump of properties at the center of the core is of magnitude enough to generate high levels of power. PAGE 63 CHAPTER 5 COUPLING OF CFD CODE WITH MCNP Modeling the shock propagation and fission power production in the Shock Wave Driven Gas Core Reactor is a complicated task for many reasons. There is no general code for coupled gasdynamic equations and neutron transport equations. The time for developing this type of software would be excessive. The solution chosen is to run separately the gasdynamic portion of the core and then create an input for a neutronic code. The gas behavior is studied with the CFD code FLUENT as extensively discussed in Chapter 4. The data extracted from the code are structured into input data files that are implemented into the general purpose Monte Carlo neutron transport code MCNP4C 1 The scope of this study is to calculate fission power neutron flux in the shock tube and the surrounding reflector. Input The gasdynamic part of the analysis is run using FLUENT in the following conditions: Real modeling of waves generation, i.e. one shock created on both sides of the tube Tube size: 1 meter long and 80 cm radius Properties of UF 4 specified in Chapter 3 Axial grid spacing: 0.167 cm Operating Temperature: 1,567 K 1 Developed by Los Alamos National Lab (LANL) 49 PAGE 64 50 Operating pressure: 10 bars Pressure Ratio of the wave: 3 Pulse time: 1.1 ms Due to the size of the grid 95,840 cells and the distance along which the shock has to propagate the computational time for this case is very long. It took 54 hours to compute 2 ms of real evolution. Time steps where in the range of 10 -7 to 5*10 -6 sec. To initiate the process, since there is a huge pressure gradient on a few cells, time steps need to be small. When the shock only propagates, bigger time steps can be used. As the two waves overlap, smaller time steps of the order of magnitude of 10 -7 seconds are used: the properties evolution obtained are refined. Density profiles were saved in a text file for many time steps of interest: there are eight relevant time steps at 2 and 70 s and 0.2, 0.55, 1.184, 1.659, 1.765 and 2.125 ms. The type of study presented here has already been led at Innovative Nuclear Space Propulsion and Power Institute [1], using an in-house MHD code coupled with MCNP. The in-house MHD code, called MHD1DSHOCK, simulates more specifically the physics of electromagnetic shock tube. An electromagnetic source term is added to the Navier-Stokes equations but a lot of simplifying assumptions are considered such as negligible thermal conductivity and viscosity, and governing single fluid MHD equations reduced to cylindrically symmetric form. FLUENT gives significantly different results than the code considered at INSPI, especially concerning the attenuation of the pressure waves propagating in Uranium Tetrafluoride. The aim of this Chapter is to present the results obtained when running MCNP with the density evolution computed by FLUENT. Files are created to input MCNP. Each file contains time step, strength of the source of neutrons, spatial grid and values of density associated to each spatial node. In the neutronic study, it is assumed that PAGE 65 51 the tube can be considered as one-dimensional. Considering a one-dimensional model of the system suppresses the physical instabilities of the 2-D or 3-D mock-up. According to the results from Chapter 4, the radial variation of physical properties can be neglected. Then 1-D mock-up is adapted to the problem. Properties in the output of FLUENT are expressed in SI units whereas the input for MCNP should be expressed in the CGS system. An example of inputs is given Figure 5.1. In the situation presented here, the mesh is discretized into 11 spatial steps. title mhdsout.dat time 2E-06 seconds delt 0.150273E-07 seconds noft 0.00000 source 0.100000E+06 leftbnd 0.000000E+00 begin z 11 0.00000 10.01670 20.03340 30.05010 40.06678 50.08347 60.10020 70.11690 80.13360 90.15020 100.00000 end z begin rho 11 2.73634E-02 2.41026E-02 2.41026E-02 2.41026E-02 2.41026E-02 2.41026E-02 2.41026E-02 2.41026E-02 2.41026E-02 2.41026E-02 2.73634E-02 end rho Figure 5.1: Example input for MCNP based on FLUENTs results. Processing of FLUENTs Output The case was run for a 1 m diameter, 1 m long, 90% enriched core; 50 cm reflector of BeO. The case considered in FLUENT was initially an 80 cm radius core. As demonstrated in Chapter 4 the radius of the core does not influence the values of the different properties; the differences of radius will not affect the results compiled by MCNP. FLUENT gives the value of density at every node of the mesh whereas MCNP PAGE 66 52 considers volumes of materials with different densities. The output created showed Figure 5.1 is collapsed by an in-house code, MHD2MCNP to adapt FLUENTs results to MCNPs input. The process of accommodation is illustrated Figure 5.2. Discretized description of density is transformed into a continuum of areas with constant averaged densities. Figure 5.2: Shock tube models used in the combined simulation codes. (a) Illustration of the shocktube model output by FLUENT with many cells. b) Illustration of the shocktube model after cell collapse by MHD2MCNP where many of the cells with nearly the same density/pressure are collapsed to simplify the MCNP model Ideally, an iterative loop should be created between FLUENT and MCNP as illustrated Figure 5.3. This logical loop should be maintained until the FLUENT code completes the interaction of the two shock waves or reaches its maximum assigned integer time step. However, this loop had been created to use the in-house program. At this stage of the study, FLUENT has not been integrated to the loop described Figure 5.3. PAGE 67 53 The CFD code is run for the total length of the phenomenon of interest and then density profiles are input in MCNP at given time steps. FLUENT mhdsout.dat Figure 5.3: Flow diagram of operating script linking FLUENT and MCNP. In the input for MCNP, the number of neutrons at a given history is the sum of the neutron source and the neutrons generated in the fission process. The neutron population is updated for each time step considering the previous k eff and prompt neutron life time according to the following formula: )1(1)()(/)1(/)1(1nnnnltknnltknneklSetNtN The initial condition is 0)0( tN and t ttnn 1 S is the background neutron source in neutrons/sec. In the case run we considered 10 5 neutrons/sec, which is a strong source. Write collapsed file cell density and right boundary location data mhdsin.dat MHD2MCNP file MCNP2MHD code outpu t file MCNP code input file code Write flux and fission heating data Yes MCNP MCNP converged? Continue Run code No PAGE 68 54 Results from the Neutronic Calculations A file is created from MCNPs output for each time step and is called mhdsin.dat. An example of output for the first time step is given Figure 5.4. The case presented corresponds to the very beginning of the shock waves propagation. Properties of Uranium Tetrafluoride are almost constant: the gas in the pipe can be described with only three areas of different densities. From MCNP analysis it is found that the core is already critical with a k eff of 1.211 at the very beginning of shock waves generation. This can be due to the high enrichment of 90% of the fissile gas assumed in this case, and also to the strength of the background source of neutrons. MCNP also outputs power density in MeV/cc and the neutron flux in each sub-volume of the core in neutrons/(cm 2 .s). The evolutions of density and power density generated in the shock tube are given in Appendix B. Figure 5.5 gives the profiles of density and power density at t=1.659 ms, time of maximum overlapping of the pressure waves. Table 5.1 compiles the maximum power generation and Figure 5.6 and 5.7 show the evolution of neutron flux and k eff with time respectively. One would expect that power density peaks when there is a maximum density. When at maximum overlapping, density in the middle of the shock tube is multiplied by 4; the corresponding power generated is only multiplied by 2. The power generation does not seem to follow the density evolution, which is of big concern. Some energy is created but not localized at the center of the core. When pressure waves are getting closer, the k effective of the core increases, but the global neutron flux decreases. The maximum k is reached at maximum overlapping and then it seems to decrease: this trend can be due to statistical error. The neutron population is greater in the shock wave where density is high. The total number of neutrons tends to decrease with time because when pressure waves travel they get further from the inlet and the BeO moderator. When PAGE 69 55 the pressure pulses overlap there is a lot of leakage and the population of neutron decreases. The prompt neutron life time calculated by MCNP is of the order of a millisecond. That is the time required for the shock waves to travel from inlet and overlap. This limits the neutron multiplication and hence the amount of power that is produced. In the present analysis, no significant power is created. To increase the power to hundreds of kilojoules per pulse the neutron life time should be decreased by at least an order of magnitude. This could be obtained by changing the configuration of the shock tube core, especially by changing the reflector from BeO to a reflector with lower moderating ratio. title mhdsout.dat time 2.00000000e-06 seconds delt 1.50273000e-08 seconds noft-1 0.00000000e+00 source 1.00000000e+05 begin cells 3 100 101 102 end cells 3 begin rho 3 2.61136000e-02 2.41073504e-02 2.73634000e-02 end rho 3 leftbnd 0.00000000e+00 begin z 3 1.66900000e-01 9.98331000e+01 1.00000000e+02 end z 3 lifetime 1.89782000e-03 keff 1.21146000e+00 sdkeff 5.34000000e-03 begin flux 3 2.99898000e-04 3.01088000e-04 3.02250000e-04 end begin q 3 1.79500447e-04 1.15830033e-04 1.85101361e-04 end noft 1.50273126e-03 Figure 5.4: Example output from MCNP. PAGE 70 56 0.0E+002.0E-024.0E-026.0E-028.0E-021.0E-011.2E-01050100Position (cm)Density (g/cc) 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04050Position (cm)Power Density (MeV/cm3) 100 Figure 5.5: Profile of density as it is input in MCNP and profile of power density as calculated by MCNP at t=1.659 ms. Table 5.1: Maximum power generation and its location at every time step. time (ms) Maximum Power Density (MeV/cm 3 ) Location (cm) 0.002 1.851E-04 0 100 0.07 3.219E-04 0 100 0.2 2.591E-04 0 100 0.55 2.076E-04 0 100 1.1184 1.898E-04 3. 97. 1.659 2.068E-04 50 1.765 1.916E-04 50 2.125 1.818E-04 50 3.1E-04 Neutron Flux (#/cm2.s) 3.0E-04 2.9E-04 2.8E-04 2.7E-04 0 0.5 1 1.5 2 2.5 Time (ms) Figure 5.6: Evolution of total neutron flux in the tube with time. PAGE 71 57 1.21.221.241.2600.511.522Time (ms)k effective .5 Figure 5.7: Evolution of k eff with time. Conclusion The analysis presented in this chapter can be considered as a preliminary coupling between FLUENT, a CFD code and MCNP, a neutronic code. The first code, FLUENT that compiles the thermal properties of the fissile gas has been run independently from MCNP, then there is no feed-back. However, this analysis shows that the power generation in the case of a pulsed gas core is not as important as expected. A possible explanation is that the time involved in shockwave propagation and interaction is of the order of the prompt neutron lifetime. Coupling between FLUENT and MCNP should give more accurate results, but this first analysis shows that some modifications should be done to the design of the Shock Wave Driven Gas Core Reactor. PAGE 72 CHAPTER 6 DESIGN ANALYSIS OF EXPANSION NOZZLE To date, the main limitation in increasing the profitability of a nuclear plant remains in the poor thermal efficiency of the design; nuclear plants thermal efficiency reaches at best 35%. In the scope of the next generation of nuclear plants, one of the objectives is to increase the efficiency of the thermal cycle, by using a magnetohydrodynamic generator (MHD) to convert thermal energy into electrical energy. Primary cycle of Gas Core Reactors will be connected to the thermodynamic cycle through an expansion nozzle that changes the direction of the flow and accelerates the working gas. The nozzle should be designed and boundary conditions chosen such as the flow remains subsonic through the nozzle, i.e. Mach number of the flow is smaller than one. If shock waves were created they would then propagate inside the nozzle and go back into the core of the reactor leading to instabilities. General Description of the Nozzle Figure 6.1 shows a sketch of the expansion nozzle. Heated gas exits the core in a large two-meter-diameter tank. The gas has an inlet velocity of approximately 50 m/s. Velocity profile is uniform. The inlet pressure is equal to 50 bars and the temperature in the tank is constant and equal to 3,000 K. Gas first goes through a converging-diverging nozzle where the throat size is 20 cm. The flow is then divided in two identical parts, consisting in two slowly converging nozzles. The exit of the pipe is 10 cm diameter and the converging part is 3 meters long. The exit back pressure is unknown; with this analysis, we should be able to find out the 58 PAGE 73 59 pressure range for which the velocity remains subsonic in the throat of the converging-diverging part of the expansion nozzle. Typically the back pressure should be of a few bars. Figure 6.1: Design and boundary conditions for the expansion nozzle. Meshing of the Nozzle The 2-D representation of the expansion nozzle is symmetric about the y-axis. To save computational time, only one half of the nozzle could be meshed and the mirror plane set as an axis. It can be noticed that the flow at inlet is in the direction of the negative y-axis. FLUENT, when dealing with one-dimensional flows usually considers a flow along the x-axis. Hence when only half of the nozzle is modeled, flow with consequent velocity comes from the axis, which is physically incorrect. The whole nozzle has to be modeled, which doubles the computational time. A good way to check that the convergence criteria chosen are adequate is to make sure that the different properties fields are symmetrical with respect to the nozzle axis. PAGE 74 60 Figure 6.2 shows the grid created with Gambit [14-15]. It is composed of 2,120 quadrilateral mesh elements. The mesh is refined in the area of the converging-diverging nozzle, as well as in the upper part of the converging pipes as can be seen Figure 6.3. Figure 6.2: Meshing of the expansion nozzle Figure 6.3: Refinement of the size of the grid in the throat of the converging-diverging nozzle. Since the initial velocity is directed vertically and down it is more appropriate to use meshes that are more or less parallel to that initial velocity and then quadrilateral PAGE 75 61 cells seem to be the best choice. The same case run using triangular cells in the rectilinear part of the nozzle the upper tank leads the solver to diverge. Specificity of the Nozzle Case Inlet Boundary Condition and Initialization Step This analysis present one difficulty: With FLUENT one can set the inlet boundary condition for either pressure or velocity, but not both of them at the same time. The inlet boundary condition sets the pressure (Define/Boundary Condition/pressure-inlet) at inlet of the tank. Setting the inlet velocity to 50m/s during the initialization step is not enough: the initialization step is only considered as a user-guess and FLUENT only takes it into account for the first iterations. After a few iterations the inlet velocity diverges from the value of 50 m/s and its direction becomes random. One method for initializing the solution is to patch values or functions for selected flow variables in selected cell zones. Once the entire flow field is initialized, different values for particular variables into different cells can be patched using the patching function in the Solve/Initialize menu. The solution to keep the magnitude and direction of velocity at inlet of the tank is to patch the fluid for the x-velocity and the y-velocity to 0 and m/s respectively. Using this method, the two inlet conditions on pressure and velocity are respected. Cases Considered and their Set-up The case is run for different back pressures, ranging from 1 to 49 bars. Contours of different properties are plotted. The function of interest is the Mach number. Solver FLUENT users manual [13] advises to use a coupled implicit solver for compressible transonic flows. A compressible gas is considered, Uranium Tetrafluoride, defined by the ideal-gas law. Residuals PAGE 76 62 At first the cases were run with default residuals, meaning that the convergence criterion is 10 -3 Once the solution has converged, the properties are not symmetrical with respect to the vertical axis. Convergence criteria should be reduced. Since it is a time-independent problem with a relatively small number of computational cells, the residuals are set low, 10 -6 for results to be very accurate. This case requires about 80,000 iterations, the number of iterations depending on the value of the back pressure. Turbulence Formulation Reynolds number at inlet is of the order of 2.*10 -7 A turbulent model is necessary. The Spalart-Allmaras turbulence model was designed specifically for aerospace applications involving wall bounded flows [13]. As discussed in Chapter 2, to treat the nozzle problem, the Spalart-Allmaras formulation is used. Operating and Boundary Conditions Operating pressure is set to zero and pressures specified at boundaries are absolute pressures. Inlet: Boundary conditions at inlet have been discussed previously. Exit: Pressure_Outlet are chosen and set to the back pressure considered. Backflow total temperature is set to 3,000 K. Walls: Adiabatic walls are considered as explained in Chapter 3. Results and Discussion With a convergence criterion of 10 -6 the profiles obtained for different properties are symmetrical with respect to the vertical axis. Profiles of absolute pressure, Mach number and static temperature are given Figures 6.7, 6.8 and 6.9 at the end of the chapter. For a better visualization, only half of the nozzle is shown. PAGE 77 63 Pressure Distribution Profiles of pressure are drawn to make sure that FLUENT is able to set the right pressure at exit. The same boundary condition is applied to two different locations, right and left exits, which could involve numerical problems. Pressure remains high in the tank but drops once it goes through the diverging-converging part of the nozzle. The pressure profile when back pressure is set to 10 bars is given Figure 6.7. It seems that the pressure is not set properly on the exit. A check is done and given Figure 6.4. Even if the pressure imposed at exit is not respected in every cell, the exit pressure tends to be 10 bars. Figure 6.4: Absolute pressure on the exit plane of the nozzle when the boundary condition sets it to 10 bars. Mach Number Profile and Repartition As expected there is an increase in the Mach number when the flow goes through the throat of the nozzle. However for reasonable back pressures greater than a tenth of a PAGE 78 64 bar the flow never reaches sonic Mach number in the converging-diverging part of the expansion nozzle. The profile of Mach number along the vertical axis of the nozzle can be seen Figure 6.5. The back pressure if smaller than 20 bars does not influence the profile of Mach number in the vertical part of the expansion nozzle. When exit pressure is increased, the Mach number at the center of the throat decreases. For the back pressures considered, flow never chokes in the vertical nozzle, which is what one looks for. The highest Mach number reached in the throat is 0.7. 00.10.20.30.40.50.60.7-0.4-0.200.20.40.60.811.21.41.6Vertical position (m)Mach Number 1 bar 10 bars 20 bars 30 bars 40 bars 49 bars Figure 6.5: Profile of Mach Numbers along vertical axis of nozzle for different exit pressures. As can be seen Figure 6.8, flow accelerates in the converging pipes. Maximum Mach number is reached at the exits of the nozzles. At the outlet of the converging pipes Mach number is close to 1 for back pressures smaller than 20 bars. This is an interesting feature: assuming that a diverging tube is added to the two-stage nozzle, velocity of the gas would reach high values [16]. PAGE 79 65 Table 6.1 presents the maximum Mach number reached by the gas in the expansion nozzle. Table 6.1: Maximum Mach number reached in the nozzle for different back pressures. Back Pressure (bars) Maximum Mach Number 0.1 1.006249 1 1.007255 2 1.007105 5 1.005936 10 1.001902 20 1.00102 21 0.989691 22 0.9543363 23 0.9201925 25 0.8534781 30 0.6993469 40 0.4182654 49 0.1225149 Non-Adiabatic Case All the problems were treated with adiabatic walls. In reality, the walls of the two-stage nozzle are cooled. An analysis has been led changing the condition on adiabatic walls into a constant temperature walls. In a second approximation walls are considered to be at a uniform temperature of 1,800 K. The temperature profile in the nozzle is given Figure 6.10 for a back pressure of 10 bars. Pressure profiles on the axis of the vertical converging-diverging nozzle are compared for adiabatic and constant temperature walls. The profiles are given Figure 6.6. The maximum Mach number reached when the walls are cooled is increased by one percent. The fact that walls were modeled adiabatic does not change the global trend of the Mach number profiles in the nozzle. PAGE 80 66 00.10.20.30.40.50.60.7-0.4-0.200.20.40.60.811.21.41.6Vertical Position (m)Mach Number adiabatic walls constant temperature 1,800 K Figure 6.6: Difference in Mach number profiles on the axis of the C-D nozzle for adiabatic walls and cooled walls. Conclusion The challenge with this analysis was to model properly the moving gas, and more specifically the inlet condition. We proved that the case is solvable with FLUENT 6.1. The scope of this part of the study was to get a first idea of whether or not the gas, Uranium Tetrafluoride, would choke in the nozzle as it has been designed. The study proved that for any exit pressure greater than 0.1 bars the flow remains subsonic in the area of the throat and that one should not expect creation and propagation of shock waves in the expansion nozzle. For back pressures smaller than 20 bars, the exit velocity of the gas is close to sonic velocity. PAGE 81 67 Figure 6.7: Contours of absolute pressure in the nozzle (only the right part is represented) when the exit pressure is set to 10 bars. PAGE 82 68 Figure 6.8: Contours of Mach number when the exit pressure is set to 10 bars PAGE 83 69 Figure 6.9: Contours of temperature when the exit pressure is set to 10 bars PAGE 84 70 Figure 6.10: In the case of cooled walls (T w =1,800 K), contours of temperature when the exit pressure is set to 10 bars PAGE 85 CHAPTER 7 SUMMARY AND CONCLUSIONS Summary of Results The main goal of this analysis is to develop a Computational Fluid Dynamics (CFD) model for both steady and pulsed Gas Core Reactors. The singularity of this type of reactors is the extreme conditions of operation: high temperatures (3,000 K), high pressure (50 bars), high magnetic field and strong gradients of pressure. The study proves that the CFD code FLUENT is an adequate tool to model GCRs. The Computational Fluid Dynamics Gas Core Reactor model is used in two different analyses. The conclusions for each study are presented separately. Pulsed Gas Core Reactor: Study of propagation of pressure pulses in the shock tube highlighted the strong attenuation of the shock waves intensity in uranium tetrafluoride. This is explained by an important viscous dissipation that is properly described by the standard k-epsilon turbulence model. The CFD analysis also shows that interaction of two approaching pressure pulses does not lead to the abrupt increase in the properties of the gas that one would have expected. The analysis shows that the best results would be obtained by operating the SWD-GCR at 10 bars and that pressure shock waves should be generated with a height of 30 bars during 10 s. In these conditions, density of the fissile gas is multiplied by 4.55 when pressure pulses overlap. 71 PAGE 86 72 Processing of the thermodynamics data obtained with FLUENT by MCNP shows that the power generated by the collision of two single pulses generates energy of the order of a few Joules. A new design for the Shock Wave Driven Gas Core Reactor should be considered since the actual design does not allow the release of enough energy when the shock waves collide. The main effort should be focused on increasing dramatically the density of the gas fuel at the center of the shock core reactor. Steady Gas Core Reactor: Analysis of the nozzle that converts enthalpy into very high kinetic energy for MHD gives good results. The design of the nozzle is validated by the study. In the converging-diverging part of the nozzle, flow remains subsonic for any back pressure. Uranium tetrafluoride gas does not choke in the throat of the two-stage nozzle. If it did, shock waves would have been able to propagate back into the core and serious safety issues would have appeared. For back pressure smaller than 20 bars, flow is close to sonicity at the exit of the converging portion. Hypothetically, a diverging system can be added after the converging portion so uranium tetrafluoride reaches high velocity at the exit of the two-stage nozzle. This high kinetic energy is needed for the magnetohydrodynamic power generator. The CFD analysis validates the global shape of the nozzle. Some additional work could be done on a mechanical engineering point of view to optimize the design of the different portions of the nozzle. Future Work Additional work should be done to this analysis that can be considered as a validation of the CFD-GCR model and a qualitative study of steady and pulsed GCRs. Some problems have not been solved by the analysis presented in this thesis. For PAGE 87 73 example, the maximum pressure ratios that have been treated by FLUENT do not exceed 4. No experiments have been led to date to figure out the typical intensity of the shock waves. If it turns out that the amplitude of the shock waves created by high magnetic fields allows pressure ratios greater than 4, a solution should be find to model it in the CFD code. Refinement of the model could be obtained by adding a thermal feed-back. In the present analysis, adiabatic walls are considered as well as uniform inlet temperature profiles. This is a rough approximation, and future studies could model more accurate temperature distributions. In the case of the SWD-GCR, overlapping of the pressure waves leads to fission. This process creates a great amount of energy and a temperature raise in the gas that has not been taken into account. A real coupling between the CFD code and MCNP should be developed to take into account these temperature changes and get more realistic results. PAGE 88 APPENDIX A EVOLUTION OF DENSITY IN THE SHOCK TUBE In this appendix, the profile of density in the shock tube is plotted for different time steps. The reference case is considered: core running at an operating pressure of 10 bars and an ambient temperature of 3,000 K with pressure pulses of 30 bars and 5s long in Uranium Tetrafluoride. The evolution can be divided into different steps: Generation and propagation of shock waves: attenuation of the intensity of the pulse Overlapping of the pressure waves, rapid increase in density to reach a peak Relaxation of the gas time = 1s 0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 4s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 5s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 6s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) 74 PAGE 89 75 time = 7s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 7.2s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 7.5s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 8s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 9s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 10s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 11s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 13s 0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) PAGE 90 76 time = 15s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 20s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 25s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) time = 30s0204060-2.5-1.5-0.50.51.52.5Position (mm)Density (kg/m3) PAGE 91 APPENDIX B POWER DENSITY EVOLUTION IN THE SHOCK TUBE This Appendix presents the results of the MCNP calculation run in the case of Chapter 5, a shock tube core of 1 m long, 1 m diameter, 50 % enrichment and surrounded by a 50 cm thick BeO moderator. For each time step, density of Uranium Tetrafluoride and power density generated by fissions are plotted. Time = 2 s 0.0E+002.0E-024.0E-026.0E-028.0E-021.0E-011.2E-01050100Position (cm)Density (g/cc) 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04050Position (cm)Power Density (MeV/cm3) 100 Time = 70 s 0.0E+002.0E-024.0E-026.0E-028.0E-021.0E-011.2E-01050100Position (cm)Density (g/cc) 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04050Position (cm)Power Density (MeV/cm3) 100 77 PAGE 92 78 Time = 0.2 ms 0.0E+002.0E-024.0E-026.0E-028.0E-021.0E-011.2E-01050100Position (cm)Density (g/cc) 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04050Position (cm) Power Density (MeV/cm3) 100 Time = 0.55 ms 0.0E+002.0E-024.0E-026.0E-028.0E-021.0E-011.2E-01050100Position (cm)Density (g/cc) 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04050Position (cm)Power Density (MeV/cm3) 100 Time = 1.184 ms 0.0E+002.0E-024.0E-026.0E-028.0E-021.0E-011.2E-01050100Position (mm)Density (g/cc) 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04050Position (mm)Power Density (g/cc) 100 PAGE 93 79 Time = 1.659 ms 0.0E+002.0E-024.0E-026.0E-028.0E-021.0E-011.2E-01050100Position (cm)Density (g/cc) 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04050Position (cm)Power Density (MeV/cm3) 100 Time = 1.765 ms 0.0E+002.0E-024.0E-026.0E-028.0E-021.0E-011.2E-01050100Position (cm)Density (g/cc) 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04050Position (cm)Power Density (MeV/cm3) 100 Time = 2.125 ms 0.0E+002.0E-024.0E-026.0E-028.0E-021.0E-011.2E-01050100Position (cm)Density (g/cc) 0.0E+005.0E-051.0E-041.5E-042.0E-042.5E-043.0E-043.5E-04050Position (cm)Power Density (MeV/cm3) 100 PAGE 94 LIST OF REFERENCES 1. Smith B, Knight T, Anghaie S: Report on Shock Wave Driven Gas Core Reactor Simulations Final Report for the PMI-GCR Concept Map Project, Innovative Nuclear Space Power and Propulsion, University of Florida, 2003. 2. Spalart P, Allmaras S: A One-Equation Turbulence Model for Aerodynamic Flows, Technical Report AIAA-92-0439, American Institute of Aeronautics and Astronautics, New York, 1992. 3. Paciorri R, Dieudonne W, Degrez G, Charbonnier J-M, Deconinck H: Validation of the Spalart-Allmaras Turbulence Model for Application in Hypersonic Flows, American Institute of Aeronautics and Astronautics, 1997. 4. Plancher J: Thermal and Fluid Design Analysis of a Square Lattice Honeycomb Nuclear Rocket Engine, Masters Thesis, University of Florida, 2002. 5. Roux V: Evaluation of RELAP5 Reactor Core Modeling Capability, Masters thesis, University of Florida, 2001. 6. Jones W and Launder B: The prediction of laminarization with a two-equation model of turbulence, Journal of Heat and Mass Transfer 15, 1972. 7. Mohammadi B, Pironneau O: Analysis of the K-Epsilon Turbulence Model (p.92), John Wiley & Sons, Chichester, NY, 1994. 8. McBride B, Gordon S, Reno M: Coefficients for Calculating Thermodynamic and Transport Properties of Individual Species, NASA TM-4513, 1993. 9. Gurvich L, Yungman V.S, Dorofeeva O, Gorokhov L, Munvez S: Thermodynamic Properties of the Uranium-Fluorine Gaseous System, Preprint 1-0018. 10. Watanabe Y, Anghaie S: Thermophysical Properties of Gas Phase Uranium Tetrafluoride, American Institute of Aeronautics and Astronautics, 28 th Thermophysics Conference, Orlando, Fl, 1993. 11. McBride B, Gordon S: Computer Program for Calculation of Complex Equilibrium Compositions and Applications I, Sec. 2.3-2.4, NASA Reference Publication 1311, 1994. 80 PAGE 95 81 12. Kim S-W: Numerical Computation of Shock Wave Turbulent Boundary Layer Interaction in Transonic Flow Over an Axisymmetric Curved Hill, NASA Technical Memorandum 101473, ICOMP-89-3, NASA, February 1989. 13. FLUENT 6.1 Users Guide, Fluent Inc., Lebanon, NH, 2003. 14. Oates G: Aerothermodynamics of Gas Turbine and Rocket Propulsion (p.71-74), American Institute of Aeronautics and Astronautics, Washington, DC, 1988. 15. Emanuel George: Gasdynamics, Theory and Applications (p.89-101), American Institute of Aeronautics and Astronautics, Washington, DC, 1986. 16. Burmeister L: Convective Heat Transfer (p.355), John Wiley & sons, Inc, Chichester, NY, 1993. PAGE 96 BIOGRAPHICAL SKETCH Anne Charmeau was born May 13, 1980, in Versailles, France. She graduated from Edouart Herriot High School in Voiron, France, in 1998. After two years of pre-engineering classes, she enrolled at the National School for Physics in Grenoble in Fall 2000 and obtained her engineer diploma in September 2003. As part of an academic exchange she joined the University of Florida in August 2002. Since then, she has been pursuing a Master of Science in the Nuclear and Radiological Engineering Department and has been working at Innovative Nuclear Space Power and Propulsion Institute as a research and teaching assistant. 82 |