Roman Thiele ModelingofDirectContactCondensationWithOpenFOAM

(1)

AlbaNova Condensation Models Master Thesis

Modeling of Direct Contact Condensation With OpenFOAM

Roman Thiele

Supervisors:

Henryk Anglart

Eric Lillberg

Division of Nuclear Reactor Technology Royal Institute of Technology Stockholm, Sweden, June 2010

(2)

(3)

ABSTRACT

Within the course of the master thesis project, two thermal phase change models for direct contact conden-sation were developed with different modeling approaches, namely inter-facial heat transfer and combustion analysis approach.

After understanding the OpenFOAM framework for two phase flow solvers with phase change capabilities, a new solver, including the two developed models for phase change, was implemented under the name of interPhaseChangeCondenseTempFoam and analyzed in a series of 18 tests in order to determine the physical behavior and robustness of the developed models. The solvers use a volume-of-fluid (VOF) approach with mixed fluid properties.

It has been shown that the approach with inter-facial heat transfer shows physical behavior, a strong time step robustness and good grid convergence properties. The solver can be used as a basis for more advanced solvers within the phase change class.

(4)

(5)

ACKNOWLEDGEMENTS

I would like to extend my gratitude all my supervisors and listeners during the course of this work. Most notably to Henryk Anglart for this guidance and the support as a supervisor and Carsten ’t Mannetje for listening and letting me bounce of ideas. Also the rest of the reactor technology department for giving valuable input during the bi-weekly meetings.

As well I would like to thank Westinghouse Electric Sweden AB and most notably Eric Lillberg for financial support and guidance in understanding the OpenFOAM source code.

(6)

(7)

LIST OF SYMBOLS

The list of symbols include bold face symbols, which indicate a vector property.

Symbol Dimensions Description

Latin Symbols

a variable Discretization coefficients

A m2 _{Inter-facial area}

cp J/kg/K = m2/s2/K Specific heat

e variable Relative error

Evl J/m3_/s _{Heat transfer between phases}

fc Hz = s−1 Combustion frequency

g m/s2 _Gravity

GCI - Grid convergence index

h W/m2_{/K = kg/s}3_/K _{Heat transfer coefficient}

i J/kg Specific enthalpy

if g J/kg = m2/s2 Latent heat

l m Representative grid size

L m Characteristic length

m kg Mass

˙

m kg/s Mass flow

N - Total number of cells

N u - Nusselt number

p P a Pressure

p - Apparent order of the solver for GCI

q W = kg · m2_/s3 _Heat

q00 W/m2_{= kg/s}3 _{Heat flux}

Q J/m3_/s _{General heat source}

r - Ratio of representative grid sizes

S variable Generic source term

T K or◦C Temperature

u m/s Velocity

Vcell m3 Volume of a cell

x m Length

Greek Symbols

α - Phase fraction of water

Γ kg/m3_/s _{Mass transfer per unit volume}

δb.l. m Boundary layer thickness

ε variable Absolute error of global variable

κ m2_/s _{Thermal diffusivity}

λ W/m/K = kg · m/s3_/K _{Thermal conductivity} Continued on next page

(10)

x CONTENTS

Continued from previous page

Symbol Dimensions Description

(11)

LIST OF FIGURES

2.1 Volume-of-Fluid method . . . 4

2.2 Energy transfer schematic . . . 6

2.3 Cell layout with interface boundary . . . 7

2.4 Graphical representation of the interface area model. . . 8

3.1 Runtime command structure interPhaseChangeCondenseTempFoam . . . 9

3.2 Class diagram . . . 15

4.1 Sketch of the simulation model . . . 20

5.1 GCI comparison . . . 24

5.2 Velocity dependency time development . . . 25

5.3 Interface distribution velocity test of the interface model . . . 25

5.4 Interface distribution velocity test of the combustion model . . . 26

5.5 Time step dependency over time development . . . 26

D.1 Condensation rate vs. time for inlet velocity 10m/s . . . 43

D.2 Time step behavior high velocity test 0.15s . . . 44

(12)

(13)

LIST OF TABLES

3.1 AlbaNovaInterfaceCoeffs . . . 17 3.2 AlbaNovaCombustionCoeffs . . . 17 4.1 Transport properties . . . 21 4.2 Grid sizes . . . 21 4.3 Test matrix . . . 22

5.1 GCI test results . . . 23

5.2 Velocity dependency test results . . . 24

5.3 Time step dependency results . . . 26

(14)

(15)

CHAPTER

1 INTRODUCTION

Energy is a huge topic in today’s society and with it comes nuclear energy. Nuclear power plants need to be efficient and above all safe in all conditions. One of the more severe conditions is a large loss of coolant accident (LOCA), which occurs if one of the main water and steam supplying pipes breaks. When this happens, large amounts of steam are released into the containment (the so called dry well) of the power plant. In order to keep the power plant safe, because the containment is not created for the high pressures as they are present in pipes, the steam must be removed from the containment. The steam is guided through pipes into so-called suppression pools (wet well), which contain large amounts of water to condensate the steam and therefore relieve the containment from high pressures.

Experimenting with large setups and going through as many possible situations as possible is expensive and time consuming. However, with modern computer tools, it is possible to simulate these conditions without experimenting. The specialization in engineering that deals with this is called computational fluid dynamics (CFD). There are many companies and software tools available for more or less complete simulations of different phenomena. One of these software packages is OpenFOAM (OF). OF is an open source softwareR package that comes with a number of predefined solvers for different tasks. However, since the package is open source, new solvers can easily be created or existing solvers be modified. The remaining question is, why should OF perform those tasks, which other software packages maybe already perform.

OF is, due to its open source nature, scalable and completely transparent when it comes to the development of new solvers. The transfer of knowledge from academic use to industrial use works without hurdles and the software can run on small and large machines (multi-core clusters) without license restrictions.

This report describes the development of a solver for direct contact condensation within the OpenFOAM framework. Two different models are developed and a sensitivity analysis is performed on both solvers. The development is meant as a basis for more sophisticated solvers and is by far not yet a complete solver. The solvers are based on the assumption of saturated steam flowing onto a sub-cooled water surface. In order to fully enjoy this report, the reader should have a basic understanding of numerical simulation, thermodynamics and fluid dynamics.

Concerning the structure of this report, chapter 2 will give an overview of the fluid dynamical and thermo dynamical theory that is used to develop the solvers. In chapter 3 the program code will be explained. Chapter 4 deals with setting up and explaining the sensitivity analysis, followed by a chapter with the results. The last chapter concludes the report and gives recommendations for further development.

(16)

(17)

CHAPTER

2 THEORY BEHIND THE SOLVERS

In this chapter the theory, which makes up the solvers and the models used in the solvers is explained. It starts with the overall fluid dynamics model, which gives the framework for the phase change models, which are subsequently explained.

2.1 Fluid Dynamics Framework

The fluid dynamics are based on the solver interPhaseChangeFoam which is provided by the OpenFOAM 1.6 software. It is a solver for two immiscible and incompressible fluids with phase change features using a volume-of-fluid approach for interface tracking. The phase change models are interchangeable, and two phase change models were developed within the scope of this work (see section 2.2). The solver is based on the original work on cavitation models by Kunz [8], without the influence of non-condensible gases. The governing equations are not based on mass continuity equations but on volume continuity as can easily be seen by the division by the density, ρ, of the present phases as shown below.

1 ρm ∂p ∂τ + ∇ · u = m˙ +_{+ ˙}_m− 1 ρl + 1 ρv , (2.1)

where one can see a derivative with respect to τ , which is a pseudo-time derivative and is used for the time iterative solution employed by Kunz. Equation 2.1 represents the volume continuity equation, where p is the pressure, ρ the density, u the velocity and ˙m+ _{and ˙}_m− _{represent the mass flow between the phases. The} positive contribution is condensation and negative contribution comes from evaporation. Within the scope of this paper, only the condensation is modeled and evaporation is set to 0.

The accompanying momentum equation is given by equation 2.2 as ∂ρmu ∂t + ∂ρmu ∂τ + ρmu · ∇u = −∇p + µm∇ 2_{u + ρ} mg , (2.2)

where µm is the mixture dynamic viscosity and g is the gravitational acceleration. In order to track the interface, a phase continuity equation is introduced.

∂α ∂t + α ρm ∂p ∂τ + ∂α ∂τ + ∇ · αu = m˙ +_{+ ˙}_m− 1 ρl (2.3) Equation 2.3 represents the phase continuity equation and introduces the phase fraction α, which is bounded by 0 and 1, representing gaseous and liquid phase, respectively. The importance and function of this function is discussed in section 2.1.1, where the volume-of-fluid approach is discussed in more depth.

The mixture density, ρm, of the governing equations is given by

ρm= αρl+ (1 − α)ρv, (2.4)

where the subscripts l and v denote liquid and vapor, respectively.

(18)

4 Theory Behind the Solvers

2.1.1 Interface Tracking with VOF

As mentioned before, the solver uses a volume-of-fluid (VOF) approach to track the interface and with this information derive the momentum exchange between the phases. A volume-of-fluid method was first devised by Hirt [6] as a method to track the interface of a fluid dispersed in another continuous fluid.

The VOF uses a scheme that tracks the interface of the dispersed phase by utilizing the phase fraction within the cell. Every cell stores a value between zero and one, which denote empty and full of the continuous phase, respectively. This function α depends on time and position, meaning α = α(t, x), can be partially differentiated and moves with the fluid, as can be seen in equation 2.3, the average phase continuity equation. A visual representation is found in figure 2.1.

α=0

α=1

n

Figure 2.1: Graphical representation of the VOF method, depicting the normal vector of the interface and the cells containing phase 1 (blue, α = 1) and cells containing phase 2 (no color, α = 0). The cells containing the interface have 0 < α < 1.

Due to the fact that the function follows the region of fluid and not the interface it is a simple and economic way to track the interface [6], since only one extra variable must be stored and can be handled as any other transport equation. With the differentiability of the function comes the convenience that the curvature and the normal vector of the interface can be computed. Those values can then be used to calculate the interface forces.

2.2 Two Phase Change Models for One Task

In this section the different models which are used for further development of the solver are introduced. There are two of them, which are based on different approaches. The first one being a combustion like approach, treating the mass transfer between the phases in analogy to chemical combustion. The second model takes an approach, which is closer to real physical behavior, using inter-facial heat transfer to calculate the mass exchange between the phases.

For both models the energy and mass upon transfer must be conserved. This means, the energy that is transfered (with the mass) away from the saturated steam must be equal to the energy being received by the liquid.

(19)

Two Phase Change Models for One Task 5

2.2.1 The Combustion-like Approach

The idea behind this approach is, that the neither heat transfer coefficient nor the heat transfer area must be known. This results in a simplified way of describing condensation.

As a starting point the energy per mass of the two sides of the interface can be expressed as

ml· cp· (Tsat− Tinf) = mv· if g, (2.5)

where ml and mv are the mass of water and steam, respectively. Tsat and Tinf represent the saturation temperature and the temperature of the water, whereas cpis the specific heat of water at saturation conditions and if gis the latent heat of water. This can be rewritten for the mass of steam transferred

mv =ml· cp· (Tsat− Tinf) if g

. (2.6)

Equation 2.6 presents how much steam is condensed given a certain amount of water.

In order to describe the mass flow across the interface a so-called combustion frequency, fc, is introduced. This frequency determines how much steam can condense within a certain amount of time. The resulting equation for the mass flow is then

˙

m+=ml· cp· (Tsat− Tinf)

if g · fc. (2.7)

The obtained mass flow rate must be considered the maximum of course. There must never be more steam condensed than the cell contains. Therefore upon programming this into the solver, one must consider this. As can be seen in equation 2.7 for the mass flow between, ˙m+, the phases, all but one value, fc, are easily obtainable from the IF97 databases [14].

2.2.2 The Inter-Facial Heat Transfer Approach

The starting point for this approach is again a balanced energy equation. The balance must be at the interface between the two fluids, as shown in figure 2.2. The energy balance is given by

q = h · (Tsat− Tinf) · Al= ˙m+· if g, (2.8) where h is the heat transfer coefficient, as determined within the next subsection, Al is the inter-facial area between water and steam, ˙m+ is the mass condensation rate and if gis the latent heat. From equation 2.8 it can be seen that two parameters, namely the inter-facial area and the heat transfer coefficient are unknowns. They have to be determined or estimated within the simulation or beforehand. As the inter-facial area mainly depends on the computational cell size, the volume ratio of steam and water and the direction of the interface within the cell, this variable must be estimated during the simulation.

As for the heat transfer coefficient, h, this can be determined before hand, knowing what the approximate conditions during the simulation are, or it can be determined at runtime, based on some heat transfer coefficient model as described in the following part of this section.

Determining the Heat Transfer Coefficient

In order to complete the inter-facial heat transfer model, the heat transfer coefficient is needed. Two different approaches are taken. Both assume a very low velocity and laminar flow.

Pure Conduction In order to determine the heat transfer coefficient, h, some kind of model must be devised. There are many models already available, which are more or less complicated. For this approach an easy model based on the way the condensation is supposed to take place is devised. Taking a look at figure 2.2, one can see that the condensation takes places on a mostly flat water surface at very slow laminar flow. The assumptions for the heat transfer is that the heat flux across the interface between the two fluids works as in conduction. The Fourier law of conduction in one direction for the water side of the interface is given by [7]

q00= −λfdT

(20)

6 Theory Behind the Solvers Interface (Tsat) }dn

Steam (Tsat)

Water (Tinf)

q

Figure 2.2: Schematic of the energy transfer at the interface of steam and water as described by equation 2.8.

where q00 is the heat flux, λf is the thermal conductivity of the material (in this case water). dT_dn is the temperature gradient. Assuming that most of the water is at bulk temperature, Tinf, and that there is a thin boundary region close to the interface equation 2.9 can be rewritten to

q00= λfTsat− Tinf δb.l.

, (2.10)

where δb.l. represents the boundary thickness. As one can notice, the minus sign has disappeared. This is due to the fact that the temperatures have been switched as well.

The heat flux due to convection is given by [7]

q00= h · (Tsat− Tinf) . (2.11)

Equating equations 2.10 and 2.11 and solving for h gives the following expression for the heat transfer coefficient

h = λf

δb.l. . (2.12)

Therefore, the total mass transfer from vapor to liquid is given by inserting 2.12 into the energy balance from equation 2.8 and solving for ˙m+_{, which results in}

˙

m+=λf· Al· (Tsat− Tinf)

δb.l.· if g . (2.13)

Nusselt Number for Flat Rectangular Plate Condensation Bejan [3] gives the following relation for condensation of steam on a flat horizontally placed plate, which is a similar configuration to the one found within the inlet channel of the condensation tank.

N u = hlL λl = 1.079 L3_{· if g}_{· g (ρl}_{− ρv)} λl· νl(Tsat− Tinf) 15 , (2.14)

where N u is the Nusselt number, hl the heat transfer coefficient, λl the thermal conductivity of the liquid, L the characteristic length of the plate, g is gravity and ρl and ρv are the densities of the water and steam, respectively.

The only unknown in this relation is the characteristic length L. However, in this case it is set to the diameter of the inlet channel. Rewriting it gives the heat transfer coefficient hlas

(21)

Thermophysical Model for Temperature Analysis 7

Inserting equation 2.15 into 2.8 the mass transfer rate across the inter-facial boundary is given as ˙

m+= hl· Al· (Tsat− Tinf) if g

, (2.16)

where hl is given by equation 2.15.

Both models should give approximately the same results, but in order to keep the model as close as possible to physical behavior, the Bejan relation will be implemented into the solver.

Interface Area Model

For this model to work, the last thing left to do is to specify the heat transfer area, Al, from equation 2.16. In order to avoid complex geometric modeling a simple model is developed. It is assumed that most cells are regular in size. Meaning that the hexahedral cells have approximately the same length dimension in all direction. It is further assumed, that the size of the inter-facial area depends on the liquid volume fraction α, which ranges between 0 and 1, as can be seen in figure 2.3. This means that at α = 0.5, the area is largest and that at 0 and at 1 no inter-facial area is present. For a good estimate of the largest inter-facial area, the largest area within a regular cube is assumed.

α

1-α

Interface

Cell boundary

Figure 2.3: Description of a regular cell with internal inter-facial boundary between the fluid (α) and the gaseous phase (1 − α).

The maximum area Amax is therefore given as

Amax=√3 · V

2 3

cell, (2.17)

where Vcell is the cell volume. And the linear function for the increase and decrease of the area is given by Al=

2 · α 0 ≤ α < 0.5

2 · (1 − α) 0.5 ≤ α ≤ 1 . (2.18)

A graphical representation is given in figure 2.4.

2.3 Thermophysical Model for Temperature Analysis

The last part that will be added to the solver is a thermo-physical model for the water. In order to achieve this, the energy equation for the water side must be solved. The general governing energy equation for the liquid phase is given by [1]

∂

(22)

8 Theory Behind the Solvers

A

max

A

α

α=0.5

Figure 2.4: Graphical representation of the interface area model.

to the water, Γcis the condensation rate from steam to water per unit volume, if g is the latent heat and Evl is the energy transfer from vapor to liquid phase from conduction.

In order to have an easily solvable expression to implement, simplifications can be performed. The first simplification is that there is no general internal heat source, meaning QL = 0. Since only the water side should be solved for, the added factor of α may be dropped. Dividing both sides by cp in order to have only temperature leaves the following expression.

∂ρT

∂t + ∇ · ρuT − ∇

2_{ρκT =} Γcif g cp

+ Evl, (2.20)

where κ is the thermal diffusivity resulting from κ = λf/ρ/cp. Equation 2.20 leaves two terms that need to be modeled, namely Evland Γc. The mass condensation rate per unit volume is modeled as described in section 2.2. For the conduction, which is also the boundary condition for the water side of the equation, a model has to be devised. The temperature on the steam side should always be kept at saturation temperature Tsat. When discretizing equation 2.20 around the point of interest (subscript p), the result is

apTp+X all

Tnan =Γcif g

cp + Evl, (2.21)

where the a gives the coefficients of all cells and the subscript n gives the nodes around node p. At the boundary, the cell temperature should be at saturation temperature, giving Tp≈ Tsat. This can be achieved by modeling Evl with large numbers in a way that

Evl= SlargeTsat− TpSlarge, (2.22)

where Slargeis a large number, much larger than the reciprocal of the precision of the computational system (for double precision, the precision is ∼ 10−8_{). Inserting equation 2.22 into equation 2.21 and solving for T}

p results in Tp= SlargeTsat ap+ Slarge − P allTnan ap+ Slarge ≈ Tsat (2.23)

(23)

CHAPTER

3 PROGRAM CODE AND USAGE OF THE SOLVER

This chapter treats the program code and gives an explanation of which parts have been edited and/or newly written to accomplish the goal of the new solver interPhaseChangeCondenseTempFoam. The programming was done in C++ with help of the libraries that are included in the OpenFOAM 1.6.x (OF) package and compiled with a gcc version higher than 4.3.3.

interPhaseChangeCondenseTempFoam – runtime commands

Reading: PISO

Time

Output: Courant number

delta t (calculated depending on CFL-condition) alphaSubCycle.H

Performs analysis of α

Creates ρ (from weighted average of mixture) Creates (flux)ϕ Runs alphaEqn.H C o rr e ct io n lo o p AlphaEqn.H

α transport equations with source terms from phase change models Solves bounded α equation

UEqn.H

Create tensor for solution of U with source from phase change model If momentum predictor is set, solves the velocity equation

P

IS

O

lo

o

p pEqn.HUses U tensor variable to create U

Sets up with ϕ U and volume flux from phase change model Sets up pressure equation with and volume flux ϕ

Solves pressure equation Updates velocity field accordingly

continuityErrs.H

Corrects the introduced continuity errors

TEqn.H (if energy equation is switched on)

Sets up energy equation with mass flux and

energy source term from phase change model for water side Solves the energy equation

ru n tim e lo o p

Figure 3.1: A graphical representation of the runtime command structure of interPhaseChangeCondenseTempFoam.

(24)

10 Program Code and Usage of the Solver

The solver is an enhanced version of the solver interPhaseChangeFoam, which is delivered with OF. The orig-inal solver includes three different phase change models that are used to simulate cavitation. The new solver omits these phase change models and only includes the two new phase change models that were developed in chapter 2. In order to comply with the programming grammar, the solvers are called AlbaNovaCombustion and AlbaNovaInterface, representing the combustion analysis approach (see section 2.2.1) and the inter-facial heat transfer approach (see section 2.2.2), respectively. A graphical representation of the code can be seen in figure 3.1. The figure omits the phase change models and only states at which part the phase change models give data back to the main code, since the phase change models are a class on their own.

The newly added features compared to the original code include a preparation the solving of the energy equation for the water side for the main class interPhaseChangeCondenseTempFoam and two new classes called AlbaNovaCombustion and AlbaNovaInterface. The energy equation is not completely implemented, see section 3.2 for more information.

3.1 The Solver Code Explained

Within the following section, the line numbers correspond to the line numbers in the solver source files, as they can be found in the solver and the appendix.

All equations that can be solved separately, e.g. pressure, velocity and phase continuity, are treated in separate .H files containing the setup of the necessary variables and the equation, as well as solving of the equation. OF handles differential equations and the solving of such, such as transport equations, in a template class called fv<type>Matrix. This class allows an easy set up and solving of the equation with different solution schemes which can be set on runtime. The code to solve a differential equation in OF is similar to the mathematical representation as it is given in the theory chapter.

3.1.1 Phase Continuity Code

As one can see in figure 3.1, the solver starts by solving the continuity equation for the liquid volume fraction. This is done in two steps. At first the alpha sub cycle is started. It first creates a field for the mass flux and sets it to zero. This is done with every start of the time loop. Meaning that a new mass flux is written at every start. 01: surfaceScalarField rhoPhi 02: ( 03: IOobject 04: ( 05: "rhoPhi", 06: runTime.timeName(), 07: mesh 08: ), 09: mesh, 10: dimensionedScalar("0",dimensionSet(1,0, -1, 0,0),0) 11: ); 12:

The next fields it creates is the flux due to interface compression phic in 38:

39: surfaceScalarFieldphic=mag(phi/mesh.magSf());

40: phic=min(interface.cAlpha()*phic,max(phic));

41:

After creating all those necessary fields for the phase continuity equation, the phase continuity equation is set up and solved in alphaEqn.H. The first step is to set up the relative flux. Since the solver relies on mixed properties and one equation models. This is achieved by multiplying the compression flux with the normalized normal vector of the interface.

04:

(25)

The Solver Code Explained 11

06:

The solver enters a loop, in which first the actual explicit flux due to phase fraction α, phiAlpha, and then the sources which result from the phase change are set up and inserted into the bounded multiphase solver that is supplied by OpenFOAM. The loop runs a predefined number of correction iterations, which is set on run time by the user.

The explicit flux of α is set up via 09: surfaceScalarField phiAlpha= 10: fvc::flux 11: ( 12: phi, 13: alpha1, 14: alphaScheme 15: ) 16: +fvc::flux 17: (

18: -fvc::flux(-phir,scalar(1) -alpha1,alpharScheme),

19: alpha1,

20: alpharScheme

21: );

22:

where the first term is due to the water phase and the second term is due to the gas phase. The “schemes” determine how to discretize the field.

During the next step, the source terms for evaporation and condensation are retrieved from the class that handles the phase change. This can be one of the models handling phase change (AlbaNovaCombustion or AlbaNovaInterface) and will explained at a later stage. The template data type Pair stands for two values, which can be in any form of a field or a scalar, depending on the definition given after tmp. The values can be read as from an array as shown below in line 25 and 26.

23: Pair<tmp<volScalarField> > vDotAlphal=

24: twoPhaseProperties->vDotAlphal();

25: const volScalarField&vDotcAlphal=vDotAlphal[0]();

26: const volScalarField&vDotvAlphal=vDotAlphal[1]();

27: 28: volScalarField Sp 29: ( 30: IOobject 31: ( 32: "Sp", 33: runTime.timeName(), 34: mesh 35: ), 36: vDotvAlphal-vDotcAlphal 37: ); 38: 39: volScalarField Su 40: ( 41: IOobject 42: ( 43: "Su", 44: runTime.timeName(), 45: mesh 46: ),

47: // Divergence term is handled explicitly to be

48: // consistent with the explicit transport solution

(26)

51: );

The two source terms volScalarField Su and volScalarField Sp represent the sources due to the differ-ence between condensation and vaporisation and due to convective transport of the water phase, respectively. After forming and establishing all fields necessary for the MULES solver that is supplied with OpenFOAM [9, MULES], it is run. The MULES code is described as a “Multidimensional universal limiter with explicit solution”, which solves convective-only transport equations with upper and lower limits, as required by the α-transport equation. The call for the MULES solver is given by

55: MULES::implicitSolve(geometricOneField(),alpha1,phi, phiAlpha,Sp,Su,1,0);

The next step is to correct the mass flux for the new found α-values. Which is done in line 57:

57: rhoPhi+=

58: (runTime.deltaT()/totalDeltaT)

59: *(phiAlpha*(rho1-rho2) + phi*rho2);

After the correction loop, the condensation rate is written into the condensation rate field, which can be used at a later stage, for example within the source term of the temperature equation.

63: 64:

65: Pair<tmp<volScalarField> >mDotAlphal=

66: twoPhaseProperties->mDotAlphal();

67: constvolScalarField&mDotcAlphal =mDotAlphal[0]();

68: mDotCondense=mDotcAlphal*(1-alpha1); // multiply by (1-alpha1) because given back that way

After leaving the alphaEqn.H file the mixture density is updated with the appropriate α values, as line 52 in the alphaEqnSubCycle.H shows.

53: rho== alpha1*rho1+ (scalar(1) -alpha1)*rho2;

3.1.2 Velocity Equation Code

After solving the phase continuity the equation for solving the velocity distribution is set up in UEqn.H. Due to the coupling of the velocity and pressure equation, the velocity equation is only set up, but not solved instantly, as can be seen below.

01: surfaceScalarField muEff

02: (

03: "muEff",

04: twoPhaseProperties->muf()

05: +fvc::interpolate(rho*turbulence->nut())

06: ); 07: 08: fvVectorMatrix UEqn 09: ( 10: fvm::ddt(rho,U) 11: +fvm::div(rhoPhi,U)

12: -fvm::Sp(fvc::ddt(rho) +fvc::div(rhoPhi),U)

13: -fvm::laplacian(muEff,U)

14: - (fvc::grad(U) &fvc::grad(muEff))

15: //- fvc::div(muEff*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))

16: );

(27)

The Solver Code Explained 13

product rule of differentiation in order to obtain only the derivative with respect to u. Consider the following derivative as it is given in line 10,

d(ρ · u) dt = u dρ dt + ρ du dt . (3.1)

In line 12, however, the time derivative of rho is multiplied with U and subtracted. Leaving only

ρdu

dt (3.2)

within the equation. The same is true for all the other derivatives that include rho (muEff contains a multiplication with rho in line 5). This avoids a problem that would arise due to the interface within the field. At the interface, the density changes abruptly, in real life this is a discontinuity. Due to this abrupt change, the gradients become very steep, and would introduce velocity sources that are not actually present across the interface, because the fluid in front of the interface moves with the same velocity as the fluid behind the interface.

3.1.3 Pressure Equation Code

With the velocity equation fully set up (note: velocity is not yet solved for), the solver enters the pressure solution step within pEqn.H. The first step is to predict the velocity, since it is needed to set up the flux for the pressure equation. This is achieved in lines 1 to 5.

01: {

02: volScalarFieldrUA =1.0/UEqn.A();

03: surfaceScalarField rUAf=fvc::interpolate(rUA);

04:

05: U =rUA*UEqn.H();

At first the central coefficient of the velocity equation matrix are retrieved and inverted (line 2), where UEqn.A() give the central coefficients [9, fvMatrix]. Line 5 then predicts the velocity by multiplying the reciprocal of the central coefficients with the operation source from the same matrix.

With the predicted velocity the flux that depends solely on the velocity is created. 07: surfaceScalarField phiU

08: (

09: "phiU",

10: (fvc::interpolate(U) &mesh.Sf())

11: +fvc::ddtPhiCorr(rUA,rho,U,phi)

12: );

The term ddtPhiCorr() introduces a correction for the divergence of the face velocity field (the velocity field was interpolated from cell central stored to cell face stored) by taking out the difference between the interpolated velocity and the real flux [9].

Within the next step, the complete flux is constructed, by taking into account not only the flux due to velocity, but also due to other forces.

16: phi=phiU+

17: (

18: fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()

19: +fvc::interpolate(rho)*(g &mesh.Sf())

20: )*rUAf;

where the first term accounts for the flux due to velocity, the second term accounts for the interfacial forces (sigmaK) and the third term account for buoyancy forces (ρ · g).

(28)

22: Pair<tmp<volScalarField> >vDotP=twoPhaseProperties->vDotP();

23: constvolScalarField&vDotcP =vDotP[0]();

24: constvolScalarField&vDotvP=vDotP[1]();

Following getting the necessary values for the source, the pressure equation is set up. This is done in a for-loop, which calls different steps depending on the number of iterations and how far the for-loop has progressed. The steps in line 38 and 47 are always performed, and are necessary to solve the problem as they represent the final solving and the correction of the flux, respectively.

26: for(intnonOrth=0;nonOrth<=nNonOrthCorr;nonOrth++)

27: {

28: fvScalarMatrix pEqn

29: (

30: fvc::div(phi) -fvm::laplacian(rUAf,p)

31: - (vDotvP-vDotcP)*pSat+fvm::Sp(vDotvP-vDotcP,p)

32: );

33:

34: pEqn.setReference(pRefCell,pRefValue);

35:

36: if(corr==nCorr-1&&nonOrth==nNonOrthCorr)

37: {

38: pEqn.solve(mesh.solver(p.name() +"Final"));

39: }

40: else

41: {

42: pEqn.solve(mesh.solver(p.name()));

43: }

44:

45: if(nonOrth==nNonOrthCorr)

46: {

47: phi+=pEqn.flux();

48: }

49: }

Finally, an updated velocity field is calculated by only adding/removing the change that was made to phi, without the contribution to phi due to the velocity field (phiU).

3.2 The Energy Equation Problem

For the energy equation, no solution could be found to implement it for one side only. The theory behind the equation assumes that a discrete interface is available within the code. This, however, is not the case with the framework provided. The solver is prepared to be compiled with an energy equation added to the source code and has the necessary switches built-in on the user side. Nevertheless, at this point the energy equation it not functional and it might be necessary to provide a two equation model in order to implement the energy equation.

If an implementation can be found, it only needs to be programmed into the TEqn.H file and the code recompiled.

3.3 The Phase Change Model Classes

(29)

The Phase Change Model Classes 15

develop the model to include evaporation due to boiling or alternative mechanisms.

Both functions return mass flow rates as a coefficient to be multiplied with a pressure term (mDotP()) or an α-term (mDotAlphal()). mDotP() is used to calculate the pressure changes that occur due to this condensation mass flow, whereas mDotAlphal() serves its purpose in the bounded α-phase transport equation. The functions were overwritten to include the temperature based phase change models described in chapter 2. These mass flow function have a non-virtual counterpart, which gives back the volume flow rates as coefficients, instead of the mass flow rates. This is necessary since, as described in the theory section, the fluid dynamical framework is based on volume continuity.

In figure 3.2 the class diagram for the two new classes AlbaNovaInterface and AlbaNovaCombustion are shown. As one can see, the models directly inherit most of their functions from the base class. There are no extra public functions implemented, but the virtual functions were overwritten, as described before. For the interface model, two new private functions have been added, which calculate the heat transfer coefficient, hCoeff(), and the heat transfer area, interfaceArea(). The decision to run separate functions for these two parameters is based on the fact, that with this method, the underlying models are easily exchangeable. The full source code of the two models can be found in appendix B and appendix C, respectively.

AlbaNovaCombustion - combFrequency_ : dimensionedScalar - Tsat_ : dimensionedScalar - Tinf_ : dimensionedScalar - ifg_ : dimensionedScalar - cp_ : dimensionedScalar - mcCoeff_ : dimensionedScalar AlbaNovaInterface - lambdaF_ : dimensionedScalar - Tsat_ : dimensionedScalar - Tinf_ : dimensionedScalar - ifg_ : dimensionedScalar - charLength_ : dimensionedScalar - g_ : dimensionedScalar - mcCoeff_ : dimensionedScalar - interfaceArea() : tmp<volScalarField> - hCoeff() : tmp<volScalarField> phaseChangeTwoPhaseMixture # phaseChangeTwoPhaseMixtureCoeffs_ : dictionary # pSat_ : dimensionedScalar + phaseChangeTwoPhaseMixture() + ~phaseChangeTwoPhaseMixture() + pSat() : dimensionedScalar

+ virtual mDotAlphal() : Pair<tmp<volScalarField> > + virtual mDotP() : Pair<tmp<volScalarField> > + vDotAlphal() : Pair<tmp<volScalarField> > + vDotP() : Pair<tmp<volScalarField> > + correct()

+ virtual read() : bool

Figure 3.2: Class diagram of the AlbaNova phase change models with inheritance from the base class phaseChangeTwoPhaseMixture. “-” denotes private attribute/operation, “#” denotes protected attribute/operation and “+” denotes public attribute/operation.

3.3.1 Functions of AlbaNovaInterface

The first function that was newly written for AlbaNovaInterface is for the heat transfer area A, as outlined below.

069: Foam::tmp<Foam::volScalarField>

070: Foam::phaseChangeTwoPhaseMixtures::AlbaNovaInterface::interfaceArea()const

071: {

072: // return the interfacial area based on model for interfacial area

073: // returns dimensions Area

074:

075: // model based on regular volume cells, taking the largest cut area

076: // as maximum for area, linear increase and decrease with alpha

077: const volScalarField&cellVolume=

078: alpha1 .db().lookupObject<volScalarField>("cellVolu");

079: volScalarFieldlimitedAlpha1=min(max(alpha1 ,scalar(0)),scalar(1));

080:

(30)

082: areaFactor("areaFactor",dimensionSet(0,2,0,0,0,0,0),0.0);

083:

084: volScalarFieldinterfaceArea=alpha1 * areaFactor;

085: volScalarFieldmaxArea=alpha1 *areaFactor;

086:

087: maxArea=sqrt(3.0)*pow(cellVolume,(2.0/3.0));

088: returntmp<volScalarField>

089: (

090: (neg(limitedAlpha1-0.5)*maxArea*2.0*limitedAlpha1) +

091: (pos(limitedAlpha1-0.5)*maxArea*(-2.0*(limitedAlpha1-1.0)))

092: );

093: 094: }

Since this class is not part of the runtime environment, it can only gain access to variables via input or via input variables, or they need to be created somewhere along the programming code. This means that the cell volume cellVolu, needs to be accessed via the passed variable alpha1. The variable contains a data base that can give access to the run time variables. The access to those is created in lines 77 and 78. After this, two new volume scalar fields are created to store the interface and maximum area of the each cell, they are initialized with 0.

The other newly introduced function is the function to calculate the heat transfer coefficient. The heat transfer coefficient is based on equation 2.15. The function is given by the following.

096: Foam::tmp<Foam::volScalarField>

097: Foam::phaseChangeTwoPhaseMixtures::AlbaNovaInterface::hCoeff()const

098: {

099: // from [Bejan, Convection heat transfer, 1995]

100:

101: return

102: (

103: 1.079* lambdaF / charLength *pow(

104: (pow3(charLength) * ifg * g * (rho1()-rho2()) ) /

105: (lambdaF * nuModel1().nu() * (Tsat -Tinf) ) ,0.2)

106: );

107: }

3.4 Usage of interPhaseChangeCondenseTempFoam

For the following section it is assumed the reader has a basic understanding of how to set up simulation cases in OF, as only the most important steps for setting up a case for the new solver are explained. For more information on how to set up cases, the reader is also referred to the user guide of OF [10, ch. 1 and ch. 2] and the test cases on the CD.

Each of the models has its own set of parameters that are to be adjusted depending on the case that is run. Three files need to be checked for consistency. The first one if the transportDict file, which contains all pa-rameters concerning water and steam properties, as well as general fluid properties and switches. Both models use the general switches and overall fluid properties. The phases need their respective density and kinematic viscosity set, where phase one is the continuous phase. The switch for phaseChangeTwoPhaseMixture can be set to either AlbaNovaInterface or to AlbaNovaCombustion, as those are the only models that are available with the solver. For the coefficient for each model tables 3.1 and 3.2 are to be consulted for the AlbaNovaInterface and AlbaNovaCombustion model, respectively. The switch for the temperature equation, temperatureEquation needs to be set to off, since the equation has not been implemented yet.

In order to have the solver work properly, the solvers for the following variables need to be specified within the fvSolution dictionary: alpha1, U, p, T (only in case the solver needs to run the temperature equation, meaning temperature is switched on), pCorr, pFinal and PISO. Also, relaxation factors may be specified for alpha1, U, p and T, but these are optional.

(31)

Usage of interPhaseChangeCondenseTempFoam 17

Table 3.1: Input coefficient for AlbaNovaInterface with sample values for ambient conditions. The dimensions correspond to the dimensions that need to be given to OF.

Label Dimensions Value Name

lambdaF [1 1 -3 -1 0 0 0] 0.6 Thermal conductivity of water Tsat [0 0 0 1 0 0 0] 373 Saturation temperature Tinf [0 0 0 1 0 0 0] 293 Bulk temperature water

ifg [0 2 -2 0 0 0 0] 2256.4e3 Latent heat charLength [0 1 0 0 0 0 0] 0.18 Characteristic length

gravity [0 1 -2 0 0 0 0] 9.81 Gravitational acceleration nuL [0 2 -1 0 0 0 0] 1.004e-6 Kinematic viscosity water

cp [0 2 -2 -1 0 0 0] 4.1813e3 Heat capacity water Pr [0 0 0 0 0 0 0] 5.5 Laminar Prandtl number

Table 3.2: Input coefficient for AlbaNovaCombustion with sample values for ambient conditions. The dimensions correspond to the dimensions that need to be given to OF.

Label Dimensions Value Name

combFrequency [0 0 -1 0 0 0 0] 0.8 Combustion frequency cp [0 2 -2 -1 0 0 0] 4.1813e3 Heat capacity water Tsat [0 0 0 1 0 0 0] 100 Saturation temperature Tinf [0 0 0 1 0 0 0] 20 Bulk temperature water

ifg [0 2 -2 0 0 0 0] 2256.4e3 Latent heat lambdaF [1 1 -3 -1 0 0 0] 0.5 Thermal conduction water

(32)

(33)

CHAPTER

4 SENSITIVITY ANALYSIS

The following chapter describes the sensitivity analysis conducted on the code that is described in chapter 3. The main input parameter that are believed to change the behavior of the code are the velocity, the cell size and the time step. The OF framework handles automatic time step adjustment on runtime, but this can also be switched off.

4.1 The Test Cases

Both models are tested and compared. As for the combustion model, the combustion frequency needs to be tuned to a value, so that the condensation rate is approximately the same as it is for the interface area model. This is done only once for the base case and kept at this value for all following cases. After the base cases have been set up, the input parameters concerning the condensation and transport behavior will not be altered. The observed quantity is the integrated condensation rate over all cells, ˙m+_.

The sensitivity analysis consists of 18 test cases. All tests are carried out on the same model with different grid sizes, different inlet velocities and different time step setups. Each of the tests is run with three different parameter values, except for the velocity. There three different velocities are chosen and a fourth velocity from the grid convergence test is added.

Except for the time step dependency tests, the so called Courant number is kept at 0.25. The Courant number gives the ratio of the time step, ∆t, to the characteristic convection time u/∆x [5]. The ratio is used to automatically determine the time step on runtime without input by the user.

In figure 4.1 the model that is used throughout the tests is shown. The model is based on the POOLEX facility set up with a reduced main tank diameter from 2.4m to 1.33m and a reduced height of the wet well compartment from 3.61 to 3.11m [11]. This was done in order to decrease the number of cells needed to describe the model, since the steam is not expected to leave the pipe which reduces the impact on the mesh surrounding the pipe.

4.1.1 Grid Convergence Index

For the grid test a parameter called the grid convergence index (GCI) is calculated based on the method given by Celik [4] as outlined below. The first step is to specify a representative grid size, l, by

l = " 1 N N X n=1 Vn #1/3 , (4.1)

(34)

20 Sensitivity Analysis Ø 1.33m 3. 11 m 2. 55 m Ø 0.22m 0.2m Water Atmosphere Inlet Pipe Steam

Figure 4.1: Simulation model used for all simulation, based on the POOLEX facility set up with reduced volume.

order, p, of the of the code can be calculated by the following three equation using a fixed point iteration.

p = 1 ln(r21)|ln |ε32/ε21| + q(p)| , (4.2) where q(p) is given by q(p) = ln r p 21− s rp₃₂− s (4.3) and s being given as

s = 1 · sign(ε32/ε21) , (4.4)

where the absolute errors between the global variables of different grids are ε32 = ˙m+₃ − ˙m+₂ and ε21 = ˙

m+₂ − ˙m+₁. The condensation rates are obtained on the different grid size test runs. The next step is to calculate the extrapolated values by

˙ m+,21_ext = r p 21m˙ + 1 − ˙m + 2 rp₂₁− 1 . (4.5)

In order to now find the GCI, the approximate relative error between the finest and second finest grid are calculated, by e21_a = ˙ m+₁ − ˙m+₂ ˙ m+₁ (4.6) and inserted into the equation giving the grid convergence index,

GCI21_{f ine}= 1.25e 21 a

r₂₁p − 1 . (4.7)

(35)

The Test Cases 21

4.1.2 Velocity and Time Step Dependency Test

In case of the velocity dependency test, the influence of the velocity on the condensation rate is examined. This is done by changing the inlet velocity. It is expected that due to higher velocities a break up of the interface can occur, which might change the condensation behavior. One lower and two higher velocity are selected. The lower velocity is set to 0.05 m/s, which simulates falling steam. The higher velocities are 0.2 m/s and 0.5 m/s, which simulate a slowly pumped insertion of steam into the water basin.

In the time step dependency test, the time step behavior is changed from automatic to fixed time steps. The influence of this change on the condensation rate is also examined. From the aforementioned tests with automatic time step control, it is known that in order to fulfill a Courant number of 0.25 on the finest grid, a time step of 0.001s is sufficient. The chosen time steps are 0.001, 0.005 and 0.01.

4.1.3 The Test Cases

All test cases run with a predetermined set of parameters for it transport properties. The transport properties were acquired from the IF97 database for water properties [14]. The pressure was chosen close to ambient and a sub-cooling of 80K was chosen. The transport properties are listed in table 4.1.

Table 4.1: Transport properties of the test cases, from the IF97 data base [14].

Property Value psat 1.1 bar Tsat 375.29 K ρl 997.71 kg/m3 ρv 0.65 kg/m3 if g 2250.4·103J/kg Tl 295.29 K λl 0.603 W/m/K cp 4183.3 J/kg/K νl 9.5·10−7 m2/s νv 1.9·10−5 m2/s

As mentioned before, the combustion model was tuned and the combustion frequency was set to 0.8 Hz. The characteristic length for the heat transfer coefficient, L, was set to the pipe diameter, 0.22m.

The chosen grids are given in table 4.2 and the test matrix is given in table 4.3. The ratios between the grids are given as r21= 1.3089 and r32= 1.6348. Where the indices represent the grid numbers which are given in table 4.2.

Table 4.2: Grid sizes for obtaining the GCI, including the representative grid size in m. Grid name Cell count Representative grid size [m]

Grid 1 626914 1.981 · 10−2

Grid 2 279808 2.593 · 10−2

(36)

22 Sensitivity Analysis

Table 4.3: Test matrix, where auto means that the time step is regulated automatically corresponding to a Courant number of 0.25

Case # Model Grid Inlet velocity Time step

(37)

CHAPTER

5 RESULTS

This chapter treats the results of the sensitivity analysis. It treats all tests separately. Firstly, the GCI test is presented, followed by the tests for velocity and time step size dependency.

5.1 GCI Test Results

The grid convergence test runs each ran for 5 seconds real time over which the global mass condensation rate, ˙m+, was recorded and averaged over time. The average is used, because as one can see from the graphs of the condensation rate against time in figure 5.1, the condensation rate changes with passing through the cells in a sinusoidal fashion. This makes it hard to compare the various time steps and therefore an averaged approach is considered.

From the graph it becomes apparent that with increasing grid density, the global condensation rate of the interface model varies about a constant value. This value is reached shortly after initiating the simulation, with a respective shorter time for a finer grid. For the combustion model however, no stable value is reached. The condensation rate rises steadily over the 5 seconds of the simulation with spikes which are similar to the sinusoidal variations seen in interface model graph in figure 5.1a.

The calculated apparent order of the interface model and the combustion model are 5.7 and 3.58, respectively. The calculated GCI for the combustion model results in 82.35% and for the interface model in 2.02%, meaning that for the finest grid with the combustion model, the acquired value is accurate to 82.35% and for the interface model to 2.02%, neglecting the computer rounding errors. The full test results can be viewed in table 5.1.

(38)

24 Results 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3x 10 −3 Total Time [s]

Total Condensation Rate, [kg/s]

Grid 3, case 1 Grid 2, case 2 Grid 1, case 3 Time smoothed case 2 Time smoothed case 3

(a) Interface Model

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5x 10 −3 Total Time [s]

Total Condensation Rate, [kg/s] Grid 3, case 4

Grid 2, case 5 Grid 1, case 6

(b) Combustion Model

Figure 5.1: Graphs showing the condensation rate vs. time on different grids used for GCI calculations. In figure (a) the time smoothed combustion rates have been added for comparison.

5.2 Velocity Dependence Test Results

As for the other test runs, the velocity test ran for 5 second, over which the average condensation rate is calculated, these results are listed in table 5.2. The three different velocities chosen are complimented by a fourth, which is taken from the GCI tests (cases 3 and 6). Figure 5.2 shows the development of the condensation rate over time.

One can see that with higher velocities, the condensation rate is larger as well. This can be explained by looking at the distribution of the interface in figures 5.3 and 5.4. Both models base their condensation on the property of 0 < α < 0, meaning with more cells being occupied by α unequal to unity or zero, condensation occurs at a larger rate.

For the interface model, a theoretical maximum condensation rate can be calculated. This maximum is 4.429·10−4kg/s. This mass flow corresponds to an inlet velocity of 1.79·10−2m/s. As one can see, the theoretical maximum is lower than than the maximum that is obtained during the lowest velocity of 0.05m/s, which results in a steady condensation.

Table 5.2: Test results of the velocity dependency test. Case # Avg. ˙m+ _[kg/s] Interface model 3 2.2574 · 10−3 7 1.8246 · 10−3 8 2.8212 · 10−3 9 3.8715 · 10−3 Combustion model 6 0.8576 · 10−3 13 0.6746 · 10−3 14 1.105 · 10−3 15 1.6147 · 10−3

5.3 Time Step Dependence Test Results

(39)

Time Step Dependence Test Results 25 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5x 10 −3 Total Time [s]

Interface

Case 7, 0.05m/s Case 3, 0.1m/s Case 8, 0.2m/s Case 9, 0.5m/s Time smoothed rate Theoretical Max

(a) Interface Model

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2x 10 −3 Total Time [s]

Total Condensation Rate, [kg/s] Case 13, 0.05m/s

Case 6, 0.1m/s Case 14, 0.2m/s Case 15, 0.5m/s

Figure 5.2: Graphs showing the global condensation rate development over time depending on different inlet veloc-ities. In figure (a) the time smoothed condensation rates as well as the theoretical maximum condensation rate have been added for comparison.

Figure 5.3: The interface distribution during the velocity dependency tests at time 2.5 s after initiation for all four cases for the interface model.

(40)

26 Results

Figure 5.4: The interface distribution during the velocity dependency tests at time 2.5 s after initiation for all four cases for the combustion model.

Table 5.3: Overview of the average condensation rates for the time step dependency test including one auto adjusted time step run.

Case # Avg. ˙m+ [kg/s] Interface model 3 2.2574 · 10−3 10 2.2593 · 10−3 11 2.2482 · 10−3 12 2.2471 · 10−3 Combustion model 6 0.8576 · 10−3 16 0.8541 · 10−3 17 0.8732 · 10−3 18 0.8809 · 10−3 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3x 10 −3 Total Time [s]

Case 10, ∆ t = 0.001s Case 11, ∆ t = 0.005s Case 12, ∆ t = 0.01s Case 3, auto

(a) Interface Model

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.2 0.4 0.6 0.8 1 1.2x 10 −3 Total Time [s]

Case 16, ∆ t = 0.001s Case 17, ∆ t = 0.005s Case 18, ∆ t = 0.01s Case 6, auto

(41)

CHAPTER

6 CONCLUSION AND RECOMMENDATIONS

This chapter presents the conclusions that have been drawn from the results with respect to the theory chapter. The latter part of this chapter presents the recommendations and how further research can be conducted using the results of this work.

6.1 Conclusions on the Phase Change Models

The first conclusion that can be drawn is that the one equation model that is implemented with the VOF method employed by OF is unsuitable for implementing the energy equation. The interface spread that has been observed in the velocity dependency tests (see section 5.2) of both models, makes the method that was supposed to be employed invalid. For this to work, a sharper interface must be present.

The results of the GCI, the velocity dependency and the time step dependency test give an indication of the physical validity of the proposed models. The range of condensation rates is close to what is expected from previous experiments, like TransAT at LUT, Finland [2]. However, due to the absence of extensive test data, the validity of the models cannot be completely confirmed.

The interface model shows great time step stability, because even with a time step increase by one magnitude, the results only differ by 0.54% from the smallest time step used, which gives a Courant number of approxi-mately 0.25. The grid convergence index shows that convergence for finer grids exists, with a GCI of 2.02% for the finest grid used within the tests carried out. However, with higher velocities, the condensation rate increases with higher velocities. This is a valid physical behavior (higher velocities give rise to turbulence, which increases the heat transfer coefficient), but should have not occurred with the developed model. The behavior stems from the fact, that with higher velocity, there is a higher numerical diffusion from the original VOF framework employed by OF. Overall, including visual inspection, it can be said that the developed interface model behaves in a physical sense and can be used for further development.

The combustion model on the other hand does not show a physical behavior. It was not possible to show a grid convergence. The GCI remained at 82.35% for the finest grid. In addition, there is no stable condensation rate present, as should be expected. This means that the current implementation of the combustion model is not capable of predicting direct contact condensation behavior.

6.2 Recommendations

There are several aspects that need to be addressed in further work on this topic. The first one being the validation of the interface model with data from experiments. This is crucial, since otherwise the model can not be used for condensation prediction. After validation or even during validation of the interface model with data from experiments reflecting the conditions it is implemented for now (low velocity steam injection), improved models for the heat transfer and the interface model should be devised.

(42)

28 Conclusion and Recommendations

advanced through the cell in order to have a better estimate of the interface area within a cell. The next step would be to write a new class for the heat transfer coefficient in order to chose the coefficient modeling before the simulation starts. This however, can also lead to an automatic selection of the heat transfer coefficient model based on velocity consideration by the phase change model itself.

(43)

BIBLIOGRAPHY

[1] Henryk Anglart. Thermal Hydraulics in Nuclear Systems (KTH, Stockholm, 2008).

[2] Henryk Anglart, Djamel Lakehal, Vesa Tanskanen, and Mikko Ilvonen. BWR thermal-hydraulic issues: dryout, core transients and steam injection to pressure-suppression pool, Dec. 2009. NURISP presenta-tion, retrieved from personal conversation with Henryk Anglart.

[3] Adrian Bejan. Convection heat transfer (Wiley, New York, 1995), 2nd ed. edn. ISBN 9780471579724. [4] Ismail B. Celik. Procedure for estimation and reporting of discretization error in CFD applications.

Journal of Fluids Engineering, 130, Jul. 2008. Editorial of JFE, Retrieved March 12, 2010.

[5] Joel Ferziger. Computational methods for fluid dynamics (Springer, Berlin ;;New York, 2002), 3rd, rev. ed. edn. ISBN 9783540420743.

[6] C Hirt and B Nichols. Volume of fluid (VOF) method for the dynamics of free boundaries1. Journal of Computational Physics, 39(1):pp. 201–225, 1981. ISSN 00219991. doi:10.1016/0021-9991(81)90145-5. [7] M Kaviany. Principles of heat transfer (Wiley, New York, 2002). ISBN 9780471434634.

[8] R Kunz. A preconditioned NavierStokes method for two-phase flows with application to cavitation predic-tion. Computers & Fluids, 29(8):pp. 849–875, 2000. ISSN 00457930. doi:10.1016/S0045-7930(99)00039-0.

[9] OpenCFD Limited. OpenFOAM programmer’s c++ documentation.

http://foam.sourceforge.net/doc/Doxygen/html/, 2010. Retrieved May 10, 2010.

[10] OpenCFD Limited. OpenFOAM user guide. http://www.openfoam.com/docs/user/, 2010. Retrieved May 05, 2010.

[11] Eija Puska. Safir2010 : the Finnish research programme on nuclear power plant safety 2007-2010 : interim report (VTT, [Espoo], 2009). ISBN 9789513872663.

[12] John W. Slater. Examining spatial (Grid) convergence. http://www.grc.nasa.gov/WWW/wind/valid/ tutorial/spatconv.html, Jul. 2008. Retrieved May 25, 2010.

[13] Bjarne Stroustrup. Programming : principles and practice using C++ (Addison-Wesley, Upper Saddle River NJ, 2009). ISBN 9780321543721.

[14] Milton Venetos. Excel engineering. http://www.x-eng.com/, 2007. Downloadable database.

(44)

(45)

APPENDIX

A

MAIN PROGRAM SOURCE CODE

001: /*--–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–--*\

002: ========= |

003: \\ / F ield | OpenFOAM: The Open Source CFD Toolbox

004: \\ / O peration |

006: \\/ M anipulation |

007: --–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–

008: License

009: This file is part of OpenFOAM.

010:

011: OpenFOAM is free software; you can redistribute it and/or modify it

012: under the terms of the GNU General Public License as published by the

013: Free Software Foundation; either version 2 of the License, or (at your

014: option) any later version.

015:

016: OpenFOAM is distributed in the hope that it will be useful, but WITHOUT

017: ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

018: FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License

019: for more details.

020:

021: You should have received a copy of the GNU General Public License

022: along with OpenFOAM; if not, write to the Free Software Foundation,

023: Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

024:

025: Application

026: interPhaseChangeCondenseTempFoam

027:

028: Description

029: Solver for 2 incompressible, immiscible fluids with condensation

030: phase-change. Uses a VOF (volume of fluid) phase-fraction based interface

031: capturing approach. It also incorporates the solution of the energy

032: equation.

033:

034: Based on interPhaseChangeCondenseFoam and the phase change models associated

035: with this solver.

036:

037: The momentum and other fluid properties are of the ”mixture” and a

038: single momentum equation is solved. The energy equation is only solved for

039: the water phase.

040:

041: Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.

(46)

32 Main Program Source Code 042: 043: \*--–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–--*/ 044: 045: #include"fvCFD.H" 046: #include"MULES.H" 047: #include"subCycle.H" 048: #include"interfaceProperties.H" 049: #include"phaseChangeTwoPhaseMixture.H" 050: #include"turbulenceModel.H" 051: 052: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 053:

054: intmain(intargc, char*argv[])

055: { 056: #include "setRootCase.H" 057: #include "createTime.H" 058: #include "createMesh.H" 059: #include "readGravitationalAcceleration.H" 060: #include "readPISOControls.H" 061: #include "initContinuityErrs.H" 062: #include "createFields.H" 063: #include "readTimeControls.H" 064: #include "correctPhi.H" 065: #include "CourantNo.H" 066: #include "setInitialDeltaT.H" 067: 068: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 069:

070: Info<<"\nStarting time loop\n"<<endl;

071:

072: while(runTime.run())

073: { 074: #include"readPISOControls.H" 075: #include"readTimeControls.H" 076: #include"CourantNo.H" 077: #include"setDeltaT.H" 078: 079: runTime++; 080:

081: Info<<"Time = "<<runTime.timeName() <<nl <<endl;

082: 083: #include"alphaEqnSubCycle.H" 084: 085: turbulence->correct(); 086: 087: // --- Outer-corrector loop

088: for(intoCorr=0;oCorr<nOuterCorr;oCorr++)

089: {

090: #include"UEqn.H"

091:

092: // --- PISO loop

093: for(intcorr=0;corr<nCorr;corr++)

(47)

33 102: { 103: #include"TEqn.H" 104: }; 105: 106: twoPhaseProperties->correct(); 107: 108: runTime.write(); 109:

110: Info<<"ExecutionTime = "<< runTime.elapsedCpuTime() <<" s"

111: <<" ClockTime = "<<runTime.elapsedClockTime() <<" s"

112: <<nl <<endl;

113: }

114:

115: Info<<"End\n"<<endl;

(48)

(49)

APPENDIX

B

SOURCE CODE COMBUSTION MODEL (ALBANOVACOMBUSTION)

001: /*--–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–--*\

002: ========= |

003: \\ / F ield | OpenFOAM: The Open Source CFD Toolbox

004: \\ / O peration |

006: \\/ M anipulation |

007: --–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–

008: License

009: This file is part of OpenFOAM.

010:

011: OpenFOAM is free software; you can redistribute it and/or modify it

012: under the terms of the GNU General Public License as published by the

013: Free Software Foundation; either version 2 of the License, or (at your

014: option) any later version.

015:

016: OpenFOAM is distributed in the hope that it will be useful, but WITHOUT

017: ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

018: FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License

019: for more details.

020:

021: You should have received a copy of the GNU General Public License

022: along with OpenFOAM; if not, write to the Free Software Foundation,

023: Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

024: 025: \*--–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–-–--*/ 026: 027: #include"AlbaNovaCombustion.H" 028: #include"addToRunTimeSelectionTable.H" 029:

030: // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

031: 032: namespace Foam 033: { 034: namespace phaseChangeTwoPhaseMixtures 035: { 036: defineTypeNameAndDebug(AlbaNovaCombustion,0);

037: addToRunTimeSelectionTable(phaseChangeTwoPhaseMixture,AlbaNovaCombustion,

038: components);

039: }

040: }

041: