Mechanistic Modeling of Annular Flows at Near-Dryout Conditions

(1)

Doctoral Thesis in Physics

Mechanistic Modeling of Annular Flows at Near-Dryout Conditions

WENYUAN FAN

Stockholm, Sweden 2020 www.kth.se

ISBN 978-91-7873-703-1 TRITA-SCI-FOU 2020:32

kth royal institute of technology

Wenyuan FanMechanistic Modeling of Annular Flows at Near-Dryout ConditionsKTH 2020

(2)

Mechanistic Modeling of Annular Flows at Near-Dryout Conditions

WENYUAN FAN

Doctoral Thesis in Physics KTH Royal Institute of Technology Stockholm, Sweden 2020

Academic Dissertation which, with due permission of the KTH Royal Institute of Technology, is submitted for public defence for the Degree of Doctor of Philosophy on Wednesday, December 2, 2020, at 10:00 a.m. in FA32, Roslagstullsbacken 21, Stockholm.

(3)

Printed by: Universitetsservice US-AB, Sweden 2020

(4)

iii

Abstract

This work aims at advancing the modeling methodology for near-dryout annular flows using two different computational fluid dynamics (CFD) approaches. One is termed as the transported film approach, where the film flow is separately modeled by a simplified liquid film model and then coupled with the gas core simulation using various closures. The other is referred to as the film-resolving approach, where flow details in the film and the gas core are directly captured.

The transported film approach is computationally efficient and able to simulate certain simplified industrial annular flows with dryout. In order to simulate more realistic scenarios, several new models have been introduced. A rewetting model was developed to make the simulation capable of modeling both the occurrence and disappearance of dryout. Carefully designed coupling schemes were used to model the conjugate heat transfer (CHT) between the annular flow and the heating structure, where the Joule heating effect can be considered. Therefore, CHT simulations were able to predict the outer wall temperature, allowing a direct comparison with experiments. In addition, thanks to the modeling of the thermal inertia of the heating structure, the complex dryout hysteresis was successfully captured by the CHT simulations.

The volume of fluid (VOF) method was employed for the film-resolving approach, where the Reynolds-averaged Navier-Stokes (RANS) approach was used for turbulence modeling. A spectrum of topics regarding the RANS-VOF approach have been discussed due to its immaturity. A long-existing implementation issue in OpenFOAM has been discovered and corrected such that self-consistent RANS- VOF formulations could be used. Tests on air-water stratified flows showed that the new implementation was able to qualitatively capture the turbulence behav- ior around the two-phase interface, while the native implementation failed catas- trophically and substantially over-predicted the turbulence level. In addition, the self-consistent implementation outperformed the native one in terms of the insen- sitivity to mesh refinement, which is crucial for numerical simulations. Still, the self-consistent implementation quantitatively over-predicted the interfacial turbulence level due to the usage of single-phase turbulence models, which was then phenomenologically corrected by a newly modified turbulence damping model.

RANS-VOF simulations were then carried out to investigate the formation and development of disturbance waves in a downward air-water annular flow. Efficient post processing methods were developed to calculate the film thickness from CFD data, allowing an intuitive comparison with the corresponding experiment. Con- trary to stratified flows where the self-consistent implementation over-predicted the interfacial turbulence level, the turbulence level was under-predicted for thin films in annular flows due to the inadequate modeling of the complex wall-interface-

(5)

iv

turbulence interaction. This unique challenge was phenomenologically resolved by further improving the turbulence damping model, with which the CFD simulations were able to reproduce the formation and development of disturbance waves. The close relation between entrainment and disturbance waves has also been captured.

(6)

v

Sammanfattning

Syftet med avhandlingen är att utveckla beräkningsmetoder för ringflöden nära tor- kokning med hjälp av två olika beräkningsströmningsdynamik (CFD)-modeller. En betecknas som filmtransportmodellen, där filmflödet modelleras separat av en för- enklad vätskefilmmodell och sedan kopplas med gaskärnsimuleringen med olika kompletterande ekvationer. Den andra kallas filmlösningsmodellen, där flödesde- taljer i filmen och gaskärnan fångas direkt.

Filmtransportmodellen är beräkningseffektiv och kan simulera vissa förenklade industriella ringflöden med torrkokning. För att simulera mer realistiska scenari- er har flera nya modeller introducerats. En återvätningsmodell utvecklades för att göra simuleringen kapabel att modellera både förekomst och försvinnande av torrkokning. Noggrant utformade kopplingsscheman användes för att modellera kon- jugerad värmeöverföringen (CHT) mellan ringflödet och uppvärmningsstrukturen, där Joule-uppvärmningseffekten kan övervägas. Därför kunde CHT-simuleringar förutsäga ytterväggstemperaturen, vilket möjliggjorde en direkt jämförelse med experimentdata. Dessutom tack vare modelleringen av värmestrukturens termis- ka tröghet fångades den komplexa torrkokningshysteresen framgångsrikt av CHT- simuleringarna.

Volume-of-fluid-metoden (VOF) användes för filmlösningsmodellen, där den Reynolds-Averaged Navier-Stokes (RANS) -metoden användes för turbulensmo- dellering. Ett spektrum av ämnen angående RANS-VOF-metoden har diskute- rats på grund av dess omogenhet. Ett länge existerande implementeringsproblem i OpenFOAM har upptäckts och korrigerats så att självkonsistenta RANS-VOF- formuleringar kan användas. Tester på stratifierade luft och vattenflödet visade att den nya implementeringen kunde kvalitativt fånga upp turbulensbeteendet runt två- fasgränssnittet, medan den ursprungliga implementeringen misslyckades katastro- falt och väsentligt felberäknade turbulensnivån. Dessutom överträffade den själv- konsistenta implementeringen den ursprungliga när det gäller okänslighet för nät- förfining, vilket är avgörande för numeriska simuleringar. Slutligen förutspådde det självkonsistenta nya modellen kvantitativt tvåfasgränsytans turbulensnivå med hjälp av enfasbaserade turbulensmodeller, som sedan korrigerades fenomenolo- giskt med en nyligen modifierad turbulensdämpningsmodell.

RANS-VOF-simuleringar genomfördes sedan för att undersöka bildandet och utvecklingen av störningsvågor i ett nedåtgående luft-vatten ringflöde. Effektiva ef- terbehandlingsmetoder utvecklades för att beräkna filmtjockleken från CFD-data, vilket möjliggjorde en intuitiv jämförelse med motsvarande experimentdata. I mot- sats till stratifierade flöden där den självkonsistenta implementeringen beräknade för hög tvåfasgränsytans turbulensnivå var turbulensnivån beräknad för låg för tun- na filmer i ringflöden på grund av den otillräckliga modelleringen av den komplexa

(7)

vi

vägg-tvåfasgränssnittet-turbulensinteraktionen. Denna unika utmaning löstes feno- menologiskt genom att ytterligare förbättra turbulensdämpningsmodellen, med vil- ken CFD-simuleringarna kunde reproducera bildandet och utvecklingen av stör- ningsvågor. Det nära sambandet mellan droppavlägsning från vattenfilmsyta och störningsvågor har också fångats.

(8)

vii

List of Publications

Publications included in this thesis

Paper 1 W. Fan, H. Li, and H. Anglart. A study of rewetting and conjugate heat transfer influence on dryout and post-dryout phenomena with a multi- domain coupled CFD approach. International Journal of Heat and Mass Transfer, 163:120503, 2020b. doi: 10.1016/j.ijheatmasstransfer.2020.

120503.

Paper 2 W. Fan and H. Anglart. varRhoTurbVOF: A new set of volume of fluid solvers for turbulent isothermal multiphase flows in OpenFOAM. Com- puter Physics Communications, 247:106876, 2020a. doi: 10.1016/j.cpc.

2019.106876.

Paper 3 W. Fan and H. Anglart. Progress in Phenomenological Modeling of Tur- bulence Damping around a Two-Phase Interface. Fluids, 4(3):136, 2019.

doi: 10.3390/fluids4030136.

Paper 4 W. Fan, H. Li, and H. Anglart. Numerical investigation of spatial and temporal structure of annular flow with disturbance waves. Interna- tional Journal of Multiphase Flow, 110:256–272, 2019b. doi: 10.1016/

j.ijmultiphaseflow.2018.10.003.

Paper 5 W. Fan, A. V. Cherdantsev, and H. Anglart. Experimental and numerical study of formation and development of disturbance waves in annular gas-liquid flow. Energy, 207:118309, 2020a. doi: 10.1016/j.energy.2020.

118309.

Author’s contribution to the included publications

Wenyuan Fan made major contribution to the first four publications and the simulation part of the fifth paper, including conceptualization, methodology, code implementation, calculation and analysis, and writing the draft.

Publications not included in this thesis

• W. Fan and H. Anglart. varRhoTurbVOF 2: Modified OpenFOAM volume of fluid solvers with advanced turbulence modeling capability. Computer Physics Communications, 256:107467, 2020b. doi: 10.1016/j.cpc.2020.

107467.

(9)

viii

• M. Bergagio, W. Fan, R. Thiele, and H. Anglart. Large eddy simulation of thermal mixing with conjugate heat transfer at BWR operating conditions. Nuclear Engineering and Design, 356:110361, 2020. doi:

10.1016/j.nucengdes.2019.110361.

• W. Fan, A. Cherdantsev, H. Li, and H. Anglart. A physical understanding of droplet entrainment in annular flow based on numerical data. In Inter- national Conference on Multiphase Flow, Rio de Janeiro, Brazil, 2019a.

• A. Cherdantsev, W. Fan, H. Li, and H. Anglart. Numerical and experimental study of formation and development of disturbance waves in annular gas-liquid flow. In International Conference on Multiphase Flow, Rio de Janeiro, Brazil, 2019.

• W. Fan, H. Li, and H. Anglart. Prediction of annular two-phase flow with heat transfer. In International Topical Meeting on Nuclear Reactor Ther- mal Hydraulics, pages 5230–5238, Portland, Oregon, USA, 2019c.

(10)

ix

Acknowledgments

I am sincerely grateful to my supervisor Prof. Henryk Anglart for his continuous support and insightful guidance, without which the current achievements would not be possible. I would like to thank my co-supervisor Dr. Haipeng Li for all the valu- able discussions. I am glad to have collaborated with Dr. Andrey V. Cherdantsev, from whom I learned the state-of-the-art knowledge of disturbance waves. I would also like to thank Dr. Jean-Marie Le Corre for his comments and suggestions on the presentation of this thesis.

I extend my thankfulness to all my colleagues for creating a nice working en- vironment. Special thanks to Dr. Mattia Bergagio for generously sharing his expe- rience with OpenFOAM, Python and Linux.

I would like to express my gratitude and love to my girlfriend Yangli for her ongoing support.

The majority of the presented simulations were carried out using resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC partially funded by the Swedish Research Council. Dr. Jing Gong at PDC is acknowledged for assistance concerning technical aspects in making code run on the PDC resources.

Financial supports from China Scholarship Council (CSC), KTH Royal Insti- tute of Technology, and the Knut and Alice Wallenberg Foundation are gratefully acknowledged.

(11)

(12)

List of Figures

1.1 Annular flow, dryout and disturbance waves. . . . 2

1.2 Multi-field approach. . . . 4

1.3 Transported film approach. . . . 5

1.4 Film-resolving approach. . . . 6

2.1 Local coordinate system for the film. . . . 12

2.2 Grids for film-gas-core coupling. . . . 14

2.3 Rewetting model for the liquid film region. . . . 16

2.4 Coupling scheme for standalone simulations. . . . . 17

2.5 Coupling scheme without thermal inertia. . . . 18

2.6 Coupling scheme with thermal inertia. . . . 19

2.7 Thermal and electrical conductivities. . . . 20

2.8 Grid sensitivity studies. . . . 21

2.9 Temperature comparison for simulations without thermal inertia. . 23

2.10 Dimensionless wall heat flux comparison. . . . 24

2.11 Heat flux during power transient. . . . 25

2.12 Film coverage at the wall during the power transient. . . . . 26

3.1 Flow configuration for the air-water stratified flow. . . . 33

3.2 Results obtained by native and new OpenFOAM VOF solvers for Run-250. . . . . 34

3.3 ν_tobtained by native and new implementations. . . . 35

3.4 ∆n calculation for a polyhedron with volume V . . . . . 39

3.5 Pressure gradients for Run-250 obtained by different δy and B values. 40 3.6 k profiles for Run-250 with δy = 0.5 mm and different B values. . 41

3.7 k profiles for Run-250 with δy = 0.5 mm and the asymmetric damping treatment. . . . 42

4.1 Experiment setup. . . . 45

xv

(17)

xvi LIST OF FIGURES

4.2 Computational domain. . . . 46

4.3 α field calculated by the VOF method. . . . 48

4.4 Glob identification. . . . 49

4.5 Identified liquid film. . . . 49

4.6 Local film thickness calculation. . . . 50

4.7 Spatiotemporal matrix for the experimental film thickness. . . . . 51

4.8 Droplet entrainment from a disturbance wave. . . . 52

4.9 δ(x, t) matrices without damping. . . . 54

4.10 δ(x, t) obtained by the simple correction approach. . . . 55

4.11 δ(x, t) obtained by the complex correction approach. . . . 56

4.12 Droplet entrainment in simulations. . . . . 57

4.13 Instantaneous film thickness distribution. . . . 58

4.14 Instantaneous axial wall shear stress. . . . 59

4.15 Instantaneous y_i⁺distribution. . . . . 60

(18)

List of Tables

1.1 Film-resolving CFD studies. . . . 7

2.1 Key parameters for simulated experiments. . . . 19

2.2 Discretized meshes for dryout simulations. . . . 20

4.1 Boundary conditions. . . . 47

4.2 Discretized meshes for disturbance wave simulations. . . . 53

xvii

(19)

(20)

List of Symbols

Acronyms

BBLIF brightness-based laser-induced fluorescence BWR boiling water reactor

CCL connected-component labeling CFD computational fluid dynamics CHF critical heat flux

CHT conjugate heat transfer CSF continuum surface force DNB departure from nucleate boiling DNS direct numerical simulation HPM hybrid particle-mesh LES large eddy simulation LS level-set

PDE partial differential equation RANS Reynolds-averaged Navier-Stokes RoI region of interrogation

SST shear stress transport VOF volume-of-fluid

WALE wall-adapting local eddy-viscosity

xix

(21)

xx LIST OF SYMBOLS

Latin symbols

f body force

g gravitational vector I identity tensor n unit normal vector

n_o unit normal vector of the outlet boundary q heat flux vector

S vector source term

Si surface vector of discretized face i

u velocity

u^turb velocity used in turbulence models x position vector

A interfacial area density

B damping factor

C the offset coefficient in the log law of velocity c_p isobaric specific heat

Dh hydraulic diameter di inner diameter do outer diameter

D_φ effective diffusion coefficient for φ E₆₀₀ entrained fraction at x = 600 mm fi interfacial friction factor

G mass flux

h specific enthalpy

h_eva heat transfer coefficient for droplet evaporation

(22)

LIST OF SYMBOLS xxi

k turbulence kinetic energy

K_ent entrainment mass transfer coefficient

p pressure

p^∗ modified pressure

prgh actual pressure solved by OpenFOAM Q volumetric energy source

q heat flux

Ql volumetric flow rate of the liquid phase Rel Reynolds number for the liquid phase S scalar source term

S_p projected area of a discretized cell onto a two-phase interface

T temperature

t time

Tsub water subcooling at inlet u⁺ dimensionless velocity u_τ friction velocity

u^k_τ friction velocity based on k V volume of a discretized cell Ve electric potential

w wet wall indicator x downstream distance y physical distance to the wall

y^∗ dimensionless wall distance based on u^k_τ y⁺ dimensionless wall distance based on uτ

y_p typical cell size in the interface normal direction

(23)

xxii LIST OF SYMBOLS

Greek symbols

α volumetric fraction

α0 threshold value for liquid phase identification β turbulence model coefficient

τ shear stress tensor τ_a, τb wall-tangential vectors τ_i interfacial shear stress τw wall shear stress

τef f effective shear stress tensor

∆n typical grid cell size across the interface

∆V_e voltage difference

∆y a typical length for turbulence damping δy size of a discretized cell

δ film thickness

∆n distance between the cell center and the boundary

turbulence dissipation rate

Γ mass source

γ threshold value used in film rewetting model

κ thermal conductivity; interface curvature; velocity log law coefficient µ dynamic viscosity

µ_t turbulent dynamic viscosity

µ⁰_t µtcalculated by single-phase turbulence models µ^σ_t µtintroduced by interfacial turbulence corrections ν kinematic viscosity

ν_t turbulent kinematic viscosity

(24)

LIST OF SYMBOLS xxiii

ω specific turbulence dissipation rate φ a generic variable

ψ symmetry factor for turbulence damping

ρ density

σ surface tension coefficient σ_e electrical conductivity

Operators

∇ nabla

∇n surface normal component of ∇

∇s surface tangential component of ∇ depth averaging

time (ensemble) averaging

˜ Favre averaging

˜˜ dynamic viscosity weighted averaging

Superscripts

0 fluctuation component in Reynolds decomposition

00 fluctuation component in Favre decomposition

∗ difference between exact and modeled forms

E exact form

G generic variable

Subscripts

σ surface tension b value at the boundary

d droplet

f film

(25)

xxiv LIST OF SYMBOLS

g gas

i value at the interface; iterator k phase indicator

l liquid

nb cell value adjacent to the boundary nw value at the near-wall cell

s solid

v vapor

w wall value

1 primary phase

2 secondary phase

(26)

Chapter 1 Introduction

A boiling water reactor (BWR) is a light water nuclear reactor that uses fission energy to heat water and produce steam in its reactor core. The steam is then used to drive a turbine that is connected to a generator. Therefore, generating more steam from the reactor core will improve the overall efficiency of the power plant.

However, a phenomenon called dryout limits our pursuit for higher efficiencies, as detailed in the following.

Annular flow is the two-phase flow regime prior to the occurrence of dryout, as depicted in Fig. 1.1. In this regime, there is a thin liquid film flowing on the channel wall and being sheared by a fast-moving gas core. The interaction between the film and the gas core produces a system of complex waves on the film surface. The dom- inant type of such wavy structures is referred to as disturbance waves, which are high-amplitude, long-spacing and fast-moving waves traveling on the residual layer of liquid that is usually referred to as substrate or base film. Disturbance waves may partially get destructed by the turbulent gas core forming dispersed droplets. This process is referred to as entrainment. When getting enough momentum towards the wall, such entrained droplets can return to the film. This process is called deposition. In a heated channel, the heat is taken away by the evaporation of the film.

This heat transfer mechanism is quite efficient thanks to the high latent heat of va- porization. Dryout occurs when the film is totally depleted, leading to the channel wall in direct contact with the gas core.

In the pre-dryout region, the wall is maintained at a temperature slightly higher than the boiling temperature. However, the wall temperature could be several hun- dred kelvins higher in the post-dryout region since the heat transfer process is dominated by the much less efficient gas core convection. For BWRs, such a dramatic temperature increase may damage the integrity of the fuel and must be avoided.

The heat flux corresponding to the occurrence of dryout is the critical heat flux

1

(27)

2 CHAPTER 1. INTRODUCTION

Figure 1.1: Annular flow and dryout in a heated flow channel. Wall temperature is proportional to the redness.

(28)

1.1. EMPIRICAL MODELING 3

(CHF) for BWRs. Therefore, a key task in BWR thermal-hydraulic calculations is to predict dryout-related phenomena accurately. Such calculations can be carried out using different levels of modeling.

1.1 Empirical modeling

The traditional approach for simulating dryout-related phenomena is based on one- dimensional empirical correlations, where macroscopic quantities, e.g., superfi- cial velocities and cross-sectional void fraction (Anglart, 2014; Nguyen and Moon, 2015; Cheng et al., 2018; Yu et al., 2018) are used. This approach is computationally efficient since there is no need to simulate flow details of the film and droplets.

Moreover, empirical correlations can reach rather high accuracies provided that there are adequate high-quality experimental data. However, this empirical modeling approach has two main disadvantages. On the one hand, empirical correlations are constructed based on certain apparatus and flow conditions. Therefore, ex- trapolating these correlations to other flow configurations may lead to inaccurate predictions. On the other hand, such correlations are usually constructed based on steady state experimental data. Therefore, their application to transient simulations should be conducted with caution.

1.2 CFD modeling

The computational fluid dynamics (CFD) method has become an attractive alter- native for the modeling of dryout-related phenomena thanks to the development in multiphase CFD algorithms and computing hardware. The reason is that CFD simulations focus on detailed and localized flow variable distributions instead of macroscopic values. This is achieved by solving basic conservation equations in a discretized 2D or 3D computational domain. Depending on how the thin liquid film is modeled, CFD simulations are categorized into three different groups, as discussed in the following.

1.2.1 Multi-field model approach

The term “field” in the multi-field model approach refers to a collection of fluid based on a certain criterion. For instance, all the liquid phase can be treated as a single liquid field, which can be further divided into a continuous liquid field and a discrete liquid field. Then each field is described by a corresponding volumetric fraction and interpenetrates with other fields, as shown in Fig. 1.2, where the real two-phase flow inside a near-wall control volume is represented by different multi- field approaches. In the multi-field model approach, each field has its own flow

(29)

variables. Due to the interpenetration of fields, the detailed interface structure is not available in the multi-field model approach. Instead, various empirical correlations are used to model the field interactions.

Figure 1.2: Multi-field representations for the real two-phase flow inside a near- wall control volume.

Kumar and Edwards (1996) developed a four-field approach for the modeling of adiabatic annular flow. A wall lateral (lubrication) force was used to keep the thin film on the wall. Empirical correlations for deposition and entrainment were used locally as closures.

Shi et al. (2016, 2017) used a three-field approach to predict the pre- and post- dryout heat transfer. Entrainment and deposition were also modeled by using empirical correlations locally. The occurrence of dryout was locally determined by an empirical correlation for the vapor quality.

Tentner et al. (2018) and Pothukuchi et al. (2019) used the most basic two-field approach to model annular flows with dryout, where entrainment and deposition were not explicitly modeled. It should be noted that the adopted method was origi- nally developed for DNB-type, instead of dryout-type, CHF (Kurul and Podowski, 1990), where DNB stands for departure from nucleate boiling. In such an approach, dryout is triggered by a predefined void fraction.

The advantage of the multi-field model approach is that it is able to cover a wide range of flow regimes ranging from single-phase liquid to superheated vapor. How- ever, the accuracy of this approach heavily depends on the quality of those empirical correlations. There is an additional challenge for this approach when considering the special nature of the liquid film as shown in Fig. 1.2. The wall-normal size of the control volume significantly affects the corresponding void fraction and vapor quality in the multi-field approach, making the dryout occurrence dependent on the near-wall mesh size. Such an issue has already been spotted by Li et al. (2019)

(30)

1.2. CFD MODELING 5

when using the two-field approach to model the DNB-type CHF. Therefore, the thin film is preferably modeled by other CFD methods.

1.2.2 Transported film approach

The transported film approach (Li and Anglart, 2015, 2016a,b, 2017; Anglart et al., 2018; Raddino et al., 2018; Li and Anglart, 2019) separates the annular flow into a liquid film region and a gas core region, as depicted in Fig. 1.3. The liquid film region is modeled based on the depth-averaged liquid film model first developed by Bai and Gosman (1996), which will be further discussed in Section 2.1.1. As shown in Fig. 1.3, the liquid film model provides the film thickness on the wall, instead of solving its detailed 3D shape.

Figure 1.3: Transported film approach separates the annular flow into a liquid film region and a gas core region. The complex film shape is represented by a film thickness in the liquid film model.

The gas core has been modeled by either the two-field approach or the Eulerian- Lagrangian method. For the former, empirical correlations were used to model the entrainment and deposition processes to couple the gas core with the liquid film (Li and Anglart, 2016a,b; Anglart et al., 2018). However, only the entrainment process was modeled empirically in the Eulerian-Lagrangian approach since the deposition process was directly captured by solving the motion of individual droplets (Li and Anglart, 2015, 2016a, 2017; Anglart et al., 2018; Raddino et al., 2018; Li and Anglart, 2019).

In this approach, the dryout occurrence is determined by the local film thickness, which is provided by the liquid film model. It should be noted that, for the Eulerian-Lagrangian gas core modeling, only the location of dryout was modeled and there was no post-dryout heat transfer modeling (Li and Anglart, 2019), while a complete pre- and post-dryout heat transfer modeling methodology was developed by Li and Anglart (2016b) using the two-field gas core modeling approach.

(31)

The transported film approach is more mechanistic than the multi-field approach in the sense that the film thickness is directly calculated and used as the dryout criterion. However, certain empirical correlations are still required, especially for the entrainment modeling. In addition, it is not possible yet for the transported film approach to simulate bubbly, slug and churn flows.

1.2.3 Film-resolving approach

The film-resolving approach aims to simulate the flow details inside the thin film, which can be achieved by using the volume-of-fluid (VOF) method (Hirt and Nichols, 1981), the level-set (LS) method (Sussman et al., 1994) and the hybrid particle- mesh (HPM) method (Liu et al., 2005). No empirical closures are required in this approach to model the deposition and entrainment processes since they are auto- matically modeled by the corresponding CFD method. As depicted in Fig. 1.4, in comparison with the multi-field and transported film approaches, much finer grids should be used for the (pure Eulerian) film-resolving approach.

Figure 1.4: Fine grids are used in the Eulerian film-resolving approach to capture flow details inside the film.

A summary of some key features of several film-resolving CFD studies is given in Table 1.1. It is clearly seen that the VOF and LS methods, especially the former, have been widely adopted in such simulations. Direct numerical simulation (DNS), large eddy simulation (LES) and the Reynolds-averaged Navier-Stokes (RANS) approach have been used for the turbulence modeling. Even though laminar simulations were carried out by Kumar et al. (2016), they stated that more correct results would be obtained by LES or similar. It should be noted that both LES and RANS simulations were carried out using turbulence models developed for single-phase flows due to the lack of their two-phase counterparts.

(32)

1.2. CFD MODELING 7

Table1.1:Film-resolvingCFDstudiesongas-liquidannularflows. SourceInterfaceTurbulence L/DhBoundarytreatmentMeshquality modelingmodelingstreamwisespanwisey+Dh/δyw JayantiandHewitt(1997)prescribed interface

laminar, RANS3.33a periodicCartesian2D16000a HanandGabriel(2007)VOFRANS52.5flatinletprofilescylindrical2DNA69 Rodriguez(2009)LSDNS1.25periodicperiodic53373 Yunetal.(2010)HPMRANS2.4periodicCartesian2DNA128 Gulatietal.(2013)LSLES10bperiodicperiodicNA100c Kumaretal.(2016)VOFlaminar6flatinletprofilesfullgeometryNA171b Saxenaetal.(2016)VOFLES49.17inletprofilesfrom precursorsimulationfullgeometry<5NA Yangetal.(2017)LSLES2πperiodicperiodic<5100c SatoandNiceno(2017)VOFLES7.97randomizedflat inletprofilesperiodicNA520 SaxenaandPrasser(2020)VOFRANS, LES2πperiodicdomain withrecyclingfullgeometry<0.43770b anominalvaluesincethegascorewasnotmodeled b valuebasedonthefinestmesh c cellnumberinahydraulicdiameter

(33)

Even though the LES approach seems to be more frequently adopted in re- cent studies, it is not guaranteed to outperform the RANS approach in the film region due to the following reason. Since extremely dense grids are required for a wall-resolving LES approach (Choi and Moin, 2012), the wall-adapting local eddy-viscosity (WALE) model (Nicoud and Ducros, 1999) has been used for the near-wall modeling by Saxena et al. (2016); Yang et al. (2017); Saxena and Prasser (2020). However, as pointed out by Sato and Niceno (2017), the applicability of this wall-modeled LES approach for vaporizing thin liquid film is unclear. In addition, since we are more interested in the thin film than the highly turbulent gas core flow, RANS simulations seem to be more efficient than the LES approach.

The computing resource requirement has not been mentioned when discussing the multi-field and transported film approaches since it is simply not an issue for both methods, which have been both applied to computational domains of several meters long. However, this requirement prevents the film-resolving approach being more practically used as evidenced by the low L/Dhvalues in Table 1.1 with L being the domain length and Dhbeing the hydraulic diameter. The reason can be intuitively understood by comparing the mesh density differences shown in Figs.

1.2-1.4. In Table 1.1, δywdenotes the height of cells adjacent to the wall and y⁺is the corresponding dimensionless wall distance. It is clearly seen that extremely fine near-wall grids are required for high-quality simulations, even for the wall-modeled LES studies. Therefore, such simulations have been forced to be carried out with much shorter computational domains that must be used together with appropriate boundary conditions. As shown in Table 1.1, the periodic boundary condition has been widely used in both the streamwise and spanwise directions. More discussions on this treatment will be given in Section 4.1.1.

1.2.4 Conjugate heat transfer

In the aforementioned studies where post-dryout heat transfer was considered (Li and Anglart, 2016b; Shi et al., 2016, 2017; Tentner et al., 2018; Anglart et al., 2018;

Pothukuchi et al., 2019), only fluid regions, i.e., film, droplets and gas core, were modeled and the solid heating effect was modeled by predefined heat fluxes im- posed on the fluid boundaries. However, it is quite challenging to specify such heat fluxes accurately for dryout-related phenomena where there is a strong heat transfer coupling between the solid and the fluid. Therefore, this fluid-solid conjugate heat transfer (CHT) problem must be properly modeled to get more precise predictions for the post-dryout heat transfer. Such simulations have been recently carried out by Shi et al. (2019) and Vegendla et al. (2019) together with the multi-field approach.

It should be noted that only steady state results were presented in both studies.

(34)

1.3. THESIS STRUCTURE 9

1.3 Thesis structure

This thesis aims at advancing the CFD modeling methodology for both the transported film and film-resolving approaches. The open-source CFD toolbox Open- FOAM is used for model implementation.

In Chapter 2, the transported film approach is further developed to model the complex fluid-solid conjugate heat transfer in pre- and post-dryout regions.

Chapter 3 describes the basic theories of general RANS-VOF simulations and corrects a long-existing implementation issue of this approach in OpenFOAM. A widely used turbulence damping model is also discussed and modified in this chapter.

Chapter 4 presents numerical simulations on disturbance waves using the RANS- VOF approach without using periodic boundary conditions. The results are com- pared with corresponding experiments with the help of newly developed post processing methods.

Chapter 5 summarizes main achievements and findings of the thesis.

(35)

(36)

Chapter 2 Transported film modeling with conjugate heat transfer

2.1 Basic theories

2.1.1 Liquid film model

The full governing equations for a film with heat and mass transfer read

∂ρ_f

∂t + ∇ · (ρ_fu_f) = S_ρ, (2.1)

∂(ρ_fu_f)

∂t + ∇ · (ρ_fu_fu_f) = −∇p + ∇ · τ + S_ρu, (2.2)

∂(ρ_fh_f)

∂t + ∇ · (ρ_fu_fh_f) = ∇ · q + S_ρh, (2.3) where Sρ, Sρu and Sρh denotes the mass, momentum and energy source terms, respectively; τ is the shear stress tensor; q is heat flux.

The liquid film model is derived using a local coordinate system shown in Fig.

2.1. Since the film is usually so thin that the flow in the wall-normal direction is negligible, consequently, the advection in this direction is negligible in comparison with that in the wall-tangential direction. In addition, the assumption for boundary layers also applies to the thin film that the spatial gradients of dependent variables are dominated by the wall-normal direction (Bai and Gosman, 1996). As a result, the key idea of the liquid film model is that advection is only considered in the wall-tangential direction and diffusion is only considered in the wall-normal direction. With this assumption, we could simplify the aforementioned governing

11

(37)

12

CHAPTER 2. TRANSPORTED FILM MODELING WITH CONJUGATE HEAT TRANSFER

Figure 2.1: A local coordinate system defined by the wall-normal vector, n, and two wall-tangential vectors, τaand τb.

equations by splitting the ∇ operator into its wall-tangential component ∇sand its wall-normal component ∇n, as shown in the following equation:

∇ · φ = ∇_n· φ + ∇_s· φ. (2.4)

Advection terms are simplified to

∇ · φadvection≈ ∇s· φadvection, (2.5) and diffusion terms are treated as

∇ · φdif f usion≈ ∇n· φdif f usion. (2.6) Then we get the following simplified governing equations

∂ρf

∂t + ∇s· (ρuf) = Sρ, (2.7)

∂(ρ_fu_f)

∂t + ∇_s· (ρfu_fu_f) = −∇_sp + ∇_n· τ + Sρu, (2.8)

∂(ρfhf)

∂t + ∇s· (ρfufhf) = ∇n· q + Sρh. (2.9) We introduce the following depth-averaged operator in the wall-normal direction

φ =1 δ

Z δ 0

φ dn, (2.10)

where δ is the film thickness. Integrating Eqs. (2.7)-(2.9) in the wall-normal direction and assuming

Z δ 0

∇s· (ρfufφf) dn = ∇s· (ρfδufφf) ≈ ∇s· (ρfδufφ_f), (2.11)

(38)

2.1. BASIC THEORIES 13

we can get the following governing equations for the liquid film model

∂(ρ_fδ)

∂t + ∇s· (ρfδuf) = Sρδ, (2.12)

∂(ρfδuf)

∂t + ∇s· (ρfδufuf) = −δ∇sp + τi− τw+ Sρuδ, (2.13)

∂(ρfδhf)

∂t + ∇s· (ρfδufhf) = qi− qw+ Sρhδ, (2.14) where Sρδ= δSρ, Sρuδ= δSρu, and Sρhδ= δSρhare the depth-integrated source terms; τiand τware interfacial shear stress and wall shear stress, respectively; qi

and qware interfacial and wall heat fluxes, respectively. These equations are solved using the one-layer finite volume method, as detailed in Appendix C in Paper 1.

2.1.2 Two-fluid two-field model

The gas core is modeled by the following governing equations, which are based on the two-fluid two-field model (Ishii and Hibiki, 2006)

∂(αkρk)

∂t + ∇ · (α_kρ_ku_k) = Γ_k, (2.15)

∂(αkρkuk)

∂t + ∇ · (α_kρ_ku_ku_k) = −α_k∇p + ∇ · (α_kτ_k) + α_kρ_kg + S_k, (2.16)

∂(αkρkhk)

∂t + ∇ · (α_kρ_ku_kh_k) = −∇ · (α_kq_k) + Q_k, (2.17) where k denotes phases; Γk is the total mass source; Sk is the total momentum source; Qkis the total energy source.

2.1.3 Solid heat transfer modeling

The following equation is used to model the heat transfer in the solid region ρscp,s

∂Ts

∂t = ∇ · (κs∇Ts) + Qs, (2.18) where ρs, cp,sand κsare the solid density, heat capacity, and thermal conductivity, respectively; Qsis the volumetric heat source in the solid region. If Joule heating is used to heat up the solid region, Qsis modeled in the following way

Qs= (σe∇Ve) · ∇Ve, (2.19)

(39)

14

CHAPTER 2. TRANSPORTED FILM MODELING WITH CONJUGATE HEAT TRANSFER

where σeis the electrical conductivity; Veis the electric potential which is obtained by solving the following equation

∇ · (σe∇Ve) = 0, (2.20)

with a prescribed voltage difference, ∆Ve, across the heated length. If ∆Veis not documented and only the total heating power is available, ∆Veis calculated using an iterative procedure before solving Eq. (2.20).

2.2 Film-gas-core coupling

In a film-gas-core-coupled simulation, the film region is created based on the boundary of the gas core mesh, as shown in Fig. 2.2. This film mesh is always present regardless of whether a discretized film surface is wet or dry. Therefore, a coupling scheme that is suitable for both the pre- and post-dryout regions should be developed.

Figure 2.2: Grids for a film-gas-core-coupled simulation. The liquid film mesh is highlighted in blue. a, b and c are used to denote the discretized cells in the film region. Gas core cells that are adjacent to the film are highlighted in yellow.

2.2.1 Pre-dryout region coupling

In the pre-dryout region, there are three heat and mass transfer mechanisms, i.e., entrainment, deposition and evaporation, to model between the gas core and the

Mechanistic Modeling of Annular Flows at Near-Dryout Conditions

Mechanistic Modeling of Annular Flows at Near-Dryout Conditions

Mechanistic Modeling of Annular Flows at Near-Dryout Conditions

Abstract

Sammanfattning

List of Publications

Acknowledgments

Contents

List of Figures

List of Tables

List of Symbols

Acronyms

Latin symbols

Greek symbols

Operators

Superscripts

Subscripts

Chapter 1 Introduction

1.1 Empirical modeling

1.2 CFD modeling

1.3 Thesis structure

Chapter 2

Transported film modeling with conjugate heat transfer

2.1 Basic theories

2.2 Film-gas-core coupling