CFD Annular Flow Modelling Based on a Three-Field Approach

(1)

CFD Annular Flow Modelling Based on a Three-Field Approach

Erik Skoog

Engineering Physics and Electrical Engineering, master's level 2020

Luleå University of Technology

Department of Engineering Sciences and Mathematics

(2)

Abstract

This master thesis aim to model the annular flow that occurs in the final section between the fuel rods inside Boiling Water Reactors, by ap- proximating the geometry to a cylindrical pipe. Simulations were performed in the software ANSYS Fluent, as a step in the development of replacing the 1D correlations currently used in the nuclear industry with CFD models in 3D. An Eulerian-Lagrangian approach was used for the three fields of steam, liquid film and liquid droplets in the model. En- trainment was modeled based on 1D correlations from Okawa [7] and deposition with the built in Discrete Phase Model in ANSYS Fluent. The work focused on making the process less time consuming, and increasing accuracy of the model by comparing the results with empirical data based on experimental values. A transverse velocity was applied on the droplets at the point of entrainment with better correlating results with the Okawa model.

(3)

Acknowledgements

First, a big thank you to my primary supervisor Tobias Strömgren. His knowledge on CFD has been extremely valuable, and has always been there to answer my questions. Also, thank you to Jean-Marie Le Corre, who has offered his great expertise on the area and provided for the main guidance on the general direction of this work. Thank you to Carl Adamsson for his constant knowledge sharing on the area and letting me use his experimental results for compari- son. Thank you to Yann Le Moigne for observing my work and providing with feedback and technical support throughout the work. Thank you to Anna-Lena Ljung, my examinator at Luleå University of Technology.

Thank you to Westinghouse Electric Sweden in Västerås for allowing me to perform my master thesis for the company. Finally, a big thank you to the SEP department at the workplace for welcoming me into the team as one of their own. It has been a great pleasure working alongside with you and getting to know you all.

(4)

List of Figures

1 Image of different flow regimes [11] . . . 8 2 Sketch of the setup of the model. Note that the image is not to

scale. . . 14 3 Cross sectional mesh of the geometry, containing 1300 cells. . . . 15 4 Injection wall . . . 16 5 Deposition rate when injecting particles in two different ways.

G = 1250kg/m²s . . . 17 6 1D calculation from onset of annular flow using the equilibrium

condition. . . 20 7 Total mass flow rate G = 750kg/m²s, compared with experimen-

tal results from the same conditions. . . 20 8 Total mass flow rate G = 1250kg/m²s, compared with experi-

mental results from the same conditions. . . 21 9 Total mass flow rate G = 1750kg/m²s, compared with experi-

mental results from the same conditions. . . 21 10 Film thickness across the perimeter of the pipe. Starting at

z=2.54 m, and slowly decreasing as the film continues further downstream. . . 22 11 Droplet mass flow for different droplet fractions of total amount

of liquid applied at the onset of annular flow. G=1250 kg/m²s and steam quality x = 51% at outlet. . . 22 12 Mass flow rates for the three fields, G = 750kg/m²s and steam

quality at outlet x = 75% . . . 23 13 Mass flow rates for the three fields, G = 1250kg/m²s and steam

quality at outlet x = 51% . . . 24 14 Mass flow rates for the three fields, G = 1750kg/m²s and steam

quality at outlet x = 43% . . . 24 15 Deposition rate divided by concentration of droplets, where the

concentration is determined as in equation (12) with the values for the droplet and steam mass flows obtained from Fluent. The solid lines are 1D calculations using the Okawa correlations and correspond color wise to the same inlet mass flow as labeled in the figure. The step wised dotted curves are values obtained from simulations in Fluent . . . 25

(6)

Nomenclature

List of Symbols

δ Liquid film thickness, m

˙

m Mass flow rate, kg/s µ Dynamic viscosity, Pa·s ρ Density, kg/m³

σ Surface tension, N/m τ Shear stress, N/m³ D Diameter, m

f Friction factor, dim-less G Total inlet mass flux, kg/m²s g Gravitational acceleration, m/s² J Volumetric flux, m/s

k Mass transfer coefficient, dim-less p Pressure, Pa

r Radius, m

S_m Mass source term, kg/m²s Smom Momentum source term, N/m² T Temperature, K

u Velocity, m/s Subscripts

d Liquid droplets dep Deposition ent Entrainment f Liquid film i Interface

(7)

k Liquid or vapor

l Liquid

v Vapor

w Wall

List of Abbreviations

BWR - Boiling Water Reactor

CFD - Computational Fluid Dynamics DPM - Discrete Phase Model

GUI - Graphical User Interface LPT - Lagrangian Particle Tracking TUI - Text User Interface

UDF - User Defined Functions 1D - One Dimension

2D - Two Dimensions 3D - Three Dimensions

(8)

1 Introduction

In Boiling Water Reactors (BWRs), liquid water is heated until it vaporizes and become steam, which in turn drives turbines to generate electricity. Heat is transferred from tall cylindrical fuel rods to the water flowing between the rods. When this heating occurs, the water passes through various flow regimes resulting in a phenomenon called annular flow at the upper part of the core.

This is the final flow regime before all liquid has been vaporized due to the heat flux from the rods. If all liquid film surrounding the rods vaporize before reaching the upper end of the core, dryout is inevitable. This must be avoided, since dryout can cause severe damage on the equipment due to the less efficient heat exchange occurring between the rods and the steam when the liquid film vanishes. Thus it is essential for companies in the nuclear industry to have well- functioning and accurate models that can predict the critical heat flux when dryout occurs.

Researchers Haipeng Li and Henryk Anglart at KTH in Stockholm has in close cooperation with Westinghouse Electric Sweden developed a computational fluid dynamics (CFD) model for this purpose using the open source software OpenFOAM. Following that, the previous master thesis student Maria Camacho worked on improving the model and implemented it in the software ANSYS Fluent [4]. The following year, the master thesis student Salvatore Raddino continued the work of Camacho and made several improvements to the model; increasing accuracy and decreasing computational time [8]. This master thesis is a continuation of the work previously performed by Raddino and Camacho on the modelling of annular flow using CFD based on a three-field approach.

During the CFD-modeling of annular flow in cylindrical pipes, three main areas are considered; the steam core occupying the center of the pipe, liquid film surrounding the pipe along the walls, and liquid droplets within the steam core. These three fields have three main ways of interacting with each other.

Entrainment is the rate at which droplets are formed and released from the film into the steam core, deposition is the rate at which liquid droplets are depositing on the liquid film and vaporization from the liquid film to the steam core. Vaporization is dependent on the power distribution along the rods. In this work, only uniform power distribution is considered.

Most of the current studies on annular flow are based on empirical correlations in one dimension along the pipe. The disadvantage in relying on these correlations mainly consist of the severe restrictions in the use of different parameters. Being able to perform meticulous CFD-simulations in 3D will allow for a more detailed understanding of all complex events occurring in the flow.

In 2003, Okawa et al. developed 1D-correlations for the rates of entrainment and deposition in perfectly cylindrical shaped pipe geometries[7]. The CFD- model in this thesis will make use of the entrainment correlation from Okawa, whereas deposition will be determined by calculating particle trajectories with a Lagrangian approach, using the Discrete Phase Model (DPM) provided in ANSYS Fluent, the software that will be used for the simulations. The purpose

(9)

with this is to obtain an understanding for how sensitive the model is to this way of modeling deposition, and results will be compared with 1D-calculations based on the Okawa correlations and experimental results.

The main objective of this thesis will be to study the sensitivity of the model to different inlet mass flows with an implemented droplet transverse velocity at the point of entrainment. At the onset of annular flow, different fractions between droplets and liquid film will be considered. Droplet- and film mass flow will then be simulated based on these different fractions at the onset of annular flow. These simulation results will then be compared with calculations based on the Okawa correlations.

Another objective will be to study the injection of droplets and how a realistic particle injection affects deposition.

2 Background

2.1 Boiling Regimes

There are a few flow regimes for the fluid to pass through before reaching annular flow and eventually completely evaporate. Bubbly flow, slug flow, churn flow and finally annular flow, all seen in figure 1. Bubbly flow occurs as soon as small vapor bubbles start to appear in the flow. As these bubbles increase in size and more liquid is vaporized, large slugs start to appear. As the fraction of steam increases, so does velocity, since pressure is kept constant in the reactor core. The slugs will later be torn apart into smaller pockets of steam as velocity is constantly increasing, causing the transition to churn flow. With further increasing velocity, the forces acting upon the pockets filled with vapor will become larger and make them unstable.

Figure 1: Image of different flow regimes [11]

Finally, annular flow with a liquid film circumfering the walls of the pipe and liquid droplets inside the steam core is formed. All simulations performed in this master thesis will only consider annular flow.

(10)

The physics behind the transitions between these flow regimes is extremely complex. Based on experimental data, G. B. Wallis provides the following correlation as a condition where the transition between churn and annular flow occurs [10],

J_v^∗= 0.4 + 0.6J_l^∗. (1)

The dimensionless volumetric fluxes for the vapor and liquid J_k^∗are each defined as

J_k^∗= Jk

r ρk

gD (ρ_l− ρ_v). (2)

Here, g is the gravitational constant, D the pipe diameter and ρ the density.

The subscript k denotes the phase in question, liquid or vapor. This method to determine the onset of annular flow was used by Okawa during the development of the 1D-correlations of deposition and entrainment [7]. It is therefore of utmost importance to use the same method in this work, since those 1D-correlations will be studied.

2.2 Properties of Annular Flow

The majority of the volume in the pipe is occupied by a steam core located in the center. The steam core mass flow is continually increasing due to the mass transfer occurring through evaporation of the liquid film. The rate of evaporation is determined by the heat flux from the wall and its distribution.

The liquid phase is divided into droplets and film. Liquid droplets inside the steam core are heavily affected by the turbulence forces in the vapor. This is particularly true for smaller droplets, since larger droplets carry more mass and thus possess a larger momentum. Droplets inside the steam core depositing on the liquid film surface are absorbed by the film, thus adding mass to it. The rate of deposition is directly proportional to the concentration of droplets in the steam core.

On the vapor-liquid interface, disturbance waves are nearly always present.

According to P. B. Whalley, droplets entrain from the film if the liquid film waves reach a certain height [12]. Experiments on this phenomenon where photos have been taken cross section wise, show that for large enough Reynolds numbers, the wave amplitudes are high enough for the liquid film to shatter into the steam core, quickly splitting up into several smaller droplets. This is called entrainment, and is directly proportional to the film thickness. The Reynolds number is in this case determined by

Ref =ρ_fJ_fD µl

(3) where µlthe liquid dynamic viscosity.

In this thesis, the mass transferring through vaporization, deposition and entrainment are considered and modeled.

(11)

3 Theory

3.1 Equations in Fluent

3.1.1 Steam Core

The steam is treated as a single-phase flow in Fluent and a two-way coupling is imposed on the steam and droplets. The continuity equation (4) is solved, containing a time derivative and a divergence term, describing the change of mass in time and space,

∂ρ_v

∂t + ∇ · (ρ_vu_v) = S_vap(z). (4) The source term on the right-hand side corresponds to the mass transfer that is constantly occurring from the film to the vapor through vaporization. Here, z denotes the axial direction of the pipe.

The source term is defined as the fraction between the power heat flux en- tering the pipe from the wall and the latent heat of vaporization described by the difference between the saturated enthalpies of vapor (hv) and liquid (hl),

Svap(z) = q(z) hv− hl

. (5)

The power heat flux q(z), can in theory possess any distribution along the pipe.

Fluent also solves the momentum conservation equation for the vapor phase, seen in equation (6),

∂(ρ_vu_v)

∂t + ∇ · (ρ_vu_vu_v) = −∇p_v+ ρ_vg + ∇ · τ_v^{ef f}+ S_mom,vap+ S_mom,drop. (6) Here, acceleration through the time derivative and divergence term is described on the left-hand side. On the right-hand side are pressure, gravity and stress forces along with two momentum source terms from vaporization of the film (Smom,vap) and droplet flow in the steam core (Smom,drop).

3.1.2 Liquid Droplets

Calculations on the particles will be performed in Fluent using the built-in Discrete Phase Model (DPM). Instead of performing calculations on each injected particle, Fluent treats the particles as part of parcels using an Eulerian- Lagrangian approach. A parcel can contain several particles with similar properties and makes the calculations in Fluent faster. Several parcel release methods can be chosen in Fluent. When injecting droplets with a User Defined Function (UDF), the parcel release method is automatically set to standard and cannot be changed. The standard release method injects a single parcel per injection stream per time step. The number of particles per parcel is determined as

N P = ˙ms

∆t md

(7)

(12)

where NP is the number of particles in a parcel, ˙msis the mass flow rate of each particle stream, ∆t is the time step and md is the particle mass. m˙s can be calculated by taking the total particle mass flow divided by the number of facets in the cross sectional mesh. With a ∆t of 0.0005 seconds, a particle diameter of 1.2mm, and a total mass flow rate of 1250kg/m²s, the number of particles per parcel is 0.0585. Since every facet in the cross sectional mesh is subject to one injection stream and thus one parcel at every time step, this number can be below one. The simulation needs to be running over a sufficient amount of time steps for enough parcels to be injected and hence create a statistically accurate result for the particles. The particle trajectories are calculated and tracked until they hit the wall and are absorbed by the film, or reach the outlet.

The calculations between the vapor and the droplets are two-way coupled;

the change in momentum between these two phases goes both ways. The momentum transfer in this model is computed as

S_mom,drop= 18µC_DRe

ρdd²_d24 (u_d− uv) + F_other

˙

m_d∆t (8)

where F_other accounts for gravity and drag forces [3] and is calculated as Fother= g

1 −ρ_v

ρd

. (9)

3.1.3 Liquid Film

For the modeling of the liquid film, the built-in function Eulerian Wall Films (EWF) is used. The function makes use of the assumption that the liquid film is thin enough for certain properties to remain constant in the film, independent of radial position. Fluent solves the conservation of film momentum,

∂ρfδuf

∂t + ∇_s· (ρ_fδu_fu_f) = −δ∇_sp_f+ ρ_fδg +3

2τ_vf−3µf

δ µ_f+

Sm,dep− Sm,ent− Sm,vap (10) which has certain similarities with the single-phase momentum equation 6. The film thickness is δ, and the third and fourth term on the right hand side of equation (10) correspond to the net viscous shear force in the vapor-film and wall-film interfaces [3]. Using the EWF, the liquid film is calculated in 2D along the wall and hence does not occupy any space in the 3D model.

3.2 1D-correlations

It is the calculation in 1D that will provide for almost all boundary conditions necessary for the simulations. One of these is the fraction of liquid film, droplets and steam at the inlet. The total mass flow rate can then be multiplied with each fraction respectively to obtain the mass flow rate for each field. The fraction between vapor and liquid at the onset of annular flow can be calculated

(13)

based on the correlations provided in equations (1) and (2). To obtain this fraction, Utsuno and Kaminaga suggest an equilibrium condition of deposition and entrainment [9]. The fraction is step wised increased from 0 to 1 and entrainment and deposition are calculated for every step. The fraction between liquid droplets and liquid film at the onset of annular flow is then found where entrainment and deposition are equal.

3.2.1 Deposition

The simulations from Fluent will be compared with 1D-calculations based on Okawa’s correlations of entrainment and deposition. For deposition, the correlation used in the 1D-calculation was presented by Okawa et al. [7], based on empirical studies:

Sm,dep= 0.0632 rCσ

D (11)

where C is the concentration of droplets in the steam core. With the assumption that the relative velocity between the droplets and the steam core is negligible, the droplet concentration becomes

C = m˙d

( ˙md/ρl) + ( ˙mv/ρv) (12) as described by Le Corre and Adamsson [6].

3.2.2 Entrainment

Entrainment is calculated as described by Okawa et al. [7], based on empirical studies. This takes into consideration that entrainment only occurs if the Reynolds number of the liquid film is above a certain critical value Ref c. The correlation provided by Okawa can be seen in equation (13) and is used in the conservation of mass calculations of the liquid film in equation (10).

S_m,ent= k_eρ_lf_iρ_vJ_v²δ σ

ρ_l ρv

ⁿ

f or Re_f > Re_{f c} (13) where the subscripts l and f denotes the liquid phase and the liquid film, respectively. ke and n are constants with recommended values of 4.79 · 10⁻⁴m/s and 0.111 respectively. The superficial vapor velocity is denoted Jv and the critical film Reynolds number is set to Ref c= 320. The interfacial friction factor fi is seen in equation (14) and is dependent on the film thickness, which is obtained directly from Fluent,

fi= 0.005

1 + 300δ D

. (14)

In the 1D correlations the film thickness is estimated from the force balance between the interfacial shear force and the wall friction force, both acting on the liquid film. It is described as

(14)

δ = 1 4

s f_wρ_l fiρv

J_f Jv

D (15)

where the wall friction factor f_wis recommended by Wallis [10] to be calculated as

f_w= max

16 Ref

, 0.005

. (16)

Equations 13 and 14 are implemented in the UDF to control the film mass source. All constants are defined in advance and the Reynolds number and film thickness are calculated by Fluent.

4 Method

4.1 Overview of previous work

This master thesis is largely based on the two previous thesis works performed by Maria Camacho and Salvatore Raddino. Camacho developed the first model in Fluent and compared it with the results obtained from an OpenFOAM model developed at KTH. Raddino continued the work of Camacho and implemented general improvements to the model, with more accurate results and vastly shorter calculation time. The model developed by Raddino contains two simulation regions. One inlet development model and one model for the main simulation.

The reason for this setup is to provide enough time for the fluid to fully develop into annular flow. Each of these models consist of two regions. First, the 1.06 meter long annular_wall, where simulation data is obtained and compared with theoretical results. This corresponds to the last section of the in total 3.65 meter tall fuel bundles. Next, a shorter injection_wall (0.1103 m) just before the previously mentioned region. From here the liquid film is injected radially and a vapor velocity is defined at the inlet. A sketch of where these two regions are found in the setup can be seen in figure 2.

(15)

Figure 2: Sketch of the setup of the model. Note that the image is not to scale.

In both previous master theses, the modeled deposition rate can be observed as either a substantial underestimation (Camacho) or an intense peak (Raddino) at the very first part of the simulated pipe. This is the result of an entrance effect, created by the defined inlet conditions for the droplets. Steam and droplets are injected at the inlet of the first model with a uniform velocity profile. This creates indeed an entrance effect, although this is irrelevant since only the outlet of this simulation is of interest. For the steam, velocity properties are saved at the outlet of the first model as a profile file and uploaded as initial condition for the steam for the main model. The particles are treated in the same way, where the properties of the particles are saved in an injection file. This file is then defining the initial conditions for the particles in the main model, where droplets are injected at the inlet of the annular_wall region, using the DPM.

One of the changes made by Raddino regarding the particles was the implementation of injection files. Instead of injecting the particles uniformly over the inlet surface of the main model, particles are injected through an injection

(16)

file containing all necessary information of the particles. Particles are first injected uniformly through the inlet surface of the InletDevelopment model and develop a more realistic distribution at the outlet during that distance. The positions and velocities of the particles are then saved in an injection file, which is implemented at the inlet of the main model. Seen in figure 3 is the cross sectional mesh used in this work. This mesh has not undergone any changes since the first OpenFOAM model constructed by Li and Anglart which, according to Camacho, also managed to confirm grid independence for this mesh. In order to better capture the small local changes that might arise in the flow closer to the wall, the mesh is more refined in this area, containing a higher number of elements in the radial direction closer to the wall.

Figure 3: Cross sectional mesh of the geometry, containing 1300 cells.

4.2 Changes and Improvements

4.2.1 Liquid Film

An overall script was constructed, with the purpose of being well synchronized with the models in Fluent. The script, provided in Appendix A, would perform 1D-calculations on various parameters of the flow, such as entrainment, deposition and film thickness based on correlations constructed by Okawa et al [7]. The results from these 1D-calculations would then be used as initial and boundary conditions for the model in ANSYS Fluent. The necessary conditions involving the three fields are flow fluxes of droplets and film, and vapor velocity.

Conditions regarding the surrounding environment include pressure, temperature, heat flux and power distribution. It is based on these latter parameters which more specifically determines the environment, constants such as density and viscosity of the two phases are defined.

The results from the simulation would then be compared with the calculations from the 1D-model, and hopefully correlate rather well. An issue with the

(17)

model was the mismatch in the manually calculated input value for the liquid film flow and the simulated film flow right at the start of annular_wall.

Figure 4: Injection wall

This problem was solved by splitting the first region seen in figure 4 into two parts (A and B) and inject the liquid film only through part B. More liquid film would then need to be injected, since the injection area would shrink with 1/9.

One of the more obscure boundary conditions is the "z-momentum flux". A force that is applied in the injection_wall, setting the liquid film in motion. In this work, it was determined through trial and error by comparing the simulated film thickness and the theoretical film thickness at the inlet. A larger force would result in a higher initial film velocity, decreasing the film thickness. A smaller force would have the opposite effect, causing an increase of the initial film thickness. For G = 1250kg/m²s, it was approximated to 55.7N/m². 4.2.2 Automation and Initialization

There are two ways of using ANSYS Fluent; by using the graphical user interface (GUI) or the text user interface (TUI) in the command window. One major progress from the previous work by Raddino, has been the automation of the simulations with the use of journal files. A journal file is a file containing instructions in text format, which are executed by ANSYS Fluent line by line.

By implementing all conditions needed for a specific simulation into a journal file, Fluent can be ran entirely without graphics, only making use of the TUI.

This modification makes it much less time consuming to run several simulations with different boundary conditions by completely avoiding using the GUI. The computational time required to perform one simulation varies greatly depending on the implemented boundary conditions, but for the three mass flow cases (G = 750kg/m²s, G = 1250kg/m²s and G = 1750kg/m²s) the computational time decreased from approximately five hours to just below one hour.

Based on studies of several simulations on the matter, it is estimated that the flow is fully developed at 0.325 seconds. That is, after 650 time steps. With this in mind, it is only desirable to export solution data after 650 time steps.

(18)

From there, solution data is exported in ASCII-files at every fifth time step until 1000 time steps (0.5 seconds) are reached. This gives 70 ASCII-files to average over, and is an effective way of reducing noise from specific time steps. The files are saved in the same folder as a Matlab script used for reading the files and plotting the result.

Along with this process, an improvement made for the particle injection was the removal of the need of reformatting the file using a separate C script.

4.2.3 Steam and Particles

The first InletDevelopment model was seen as a redundant element of the simulation process. With the existing setup, in order to run a simulation with certain conditions, the first InletDevelopment model had to be initialized and completed first. For its removal, the steam and particles had to be injected in a sufficiently realistic way, without relying on the first model to generate reasonable particle and velocity profiles.

For the particles, two types of injections are implemented. One at the inlet using a particle injection file, and particles injected into the core from the wall through entrainment using a surface injection. For the injection at the inlet, various injections were tried and compared with the results from the Okawa 1D-model. Keeping one set of original particle velocities (first obtained by running the InletDevelopment model), and multiplying their transverse velocity components in the XY-plane with 0.55 significantly reduced the first peak in deposition as seen in figure 5.

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0 0.2 0.4 0.6 0.8 1 1.2 1.4

Mass flux [kg/m2 s]

Deposition

New injection Old injection 1D Okawa model

Figure 5: Deposition rate when injecting particles in two different ways. G = 1250kg/m²s

(19)

As for their axial velocities, keeping an original distribution but varying the mean axial velocity proportional to 95% of the injected steam velocity provided the most accurate results. The droplet mass flow was calculated manually and multiplied with each cell area in the injection surface to create a uniform injection of droplet mass flow. This also proved to be a more accurate method, rather than applying a droplet mass flow with a certain distribution. Looking at figure 5, it is interesting to notice the slope of the curve from the new injection at z=3.1. Although, it is slightly overpredicting the Okawa correlation, it seems to decrease at a rather similar rate as the correlation under these conditions.

In the previous model, newly formed droplets detached from the film were given the same velocities in all directions as the steam in that specific position of entrainment. A more realistic way of modeling this would involve assigning these droplets a transverse velocity at the point of entrainment. P. Andreussi and B. J. Azzopardi provide a correlation for this purpose [2], that has been used by Haipeng Li in previous research. The particle transverse velocity correlates with

ud= 11.1uτ

r ρv

ρl

(17) where

u_τ= s

f f · u²_{f ree}

2 (18)

f f is a friction factor determined by

f f = 0.045 u_{f ree}δ ν

−1/4

(19) and

u_{f ree}= u_v8

7 (20)

where uv is the mean vapor velocity at the current the cross section. The transverse velocity would be calculated and applied in the UDF with X- and Y-components, assuming a radial direction towards the center of the pipe of this velocity. In this work, only 60% of this calculated transverse particle velocity u_d was applied, due to convergence issues for larger velocity magnitudes. All UDFs used in this work are provided in Appendix B.

In the previous model, steam was injected as a velocity profile obtained at the outlet of the InletDevelopment model. This made the profile fixed, impos- sible to adjust the velocity magnitude for different conditions. Injecting steam uniformly, only possessing a velocity in the axial direction proved to not change the final results noticeably. The steam would be injected at the first inlet of the pipe, passing through the injection_wall region.

(20)

4.3 Post-Processing

Much of the focus in this work has been put on making the results obtained through post-processing more reliable. Deposition rate is now obtained directly from Fluent and not manually calculated using a mass balance equation as in the previous work. During the simulation, ASCII-files containing the desired data are saved at every fifth time step. The data is then summed and averaged over the time steps where the flow is fully developed.

Since the liquid film is calculated in 2D on the walls, the entire cross sectional area can be assumed to be filled with droplets and vapor. By creating cross sectional planes through the pipe, the vapor velocities at these locations can be used to calculate the vapor flow rate. This is done by multiplying the mean vapor velocity with the vapor density and the cross sectional area. The droplet mass flow rate is calculated as a mass balance, by subtracting the vapor and film flow rates from the total injected mass flow rate.

4.4 Assumptions

Assumptions have been made based on assumptions documented in previous work regarding annular flow and the accuracy of the model used for this study.

No vaporization is assumed to occur of the liquid droplets to the steam core.

That is, all droplets are tracked by the Lagrangian Particle Tracker (LPT) until they deposit on the film or reach the outlet. All droplets are also assumed to be perfectly spherical, and possess a diameter of 1.2 mm. In the previous master thesis work, Raddino performed some studies on the droplet diameter and 1.2 mm was considered to be rather optimal. Studies on droplet diameter have also been performed by Le Corre et al., who concluded that particles in annular flow regimes like these with a diameter of 1.2 mm contain the majority of the total droplet mass [5]. For the simulation, droplets are also assumed to not merge with each other and form bigger droplets, as well as not splitting up into smaller droplets. This phenomenon is ignored because of the current absence of such a calculation model within the field of CFD.

Saturation conditions are assumed to be present throughout the entire flow region, and all heat flux on the wall is used for the evaporation of the liquid film. All simulations are performed under the condition of 10 K subcooling and a pressure of 70 bar. Based on these assumptions, a uniform pressure and temperature can be considered across the entire pipe.

5 Results

Before initiating a simulation the 1D-script containing the correlations by Okawa performs calculations on the specific case. First, the position where the onset of annular flow occurs is determined. Based on that result, the boundary conditions can be obtained for the simulation model. Figures 6a and 6b show the results from the script using a total flow rate of G = 1250kg/m²s with a steam quality of x = 51% at the outlet and uniform power distribution. Under these

(21)

conditions, the onset of annular flow appear at z = 1.13 m and the three fields along with the rates of deposition and entrainment can be calculated and plotted onward from that point.

1 1.5 2 2.5 3 3.5 4

z [m]

0 0.02 0.04 0.06 0.08 0.1 0.12

Mass flow rates [kg/s]

Mass flow rates

Liquid film Vapor Droplets

(a) Mass flow rates for the three fields

1 1.5 2 2.5 3 3.5 4

z [m]

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Mass flux [kg/m2s]

Deposition and Entrainment rates Deposition Entrainment

(b) Deposition and entrainment rates

Figure 6: 1D calculation from onset of annular flow using the equilibrium condition.

Experimental results in figures 7-9 have been taken from Adamsson and Anglart [1]. For all three cases, the deposition rate recover from the under- prediction at z = 3.1 and transitions into an overprediction in figures 8 and 9.

However disregarding that, the slope of the deposition rate for all cases coincide fairly well with the 1D Okawa model after about 3.1 m. This pattern is observable for the film flow rates from Fluent as well, most notably in figures 8b and 9b where flow rate initially decreases but start to recover at z = 3.1.

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0.05 0.1 0.15 0.2 0.25 0.3 0.35

Mass flux [kg/m2s]

Deposition, G=750 kg/m²s From Fluent 1D Okawa model

(a) Deposition rate

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0 0.005 0.01 0.015 0.02

Mass flow [kg/s]

Film Flow Rate, G=750 kg/m²s From Fluent 1D Okawa model Experimental results

(b) Film flow rate

Figure 7: Total mass flow rate G = 750kg/m²s, compared with experimental results from the same conditions.

(22)

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 z [m]

0.1 0.2 0.3 0.4 0.5 0.6 0.7

Mass flux [kg/m2s]

(a) Deposition rate

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0 0.005 0.01 0.015 0.02 0.025 0.03

Mass flow [kg/s]

(b) Film flow rate

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0.2 0.3 0.4 0.5 0.6 0.7 0.8

Mass flux [kg/m2s]

(a) Deposition rate

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0 0.005 0.01 0.015 0.02 0.025

Mass flow [kg/s]

(b) Film flow rate

Figure 10 shows the film thickness along the wall of the pipe, averaged over several time steps where the flow is fully developed. Small variations are expected because of the constant deposition of droplets causing changes in the film thickness on a local level. That is probably the explanation for the sudden peak observable at approximately 25 degrees. Droplets are injected just at z=2.54 but do not collide with the wall instantly. Hence, variations can be seen 6.1 cm further downstream when the first set of droplets has deposited and not on the blue curve where the droplets are first injected.

(23)

-200 -150 -100 -50 0 50 100 150 200 degrees

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

Film thickness [m]

10^-4 Film thickness across perimeter, 6.1 cm spaced

Figure 10: Film thickness across the perimeter of the pipe. Starting at z=2.54 m, and slowly decreasing as the film continues further downstream.

Figure 11 shows the development of droplet mass flow for various conditions applied at onset of annular flow. Instead of imposing an equilibrium condition between entrainment and deposition, the fraction between liquid film and liquid droplets is controlled manually, while the onset of annular flow is still assumed to occur at the same position (z=1.13). When the equilibrium condition is active, the fraction of droplets is 60.3%. In this case, the power is uniformly distributed, G=1250 kg/m²s and steam quality x = 51% at outlet.

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0.085 0.09 0.095 0.1 0.105 0.11 0.115 0.12

Mass flow [kg/s]

Mass flow rate

30% droplets 40% droplets 50% droplets 60% droplets 70% droplets 80% droplets 90% droplets

(a) Simulated results

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0.085 0.09 0.095 0.1 0.105 0.11 0.115 0.12

Mass flow [kg/s]

1D Okawa model

30% droplets 40% droplets 50% droplets 60% droplets 70% droplets 80% droplets 90% droplets

(b) 1D-Okawa model

Figure 11: Droplet mass flow for different droplet fractions of total amount of liquid applied at the onset of annular flow. G=1250 kg/m²s and steam quality x = 51% at outlet.

(24)

Figures 12 to 14 show the mass flow for all three fields present in the pipe for different total inlet mass flows compared with calculations in 1D based on the Okawa correlations. The field of each 1D curve is corresponding to the curve from Fluent it is closest to. The vapor seems to coincide with the Okawa correlation fairly well for all cases, but for the liquid film and droplets slight differences are noticed, specifically in figure 12. The droplet mass flow is calculated as a mass balance from the vapor and film flow rates, so it is not unexpected with a larger difference in droplet flow rate when such is noticed in the film flow rate, as in figure 12.

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

Mass flow rate [kg/s]

Mass flow rates, G=750 kg/m²s

Droplets from Fluent Liquid film from Fluent Vapor from Fluent 1D Okawa

Figure 12: Mass flow rates for the three fields, G = 750kg/m²s and steam quality at outlet x = 75%

(25)

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 z [m]

0 0.02 0.04 0.06 0.08 0.1 0.12

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

Since deposition is strictly dependent on the concentration of droplets, di- viding the deposition rate with the concentration of droplets can provide some

(26)

indications for the deposition sensitivity. Seen in figure 15 is the fraction between deposition rate and concentration of droplets. The cases where G = 750kg/m²s and G = 1250kg/m²s follow their corresponding 1D Okawa correlation rather well. For the case where G = 750kg/m²s, the liquid film has completely vaporized at z = 3.4 as can be seen in figure 7. Since the CFD-model is only valid for as long as the liquid film is still present, the results beyond this point are not trustworthy and have thus been removed. The case where G = 1750kg/m²s is clearly below its corresponding 1D correlation, indicating that the model has difficulties in performing calculations on the droplets for this inlet mass flow.

As can be seen in figure 9a, deposition differed the most from the Okawa correlation for this inlet mass flux, compared to the other cases in figures 7a and 8a.

2.4 2.6 2.8 3 3.2 3.4 3.6 3.8

z [m]

0 0.005 0.01 0.015 0.02 0.025

[m/s]

Deposition over droplet concentration

G=750 kg/m²s G=1250 kg/m²s G=1750 kg/m²s

Figure 15: Deposition rate divided by concentration of droplets, where the concentration is determined as in equation (12) with the values for the droplet and steam mass flows obtained from Fluent. The solid lines are 1D calculations using the Okawa correlations and correspond color wise to the same inlet mass flow as labeled in the figure. The step wised dotted curves are values obtained from simulations in Fluent

(27)

6 Conclusions

It is clear that the LPT in the CFD-model is highly sensitive to changes in properties of the particles. Many different properties of the injected particles were tried and rate of deposition was always highly unpredictable. In order to obtain an initial velocity distribution of the droplets, particles were first injected through a surface injection with uniform velocity in axial direction at the inlet.

The positions and velocities of all particles in every facet at the outlet were then saved. This way, the particles would already have an initial realistic velocity distribution axially as well as radially.

Keeping the velocity distribution in axial direction, although multiplying the magnitude with 90% of the inlet steam velocity provided satisfying results for deposition rates. The particles would also have an initial transverse velocity in the (cross sectional) XY-plane. Keeping this velocity unchanged resulted in a large initial spike in deposition. This can be explained by droplets having a large transverse velocity depositing on the walls right after the droplets are injected.

With a reduced initial transverse velocity, this spike could be reduced. The initial transverse velocity was set to 55% of the transverse velocity the particles originally obtained at the outlet of the InletDevelopment model.

The implementation of a transverse velocity for entrained droplets constructed by [2] proved to generate results more accurate with the Okawa correlations for the deposition rate, despite only using a 60% fraction of the suggested velocity. Even though it could not be noticed for the first part of the simulation, the deposition rate would be much closer to the Okawa correlation during the second half of the geometry. Figure 11 shows an interesting pattern, where the simulation results seem to follow a somewhat similar shape as the calculated results. Seen in figure 11a, the different fractions seem to merge together further downstream, just as in the 1D-calculation in figure 11b. This is yet another indication that the model treats the droplet mass flow moderately accurate.

This master thesis has achieved in making the simulation process and post- process more efficient and less time consuming. An overall 1D script with attached functions was constructed based on the Okawa correlations for calculating the boundary conditions and compare the simulation results with. The possibility to automatically initiate simulations in Fluent through journal files was implemented, using all necessary boundary conditions obtained from the 1D script. The "InletDevelopment"-model previously used for obtaining initial velocity distributions for the droplets and vapor was removed, without decreasing accuracy in the final results. By using the TUI through journal files instead of the GUI, computational time decreased in many cases from several hours to just below one hour. Based on this, it can be determined that merely constantly updating the GUI requires a considerable amount of time for Fluent.

The error in the amount of initial film injected not matching the amount manually calculated was fixed. The post-processing has been made more efficient and more trustworthy since the data from Fluent is now averaged over several time steps, reducing noise from individual time steps.

(28)

7 Future Work

Since this model has been constructed based on a very basic geometry, there is much to unravel by developing the geometry to more complex examples. No doubt, it would be desirable to reform the geometry to perfectly match flow channels between the fuel rods in the BWR core. In order to do so, the simulations cannot rely on correlations in 1D, performed on other geometries. In this work, entrainment is based on the Okawa correlation whereas deposition is modeled using LPT. Replacing the current entrainment correlation is therefore absolutely necessary for further studies on other geometries than perfectly shaped cylindrical pipes. Reducing the sensitivity for deposition would also be a major step towards an accurate 3D model.

As previously mentioned, the transverse velocity for entrained droplets provided by [2] was a success for the cases tried in this work. Further studies on this matter is needed in order to validate its functionality in the model. Future work on the matter should also investigate what might cause the convergence errors appearing for transverse velocities above 60% of the calculated value.

Knowing the relationship between the boundary condition "z-momentum flux" and the initial film thickness would be advantageous for future simulations with this model. Since entrainment is directly proportional to film thickness, it would allow for a more accurate initial prediction of entrainment as well.

The suggestion of implementing an unstructured mesh was first proposed by Raddino. This would require major changes in the UDF, since its construction is based on the structured mesh currently used. It would also work as another tool for investigating the mesh dependency of the model. Making a transition from a structured to unstructured mesh paves the way for easier future development of more complex geometries. A first change in geometry would be an extension of the pipe in order to perform simulations from the onset of annular flow.

(29)

8 Appendix A

In order to implement all boundary conditions correctly, this 1D script with attached functions based on the Okawa correlations were created. To calculate all properties of the saturated water, steam table functions in "XSteam" were used, provided by Magnus Holmgren. The file "bc.inp" contains most if the inputs and are read in by the script at the start, defining several of the present conditions. The first script in this appendix is the major 1D script used for calculating the boundary conditions. Below, is first the function for determining the heat flux based on the outlet steam quality. Secondly, the function to determine the fraction between droplets and liquid film at the onset of annular flow based on the equilibrium conditions. Next is the function to determine the film thickness and finally, the functions for determining entrainment and deposition according to the correlations developed by Okawa.

1 close all 2

3 inputf=’bc.inp’;

4 bc=readInputFile(inputf);

5

6 % List of available power distributions

7 inlet = @(z) 0.01811∗z.^6−0.23174∗z.^5+1.179∗z.^4−2.878∗z.^3+2.8861∗z .^2−0.11836∗z+0.4386179;

8 middle= @(z) −0.041572∗z.^6+0.43222∗z.^5−1.5111∗z.^4+1.7822∗z .^3−0.036178∗z.^2+0.094166∗z+0.33886014;

9 outlet= @(z) −0.012984∗z.^6+0.13411∗z.^5−0.51332∗z.^4+0.78012∗z .^3−0.28829∗z.^2+0.20122∗z+0.5568055;

10 uniform=@(z) 1+0∗z;

11 pow_distrib=uniform;

12

13 %input variables

14 r=bc.D/2;%pipe radius [m]

15 length=bc.L;% [m]

16 pressure=bc.POUTL;%[N/m2]

17 totalflow=bc.FLOWFUE;%total mass flow rate [kg/s]

18 total_power=bc.ACTPOW;%total power [J/s]

19 segments=bc.CZMESH;%power segments 20 distrib=bc.PINPOW;%power distribution 21

22 N=1000;%number of elements 23

24 % Properties

25 g=9.81;%gravitational acceleration

26 rho_liq=XSteam(’rhoL_p’,pressure/1E5);%density of liquid 27 rho_vap=XSteam(’rhoV_p’,pressure/1E5);%density of vapor

28 t_sat=XSteam(’tsat_p’,pressure/1E5)−0.0001;%saturated temperature 29 temp=t_sat−10;% 10K subcooling

30 mu_liq=XSteam(’my_pT’,pressure/1E5,t_sat);%saturated liquid viscocity 31 hin=XSteam(’h_pT’,pressure/1E5,temp)∗1E3;%enthalpy at inlet [J/kg]

(30)

32 hl=XSteam(’hL_p’,pressure/1E5)∗1E3;%saturated liquid enthalpy [J/kg]

33 hv=XSteam(’hV_p’,pressure/1E5)∗1E3;%saturated vapor enthalpy [J/kg]

34 sigma=XSteam(’st_p’,pressure/1E5);%surface tension [N/m]

35 quality=@(h) (h−hl)./(hv−hl);

36 H=hv−hl;%latent heat of vaporization 37

38 x_outlet=0.43;%steam quality at outlet

39 q=get_heat_flux(pressure,temp,totalflow,x_outlet,r,length);%function to determine heat flux

40

41 % Enthalpy and quality distribution 42 z=linspace(0,length,N);

43 PINPOW=distrib./sum(distrib);% Normalization 44 h=zeros(1,N);

45

46 fori=1:1000

47 h(i)=hin+2∗pi∗r∗q∗integral(pow_distrib,0,z(i))/totalflow;

48 end

49 x=quality(h);%steam quality as function of position 50

51 % Onset of annular flow

52 vap_factor=sqrt(rho_vap/(g∗2∗r∗(rho_liq−rho_vap)));%factor multiplied with J to gain dimensionless J_vapstar

53 liq_factor=sqrt(rho_liq/(g∗2∗r∗(rho_liq−rho_vap)));%factor multiplied with J to gain dimensionless J_liqstar

54 x_onset_from_Wallis=(0.4∗pi∗r∗r/totalflow+0.6∗liq_factor/rho_liq)/(

vap_factor/rho_vap+0.6∗liq_factor/rho_liq);

55 z_onset_from_Wallis=interp1(x,z,x_onset_from_Wallis);% [m]

56

57 %1D model

58 step_size=length/(N−1);

59 d_frac_answer=droplet_fraction(x_onset_from_Wallis,mu_liq,r, totalflow,rho_liq,rho_vap,sigma);%call droplet fraction function 60 Wf_onset=(1−x_onset_from_Wallis)∗(1−d_frac_answer)∗totalflow;%

film flow rate at onset [kg/s]

61 Wd_onset=(1−x_onset_from_Wallis)∗d_frac_answer∗totalflow;

62

63 % check where onset occurs 64 k=0;

65 onset_value=0;

66 whileonset_value<z_onset_from_Wallis 67 onset_value=onset_value+step_size;

68 k=k+1;%onset occurs at k−th element in z 69 end

70

71 %create vectors that goes from onset to end 72 dep_okawa=zeros(1,N−k);

73 ent_okawa=zeros(1,N−k);

74 Wg=zeros(1,N−k+1);

75 Wf=zeros(1,N−k+1);

(31)

76 Wd=zeros(1,N−k);

77 film_thickness=zeros(1,N−k+1);

78 steamvel=zeros(1,N−k+1);

79 film_area=zeros(1,N−k);

80 film_velocity=zeros(1,N−k);

81

82 %entrainment and deposition are first set to zero 83 Ment=0;

84 Mdep=0;

85

86 %starting values for vapor and film flow rate 87 Wg(1)=x_onset_from_Wallis∗totalflow;

88 Wf(1)=Wf_onset;

89 film_thickness(1)=r/2;%guess starting value for film thickness 90 fori=1:N−k%start loop from onset and onwards

91 Wg(i+1)=Wg(1)+(q/H)∗integral(pow_distrib,z(k),z(k+i))∗pi∗2∗r;%vapor flow rate [kg/s]

92 steamvel(i+1)=Wg(i+1)/(rho_vap∗pi∗r∗r);%superfiscial steam velocity [m/s]

93

94 %call film thickness function

95 film_thickness(i+1)=get_film_thickness(Wf(i),film_thickness(i),mu_liq,r, Wg(i),rho_liq,rho_vap);% [m]

96 film_area(i)=pi∗(2∗r∗film_thickness(i+1)−film_thickness(i+1)∗

film_thickness(i+1));% [m2]

97 film_velocity(i)=Wf(i)/(rho_liq∗film_area(i));% [m/s]

98

99 Wd(i)=totalflow−Wg(i)−Wf(i);%droplet mass flow [kg/s]

100 ent_okawa(i)=entrainment(film_thickness(i+1),Wg(i),rho_liq,rho_vap, sigma,r);%call entrainment function [kg/m2s]

101 dep_okawa(i)=deposition(Wg(i),Wd(i),rho_liq,rho_vap,sigma,r);%call deposition function [kg/m2s]

102 Mdep=Mdep+dep_okawa(i)∗step_size∗pi∗2∗r;%add and multiply over area 103 Ment=Ment+ent_okawa(i)∗step_size∗pi∗2∗r;%add and multiply over area 104

105 Wf(i+1)=Wf(1)−Ment+Mdep−Wg(i+1)+Wg(1);%calculate new film flow rate from mass balance [kg/s]

106 ifWf(i+1)<0

107 Wf(i+1)=0;%since film flow rate cannot be negative 108 end

109 end%loop from onset to end of pipe 110

111 % DETERMINE INITIAL CONDITIONS 112 % liquid phase in Main model injected at z=2.59 113 % vapor phase in Main model injected at z=2.48

114 z_Main_liq=round(interp1(z,(1:N),2.59));%obtain element in z where z=2.59 115 z_Main_vap=round(interp1(z,(1:N),2.48));%obtain element in z where z=2.48 116

117 % initial conditions for Main model

118 Main_first_inlet_steamvel=steamvel(z_Main_vap−k+1);%vapor velocity at first Main inlet (z=2.48)

(32)

119 Main_second_inlet_steamvel=steamvel(z_Main_liq−k+1);%vapor velocity at second Main inlet (z=2.59)

120 Main_inlet_dropletmassflux=Wd(z_Main_liq−k)/(pi∗r∗r);%droplet mass flow [kg/m2s]

121

122 % initial film mass flux boundary condition

123 injection_wall_area=0.0042847384;% area of injection wall obtained from Fluent

124 %injection_wall_area=0.020812442; % area of of injection wall for "longpipe"−

geometry

125 filmBC=Wf(z_Main_liq−k)/injection_wall_area;%film mass flux [kg/m2s]

126

127 %plot section of pipe from onset to end 128 z_annular=(z(k+1):step_size:z(N));

129 figure(2)

130 plot(z_annular,dep_okawa,’b’,z_annular,ent_okawa,’r’) 131 title(’Deposition and entrainment’)

132 ylabel(’Deposition and entrainment [kg/m2s]’) 133 xlabel(’z [m]’)

134 legend(’Deposition’,’Entrainment’) 135 gridon

136

137 %creating vectors for final section of pipe (Main model) 138 dep_main=zeros(1,N−z_Main_liq+1);

139 ent_main=zeros(1,N−z_Main_liq+1);

140 filmthickness_main=zeros(1,N−z_Main_liq+1);

141 steamvel_main=zeros(1,N−z_Main_liq+1);

142 filmvel_main=zeros(1,N−z_Main_liq+1);

143 film_area_main=zeros(1,N−z_Main_liq+1);

144 Wf_main=zeros(1,N−z_Main_liq+1);

145 Wg_main=zeros(1,N−z_Main_liq+1);

146 Wd_main=zeros(1,N−z_Main_liq+1);

147 z_main=(z(z_Main_liq):step_size:z(N));

148

149 %calculating values for final section (last 1.06 m) 150 fori=z_Main_liq:N

151 dep_main(i−z_Main_liq+1)=dep_okawa(i−k);

152 ent_main(i−z_Main_liq+1)=ent_okawa(i−k);

153 filmthickness_main(i−z_Main_liq+1)=film_thickness(i−k);

154 steamvel_main(i−z_Main_liq+1)=steamvel(i−k+1);

155 filmvel_main(i−z_Main_liq+1)=film_velocity(i−k);

156 film_area_main(i−z_Main_liq+1)=pi∗(2∗r∗film_thickness(i−k)−

film_thickness(i−k)∗film_thickness(i−k));

157 Wf_main(i−z_Main_liq+1)=rho_liq∗filmvel_main(i−z_Main_liq+1)∗

film_area_main(i−z_Main_liq+1);

158 Wg_main(i−z_Main_liq+1)=Wg(i−k);

159 Wd_main(i−z_Main_liq+1)=Wd(i−k);

160 end 161

162 %plot final section

(33)

163 figure(3)

164 plot(z_main,dep_main,’b’,z_main,ent_main,’r’) 165 title(’Deposition and entrainment’)

166 ylabel(’Deposition and entrainment [kg/m2s]’) 167 xlabel(’z [m]’)

168 legend(’Deposition’,’Entrainment’) 169 gridon

170

171 Wf(end)=[];

172 Wg(end)=[];

173

174 figure(4)

175 plot(z_annular,Wf,’r’,z_annular,Wg,’k’,z_annular,Wd,’b’) 176 title(’Film flow rates’)

177 xlabel(’z [m]’)

178 ylabel(’Mass flow rates [kg/s]’) 179 legend(’Liquid film’,’Vapor’,’Droplets’) 180 gridon

1 function[ heat_flux ] = get_heat_flux( pressure,temp,totalflow,x_outlet,r, length)

2 %Input: pressure [Pa], temperature [K], massflow [kg/s], 3 %steam quality at outlet [dim−less (%)]

4 %Output: heat flux [W/m2]

5

6 hin=XSteam(’h_pT’,pressure/1E5,temp)∗1E3;%enthalpy at inlet [J/kg]

7 hf=XSteam(’hL_p’,pressure/1E5)∗1E3;%saturated liquid enthalpy [J/kg]

8 hg=XSteam(’hV_p’,pressure/1E5)∗1E3;%saturated vapor enthalpy [J/kg]

9 heat_flux=(x_outlet∗(hg−hf)+hf−hin)∗totalflow/(2∗pi∗r∗length);

10 11 end

1 function[ d_frac_answer] =droplet_fraction( x_onset_from_Wallis, mu_liq,r,totalflow,rho_liq,rho_vap,sigma )

2 %Output: droplet/film fraction [dim−less]

3 %Input: steam quality at onset [dim−less], viscocity at onset [Pa s]

4

5 film_thickness=r/2;%starting value for film thickness (r/2) 6 d_step_size=0.001;

7 d_frac=d_step_size;%starting value for droplet/film fraction 8 dep=zeros(1,1/d_step_size);

9 ent=zeros(1,1/d_step_size);

10 Wg=x_onset_from_Wallis∗totalflow;%steam flow rate [kg/s]

11 fori=1:1/d_step_size%loop to find droplet fraction

12 Wf=(1−x_onset_from_Wallis)∗(1−d_frac)∗totalflow;%film flow rate [kg/s]

13 Wd=(1−x_onset_from_Wallis)∗d_frac∗totalflow;%droplet flow rate [kg/s]

14 film_thickness=get_film_thickness( Wf,film_thickness,mu_liq,r,Wg, rho_liq,rho_vap);%call film thickness function

15

CFD Annular Flow Modelling Based on a Three-Field Approach