Document [original]

A three-dimensional two-phase model for flow,

transport and mass transfer processes in sewers

vorgelegt von

Master of Science

Katharina Teuber

bei der Fakultät VI - Planen Bauen Umwelt

der Technischen Universität Berlin

zur Erlangung des akademischen Grades

Doktorin der Ingenieurwissenschaften

- Dr.-Ing. -

genehmigte Dissertation

Promotionsausschuss:

Vorsitzender: Prof. Dr.-Ing. Matthias Barjenbruch

Gutachter: Prof. Dr.-Ing. Reinhard Hinkelmann

Gutachter: Prof. Dr. rer. nat. Gunnar Nützmann

Gutachter: Prof. Dr.-Ing. Dirk Muschalla

Gutachter: Prof. Dr.-Ing. Jann Strybny

Tag der wissenschaftlichen Aussprache: 01.10.2019

Berlin 2020

Preface

I am very grateful for having had the possibility of conducting this doctoral thesis at the

Chair of Water Resources Management and Modeling of Hydrosystems at TU Berlin within

the DFG funded Research Training Group “Urban Water Interfaces" (UWI).

I would like to thank all current and former members of the Chair for creating such a

great working environment and conducting high-level research at the same time: my su-

pervisor Reinhard (Phillip) Hinkelmann, Ralf Duda for both his technical support and his

kindness in every respect. My friends and colleagues Ilhan Özgen-Xian and Elena Matta

who often took the part of being my mentors. And Tabea Broecker for being a good friend,

sharing an office with me and being a peer in using OpenFOAM.

My thanks also go to the students who supported me during the course of the disserta-

tion: Nicholas R. Bowsher, who was such a great support in the end of this thesis, as well

as Waldemar Elsesser, Gerie Goitom and Shibashish D. Jaydev.

I am also very grateful for the members of the Research Training Group: Gunnar Nütz-

mann and Matthias Barjenbruch for valuable feedback and always offering a new perspec-

tive on my work; to Tosca Piotrowski and Gwendolin Porst for the great support in admin-

istrational and interdisciplinary work. My fellow doctoral students in UWI: Clara Romero,

Mikael Gillefalk, Robert Ladwig and Daneish Despot, thanks for the support and so many

nice social activities.

I am very thankful for the great support I received from my parents and brothers, not

only during the time of my dissertation but also all the way along my education. Thank

you for supporting me in every decision.

With all my love I am thankful to my partner Manuel and our son Leander for all the

support and understanding. Thank you for showing me what is truly important in life.

Berlin, April 16, 2019

Publications of cumulative doctoral thesis

First author publications (in chronological order):

1. Teuber, K., Broecker, T., Bayón, A., Nützmann, G. & Hinkelmann, R.: CFD-modelling

of free-surface flows in closed conduits, Progress in Computational Fluid Dynamics,

19 (6), 368-380, 2019; doi:10.1504/PCFD.2019.103266. Postprint.

2. Teuber, K., Broecker, T., Bentzen, T.R., Stephan, D., Nützmann, G. & Hinkelmann,

R.: Using computational fluid dynamics to describe H2S mass transfer across the

water-air interface in sewers, Water Science and Technology, 79 (10), 1934-1946, 2019;

doi:10.2166/wst.2019.193. Postprint.

3. Teuber, K., Broecker, T., Nützmann, G. & Hinkelmann, R.: CFD simulation of H2S

mass transfer under turbulent conditions in a stirring tank. In preparation for submis-

sion to Water Science and Technology (IWA Publishing).

Book chapter:

1. Teuber, K., Broecker, T., Jaydev, S.D., Goitom, G.M., Sielaff (née Grüneberger), M.,

Despot, D., Stephan, D., Barjenbruch, M. & Hinkelmann, R.: Multiphase CFD-Simulation

of Transport Phenomena in Sewer Systems, in: Mannina G. (ed) New Trends in Urban

Drainage Modelling. UDM 2018. Green Energy and Technology, 584-853, Springer,

Cham., 2019; 10.1007/978-3-319-99867-1_146. Postprint.

Conference papers (in chronological order):

1. Teuber, K., Broecker, T., Barjenbruch, M. & Hinkelmann, R.: High-resolution numeri-

cal analysis of flow over a ground sill using OpenFOAM, Proceedings of the 12th In-

ternational Conference on Hydroscience & Engineering (ICHE), Tainan, Taiwan, 2016.

Postprint.

2. Teuber, K., Sielaff (née Grüneberger), M., Despot, D., Stephan, D., Barjenbruch, M. &

Hinkelmann, R.: Modeling and Measuring of Interfaces in Sewer Systems, Proceedings

of the 37th IAHR (International Association for Hydro-Environment Engineering and

Research) World Congress, Kuala Lumpur, Malaysia, 2017. Postprint.

Supplementary contributions

Co-authored publication:

1. Broecker, T., Elsesser, W., Teuber, K., Özgen, I., Nützmann, G. & Hinkelmann, R.:

High-resolution simulation of free-surface flow and tracer transport over streambeds

with ripples, Limnologica, 68: 46-58, 2018.

Book chapter:

1. Broecker, T., Teuber, K., Elsesser, W. & Hinkelmann, R.: Multiphase Modeling of

Hydrosystems Using OpenFOAM, in: Gourbesville P., Cunge J., Caignaert G. (eds)

Advances in Hydroinformatics, Springer Water, 1013-1029, Springer, Singapore, 2018.

Conference poster:

1. Dixit, A., Teuber, K., Barjenbruch, M., Stephan, D. & Hinkelmann, R.: Extension of a

3D two-phase flow model to multicomponent reactive transport for odour and corro-

sion control in sewer systems, Workshop on Applications of Multi-scale Approaches

in Environmental Chemistry (AMARE), Rennes, France, 2019.

Presentation:

1. Teuber, K., Broecker, T., Elsesser, W. & Hinkelmann, R.: Beyond shallow water flow:

Navier-Stokes simulations with OpenFOAM, BIMoS Day: Shallow Water Flow Simula-

tions, Berlin International Graduate School in Model and Simulation based Research,

TU Berlin, Berlin, Germany, 2017.

An overview of all supplementary scientific work is given in Chapter 8.

This thesis was carried out as project T3 of the Research Training Group “Urban Water In-

terfaces (UWI)" (GRK 2032/1), which is funded by the German Research Foundation (DFG).

Abstract

Sewer networks are one major pillar of modern cities’ infrastructure. Their functionality

ensures the transport of wastewater to the sewage treatment plant and the transport of rain-

water from residential areas. Damages to sewers cause infiltration and exfiltration and at

the same time high costs for rehabilitation. The formation of hydrogen sulphide (H2S) rep-

resents a risk factor for the conditions of concrete channels. Its emission cannot only cause

the destruction of sewer walls by concrete corrosion, but can also represent a safety risk

for sewer workers. Within the last decades, the characteristics of H2S emissions were inten-

sively investigated and various models for predicting odour and corrosion were developed.

The current state of the art are one-dimensional model approaches. At the same time, some

predominant processes, e.g. the flow velocities in the air phase, are three-dimensional, and

H2S emissions are very relevant on locations with high turbulence and complex flow fields

(e.g. drops).

This work continues at this point. It investigates and extends a three-dimensional two-

phase model with regard to different aspects. For this purpose the two-phase solver in-

terFoam of the software OpenFOAM is used. Initially, the hydrodynamic properties for

different models in closed conduits are investigated by analysing hydrodynamic proper-

ties for different models in closed cross sections. The analysis begins with the simulation

of a simple single-phase water flow over a ground sill and is then extended to a highly

complex sewer geometry. The complex sewer network geometry is compared with results

of a 1:20 scale model and existing CFD simulations for an open geometry. The results

show a good agreement. Extensions are based on the description of mass transfer using

the Henry coefficient. Furthermore, adjustments are made to improve the specifics of H2S

emissions in sewers. These include the description of the temperature dependency of the

Henry coefficient, the equilibrium between H2S and the bisulphide ion (HS-) in the water

phase and the influence of the pH value on this equilibrium. An additional extension de-

scribes the concentration of H2S in the air phase as partial pressure. The extensions and

adaptations are validated using different analytical examples and the advantages of using a

three-dimensional model over a one-dimensional approach are demonstrated using the ex-

ample of the complex sewer geometry. Finally, the extended solver is coupled with a solver

for dynamic geometries to validate the simulated mass transfer under turbulent conditions.

The comparison of simulation results for mass transfer in a stirring tank with different stir-

ring rates leads to a good agreement with experimental results from laboratory experiments.

This work results in two new solvers, the difference of which lies in the geometry to

be described. The first solver can be applied to static meshes, while the second solver can

describe dynamic meshes, such as rotating geometries.

Zusammenfassung

Kanalnetze stellen eine wichtige Säule in der Infrastruktur moderner Städte dar. Ihre Funk-

tionsfähigkeit sichert den Transport von Abwasser zur Kläranlage und den Transport von

Regenwasser aus Siedlungsgebieten. Schäden an Kanälen können zu In- und Exfiltrationen

und gleichzeitig hohen Instandsetzungskosten führen. Einen Risikofaktor für den Zustand

von Betonkanälen stellen Umwandlungsprozesse von Schwefelwasserstoff (H2S) dar. Ihre

Emission kann nicht nur Kanalwände durch Betonkorrosion zerstören, sondern auch ein

Sicherheitsrisiko für Kanalarbeiter*innen darstellen. Innerhalb der letzten Dekaden wurden

H2S Emissionen intensiv erforscht und verschiedenste Modelle zur Vorhersage von Geruch

und Korrosion entwickelt. Der aktuelle Stand der Technik sind eindimensionalle Modellan-

sätze. Gleichzeitig sind einige vorherrschende Prozesse, beispielsweise die Fließgeschwin-

digkeiten in der Luftphase, dreidimensional und H2S Emissionen spielen eine besondere

Rolle an Stellen mit hoher Turbulenz und komplexen Strömungsfeldern (z.B. Abstürze).

Diese Arbeit setzt an diesem Punkt an und untersucht und erweitert ein dreidimensiona-

les Zweiphasenmodell hinsichtlich unterschiedlicher Aspekte. Hierfür wird der Zweiphasen-

Löser interFoam der Software OpenFOAM verwendet. In einem ersten Schritt werden die

hydrodynamischen Eigenschaften für unterschiedliche Modelle in geschlossenen Querschnit-

ten untersucht. Die Testfälle beginnen bei der einfachen Simulation einer Einphasen-Wasser-

strömung über eine Schwelle und werden dann hin zu einer hochkomplexen Kanalnetz-

geometrie erweitert. Die Ergebnisse der Simulationen werden mit experimentellen Ergeb-

nissen und analytischen Lösungen verglichen. Die komplexe Kanalnetzgeometrie wird mit

Ergebnissen eines im Verhältnis 1:20 gebauten Modellversuchs und mit existierenden CFD-

Simulationen für eine offene Geometrie verglichen, und die Ergebnisse zeigen eine gute

Übereinstimmung. Es werden Transport-und Massentransferprozesse anhand verschiede-

ner Beispiele untersucht, wozu vorhandene Erweiterungen des Lösers verwendet werden.

Diese Erweiterungen basieren auf der Beschreibung von Massentransfer mittels des Henry-

Koeffizienten. Es werden Anpassungen vorgenommen, um H2S Emissionen im Kanal bes-

ser beschreiben zu können. Diese umfassen die Beschreibung der Temperaturabhängigkeit

des Henry-Koeffizienten sowie das Gleichgewicht zwischen H2S und dem Hydrogensulfid-

Anion (HS-) in der Wasserphase und den Einfluss des pH Wertes auf dieses Gleichgewicht.

Eine zusätzliche Erweiterung beschreibt die Konzentration von H2S in der Luftphase als

Partialdruck. Die Erweiterungen und Anpassungen werden anhand unterschiedlicher ana-

lytischer Beispiele validiert, und es werden die Vorteile der Anwendung eines dreidimensio-

nalen Modells gegenüber eines eindimensionalen Ansatzes wird am Beispiel der komplexen

Kanalgeometrie gezeigt. Zuletzt wird der erweiterte Löser mit einem Löser für dynamische

Geometrien gekoppelt, um den simulierten Massentransfer unter turbulenten Bedingun-

gen zu validieren. Der Vergleich von Simulationsergebnissen für Massentransfer in einem

Rührbehälter mit unterschiedlichen Rührgeschwindigkeiten führt zu einer guten Überein-

stimmung mit experimentellen Ergebnissen aus Laborversuchen.

Insgesamt resultieren aus der Arbeit zwei neue Löser, deren Unterschied in der zu be-

schreibenden Geometrie liegt. Der erste Löser kann im Bereich von statischen Netzen an-

gewandt werden, während der zweite Löser dynamische Geometrien, wie zum Beispiel die

Bewegung von Rotoren, beschreiben kann.

Contents

1 Introduction 1

1.1 Odour and corrosion in sewers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Scientificbackground.................................. 5

1.2.1 Preliminary model developments for odour and corrosion . . . . . . . . 5

1.2.2 Three-dimensional two-phase modelling using OpenFOAM . . . . . . . 8

1.2.3 Transport and mass transfer modelling . . . . . . . . . . . . . . . . . . . 12

1.2.3.1 Transport............................... 12

1.2.3.2 Masstransfer............................. 14

1.3 Scopeofthisthesis ................................... 16

2 Research on interfaces in sewer systems within the DFG Research Training Group

“Urban Water Interfaces" 19

2.1 Introduction....................................... 20

2.2 Pilotplantinvestigations................................ 22

2.2.1 Sulfide formation in pressure or rising mains . . . . . . . . . . . . . . . 22

2.2.2 Elevated H2S emissions in gravity sewers due to flow interruptions

caused by the accumulation of gross solids (debris / solid materials) . 23

2.2.3 Investigation of H2S emissions at drop structures of sewer systems . . . 23

2.2.4 Impact of countermeasures at hotspots . . . . . . . . . . . . . . . . . . . 24

2.3 Investigations of corrosion processes in sewer and laboratory tests . . . . . . . 25

2.3.1 Concrete and mortar types . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.3.2 Pilot plant acid resistance tests . . . . . . . . . . . . . . . . . . . . . . . . 26

2.3.3 Laboratory performance tests . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.3.4 Analysis of concrete and mortar specimens . . . . . . . . . . . . . . . . . 26

2.4 Development of a three-phase CFD simulation model . . . . . . . . . . . . . . . 27

2.4.1 Numericalmodel................................ 27

2.4.2 Validation of the water phase behavior . . . . . . . . . . . . . . . . . . . 29

2.5 Conclusions ....................................... 31

3 Validation of hydrodynamic two-phase simulations 33

3.1 Introduction....................................... 33

3.2 Methodsandmaterials................................. 35

3.2.1 Geometryandmesh .............................. 35

3.2.2 Numericalmodel................................ 35

3.2.3 Boundaryconditions.............................. 37

3.3 Applications....................................... 38

3.3.1 Single phase water flow over a ground sill . . . . . . . . . . . . . . . . . 38

3.3.2 Free surface flow over a hill . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.3.3 Complex free surface flow in a sewer model . . . . . . . . . . . . . . . . 46

3.4 Conclusions ....................................... 51

4 Parameter study on hydrodynamic two-phase simulations of flow over a ground

sill 52

4.1 Introduction....................................... 52

4.2 Computational framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.2.1 Numericalmodel................................ 53

4.2.2 Turbulence modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

4.2.3 Boundaryconditions.............................. 54

4.2.4 Geometryandmesh .............................. 54

4.3 Single-phaseflow.................................... 55

4.4 Two-phaseflow..................................... 57

4.5 Conclusions ....................................... 60

5 Single-phase transport simulations 61

5.1 Introduction....................................... 61

5.2 Materialandmethods ................................. 62

5.2.1 Hydrodynamic simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.2.2 Transport simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

5.2.3 Masstransfer .................................. 63

5.3 Resultsanddiscussion................................. 63

5.3.1 Spreading of hydrogen sulphide in headspace of sewer pilot plant (air

phase) ...................................... 63

5.3.2 Tracer transport in sewer pipe (water phase) . . . . . . . . . . . . . . . . 66

5.4 Conclusions ....................................... 68

6 Mass transfer simulations under equilibrium conditions: Validation and solver

extension 69

6.1 Abstract ......................................... 69

6.2 Introduction....................................... 70

6.3 Methods ......................................... 71

6.3.1 Numericalmodel................................ 71

6.3.2 Henrycoefficient................................ 74

6.3.3 Extensions.................................... 74

6.3.4 Casestudies................................... 77

6.4 Resultsanddiscussion................................. 78

6.4.1 Saturation of H2Sinatank .......................... 78

6.4.2 Mass transfer in a rectangular channel . . . . . . . . . . . . . . . . . . . 82

6.4.3 Complex sewer geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

6.5 Conclusions ....................................... 86

6.6 Acknowledgements................................... 86

7 Mass transfer under turbulent conditions 87

7.1 Abstract ......................................... 87

7.2 Introduction....................................... 87

7.3 Methods ......................................... 89

7.3.1 Geometryandmesh .............................. 89

7.3.2 Numericalmodel................................ 89

7.3.3 Turbulence.................................... 92

7.3.4 Casestudy.................................... 93

7.3.5 Initial and boundary conditions . . . . . . . . . . . . . . . . . . . . . . . 94

7.4 Resultsanddiscussion................................. 95

7.4.1 Graphicanalysis ................................ 95

7.4.2 Sensitivityanalysis............................... 96

7.4.3 Quantitative analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

7.5 Conclusions .......................................103

8 Supplementary contributions 104

8.1 Advantages of three-dimensional flow simulations . . . . . . . . . . . . . . . . 104

8.2 Possibilities extending the interH2SFoam solver to multicomponent reactive

transport.........................................105

8.3 Research on further multiphase CFD applications carried out at the Chair of

Water Resources Management and Modeling of Hydrosystems, TU Berlin . . . 106

8.4 Application of free-surface flow and transport simulations to natural river beds107

9 Synthesis 108

9.1 Conclusions .......................................108

9.1.1 Generaloutcomes ...............................108

9.1.2 Outcomes of hydrodynamic modelling . . . . . . . . . . . . . . . . . . . 109

9.1.3 Outcomes of single-phase transport modelling . . . . . . . . . . . . . . . 110

9.1.4 Outcomes of mass transfer modelling . . . . . . . . . . . . . . . . . . . . 111

9.1.5 Finalnotes....................................113

9.2 Outlook..........................................113

Appendix A Flow charts of solver extensions 116

A.1 Integration into the solution procedure . . . . . . . . . . . . . . . . . . . . . . . 117

A.1.1 interH2SFoam solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

A.1.2 interDyMH2SFoam solver . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

A.2 Code extensions and integration into the object-oriented framework . . . . . . 121

A.2.1 Filestructure ..................................121

A.2.2 Make/files....................................126

A.2.3 createFields.H..................................127

A.2.4 CEqn.H .....................................129

A.2.5 interH2SFoam.C ................................130

A.2.6 Integration of interDyMH2SFoam . . . . . . . . . . . . . . . . . . . . . . 133

A.2.7 Passive scalar transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

A.3 Function objects and boundary conditions . . . . . . . . . . . . . . . . . . . . . 135

A.3.1 Equilibrium conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

A.3.2 Calculation of partial pressure in ppm . . . . . . . . . . . . . . . . . . . 136

Appendix B Test case overview 137

B.1 Single-phase flow over a ground sill . . . . . . . . . . . . . . . . . . . . . . . . . 138

B.2 Two-phase flow over a two-dimensional hill . . . . . . . . . . . . . . . . . . . . 140

B.3 Complexsewergeometry ...............................142

B.4 Single-phase transport in rectangular pipe . . . . . . . . . . . . . . . . . . . . . 144

B.5 Concrete probes in sewer pilot plant . . . . . . . . . . . . . . . . . . . . . . . . . 146

B.6 Quasi one-dimensional cubic tank . . . . . . . . . . . . . . . . . . . . . . . . . . 148

B.7 Mass transfer in rectangular duct . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

B.8 Stirringtank.......................................152

Bibliography 154

List of Figures

1.1 Sewer processes (following Hvitved-Jacobsen et al. (2013)) . . . . . . . . . . . . 3

1.2 Connection between rising and gravity main (Jiang et al., 2015) . . . . . . . . . 4

1.3 Previous model developments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.4 Three-dimensional air phase velocities [m/s] in circular sewer headspace (Edwini-

BonsuandSteffler,2004)................................ 7

1.5 Transport processes, Schankat (2010) after Barlag et al. (1998) . . . . . . . . . . 13

1.6 Overview over different chapters . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.1 Interface processes in sewer systems. . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.2 Researchpilotplant. .................................. 22

2.3 Schematic drawing of the research pilot plant. . . . . . . . . . . . . . . . . . . . 22

2.4 Experiment 2 – Elevated H2S emissions due to gross solids in wastewater. . . . 23

2.5 Experiment 3 – Pilot plant scale investigation of H2S emissions at drop struc-

tures. ........................................... 24

2.6 Model domain of water surface validation case. . . . . . . . . . . . . . . . . . . 30

3.1 Domain of single-phase, two-dimensional hill flow including measurement

locations, three selected indicator locations for grid convergence study and

computational mesh (Davroux et al., 1995). . . . . . . . . . . . . . . . . . . . . . 39

3.2 Mesh sensitivity analysis, (a) Comparison of simulated indicator values φ

(flow velocities) for different grid sizes in three points within the domain; (b)

RMSE between different refinement steps over all 15 indicator values . . . . . 41

3.3 Comparison of simulated velocities Uxusing different turbulence models

with experimental data for flow over a two-dimensional hill, (a) Uxveloci-

ties at x−03 = -0.30 m, (b) Uxvelocities at x00 = 0.00 m, (c) Uxvelocities at x01

= 0.03 m (d) Uxvelocities at x02 =0.05m ...................... 42

3.4 Comparison of simulated velocities Uyusing different turbulence models with

experimental data for flow over a two-dimensional hill, (a) Uyvelocities at x00

= 0.00 m, (b) Uyvelocities at x03 =0.03m ...................... 43

3.5 Free surface flow over a hill - Model domain of case 1 . . . . . . . . . . . . . . 45

3.6 Complex sewer geometry - Geometry of the sewer stretch (Bayón et al., 2015) 48

3.7 Complex sewer geometry - Overview of monitoring points for water surface

profiles.......................................... 49

3.8 Complex sewer geometry - Comparison of water surface profile at different

locations (existing CFD model - results obtained by Bayón et al. (2015); new

CFD model - results obtained in this paper) . . . . . . . . . . . . . . . . . . . . . 50

4.1 Model setup and measurement locations of single-phase case. . . . . . . . . . . 54

4.2 Velocity profiles at x−03 =−0.30m(inlet profile). . . . . . . . . . . . . . . . . . . 56

4.3 Velocity profiles at x01 =0.03m(endofsill). .................... 56

4.4 Filling of the domain, segment of the computational domain for different time

steps. ........................................... 59

4.5 Comparison of water level drawdown for different sill structures . . . . . . . . 59

4.6 Comparison of water level drawdown for two- and three-dimensional model

setup ........................................... 59

5.1 Tracer transport in water phase of rectangular pipe with non-physical spread-

ingacrosswatersurface................................. 62

5.2 Top view of the model domain (pipeline containing the concrete probes). . . . 64

5.3 Resulting tracer distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.4 Movement and spreading of tracer along sewer. . . . . . . . . . . . . . . . . . . 66

5.5 Movement and spreading of tracer along sewer. . . . . . . . . . . . . . . . . . . 67

6.1 H2S saturation in a tank (left: phase fraction value, right: concentration pro-

files along the vertical axis over time). . . . . . . . . . . . . . . . . . . . . . . . . 79

6.2 H2S saturation in a tank for different temperatures (298.15K (cp. Figure 6.1,

right) and 288.15K) (left: phase fraction value, right: concentration profiles

alongtheverticalaxis).................................. 80

6.3 Application example for HS-and H2S equilibrium and partial pressure of air

phase concentration (left: Phase fraction value profile over domain height,

middle: concentration profile, right: partial pressure of gas-phase concentra-

tion). ........................................... 81

6.4 Mass transfer in rectangular channel for test 7 (see Table 6.1) (left: phase

fraction value, middle: velocity, right: concentration profiles). . . . . . . . . . . 82

6.5 Mass transfer in rectangular channel for test 21(see Table 6.1) (left: phase

fraction value, middle: velocity, right: concentration profiles). . . . . . . . . . . 83

6.6 Mass transfer simulations in complex sewer geometry at t=10s (a) overview

of the domain filled with water under steady-state conditions, b) flow veloc-

ities and water surface behaviour in hydraulic jump, c) tracer distribution in

hydraulic jump, d) top view on tracer distribution). . . . . . . . . . . . . . . . . 85

7.1 System sketch of stirring tank including relevant variables (Wu (1995), modi-

fied). ........................................... 88

7.2 Computational grid used for the simulations including boundaries (red: stir-

rer, grey: shaft, blue (outer cylinder): cylinder walls, blue (inner cylinder):

AMI, blue (upper plane): top). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

7.3 Vortex shapes for different stirring rates (left: N=0.8333 1/s, middle: N=1.667

1/s, right: N=2.333 1/s). ................................ 95

7.4 Influence of different turbulence models on the resulting mass transfer coef-

ficient (inundation ratio of h/R = 0.7). ........................ 97

7.5 Influence of stirrer inundation on the mass transfer coefficient. . . . . . . . . . 98

7.6 Influence of stirring rate and inundation ratio on mass transfer (Wu, 1995). . . 100

7.7 Influence of stirring rate on mass transfer (Carrera et al., 2017). . . . . . . . . . 101

7.8 Comparison of relation between Reynolds number and volumetric mass trans-

fer coefficient (Wu, 1995, Carrera et al., 2017). . . . . . . . . . . . . . . . . . . . . 102

9.1 Sewer processes (following Hvitved-Jacobsen et al. (2013)) . . . . . . . . . . . . 115

A.1 Flow chart of interFoam solver (following Devolder et al. (2015), Lopes et al.

(2017)) ..........................................117

A.2 Flow chart of interH2SFoam solver (following Devolder et al. (2015), Lopes

et al. (2017)), solver extensions compared to interFoam solver in green . . . . . 118

A.3 Flow chart of interDyMFoam solver (following Devolder et al. (2015), Lopes

etal.(2017)) .......................................119

A.4 Flow chart of interDyMH2SFoam solver (following Devolder et al. (2015),

Lopes et al. (2017)), solver extensions compared to interFoam solver in green . 120

A.5 File structure of interFoam solver . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

A.6 File structure of interH2SFoam solver . . . . . . . . . . . . . . . . . . . . . . . . 123

A.7 File structure of interDyMH2SFoam solver . . . . . . . . . . . . . . . . . . . . . 124

A.8 File structure of passiveScalarInterFoam solver . . . . . . . . . . . . . . . . . . . 125

List of Tables

2.1 Composition and assumed acid resistances of concrete types used for acid

resistance tests. Cement types are classified after DIN EN 197-1:2011-11. MS:

microsilica, FA: hard coal fly ash. . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.2 Testing parameters of the four laboratory performance tests. . . . . . . . . . . . 28

3.1 Water phase validation - RMSE/U0[-] between experimental results and sim-

ulation .......................................... 40

3.2 Water-air-interface - Properties of the different test cases . . . . . . . . . . . . . 45

4.1 Overview over different sill structures. . . . . . . . . . . . . . . . . . . . . . . . 57

4.2 Subcritical two-phase flow: Properties of different testcases. . . . . . . . . . . . 58

6.1 Mass transfer in rectangular channel: flow properties of analysed test cases. . 78

7.1 Fluid properties as defined in the CFD simulations. . . . . . . . . . . . . . . . . 93

7.2 Setup of stirring tank in CFD simulations and known metrics of experimental

setup by Carrera et al. (2017) and Wu (1995). . . . . . . . . . . . . . . . . . . . . 99

B.1 Model setup of ERCOFTAC test case . . . . . . . . . . . . . . . . . . . . . . . . . 138

B.2 Boundary conditions for hydrodynamic simulations of single-phase flow over

agroundsill.......................................139

B.3 Model setup, two-phase flow over a two-dimensional hill . . . . . . . . . . . . 140

B.4 Boundary conditions for hydrodynamic simulations of two-phase flow over a

hill ............................................141

B.5 Model setup of complex sewer geometry . . . . . . . . . . . . . . . . . . . . . . 142

B.6 Boundary conditions for simulations of complex sewer geometry . . . . . . . . 143

B.7 Model setup of single-phase transport in rectangular pipe test case . . . . . . . 144

B.8 Boundary conditions for simulations of single-phase transport in rectangular

pipe............................................145

B.9 Concrete probes in sewer pilot plant . . . . . . . . . . . . . . . . . . . . . . . . . 146

B.10 Boundary conditions for transport simulations around concrete probes . . . . 147

B.11 Model setup of quasi one-dimensional cubic tank . . . . . . . . . . . . . . . . . 148

B.12 Boundary conditions for simulations of quasi one-dimensional cubic tank . . . 149

B.13 Model setup of mass transfer in rectangular duct . . . . . . . . . . . . . . . . . 150

B.14 Boundary conditions for mass transfer in rectangular duct . . . . . . . . . . . . 151

B.15 Model setup of stirring tank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

B.16 Boundary conditions for stirring tank . . . . . . . . . . . . . . . . . . . . . . . . 153

xii

Chapter 1

Introduction

1.1 Odour and corrosion in sewers

Every year damages of sewers due to concrete corrosion caused by hydrogen sulphide (H2S)

result in high costs for sewer maintenance. In 1998, the restoration costs for corroded sew-

ers in Germany were estimated to be in a range of billions of US $ (Kaempfer and Berndt,

1999). Barjenbruch et al. (2008) state that Veolia Eau Germany employs 800 t/a for the coun-

termeasure Nutriox which is surpressing the production of H2S for around 200,000 €/a. At

the same time, high H2S concentrations in the sewer atmosphere are a health risk for sewer

workers. This is why three research projects within the DFG Research Training Group “Ur-

ban Water Interfaces" (UWI) focus on the enhanced understanding of odour and corrosion

processes as well as on the development of new prediction models. A detailed explanation

of the ongoing research will be given in the course of this thesis.

In general, odourous substances from sewers stem from multiple substances such as

sulphides, nitrogen containing compounds, acids, aldehydes and ketones. Depending on

the pH value in the wastewater, different compounds are emmitted from the sewer due to

different chemical states with varying pH values. Low pH value emissions are considered

as more intense which is why H2S is often considered as key parameter for odour nuisance

(Matsché et al., 2005).

There are several factors, why odour and corrosion due to hydrogen sulphide have

become more relevant in the last years and will strongly gain importance in future:

• Decreasing water consumption: Due to demographic change and an increased water

efficiency the water consumption decreases (Barjenbruch et al., 2008);

• Higher detention times: More centralized sewer systems cause higher detention times

of wastewater in the sewers (Barjenbruch et al., 2008). The current discussion on de-

centralized urban water systems could help decreasing the problem of long detention

times in the future (Leigh and Lee, 2019);

• Climate change: Warmer climates reduce the concentration of dissolved oxygen (DO),

which enhances more rapid depletion of oxygen in the sewers and stimulates the

generation of H2S by the sulfate reducing bacteria (SRB) (Rootsey and Yuan, 2010).

These aspects emphasize why it is important to better understand and to be able to estimate

transformation processes in sewers. In the past 70 years several models have been developed

in order to predict the occurrence of odour and corrosion. They will be outlined in Chapter

Chapter 1: Introduction

1.2.1 but it is important to note that the state-of-the-art models are one-dimensional models,

neglecting three-dimensional effects of flow velocities, non-uniform flow and highly turbu-

lent effects (e.g. most hydraulic structures). This is why in this thesis, a three-dimensional

model was developed to better describe local effects of H2S emissions, to provide advice

for the optimization of H2S hotspots and -as a long-term goal- to provide parametrizations

for the improvement of existing model approaches. To be able to better understand the im-

portance of model developments in the field of H2S emissions, first a general introduction

regarding relevant processes, locations with high H2S emissions and mitigation strategies is

given.

Figure 1.1 shows an overview of the relevant processes. For the work of this thesis and

the different projects working on interfaces in sewer systems in UWI, three main interfaces

are highly relevant and of main interest. These are namely:

• The biofilm - (waste)water interface;

• The (waste)water - (sewer) air interface;

• The (sewer) air - biofilm - (concrete) wall interface.

At the biofilm - (waste)water interface sulphate-reducing bacteria (SRB), which are located

in the biofilm, can reduce sulphate in the wastewater. This process can take place under

anaerobic conditions in the wastewater. In the water phase, an equilibrium between H2S

and bisulphide ion (HS-) exist, together they are described as the total dissolved sulphide.

The pH value and the temperature impact this equilibrium. At the (waste)water - (sewer)

air interface, H2S can be emitted from the water phase into the air phase, depending on

factors such as the pH value, temperature and the concentration of oxygen and nitrate.

H2S is classified to be a volatile compound, which defines a compound that has a high

vapor pressure at normal temperatures causing molecules to change phases, in this case

considered as emission from the water into the air phase. For these compounds, the water-

air equilibrium can be described by the Henry coefficient, which is a key parameter in this

thesis. It describes the state at which the transfer rates between the two phases are equal

and thus no net transfer between the phases occurs (Hvitved-Jacobsen et al., 2013). The

Henry coefficient itself depends on the overall temperature in the domain. In the literature,

many different ways to describe this equilibrium exist. An overview and a conversion table

are given in Sander (2015). One way to describe the temperature dependency of the Henry

coefficient is by using the van’t Hoff equation (Hvitved-Jacobsen et al., 2013).

The mass transfer between the phases depends on factors, such as the concentration of

H2S in the water. The concentration is influenced by the pH value and hydrodynamic pa-

rameters, which enhance the transport of H2S in the sewer. The transport is furthermore

influenced by the turbulence of the phases. The higher the turbulence in a system, the

higher the turbulent diffusion, the faster the mass transfer and therefore the faster equilib-

rium conditions develop. Due to the strong effect of turbulence on the local mass transfer,

it is important to account for the effects properly.

At the (waste)water - (sewer) air interface, the transer of oxygen, called reaeration, is

another relevant process. Reaeration is the only way to aerate the water phase and thus

influences the potential of aerobic and anaerobic processes in the wastewater (Hvitved-

Jacobsen et al., 2013).

Chapter 1: Introduction

Moving on to the third interface, the (sewer) air - biofilm - (concrete) wall interface, H2S

can cause the corrosion of concrete or cement-bound materials depending on the humid-

ity of the pipe walls. All processes described above highly depend on further factors, i.e.

the activity of SRB strongly depends on factors such as the availability of oxygen, redox

potentials, temperature, retention time etc., which cannot be addressed as a whole in the

framework of this thesis. A complete overview can be found in Barjenbruch et al. (2008).

Figure 1.1: Sewer processes (following Hvitved-Jacobsen et al. (2013))

Apart from the underlying processes, it is important to fully understand the different

designs of sewers and the influence of the design on H2S formation and emission. The

design and age of the city, as well as its topography, determine, which type of sewer is used.

These sewer networks can be designed in very different ways. Hvitved-Jacobsen et al. (2013)

propose three different ways to classify sewers:

1. Based on the type of sewage collected (sanitary sewers, storm sewers and combined

sewers);

2. Based on the transport model applied (gravity sewers and pressure pipes);

3. Based on the size and function.

Wastewater is generally collected in sanitary and combined sewers. In a gravity sewer, the

flow is driven by the slope of the pipe and has a free water surface whereas in a pressure

sewer the flow is mainly driven by a pump and only water is present in the pipe. In a

Chapter 1: Introduction

gravity sewer, mass transfer processes across the wastewater-air interface (H2S emissions,

reaeration) are enabled due to the exposure to the sewer atmosphere. In pressure pipes,

oxygen within the pipe can be consumed leading to anaerobic conditions which enable the

reduction of sulphate.

All of these design classes influence the processes in the sewer and therefore the possible

occurrence of odour and corrosion. Due to the variety of influence factors, odour and corro-

sion can develop in many different places in a sewer network, however, there are some areas

which are more likely to suffer these problems. As reported in Barjenbruch (2003), one par-

ticularly problematic location are connection shafts between rising mains and gravity sewers

with long detention times. These structures enable a direct contact between wastewater from

anaerobic conditions and the surrounding air and are in general locations of higher turbu-

lence. Figure 1.2 illustrates such a connection. Nevertheless, it should be pointed out that

also pure gravity sewers could suffer from odour and corrosion when anaerobic conditions

prevail due to long retention times, low slopes and therefore low flow velocities.

As this short overview shows, the design of a sewer is a major factor in the possible

occurrence of odour and corrosion problems. One main aspect is the occurrence of anaerobic

conditions, the other is the level of turbulence further downstream which enhances H2S

emissions.

Figure 1.2: Connection between rising and gravity main (Jiang et al., 2015)

A number of countermeasures are available to mitigate the development of odour and

corrosion. There are chemical solutions on the one hand, which include the enrichment

with oxidants, precipitation of sulphide and pH-regulation techniques. On the other hand,

physical measures, such as gas treatment with biofilters, covering systems and constructive

solutions, exist (Barjenbruch et al., 2008). The range of these countermeasures shows that

a model should be able to account for the countermeasures in order to become a good

planning tool. The physical countermeasures, such as covering systems, could be accounted

for by a purely physical model but chemical countermeasures would need to be modelled

using reactive transport modelling. The long-term goal of the model developed in this thesis

is to be able to model countermeasures such as ventilation, dosages and flushing.

Chapter 1: Introduction

1.2 Scientific background

1.2.1 Preliminary model developments for odour and corrosion

In this Chapter, a short overview over existing model approaches is given. Their strenghts

and weaknesses are outlined and the possible advantages of a new model are highlighted.

Within the last 50 years, odour and corrosion in sewer systems caused by H2S has been

intensively analysed by different researchers. The main milestones are illustrated in Figure

1.3. Vollertsen et al. (2005) give a short overview about the history of research on hydrogen

sulphide:

Figure 1.3: Previous model developments

In the 1970s, research on H2S was intensified in the USA and Australia, and empirical

models were developed by Pomeroy and Parkhurst (1978) as well as Thistlethwayte (1972).

From the beginning of the 1990s to the late 2000s, the Danish research group around

Thorkhild Hvitved-Jacobsen has made major contributions to developing a conceptual un-

derstanding of major processes during H2S transformations: In the 1990s, research focussed

on an conceptual understanding of organic matter transformations in sewers. A concep-

tual model inspired by the Activated Sludge Model (ASM) called WATS (Wastewater Aer-

obic/Anaerobic Transformations in Sewers) has been developed by Bjerre et al. (1998). The

ASM is a one-dimensional model to describe biochemical transformations solving transport-

reaction or mass-balance equation for different components based on the description of dif-

ferent interacting processes. This model is mainly used for the description of processes in

wastewater treatment plants. Due to several differences between a sewer and a wastewa-

ter treatment plant, the applicability of the model is limited (Barjenbruch et al., 2008). In

the late 1990s, research on anaerobic transformations was performed with a focus on trans-

formations in pressure mains and gravity sewers under alternating aerobic and anaerobic

conditions. Especially H2S and anaerobic transformations of organic matter were investi-

gated. In the early 2000s, the release of H2S from bulk water into the sewer gas phase, the

Chapter 1: Introduction

oxidation of hydrogen sulphide from the biofilm and bulk water and air injections into pres-

sure pipes were investigated. The latest version of the WATS model includes the following

functionalities (HV-Consult, 2019):

• Hydrodynamics: Description of gas and water flow along the sewer, ventilation of

sewer gas to atmosphere;

• Water-air interface: Transfer of oxygen, H2S, carbon dioxide, mercaptanes, others;

• Air phase: Oxidation of H2S on moist concrete walls, concrete corrosion;

• Water phase: Tranformation in bulk water, biofilms, sediments of organic matter, sul-

furous compounds, specific organic compounds, oxygen, nitrate, nitrite; precipitation

of H2S by iron; chemical oxidation of H2S by strong oxidizing agents; pH and buffer

strength of wastewater;

• Focus on dry weather conditions.

The model’s focus lies on the description of dry weather conditions: Furthermore, uniform

flow is being considered under steady-state conditions. The gas phase velocity is estimated

as a proportion of the water phase velocity (less than 35 to 50 %), dispersion of transported

substances is not considered (Hvitved-Jacobsen et al., 2013). Carrera et al. (2016) give an

overview over the capabilities of existing model approaches and discuss the ability of the

WATS model to describe turbulent effects. Different simplified expressions have been used

in different applications to link H2S emissions to O2emissions or to flow properties in the

pipe. Furthermore, an empirical approach exists to quantify H2S emissions in drop struc-

tures (Matias et al., 2017).

In late 2008, a research project was launched in Australia. It was funded by the major wa-

ter utilities in Australia in cooperation with the Australian Research Council. The so-called

Sewer & Corrosion & Odour Research (SCORe) project ran for five years with a total budget

around $20 million. The project focussed on different research areas, which are namely cor-

rosion processes, gas phase technologies, liquid phase control and knowledge management

(Rootsey et al., 2012). A new model named SeweX has been developed. In this model the

mathematical model of Hvitved-Jacobsen for predicting sulphide generation in sewers has

been improved in order to predict spatial and temporal variations in sulphide concentrations

as well as other parameters. The model has also been linked to sewer hydraulic models to

predict dynamic changes in sulphur compounds within the sewer system due to changing

sewer characteristics such as diurnal variations (Rootsey et al., 2012). In SeweX, anaerobic

and aerobic carbon and sulphur transformation processes in a rising main are described

using the model developed by Freudenthal et al. (2005): instead of the ASM, the Anaerobic

Digestion Model No. 1 (ADM1) is used in order to account for anaerobic fermentation pro-

cesses happening downstream of a sewer pipe (Sharma et al., 2008b). Regarding corrosion

and odour control methods, the SeweX model can be used to predict the performance of

different chemicals and various dosing locations (Rootsey et al., 2012, Sharma et al., 2008a).

Ventilation is an important factor when analysing odour and corrosion processes in sewer

systems. It can be used as a measure against unwanted processes in two different ways. The

most common approach is to change the air sufficiently to maintain dry sewer structures at

all times and to minimise hydrogen sulphide buildup in the sewer air. But it is also pos-

sible to maintain zero relative velocity between wastewater and ventilating air to minimise

the rate of H2S emission and evaporation from the wastewater surface. However, it has to

be emphasized that the second approach is not used as a practical measure against H2S

Chapter 1: Introduction

emissions. Within the SeweX project, a new algorithm has been developed to describe air

movement (Rootsey et al., 2012). The algorithm describes ventilation from a viewpoint of

conservation of momentum. A force balance is used to describe the influences bearing on

ventilation. The influences accounted for are pressure differences, gravitational forces, drag

forces at the air/water interface and friction forces at the air/pipe interface (Ward et al.,

2011).

After analysing the two main models that have been developed in recent years, different

deficits can be identified. One main disadvantage of the WATS model has been removed in

the SeweX model: the fact that it only accounts for stationary flow conditions. Nevertheless,

both models are one-dimensional models. Further deficits are the lack of possibilities to

describe non-uniform flow, turbulent flow, and the simplifications when describing highly

turbulent locations (e.g. drops), gas flow velocities and transport phenomena (Barjenbruch

et al., 2008). Another deficit of the models is their availability. Both the WATS and the

SeweX model are not public domain. These advantages and disadvantages show that the

existing models are a suitable tool to analyse long sewer networks under the assumption of

steady flow conditions. They can account for typical H2S hotspots such as drop structures

by parametrizations but lack a detailled view on specific flow and emission characteristics

on a local scale. A three-dimensional model accounting for non-uniform flow conditions,

transient flows as well as highly turbulent flow characteristics can help improving the de-

sign of H2S hotspots such as connecting shafts between rising mains and gravity sewers.

Furthermore, different researchers concerned with sewer processes are describing three-

dimensional effects in sewers. Vollertsen et al. (2008) mention, that corrosion processes start

close to the water surface and develops further in direction of the sewer crown. Edwini-

Bonsu and Steffler (2006) show secondary flows in the air phase that do not move in the main

flow direction (Figure 1.4). These observations show that processes occurring in a circular

sewer are three-dimensional. The error caused by simplifications due to one-dimensional

models have not been analysed in publications yet. A three-dimensional model could fur-

thermore analyse the effects of physical or chemical countermeasures in detail. To begin

with, physical countermeasures could be analysed with an implementation without reactive

transport. With reactive transport included, the model could analyse the effect of chemical

countermeasures.

Figure 1.4: Three-dimensional air phase velocities [m/s] in circular sewer headspace

(Edwini-Bonsu and Steffler, 2004)

Chapter 1: Introduction

1.2.2 Three-dimensional two-phase modelling using OpenFOAM

There are many different factors which influence H2S formations in a sewer, ranging from

hydrodynamic conditions, such as the flow velocities and turbulence rate, to biochemi-

cal factors, such as the pH value. Numerical models and Computational Fluid Dynamics

(CFD) software products offer the possibility of describing even complex hydrodynamics

for single- or multiphase systems as well as transport processes under various conditions.

Their application ranges from natural systems, such as reservoirs with a scale of kilometers,

to technical systems, such as combined sewer overflows or even biofilters with a scale of

millimeters. Numerical modelling can therefore be considered as the simulation of single-

or multiphase systems in response to certain forcing conditions. They can help in identify-

ing the influence of different factors and future needs (Matta, 2018). In the context of sewer

processes, models can help estimate pollutant load as well as the influence of hydrodynam-

ics on the transport and mass transfer of these loads. They can support deciding for correct

countermeasures ranging from construction measures to dosing strategies. “The processes

occurring in a water body are described by mathematical models, often a set of coupled,

non-linear, partial differential equations, which are known as the Navier-Stokes equations

for momentum and the continuity equation (Ji, 2008, Hinkelmann, 2005). In order to deter-

mine flow and transport in complex natural hydrosystems, time and space are discretized

through various methods, leading to an approximated solution, which must be fairly close

to the reality and should be reached sufficiently fast"(Matta, 2018). Depending on the “key

hydrodynamic processes, the water quality concerns, as well as the proper scaling (time,

space) and the most suitable simulation scenarios"(Matta, 2018) the right model for the ap-

plication has to be chosen. In the case of sewer systems, the choice has to be made between

a horizontal one-dimensional (1D) or a three-dimensional (3D) approach. A 1D approach

is commonly used when the effects in two of the three dimensions can be neglected. This

is common for sewer network modelling, when filling ratios of pipes and main flow veloc-

ities are of main importance. When specifics of the detailled built environment - as in this

thesis - are of interest, i.e. the effects of a combined sewer overflow on the pollutant load,

3D modelling is preferred. This highlights one main advantage of the model developed in

this thesis: Instead of simplified 1D simulations for long sewer networks, local small-scale

effects can be analysed using a three-dimensional approach. The findings can lead to an

improved design of H2S hotpots.

As a basis for a three-dimensional model for odour and corrosion, the open source soft-

ware OpenFOAM is chosen. OpenFOAM offers the possibility of solving complex partial

differential equations and contains multiple solvers and utilities addressing different nu-

merical problems. It contains pre- and postprocessing environments and can describe two-

and three-dimensional problems (Greenshields, 2015). The partial differential equations are

discretized based on the Finite Volume Method in space and the Finite Differences Method

in time (Schulze and Thorenz, 2014). Due to the nature of the problem addressed in this

thesis, a two-phase flow model is needed in order to be able to account for water and air

flow. In order to introduce the software, a short overview over the model setup in Open-

FOAM is given and relevant solvers for multiphase problems are introduced. Setting up a

model in OpenFOAM can be divided into three main steps:

1. Preprocessing: Mesh generation, definition of boundary and initial conditions;

2. Solving;

3. Postprocessing.

Chapter 1: Introduction

The preprocessing step consists of the generation of the computational mesh and the def-

inition of boundary and initial conditions. A computational mesh is needed in order to

solve the relevant equations. By generating the mesh, the computational domain is divided

into a finite number of volumes where the solution is calculated. The number of volumes

influences the calculation effort (Schulze and Thorenz, 2014). Structured meshes can be cre-

ated with the built-in utility blockMesh. Several other mesh generators such as gmsh and

snappyHexMesh exist and come with different advantages and disadvantages (Kortelainen,

2009). Importing meshes from other CFD programs is also possible. The development of

a high-quality three-dimensional mesh is a complex and time consuming process, but the

effect of the mesh’s quality on the simulation results is substantial and should not be un-

derestimated.

For the second part of preprocessing, boundary and initial conditions need to be defined

for each variable that is being solved. Good initial conditions can lead to faster convergence.

For boundary conditions, the user has the possibility to define Dirichlet and Neumann con-

ditions but also variations combining both types exist.

For multiphase flow simulations, a number of different solvers are available in Open-

FOAM. Multiphase flows can be classified in different ways and depending on their classifi-

cation, a different kind of solver can be chosen (Schulze and Thorenz, 2015). A phase in this

context is a defined material class with homogeneous features, which can also be different

physical states of one fluid. Multiphase flows can then be classified by the physical state

of the different fluids as well as the topology of the interface between phases. According

to the physical state, one can therefore distinguish between gas-solid, gas-fluid, fluid-solid

and fluid-fluid flows. Gaseous phases can additionally be considered as compressible or in-

compressible flows depending on the flow velocity and speed of sound. When considering

the properties of the interface, one can differentiate between disperse and separated flows.

Separated flows are not mixed and share a continuous interface at which the phases gener-

ally share the same velocity. Disperse flows consist of one continuous phase and a disperse

phase. The phases penetrate each other and the interface between phases is not continuous

(Schulze and Thorenz, 2015). For the free-surface flows considered in this thesis, incom-

pressible water-air flows are of interest. The compressibility of air is only important for

very high velocities or pressure differences (Schulze and Thorenz, 2015). The free-surface

flows considered in this thesis are very closely related to the field of hydraulic engineer-

ing, where the Volume of Fluid (VOF) method is the most popular method used in CFD

applications. This method is able to describe non-mixing fluids but exchange of mass and

momentum are not accounted for in their original form. Disperse particles are only de-

scribed correctly when their size is bigger than multiple mesh cells (Schulze and Thorenz,

2015). Since the incompressible, separated two-phase flows are best described by using the

VOF method, this method will be used in this thesis.

In OpenFOAM, the VOF method is implemented in a solver called interFoam. The

governing equations that have to be solved simultaneously are the continuity equation for

compressible flow, the momentum equations and the αtransport equation together with the

constitutive relations for the density and dynamic viscosity (Ubbink, 1997).

The continuity equation is solved in the non-conservative form, since for two-fluid systems

with high density differences it is much more suitable for numerical solution, because the

velocity field Uis in contrast to the momentum ρUcontinuous at the interface (Ubbink,

1997). Additionally, incompressible flow is considered here, and therefore density changes

are negligible. As a result, the non-conservative form for incompressible flow can be written

Chapter 1: Introduction

as:

∇ · U=0 (1.1)

Here, Uis the velocity field [m/s].

As for the momentum equation, the basic form is:

∂ρU

∂t+∇ · (ρUU)=−∇p+∇ · T+ρfb(1.2)

Where ρis the density [kg/m3]; t is time [s]; p is the pressure [Pa]; Tis the deviatoric

viscous stress tensor [N/m2] and fbare body forces per unit mass [m/s2].

Since both fluids are considered to be Newtonian and incompressible, the stress tensor can

be decomposed into a more convenient form for discretization (Berberovi´c et al., 2009):

∇ · T=µ[∇U+ (∇U)T] = ∇ · (µ∇U) + (∇U)· ∇µ(1.3)

Where µis the dynamic viscosity [Ns/m2]. A detailed derivation and explanation can be

found in Ubbink (1997).

The VOF method used in the interFoam solver considers a single pressure system. prgh is a

modified pressure which is used to avoid the occurrence of steep pressure gradients caused

by hydrostatic effects and in order to simplify the definition of boundary conditions, thus

describing the static pressure minus hydrostatic pressure component:

prgh =p−ρ·g·x(1.4)

Here, gis the acceleration vector due to gravity [m/s2] and xis a spatial position vector [m].

This simplification has been first defined by Rusche (2003).

Taking these assumptions into account, the resulting momentum equation is:

∂ρU

∂t+∇ · (ρUU)=−∇prgh +∇ · (µ∇U)

+(∇U)· ∇µ−g·x∇ρ

(1.5)

Since the VOF method considers the two immiscible fluids liquid and gas as one fluid, the

fluid properties density and viscosity in the domain are described by a function based on

the volume fraction α:

ρ=αρL+ρG(1−α)(1.6)

µ=αµL+µG(1−α)(1.7)

(1.8)

Here, αis a volume fraction or indicator function [-] and the subscripts L and G denote the

fluids water (L - liquid) and air (G - gas). The viscosity divides into the physical viscosity

µphys as well as the turbulent viscosity µturb, where µturb is calculated by a turbulence model:

µi=µi,phys +µi,turb (1.9)

Where i is the subscript for the respective phase L or G.

The αtransport equation - also known as VOF equation - is used with an advanced formu-

lation. In general, it is implemented as an additional transport equation, which is used as

a marker to describe the distribution of the phases through the domain and therefore the

position of the free surface (Morgan, 2013). These marker variables are known as volume

Chapter 1: Introduction

fractions α, that are strictly bounded: 0 ≤α≤1 (Morgan, 2013). Both phases share an

interface at which both phases have the same velocity (Schulze and Thorenz, 2015). “In this

model an additional convective term originating from modelling the velocity in terms of a

weighted average of the corresponding liquid and gas velocities is introduced into the trans-

port equation for phase fraction, providing a sharper interface resolution. The model makes

use of the two-fluid Eulerian model for two-phase flow, where phase fraction equations are

solved separately for each individual phase" (Berberovi´c et al., 2009):

∂α

∂t+∇(αUL)=0 (1.10)

∂(1−α)

∂t+∇((1−α)UG)=0 (1.11)

“Assuming that the contributions of the liquid and gas velocities to the evolution of the free

surface are proportional to the corresponding phase fraction, and defining the velocity of

the effective fluid in a VOF model as a weighted average"(Berberovi´c et al., 2009) leads to:

U=αUL+ (1−α)UG(1.12)

Equation 1.10 can be rearranged and used as an evolution equation for the phase fraction α:

∂α

∂t+∇ · (αU)+∇ · ((1−α)Urα)=0 (1.13)

Where “Ur=UG−ULis the vector of relative velocity, designated as the “compression

velocity""(Berberovi´c et al., 2009). A detailled derivation of the equations can be found in

Damián (2012).

The values of αcan be interpreted for a water-air flow as follows:

α=









1 fluid L - water

0<α<1 transitional region

0 fluid G - air

(1.14)

In the transitional region, the water surface is considered to be the area where α=0.5.

In the range of high Reynolds numbers, turbulence models become relevant to correctly

solve CFD problems. OpenFOAM enables the use of a wide number of turbulence models,

ranging from Reynolds averaging to Direct Numerical Simulation (DNS).

DNS leads to the highest accuracy of the simulations, because the turbulence is discretized

by using small enough cell sizes and time steps to display all vortices (Maric et al., 2014).

This leads to a computational effort, which is in most cases still not practically applicable

today.

Large Eddy Simulation (LES) offers a compromise between DNS and RANS models in terms

of complexity. Small enough cell sizes enable the model to resolve large eddies spatially

and small eddies are modelled using a subgrid scale model. The necessary grid resolution

is therefore not as small as for DNS and the computational effort smaller (Maric et al., 2014).

In this thesis, the majority of cases are simulated using Reynolds averaging. Here, turbulent

flow is described by dividing the mean velocity into an average velocity component and a

fluctuating velocity component, resulting in a Reynolds stress tensor in the Navier-Stokes

equations (RANS). This tensor describes an additional eddy viscosity and leads to new un-

knowns in the equations, which can be either k (the turbulent kinetic energy) and e(the

Chapter 1: Introduction

turbulent dissipation) for k-eturbulence models or k and ω(the specific dissipation) for

k-ωturbulence models. These unknowns can be solved by using a two-equation model,

that adds two coupled transport equations in order to describe convection and diffusion of

turbulent energy.

The partial differential equations that describe the flow need to be transformed to alge-

braic equations by integrating them over a certain time step and a control volume. The user

has to define spatial and temporal discretization schemes for the above mentioned system

of equations (Schulze and Thorenz, 2014). Therefore, every differential operator as well

as the interpolation between fields have to be defined. OpenFOAM comes with the option

of defining one default parameter or one parameter for every term in the solution procedure.

The discretized equations then result in a sparse matrix which can be solved using it-

erative solution techniques (Schulze and Thorenz, 2014). The parameters of the solution

procedure also need to be defined by the user. For every variable, the procedure to solve

the systems of equations and possible preconditioners and tolerances for the exact solution

need to be defined.

For the VOF method, a special solution procedure is required for the four variables of

the flow field, the velocities in three directions and the pressure. The velocities can be cal-

culated with the Navier-Stokes equations but the mass conservation equation is not capable

to solve for the pressure. This issue can be solved by a pressure linking equation. In the

interFoam solver, this is done with an extended algorithm called PIMPLE (PISO-SIMPLE)

(Greenshields, 2015) which is based on the PISO (Pressure-Implicit with Splitting of Op-

erators) and the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm

(Caretto et al., 1973).

Due to the high number of cells in the computational domains, the OpenFOAM sim-

ulations presented in this thesis were mostly computed on supercomputers of the Nord-

deutscher Verbund für Hoch- und Höchstleistungsrechnen (hlrn) in Berlin or on high per-

formance computing (HPC) clusters of TU Berlin.

1.2.3 Transport and mass transfer modelling

1.2.3.1 Transport

When developing a model to predict the occurrence of odour and corrosion, it becomes clear

that the model needs to include transport processes in order to account for the movement of

H2S in the sewer. The previously outlined interFoam solver in OpenFOAM, which uses the

VOF method, is able to describe two-phase hydrodynamics but the transport of substances

need to be included separately. This is why in this Chapter, the basics of transport processes

are outlined.

Transport processes can be divided in conservative and reactive transport, where reactive

transport considers the formation of substances due to chemical reactions. In this thesis only

conservative transport will be considered.

Conservative transport can generally be divided into three driving processes which are

illustrated in Figure 1.5:

• Advection / Convection;

Chapter 1: Introduction

• Molecular Diffusion;

• Dispersion / Turbulent diffusion.

Advection is the linear transport of a particle along the flow direction without any spread-

ing taking place (Hvitved-Jacobsen et al., 2013) (see Figure 1.5).

Diffusion can be divided into molecular and turbulent diffusion. Molecular diffusion

describes Brown’s molecular movement and is caused by temperature-induced impact of

the molecules and is a disordered movement. On a macroscopic scale it can be expressed by

Fick’s first law of diffusion. The diffusivity of a substance depends on the substance itself

and on the transport medium. It is also temperature dependent (Hvitved-Jacobsen et al.,

2013).

Dispersion or turbulent diffusion is caused by flow velocity variations in space and time

and is related to turbulent conditions. Depending on the dimensionality of the computa-

tional model (1D, 2D or 3D), dispersion is included in different ways. When using a 1D or

2D model, dispersion is accounted for by using a dispersion term. In 3D models, it can be

accounted for by including a turbulent diffusion by employing a turbulent Schmidt number

(Bates et al., 2005).

Figure 1.5: Transport processes, Schankat (2010) after Barlag et al. (1998)

OpenFOAM’s interFoam solver does not have a built-in solver for transport equations

but the relevant equations were implemented into the solver. An advection-diffusion equa-

tion to describe transport of a passive tracer with a concentration C is shown in Equation

1.15). Here, the diffusivity is divided into the physical diffusivity Dphys and the turbulent

diffusivity Dturb which is calculated using Equation 1.16. The physical diffusivity Dphys and

the turbulent Schmidt number Scturb have to be defined by the user.

Advection-diffusion equation:

∂C

∂t+∇ · (UC) = ∇ · ((Dphys +Dturb)∇C) + W(1.15)

Chapter 1: Introduction

with

Dturb =µturb/ρ

Scturb

(1.16)

Where W can be a sink or source term or the production term related to chemical reaction.

1.2.3.2 Mass transfer

The description of relevant processes in Chapter 1.1 shows the relevance of accurately in-

cluding H2S mass transfer into a three-dimensional model for odour and corrosion. Mass

transfer between liquid and gas phase depends on many influence factors. The exchange

of mass itself is defined by mass transfer rates and equilibrium conditions. The mass flow

between the phases is in general decribed by the mass transfer coefficient. The equilibrium

between the two phases can be described by the Henry coefficient.

In general, three main approaches exist to describe mass transfer across the air-water inter-

face in a model. These are namely:

1. Two-film theory,

2. Penetration theories, and

3. Surface renewal theory.

The two-film theory assumes that two static films exist at the interface between the two

fluids. Through these films the volatile substance diffuses based on the molecular diffusion.

During this process, Fick’s first law is assumed to be consequently valid. This theory has

been applied in the WATS model. The advantage of this theory is the fact that it leads

to simple expressions for gas-liquid mass transfer that are useful for basic practical appli-

cations (Hvitved-Jacobsen et al., 2013). A formulation has been integrated into the WATS

model to account for mass transfer in drop structures and its formulation has been further

developed by the work by Matias et al. (2017). However, the simplicity comes with a price

of certain assumptions:

1. The species concentrations in any location is constant over time.

2. The films at the interface of both phases are laminar.

3. Equilibrium conditions at the interface are achieved immediately.

Wang et al. (2018) give an overview over multiple criteria, why these assumptions are com-

plicated to be applied to real-life applications. Local flow regimes could vary between

laminar, transitional or turbulent. Furthermore, variations of fluid properties could occur.

These variations cannot be accounted for when applying the two-film theory. Flow insta-

bilities can occur and change the liquid film. Most importantly, the two-film theory is a

one-dimensional approach.

Other approaches, such as the penetration theory and the surface renewal theory offer some

improvements compared to the two-film theory. Both account for the variability of the flux

over time. But still they do not account for local variations of fluid velocities, physical prop-

erties and different flow regimes (Wang et al., 2018).

This is why in recent years, more and more advanced CFD models have been developed

to analyse mass transfer phenomena, many of them using the VOF method. A one-fluid for-

mulation accounting for advection, diffusion and possibly chemical reactions can be used

Chapter 1: Introduction

to describe flow and chemical kinetics (Wang et al., 2018).

Using a CFD approach to describe mass transfer phenomena offers a number of advan-

tages, such as the possibility to account for three-dimensional effects. As outlined above,

H2S mass transfer processes in sewers can be influenced by high turbulent effects and three-

dimensionalities, even in simple structures such as circular pipes. It is therefore necessary

to apply a three-dimensional two-phase model which can account for mass transfer to the

environment of a sewer system.

In OpenFOAM, mass transfer can be simulated with different solvers. But as mentioned

above, the VOF method offers different advantages to describe the systems considered in

this thesis. Therefore, mass transfer will be considered using the one-fluid approach for the

VOF method as defined by Haroun et al. (2010a,b) as it has been implemented by Nieves-

Remacha et al. (2015) and Severin (2017). The approach is based on the interFoam solver

and considers one additional transport equation for both phases as outlined in Equations

1.15 and 1.16.

Haroun et al. (2010a) introduced a one-fluid formulation for transport and mass transfer

using Henry’s law as a constant coefficient. This approach has been implemented into the

JADIM code and compared to the penetration theory and a good agreement has been ob-

served. Nieves-Remacha et al. (2015) implemented the formulation into OpenFOAM. Note

that the formulation by Haroun et al. (2010a) considers D as the molecular diffusivity and

thus does not account for the turbulent diffusivity Dturb. This is due to the fact, that DNS

simulations have been performed, and no turbulent viscosity µturb has been computed by

the solver. The work shown in Chapter 7 will highlight the importance of including Dturb

into the mass transfer formulation presented in this thesis. The following equations contain

the turbulent diffusivity Dturb:

∂C

∂t+∇ · (UC) = ∇ · ((Dphys +Dturb)∇C+φ) + W(1.17)

With φbeing a concentration flux expression at the interface resulting in the following:

φ=−(Dphys +Dturb)C(1−He)

α+He(1−α)∇α(1.18)

The production term will not be considered in this thesis. As shown in Hvitved-Jacobsen

et al. (2013), the Henry coefficient is highly dependent on the temperature. The mass trans-

fer rate depends on the Reynolds number and pH value (Yongsiri et al., 2004). These driving

factors are not accounted for in the original form of the approach.

At the interface, Henry’s law must be fulfilled and the molecular diffusive transport

expressed by Fick’s law is supposed to be the same in both phases:

He =CL

(1.19)

(Dphys,L+Dturb,L)∇CL= (Dphys,G+Dturb,G)∇CG(1.20)

Like density and viscosity, the concentrations and diffusion coefficients are considered as

single-phase properties depending on the phase fraction value α. While the concentrations

Chapter 1: Introduction

are calculated by linear averaging, the diffusion coefficient is calculated by using a harmonic

average:

C=αCL+CG(1−α)(1.21)

Dphys =Dphys,L·Dphys,G

αDphys,L+ (1−α)Dphys,G(1.22)

The diffusion coefficients for Dphys,Land Dphys,Gare defined by the user. Note that these

coefficients are also temperature dependent and that the temperature dependency of the

coefficients has to be taken into account by the user when defining the values.

In order to validate mass transfer rates, the local mass transfer can be computed as

follows:

KL,local =−((Dphys,L+Dphys,G)∇C+φ)·nL

∆CL,local

(1.23)

with nLbeing the normal to the interface pointing into the liquid which is obtained by

nL=∇α

|∇α|,∆CL,local is the difference between the free surface region and the bulk liquid

concentration. The overall mass transfer across the interface can be obtained by integrating

over the interface

KL=1

λZλ

0KL,localdξ(1.24)

where λis the interface length and ξthe curvilinear coordinate associated to the interface.

Following Carrera et al. (2017), the mass transfer can also be obtained by fitting the decrease

of concentration in the liquid phase to the following function:

CL,H2S=CL,0e−KL,H2Sa(t−t0)(1.25)

The approach by Haroun et al. (2010a) offers another advantage when modelling transport

processes using the interFoam solver. Because one set of Navier-Stokes equations is solved,

a standard advection-diffusion equation cannot account for the interface between phases.

If a tracer is modelled in the nearfield of the interface between two phases, an unphysical

spreading can occur. Using the approach by Haroun et al. (2010a) with a small value for

the Henry coefficient can avoid the unphysical spreading to a certain extent and therefore

enable the user to simulate single-phase transport.

1.3 Scope of this thesis

In this thesis, a model for the correct description of processes occurring at the water-air inter-

face was developed. The main focus lies in the use and extension of the three-dimensional

two-phase solver interFoam, which is implemented in OpenFOAM, to describe water-air

flow and H2S emissions across the air-water interface in closed conduits.

Therefore, existing approaches are applied and extended to describe three-dimensional

water (single phase) and water-air flow (two phases) and conservative tracer transport of

H2S. First, a validation of the hydrodynamics of the interFoam solver is performed. Then, an

existing mass transfer solver (interHarounFoam) which accounts for mass transfer depend-

ing on the Henry coefficient as described by Haroun et al. (2010a) is extended to describe

the specifics of H2S transfer in sewers. As a result, a three-dimensional implementation of

parts of an in-sewer model exists and this solver is then applied and tested under turbulent

Chapter 1: Introduction

conditions in a stirring tank.

This thesis is structured into nine chapters consisting of the introduction (the current

Chapter 1), three peer-reviewed journal articles (one in press, two submitted), one book

chapter, two peer-reviewed conference contributions, further related work and a synthesis.

An overview over the different topics covered in each Chapter is given in Figure 1.6.

The thesis has several connection points to the two other projects in UWI which are

working on interfaces in sewer systems. The relevance of the collaboration and the mutual

influence will be outlined in Chapter 2. The projects focus on two main aims: The enhanced

understanding of odour and corrosion mechanisms and the development of a CFD simula-

tion model. The research is carried out at the three different interfaces biofilm - (waste)water

interface, the (waste)water - (sewer) air interface and the (sewer) air - biofilm - (concrete)

wall interface. This Chapter aims to give an overview over the integration of the work into

the Research Training Group.

After this overview, Chapters 3 and 4 show results of a hydrodynamic study with the

interFoam solver for two-phase flow in closed conduits. Chapter 3 is a validation of the

water phase described by the numerical model concerning different hydraulic conditions,

the consideration especially focuses on flow in closed conduits. The validation is carried

out by analysing the simulation results of three different test cases: single-phase flow over

a hill under consideration of different turbulence models, free-surface flow over a ground

sill and the analysis of a complex sewer geometry. The latter case has been analysed to in-

vestigate the stability of the simulations using the setup of a closed conduit and the overall

accuracy in such a complex model domain. The results of the complex sewer are compared

with results of an existing CFD model as well as measured results from a 1:20 scale model.

Chapter 4 deals with a further parameter study of a test case analysed in Chapter 3. A

two-phase free-surface flow over a ground sill has been simulated and different parameters

such as the structure of the ground sill, discharge, water level and flow regime have been

evaluated concerning their influence on the water level drawdown.

Chapter 5 introduces two different ways to analyse single-phase transport in a mul-

tiphase system. The implementation of a standard advection-diffusion equation into the

interFoam solver can lead to an unphysical spreading across the water-air interface. In a

first example, tracer transport around concrete probes in the headspace of the sewer pilot

plant of the Berliner Wasserbetriebe (BWB) is simulated. The model setup considers the

two-phase system of the sewer as a single-phase system by describing the water surface as

a boundary condition. Then, a standard advection-diffusion tracer transport equation is ap-

plied. In a second case, tracer transport in the complex sewer stretch described in Chapter

3is simulated. Instead of a passive tracer, the interHarounFoam solver by Haroun et al.

(2010a) is used to achieve phase-constrained transport.

Chapter 6 explores the capabilities of the interHarounFoam solver for mass transfer

which is based on the two-phase flow approach used in Chapters 3 and 4 and extends the

solver for specifics of H2S mass transfer in sewers, namely the temperature dependency

of the Henry coefficient, a conversion from pH value and total dissolved sulphide to the

equilibrium between H2S and bisulphide ion (HS-) and a conversion of the gas phase con-

centration H2SGto its respective partial pressure. The latter two are especially important

when considering the experimental work carried out at the sewer pilot plant of the BWB

as outlined in Chapter 2: Being able to insert the measured pH value and the total dis-

Chapter 1: Introduction

solved sulphide, enables the user to directly feed measurements into the model. The partial

pressure of H2SGis another parameter that is measured in the pilot plant and supports the

user in the direct interpretation of modelling results in the close future. The advantages

of three-dimensional mass transfer simulations in complex geometries are highlighted by

performing mass transfer simulations in the complex sewer geometry presented in Chapter

Chapter 7 uses the functionalities of the solver introduced in Chapter 6 to analyse and

validate the effect of turbulence on the H2S mass transfer coefficient in a stirring tank.

Therefore, the interH2SFoam solver is connected to a dynamic meshing functionality and

extended to compute the occurring mass transfer coefficient. Being able to describe the mass

transfer in a three-dimensional model can help to improve existing model concepts as well

as to better examine hotspots of H2S emissions in order to give recommendations for design

improvement. The analysis of turbulence effects of drop structures on H2S mass transfer, as

it is introduced in Chapter 2, is one example of the importance of the correct description of

turbulent effects.

Chapter 8 gives an overview of supplementary scientific work and Chapter 9 synthesizes

the outcomes of this thesis and gives an outlook for future research.

Figure 1.6: Overview over different chapters

Chapter 2

Research on interfaces in sewer

systems within the DFG Research

Training Group “Urban Water

Interfaces"

This study was published as:

Teuber, K., Sielaff (née Grüneberger), M., Despot, D., Stephan, D., Barjenbruch, M. &

Hinkelmann, R.: Modeling and Measuring of Interfaces in Sewer Systems, Proceedings

of the 37th IAHR (International Association for Hydro-Environment Engineering and Re-

search) World Congress, Kuala Lumpur, Malaysia, 2017.

This is the postprint version of the article.

The test cases’ setups are listed in Appendix B (single-phase flow over ground sill: B.1,

free-surface flow over ground hill B.2, complex sewer: B.3).

Abstract

This paper presents the research focusing on interfaces in sewer systems as it is carried out

within the DFG Research Training Group Urban Water Interfaces. During its transporta-

tion, wastewater in sewer systems undergoes a number of physical, biological and chemical

processes and transformations. Under certain conditions such as high detention times, the

formation of hydrogen sulfide (H2S) in sewer systems leads to odor in the sewer atmosphere,

and after a biogenic oxidation to sulfuric acid (H2SO4) to corrosion at sewer walls. In 1998,

the restoration costs for corroded sewers in Germany were estimated to be in a range of bil-

lions of US $ (Kaempfer and Berndt, 1999). High concentrations of odorous substances in the

atmosphere can even lead to death of sewer workers. Within the research training group,

three dissertation projects focus on two driving aspects: The enhanced understanding of

odor and corrosion mechanisms and the development of a CFD simulation model. In order

to gain a deeper understanding, a research pilot plant owned by the Berliner Wasserbetriebe

(BWB) is operated. The work of this group enables a deeper in-detail understanding of H2S

formation in sewer systems and can therefore lead to an improvement of existing models as

well as to developments that enable a detailed analysis of odor and corrosion problems in

sewers.

Chapter 2: Research on interfaces in sewer systems

2.1 Introduction

Every year damages of sewers due to corrosion cause high costs for sewer maintenance.

At the same time, high hydrogen sulfide (H2S) concentrations in the sewer atmosphere are

a health risk for sewer workers. The processes leading to odor and corrosion as well as

the empirical and conceptual description of these processes have been investigated since

more than 70 years (e.g. Parker (1945), Gilchrist (1953), Thistlethwayte (1972)). However,

within the last 20 years a deeper understanding has been gained thanks to the efforts of

research groups in Denmark and Australia (Rootsey and Yuan (2010), Rootsey et al. (2012),

Hvitved-Jacobsen et al. (2013)). Conceptual model approaches have been developed in or-

der to estimate corrosion risks. The aim of our research is to gain a deeper insight into

the influence of small-scale structures such as obstacles and drop structures in sewers, to

evaluate the comparability of different laboratory tests for sewer corrosion and to develop

a Computational Fluid Dynamics (CFD) simulation model that is able to account for three-

dimensional effects.

In order to better understand the processes investigated, the main processes leading to

odor and corrosion shall be outlined in the following (see Figure 2.1). Under anaerobic

conditions in sewage, sulfate present in the wastewater can be reduced to sulfide by sulfate-

reducing bacteria (SRB) residing in the biofilms on the walls of sewer pipes (Sharma et al.,

2008a). Sulfide is diffused from the biofilm into the water phase as H2S. Influenced by the

pH value and the temperature, different amounts of sulfide as H2S and bisulfide ion (HS-)

are present in the water phase. Described by the air-water equilibrium, emission of H2S

from the water into the air phase can occur. This process depends on factors such as the

air and water phase velocities, pH value, temperature and the concentration of oxygen and

nitrate. The air-water equilibrium for a volatile compound such as H2S can be described by

Henry’s law which describes the relative amount of a volatile compound in the gas phase

as a function of its relative occurrence in the water phase under equilibrium conditions and

at constant temperature. The temperature dependency of Henry’s law can be described by

different equations for example, the van’t Hoff equation. The concentration of H2S in the air

phase defines the intensity of odor. Another process taking place at the air - water interface

is reaeration which is the transfer of oxygen across the air - water interface. This process is

the only way to supply oxygen to the water phase and therefore influences the potential of

aerobic and anaerobic processes in the wastewater (Hvitved-Jacobsen et al., 2013).

At moist concrete pipe walls, H2S can be oxidized to sulfuric acid (H2SO4) by aerobic

microbial reactions (Hvitved-Jacobsen et al., 2013). Sulfuric acid can lead to the corrosion of

the building material of the sewer system, especially of concrete and other cement-bound

materials. A combined attack of the cement stone mineral’s chemical dissolution and of

expansive mineral formation (e.g. gypsum) can result in corrosion rates up to 1 cm yr-1 or

higher, even for concretes developed for such conditions (Grengg et al., 2015).

The research is carried out in three different sub-projects focusing on different interfaces

in sewer systems as well as applying different principles. The projects collaborate in diverse

states and exchange results of their investigations. Measurements of H2S emissions as well

as corrosion rates are performed in a research pilot plant. The findings of this research pilot

plant are used in order to validate a CFD model describing H2S emissions across the air -

water interface. Corrosion rates are measured on a laboratory scale, a pilot plant scale and

on a sewer scale. Different concrete samples are compared and different laboratory tests are

compared among each other in order to generalize the test procedure for corrosion at the

Chapter 2: Research on interfaces in sewer systems

sewer atmosphere - biofilm - sewer wall (concrete) interface.

In the following, the respective projects will be presented. In Section 2.2, the experi-

ments planned in the research pilot plant will be outlined. Section 2.3 presents the different

tests that will be performed in order to investigate corrosion rates. An overview of the

development of the CFD model will be given in Section 2.4.

Figure 2.1: Interface processes in sewer systems.

Chapter 2: Research on interfaces in sewer systems

2.2 Pilot plant investigations

The research pilot plant provided by Berliner Wasserbetriebe (BWB) facilitates a unique op-

portunity for performing experiments under realistic but controlled conditions (Figure 2.2).

Here, wastewater from a wastewater pumping station in Berlin is directly fed into the ex-

perimental sewer pipes of the research pilot plant. The facility comprises of two pressure

sewers with diameters of 100 mm that act as a fermenter, providing anaerobic conditions

for the wastewater to reach the required levels of septicity and thus adequate H2S forma-

tion (Figure 2.3). The wastewater is then conveyed into two corresponding gravity sewers

having a length of 25 m each and a diameter of 400mm. The experimental sewer lines are

duplicated so that one is used for experiments dealing with containment strategies and the

other is used as a reference line for comparison. The concrete samples outlined in Section

2.3.2 will be placed in both the experimental line and the reference line. Sensor probes are

used for continuous measurements in both the liquid and gas phases which can be accessed

digitally by storage on a database. Measurement of pH values, temperature, dissolved oxy-

gen (DO), equivalent chemical oxygen demand (CODeq) and equivalent suspended solids

(TSeq) in the water phase can be made at various sampling points of the pilot plant setup.

The s::can spectro::lyserTM UV-VIS spectrometer probes (s::can Messtechnik GmbH, Austria)

are used for continuous measurement of total dissolved sulfides and nitrate concentration

in the liquid phase. As for the gas phase, it is possible to measure the H2S, temperature and

relative humidity. The H2S concentration in the gas phase is measured using the Kemira

H2S-GuardTM.

The experimental devises are used to further understand the key physicochemical and

biochemical processes responsible for odor and corrosion in four different experiments. The

experiments will be outlined in the following.

Figure 2.2: Research pilot plant. Figure 2.3: Schematic drawing of the research pilot

plant.

2.2.1 Sulfide formation in pressure or rising mains

Severe problems are associated with the microbial reduction of sulfate to sulfides in sewer

systems, especially in pressurized flows (Hvitved-Jacobsen et al., 1995). Previous studies

on sulfide production in pressure mains have resulted in a number of empirical formulas

Chapter 2: Research on interfaces in sewer systems

Figure 2.4: Experiment 2 – Elevated H2S emissions due to gross solids in wastewater.

which are still used in praxis and are based on parameters such as biochemical oxygen

demand (BOD), chemical oxygen demand (COD), sulfate concentration and temperature

etc. (Thistlethwayte (1972), Boon and Lister (1975), Pomeroy and Parkhurst (1978), Hvitved-

Jacobsen et al. (1988), Nielsen et al. (1998)). In this investigation, the mostly used empirical

formulas are to be evaluated using statistical analysis to determine the relationships of the

key parameters promoting sulfide formation in pressure sewers, especially sulfide produc-

tion at the biofilm - wastewater interface. The main goal of the experiments carried out in

the first stage is to gain an improved understanding of the conversion processes within the

sewer sulfur cycle under anaerobic conditions.

2.2.2 Elevated H2S emissions in gravity sewers due to flow interruptions caused

by the accumulation of gross solids (debris / solid materials)

Gross solids are of particular concern for sewer systems since they can cause maintenance

problems such as blockages and their sedimentation can increase the formation of toxic

gases (e.g. sulfide and methane) (Rutz, 2016). The sewer blockage formation process is

attributed by materials such as plastics, wet-wipes etc. that enter into gravity sewers. The

accumulation of gross solids in gravity sewers can contribute to the emission of H2S in

two ways: Due to a disruption of the flow regime, turbulence effects can occur at the sur-

rounding regions of the obstacle and a higher amount of pollutants develops due to the

accumulation of sediments (Ashley, 2004).

Until now, the influence of these small-scale structures in the water phase on the emis-

sions has not been quantified. Investigating this relation can lead to better predictions in

simulation models as well as better recommendations for countermeasures. The aim of

this experiment is to quantify the levels of H2S emitted into the sewer atmosphere at these

hotspots with respect to different hydraulic conditions (Figure 2.4). The highlighted inter-

face for this experiment is the wastewater - sewer atmosphere interface in gravity sewers.

2.2.3 Investigation of H2S emissions at drop structures of sewer systems

The release or emission of H2S to the sewer atmosphere is known to be related to turbulence,

pH, temperature and wastewater constituents (Matias et al., 2016). Drop structures found at

manholes and transfer zones of sewer networks in sewers lead - similar to small-scale obsta-

cles - to changed flow conditions and therefore facilitate increased H2S emissions. Previous

efforts in understanding the H2S emission processes at drop structures of sewer systems

have been only elaborated in laboratory experiments (i.e. with artificial wastewater) (Matias

Chapter 2: Research on interfaces in sewer systems

Figure 2.5: Experiment 3 – Pilot plant scale investigation of H2S emissions at drop structures.

et al., 2016). Extending these investigations to a pilot scale study can help to verify the

laboratory experiments to more realistic conditions (sewer geometry and real wastewater).

The knowledge obtained from this investigation provides essential information to predict

the effects of sulfides in wastewater infrastructures, and has the potential to substantially

improve the design and management approach of sewer systems (Matias et al., 2016).

A drop structure as displayed in Figure 2.5 will be added to the pilot plant and different

parameters will be varied. Again, the H2S emissions into the sewer atmosphere will be

measured. As in the previous experiment, the wastewater - sewer atmosphere interface is

the main focus. The tailwater depth will be calculated by using the validated CFD model

which will be described in Section 2.4.

2.2.4 Impact of countermeasures at hotspots

In a fourth experiment, the impact of different countermeasures at hotspots is being tested.

Problems of odor and corrosion in sewers are solved by several methods which have been

intensively investigated. Today, most water companies and operators of sewer networks

rely on several chemicals which have been proven to be effective in controlling sulfides in

sewers. Chemical dosing of oxygen, nitrate, iron salts caustic, free nitrous acids and mag-

nesium hydroxide (Mg(OH)2) are all examples of effective chemicals that are most widely

used (Ganigué et al., 2016). Recent research on the application of countermeasures for

controlling and reducing sulfides has been carried out on the optimization of the dosage

strategies. Two optimization techniques that can be used are: (1) to consider the dosage lo-

cation and response time for the chemicals to take effect, (2) online dosing control strategies.

The countermeasures to be used in this research project were selected with respect to the

optimization techniques mentioned above. The first countermeasure application is based on

a new method proposed by Auguet et al. (2015) which studied the effectiveness of down-

stream nitrite dosage. This study was conducted using a lab-scale sewer system and until

now there is no report on its application to field studies. Evaluation of this method under

Chapter 2: Research on interfaces in sewer systems

more realistic sewer network conditions can be studied in the pilot plant setup. A further

step in this investigation will consider an online dosing control strategy.

The second countermeasure application will be based on the online control of Mg(OH)2

dosing. This investigation follows Ganigué et al. (2016) which developed an online control

algorithm for the optimized dosing of Mg(OH)2for sulfide mitigation in sewers. The use of

Mg(OH)2for sulfide control in sewer systems in Germany is undocumented. Therefore, an

assessment of this chemical dosing usage will also be made during this investigation.

2.3 Investigations of corrosion processes in sewer and laboratory

tests

In the second research topic, corrosion processes in sewers at the interfaces of sewer atmo-

sphere, biofilm and concrete sewer is being investigated. The use of concrete with a higher

resistance against biogenic sulfuric acid corrosion (BSC) can increase the lifetime of sewer

system components and reduce rehabilitation and replacement costs. In the past much re-

search on the development and testing of high resistance concrete has been done and a

large variety of laboratory performance tests for acid resistance of concrete has been devel-

oped (e.g. De Belie et al. (2002), Belgium; Petersen and Lohaus (2006), Germany; Fourie

and Alexander (2008), South Africa; OENORM B 4710 appendix K, Austria; Hüttl (2008),

Germany). Those performance tests have not been experimentally compared. Also, most

laboratory performance tests have not been directly compared to BSC in sewer systems by

using the same concrete mixtures. Therefore, concrete samples in a sewer pilot plant will

be compared with the same concrete mixtures tested in four laboratory performance tests

developed in Germany and Belgium. The test developed by the Materialprüfungsanstalt

(MPA) Berlin- Brandenburg (hereafter MPA Berlin-Brandenburg test) is performed in the

original facility of MPA-Berlin- Brandenburg (now Kiwa Berlin) at TU Berlin. The facilities

according to the LPI test and the E DIN 19573 appendix B test have been installed at TU

Berlin for this project. The TAP test will be performed in the original facility at the Mag-

nel Laboratory for Concrete Research at Ghent University (Belgium). A special focus lies

on the time laps effect of the laboratory performance tests and the corrosion mechanisms

and sulfate containing components on the interface between sewer atmosphere, biofilm and

concrete.

2.3.1 Concrete and mortar types

Five concrete types with different binders and water/binder value (w/b) are being pro-

duced to be tested in the pilot plant and in three laboratory performance tests. The different

compositions result in varying acid resistances. The use of concrete with low, medium and

high acid resistance will allow a comprehensive comparison of the pilot plant processes and

the laboratory performance tests. Table 2.1 presents the composition of the five concrete

types and their assumed acid resistances. A maximum grain size of 16 mm and the grading

curve between A16 and B16 according to DIN 1045-2:2008-08 appendix L for all concrete

types were chosen.

The laboratory performance test after E DIN 19573 appendix B (E DIN 19573:2013-06)

requires mortar samples. The binder composition and w/b values of the five mortar types

Chapter 2: Research on interfaces in sewer systems

Table 2.1: Composition and assumed acid resistances of concrete types used for acid resis-

tance tests. Cement types are classified after DIN EN 197-1:2011-11. MS: microsilica, FA:

hard coal fly ash.

Concrete type 1 2 3 4 5

Type of binder white CEM I CEM I 42.5 R CEM III B 42.5 CEM III B 42.5 CEM I 42.5 R

42.5 R (dw) N-NA N-NA - MS - FA

mbinder/Vconcrete 360 360 360 373 350*

[kg/m3]

w/b value 0.45 0.45 0.45 0.35 0.42

assumed acid very low low medium high high

resistance

*cement 270 kg/m3; microsilica 27 kg/m3; hard coal fly ash 53 kg/m3.

are equivalent to the concrete types. For each mortar type 450 g binder and 1350 g CEN-

Reference-Sand DIN EN 196-1:2005-05 are being used.

2.3.2 Pilot plant acid resistance tests

The concrete samples are positioned in the concrete sample point 1 (see Figure 2.4) of both

gravity pipes in the pilot plant. Six cuboids per concrete type are stored in each gravity

pipe. The cuboids’ edge lengths are 150 mm x 100 mm x 40 mm. As a reference sample,

another six cuboids per concrete type are stored over water in a closed container at 20◦C.

One cuboid per concrete type is being taken from each gravity pipe of the pilot plant every

six months. Also one cuboid per concrete mixture of the reference samples is being taken

for analysis every six months. Results will be evaluated using the parameters monitored in

the pilot plant, such as H2S concentration in the sewer gas phase.

2.3.3 Laboratory performance tests

All concrete types are analyzed using the MPA Berlin-Brandenburg (Hüttl, 2008), the LPI

(Petersen and Lohaus, 2006) and the TAP test (De Belie et al., 2002). All mortar types are

examined according to E DIN 19573:2013-06. An overview of the four tests is given in Table

2.2.

2.3.4 Analysis of concrete and mortar specimens

Concrete and mortar specimens from the pilot plant and the laboratory performance tests

are being analyzed with reflected-light microscopy of polished samples and polarizing light

microscopy of thin sections to measure the damage depth. In addition scanning electron

microscopy with energy-dispersive X-ray spectroscopy will be performed. Powder X-ray

diffraction will be done to characterize minerals formed during the experiments. Micro

X-ray fluorescence spectroscopy will be performed to produce element maps of polished

samples and to measure the damage depth. The proton consumption is calculated using the

titrated acid in E DIN 19573 appendix B test. The concrete degradation and surface rough-

ness is measured by laser measurements during the TAP test (see De Belie et al. (2002)).

Because for all tests the same concrete types are being used, the corrosion in the pilot

plant and in the laboratory performance tests can be directly compared. New information

Chapter 2: Research on interfaces in sewer systems

is being gained on the applicability of laboratory acid resistance tests concerning sewer sys-

tems. The time laps effect of the laboratory performance test will be estimated. Thereby

more detailed insights are been gained concerning corrosion mechanisms caused by BSC.

2.4 Development of a three-phase CFD simulation model

In the third project, a CFD simulation model is being developed in order to numerically de-

scribe formations across the wastewater - sewer atmosphere interface. So far, existing model

approaches (Sharma et al. (2008a); Hvitved-Jacobsen et al. (2013)) are one-dimensional ap-

proaches, neglecting three-dimensional flow effects in sewage and air. The three-dimensional

model approach will be able to verify these assumptions and can be an extension of these

models in respect to a more detailed analysis of hydraulic aspects.

The work carried out in this research divides in three different parts. The first step

is a validation of the water phase described by the numerical model concerning different

hydraulic conditions. Results obtained within this working package will be presented in

Section 2.4.2. Because of the fact that sewer systems are focused on in the work of this

group of researchers, the consideration especially focuses on flow in closed conduits. After

the first validation step is completed, the model can be used to support the experimental

work in the pilot plant in hydraulic questions. Possible applications are the estimation of

the tailwater depth in the experiment outlined in Section 2.2.3. The two remaining working

steps are the validation of the air phase behavior in closed systems as well as the implemen-

tation of transport as well as mass transfer processes depending on factors such as Henry’s

law in order to describe H2S formations. At a later stage, the model could even be extended

to describe corrosion effects at the atmosphere - biofilm - sewer wall (concrete) interface and

interact with the research project outlined in Section 2.3.

2.4.1 Numerical model

Surface water flow is calculated by using the two-phase flow solver interFoam based on

a volume of fluid (VOF) approach for one- and two-phase flows as it is implemented in

the open source model OpenFOAM. Both phases are considered as one fluid with rapidly

changing fluid properties, therefore one set of Navier-Stokes equations is solved. The phases

are distinguished by an additional transport equation for the volume fraction which is used

as a marker to describe the distribution of the phases throughout the domain. The equations

can be formulated as follows Rusche (2003):

Mass conservation equation:

∇ · −→

U=0 (2.1)

Momentum conservation equation:

∂ρ−→

∂t+∇ · ρ−→

U−→

U=−∇prgh +∇ · µ∆−→

U+∇−→

U∇µ−−→

g·−→

x∇ρ(2.2)

Where prgh is the static pressure minus hydrostatic pressure:

prgh =p−ρ·g·h(2.3)

Chapter 2: Research on interfaces in sewer systems

Table 2.2: Testing parameters of the four laboratory performance tests.

Testing

parameter

MPA Berlin-

Brandenburg test

LPI-test E DIN 19537

appendix B test

(bath test)

TAP test

Test specimen:

Type concrete concrete or mortar mortar concrete

Geometry 4x15 cuboids:

150 mm x 100mm x 40

2x4 concrete cuboids:

70 mm x 35 mm x

150 mm

2x6 mortar cuboids:

40 mm x 20 mm x

160 mm

5 cuboids:

40 mm x 40 mm x

80 mm

cylinder:

d = 270 mm,

h= 70 mm

Testing age 28 days* 28 days* 28 days* 28 days*

Testing medium:

Type sulfuric acid sulfuric acid* sulfuric acid sulfuric acid or

acetic/lactic acid mix

Concentration pH 3.5 pH 3.0* pH 4.0 sulfuric acid:

pH 0.8 - 1.0

acetic/lactic acid mix:

pH 2.0 - 2.2

Performance:

Test duration 12 weeks 12 weeks* 4,000 h dependent on

testing procedure

Testing facility 5 tanks (45 l),

1 reservoir tank (80 l)

tank (13 l) tank (4 l) tank (2 l)

Acid renewal every 2 weeks every week* every 1,000 h after every attack cy-

cle

Constant pH yes

(automatic titration)

yes

(automatic titration)

yes

(automatic titration)

Mixing of

medium

circulation of

medium

rotation of sample in

acid tank

(150 s per revolution)

magnetic stirrer rotation of sample in

acid tank (1.04

revolutions per h)

Abrasion manual brushing of a

part of the samples

once per week

automatic brushing of

a part of the samples

every 150 s

none automatic brushing

after every attack

cycle

Evaluation:

Analysis thin-section polariz-

ing light microscopy

and mass loss

thin-section polariz-

ing light microscopy

proton consumption

during the test

laser measurements

Assessment

criterion

measured damage

depth in comparison

with reference

concrete

measured damage

depth

calculated damage

depth

concrete degradation,

surface roughness

*Parameters can be changed according to the behavior and application of the building material.

Chapter 2: Research on interfaces in sewer systems

Volume of Fluid equation:

∂α

∂t+∇ · α−→

U+∇ · ((1−α)Urα)=0 (2.4)

with the following parameters:

ρ=αρw+ρa(1−α)(2.5)

µ=αµw+µa(1−α)(2.6)

µ=µw+µa(2.7)

where −→

Uis the velocity field [m/s]; ρis the density [kg/m3]; t is time [s]; p is pressure

[Pa]; µis dynamic viscosity [Ns/m2]; g is acceleration vector due to gravity [m/s2]; x is a

spatial position vector [m]; αis volume fraction or indicator function [-]; Uris the relative ve-

locity between the phases [m/s]; the subscripts a and w denote different fluids air and water.

The indicator function αis defined as:

α=









1 fluid water

0<α<1 transitional region

0 fluid air

The water surface is considered as the area described by a volume fraction α= 0.5. For

single-phase test ases the volume fraction αis 1 and constant over the whole domain and

during the simulation time. A special definition of boundary conditions for two-phase flows

is necessary in order to obtain stability of simulations in closed ducts. The upper and lower

walls are defined as no-slip conditions. Depending on the fact whether the simulations

are two-or three-dimensional, the sidewalls are defined as no-slip conditions for three- di-

mensional cases and as empty boundary conditions for two-dimensional cases. The empty

boundary condition in OpenFOAM is a special boundary condition for two-dimensional

model setups. The definition of inlet and outlet is displayed in Figure 2.6. The inlet is di-

vided in two parts, an inlet for the air phase, where a total pressure is defined. The second

part is the inlet for the water phase, defined by a fixed discharge or flow velocity. At the out-

let, a constant pressure is defined leading to a free outflow of the water without definition

of a certain upstream water level. This upstream water level is achieved by adding a weir

structure in the proximity of the outlet to the model geometry. This type of outlet boundary

condition has originally been outlined for open systems with top atmospheric boundaries

by Bayón and Lopez-Jimenez (2015). For a more detailed description of further properties,

such as the chosen boundary conditions for each test case and the turbulence models, the

reader is referred to Teuber et al. (2019a).

2.4.2 Validation of the water phase behavior

In order to secure the applicability of the chosen numerical model, flow properties such as

the water flow velocity and the water surface behavior described by the VOF approach have

been validated (Teuber et al., 2016, 2019a).

First, a validation of the water velocity has been performed by using experimental re-

sults of single-phase water flow over a ground sill in a closed duct by (Almeida et al., 1993).

Different Reynolds averaged (RANS) turbulence models and Large Eddy Simulations (LES)

Chapter 2: Research on interfaces in sewer systems

Figure 2.6: Model domain of water surface validation case.

were compared as well. The RANS models investigated were the Standard k-e(Launder

and Sharma, 1974), k-ω(Wilcox, 1988) and k-ωShear Stress Transport (SST) model (Menter,

1993, 1994). As the subgrid-scale model for the LES simulations, the Smagorinsky model

(Smagorinsky, 1963) was chosen. The results led to the conclusion that the chosen RANS

models as well as the LES simulations lead to a good approximation of the experimental re-

sults. Due to the nature of the different turbulence models, the LES simulations were able to

capture fluctuations in the flow velocity which can be of interest for some application areas.

With this advantage of LES comes the disadvantage of a higher necessary mesh resolution

which leads to significantly increased computation times.

In order to validate the behavior of the water surface, two-phase flow over a ground

sill was simulated. The basic setup of the model domain is displayed in Figure 2.6. Simula-

tions were performed for different two- and three-dimensional model setups with variations

concerning the structure of the sill, discharge, water level and the flow regime. Since fluctu-

ations of the flow velocity were of minor interest, the k-eturbulence model was used. The

water level drawdown was compared to analytical solutions obtained by using continuity

and Bernoulli’s equation. A detailed discussion of the results can be found in Teuber et al.

(2016). Overall, the results showed a good agreement of the simulations with measured wa-

ter velocities and a reasonable behavior of the water - air interface. However, the numerically

computed water level drawdown was smaller than the analytically calculated drawdown.

Explanation can be found in the treatment of single losses. The three-dimensional CFD

model is able to account for single losses caused by the structure of the sill whereas the

analytical solution is a one-dimensional solution that neglects single losses.

In a last step, simulations have been performed in a complex sewer geometry described

in Bayón et al. (2015), proving the applicability of the simulations under demanding hy-

draulic conditions. This geometry contains different hydraulic structures such as weirs,

different flow regimes and a hydraulic jump. The simulation results were compared to ex-

perimental measurements carried out on a 1:20 scale model as well as to CFD simulations

previously performed by Bayón et al. (2015) by using a different model setup (open system,

without upper wall). A detailed analysis will be published in Teuber et al. (2019a). The

results showed a good agreement between both the existing simulation model as well as

the experimental results and lead to the conclusion, that the solver is applicable to complex

hydraulic cases in closed conduits such as sewer systems.

Chapter 2: Research on interfaces in sewer systems

Overall, the simulations performed show that the solver is able to simulate different hy-

draulic cases. The chosen turbulence models were comparable in their accuracy, however,

the LES simulations were able to account for fluctuations in the velocity under increased

computational effort. A second validation step will be to analyze the accuracy of the air-

phase behavior.

2.5 Conclusions

In this paper we presented the research as it is carried out in the sewer interfaces group

within the DFG Research Training Group Urban Water Interfaces. Different sub-projects fo-

cus on a deeper understanding of H2S formations at and across the wastewater - biofilm, the

wastewater - sewer atmosphere and the sewer atmosphere - biofilm - sewer wall (concrete)

interface.

The first project conducts experimental investigations on a pilot plant in order to gain an

enhanced understanding of sulfide formations in pressure mains as well as on investigating

the influence of flow interruptions caused by gross solids and drop structures. Further, the

effects of different countermeasures at hotspots are being analyzed. In another project, this

research pilot plant is used in order to investigate BSC in sewers. The main goal is to in-

vestigate different laboratory performance tests and to compare them among each other as

well as to concrete probes installed in the pilot plant. In the third project, a CFD simulation

model is being developed in order to describe H2S emissions across the wastewater - sewer

atmosphere interface. First results show, that the model is able to accurately describe the

behavior of the water phase. Therefore, the model can be used in order to estimate different

hydraulic properties within the sewer pilot plant. On the other hand, results concerning

H2S formations obtained by the sewer pilot plant can be used to validate the extended CFD

model.

The added value of the work carried out by this group lies in the detailed insight into

H2S formation processes and sensitivity analyses that can be carried out. The influence

of gross solids or drop structures on the H2S release can be investigated on a pilot plant

scale. The validated CFD model can be used to investigate the effect of measures such as

i.e. changed discharges or ventilation measures on H2S emissions. These findings can be

integrated into existing model approaches in order to improve predictions. The develop-

ment of a CFD model can generalize the findings of the investigations in the pilot plant to

different conditions. The detailed analysis of different laboratory performance tests helps

to generalize the development of acid-resistant concrete.

Acknowledgements

The funding provided by the German Research Foundation (DFG) within the Research

Training Group Urban Water Interfaces is gratefully acknowledged.

Parts of the simulations were computed on the supercomputers of Norddeutscher Verbund

für Hoch- und Höchstleistungsrechnen in Berlin.

The sewer pilot plant is provided by the Berliner Wasserbetriebe (BWB). We also thank the

Chapter 2: Research on interfaces in sewer systems

Kompetenzzentrum Wasser Berlin for providing s::can probes, the Kiwa (Berlin) for provid-

ing the apparatus for the MPA Berlin-Brandenburg test, and the Magnel Laboratory of the

University Ghent for providing the TAP apparatus.

Chapter 3

Validation of hydrodynamic

two-phase simulations

This study was published in Progress in Computational Fluid Dynamics (Inderscience) as:

Teuber, K., Broecker, T., Bayón, A., Nützmann, G. & Hinkelmann, R.: CFD-modelling

of free-surface flows in closed conduits, Progress in Computational Fluid Dynamics, 19 (6),

368-380, 2019; doi:10.1504/PCFD.2019.103266.

©Inderscience 2019. The definitive peer-reviewed and edited version of this article is pub-

lished in doi:10.1504/PCFD.2019.103266 and available at www.inderscienceonline.com.

Reproduced from Teuber et al. 2019 Progress in Computational Fluid Dynamics 19(6)

368-380, with permission from the copyright holders, Inderscience.

The test cases’ setups are listed in Appendix B (single-phase flow over ground sill: B.1,

free-surface flow over ground hill B.2, complex sewer: B.3).

Abstract

Computational Fluid Dynamics (CFD) is gaining an increasing importance in the field of hy-

draulic engineering. This publication presents different application examples of a two-phase

approach as implemented in the open source software OpenFOAM. The chosen approach

is based on the volume of fluid method focusing on the simulation of flow in closed con-

duits. Three examples are presented: single-phase flow over a ground sill and free surface

flow over a hill as well as complex free surface flow in a sewer model. The first example

compares the results of different RANS turbulence models with experimental results. The

results of the second example are compared with an analytical solution. In the last exam-

ple the behaviour of the free surface flow is compared with the results of a model test and

existing simulations using a simplified, open channel geometry for the closed conduit. For

the examples analysed, the two-phase approach provides stable and reliable results.

3.1 Introduction

Numerical multiphase solvers have gained an increasing importance in the field of hydraulic

research. In this paper, the open field operation and manipulation (OpenFOAM) solver in-

terFoam (Rusche, 2003) which is based on a VOF interface capturing approach (Section 3.2.2)

will be analysed.

Chapter 3: Validation of hydrodynamic two-phase simulations

In a number of recent publications, such as Bayón et al. (2015), Bayón and Lopez-Jimenez

(2015), Schulze and Thorenz (2014) and Thorenz and Strybny (2012) the interFoam solver

was used to analyse complex hydraulic cases in open channels. Possible application areas

are inland waters (Thorenz and Strybny, 2012), coastal areas (Higuera et al., 2013) as well

as sewer systems (Bayón et al., 2015). The previously mentioned publications simulated the

water phase behaviour and used the two-phase approach in order to display the movement

of the free water surface.

Bayón and Lopez-Jimenez (2015) used the interFoam solver to analyse a hydraulic jump

in a rectangular channel with smooth walls. A comparison of different variables with ex-

perimental results led to the conclusion that the model can be applied to real-life cases of

designing hydraulic structures. In Bayón et al. (2015), the test case of the hydraulic jump

was extended to a more complicated surrounding consisting of an existing sewer stretch

with different hydraulic structures such as weirs, quickly-varying shapes, macro-roughness

elements, fast and slow flow regimes as well as hydraulic jumps. The results were compared

with experimental results from a 1:20 scale model and also showed a good agreement. In

our last example we will use this sewer geometry. Due to stability problems, the sewer has

been simulated as an open channel with an atmospheric top boundary.

The novelty applied to the existing setup is the closed system setup by defining a closed

(wall) top boundary instead of an open atmospheric top boundary as it has been used in

Bayón et al. (2015) and a different definition of the outlet boundary condition.

The importance of the correct choice of a turbulence model was shown in Bayón et al.

(2015) and Bayón and Lopez-Jimenez (2015). A sensitivity analysis in Bayón and Lopez-

Jimenez (2015) showed a good performance of the Standard k-eturbulence model.

A number of pipe flow simulations have already been performed applying the VOF ap-

proach in closed ducts. A few of them shall be outlined in the following and their certain

characteristics shall be pointed out.

Shuard et al. (2016) compared simulations of two-phase flow in a circular pipe using the

interFoam solver with results of a mechanistic model. Similar simulations have been carried

out by Thaker and Banerjee (2013) analysing the transition between different flow regimes

as well as the development of flow regimes. The results were compared to experimental

measurements using similar boundary conditions. The boundary conditions that have been

used by Shuard et al. (2016) have been used for a different area of application by Kinyua

et al. (2016). Here, tubular anaerobic digesters have been modelled by including tracer simu-

lations. The simulations carried out by Shuard et al. (2016), Thaker and Banerjee (2013) and

Kinyua et al. (2016) analysed flow behaviour in closed ducts. By applying a constant value

as an outlet pressure boundary condition, a free outflow out of the domain was achieved,

assuming that the water level is not impounded by any downstream water level.

This short literature review shows that CFD models are an upcoming topic in hydraulic

engineering but the simulation results strongly depend on factors such as boundary condi-

tions, mesh quality and turbulence models. Existing validation shows a good accuracy of

the interFoam solver with experimental data for open systems with atmospheric top bound-

aries and for closed ducts with free outflows, however, there is a lack of research on the

behaviour of free surface flows in closed ducts with a defined water level.

In this paper we want to investigate simulations that describe free surface flow in closed

Chapter 3: Validation of hydrodynamic two-phase simulations

pipes where at the same time a certain outlet water level is desired (see also Bayón et al.

(2015)). This problem occurs in cases where the area of interest comprises a large system,

i.e. a complete sewer system. Due to limitations of computational resources it is usually not

feasible to simulate the whole system, however, it is possible that the water level at the end

of the stretch that is being considered in the simulation influences the water level within

the model domain. Therefore it is not possible to use a free outflow as it has been used in

previous studies. The literature analysis shows that these systems have not been addressed

in previous publications.

The simulation of closed systems is interesting for different areas of application, i.e. for

the modelling of in-sewer processes (Edwini-Bonsu and Steffler, 2004, Gessner et al., 2014,

Hvitved-Jacobsen et al., 2013, Rootsey et al., 2012). A detailed analysis of the air phase be-

haviour will be subject to future research. To ensure a correct behaviour of the water phase,

a thorough validation of flow simulations in closed conduits using the solver interFoam is

performed in three steps. First, the behaviour of the water phase is analysed by comparing

the numerical data of flow over a hill for single-phase flow with experimental results by

Almeida et al. (1993). The results of the first test case show us which turbulence model is

most suitable for the problems analysed. Then, the different outlet boundary conditions

are compared using a simple two-phase simulation. The interface behaviour is analysed by

comparing the simulated results of water level drawdown due to a ground sill with ana-

lytical results based on continuity and Bernoulli’s equation. In a last step, the stability and

the accuracy in describing a practical test case is checked by simulating a complex sewer

geometry (Bayón et al., 2015).

The paper is organized as follows: Section 3.2 gives an overview over the methods and

materials. The computational test cases are presented in Section 3.3. The validation of the

different parts is done as follows: single phase water flow over a ground sill, free surface

flow over a hill, complex free surface flow in a sewer model. In Section 3.4, the results of

the validation are summarized.

3.2 Methods and materials

3.2.1 Geometry and mesh

Unstructured meshes created in the open source mesh generation tool gmsh were used for

the single phase validation case (see Section 3.3.1) and for the free surface flow over a hill

(see Section 3.3.2). The complex sewer case of Section 3.3.3 was discretized with a structured

grid using the OpenFOAM utility snappyHexMesh. The detailed geometry of the different

test cases will be outlined in the subsequent sections.

3.2.2 Numerical model

The Open Source CFD software OpenFOAM (Open Field Operation and Manipulation) ver-

sion 2.4.0 was used to simulate different test cases. Single and two-phase flow is calculated

by using the solver interFoam based on a VOF approach. interFoam is a multiphase solver

for immiscible and isothermal fluids that solves the three-dimensional Navier-Stokes equa-

tions using the finite-volume-method in space and the finite-differences-method in time.

OpenFOAM allows parallel computations on a theoretically unlimited number of processor

cores (Keough, 2014).

Chapter 3: Validation of hydrodynamic two-phase simulations

The conservation of mass (Equation 3.1) and momentum (Equation 3.2) for incompressible

flow can be written as:

∇ · U=0 (3.1)

∂ρU

∂t+U· ∇ρU=−∇p+µ∆U+ρg(3.2)

The viscosity term µreferred to in Equation 3.2 contains the physical viscosity µphys as

well as the turbulent viscosity µturb which will be obtained by a turbulence model:

µ=µphys +µturb (3.3)

The VOF method used in the interFoam solver uses a specific pressure formulation where

prgh is a modified pressure which is used in order to avoid the occurrence of steep pressure

gradients caused by hydrostatic effects. In the following, prgh will be referred to as:

prgh =p−ρ·g·h(3.4)

The two immiscible fluids liquid and gas “are considered as one effective fluid through-

out the domain, the physical properties of which are calculated as weighted averages based

on the distribution of the liquid volume fraction, thus being equal to the properties of each

fluid in their corresponding occupied regions and varying only across the interface"(Berberovi´c

et al., 2009), leading to a definition of ρand µas follows:

ρ=αρwater +ρair (1−α)(3.5)

µ=αµwater +µair (1−α)(3.6)

The αtransport equation - also known as VOF equation - is used with an advanced formu-

lation which can be considered as an evolution equation for the phase fraction α(Berberovi´c

et al., 2009):

∂α

∂t+∇ · (αU)+∇ · ((1−α)Urα)=0 (3.7)

“where Ur=Ug−Ulis the vector of relative velocity, designated as the “compression

velocity""(Berberovi´c et al., 2009). A detailed derivation of the equations can be found in

Damián (2012).

In the equations above, Urepresents the ensemble averaged velocity field shared by the two

fluids throughout the flow domain [m/s]; ρis the density [kg/m3]; t is time [s]; p is pressure

[Pa]; µphys and µturb are the physical and turbulent viscosity [Ns/m2]; gis the acceleration

vector due to gravity [m/s2]; his a spatial position vector [m]; αis the volume fraction or

indicator function [-]; Uris the relative velocity between the phases [m/s]; the subscripts l

and gdenote different fluids liquid (water) and gas (air). The indicator function αis defined

as:

α=









1 for a point inside fluid water

0<α<1 for a point in the transitional region

0 for a point inside fluid air

The water surface is defined as the transition area where α=0.5. The solver can be used

as well for single-phase flow simulations. The volume fraction αis then 1 and constant over

the whole domain.

Chapter 3: Validation of hydrodynamic two-phase simulations

Since the relative velocity cannot be computed directly from the one-fluid formulation

in OpenFOAM, the numerical implementation of the relative velocity is as follows (Cifani

et al., 2016):

Ur=nαmin[Cα

|φ|

|Sα|,max(|φ|

|Sα|)] (3.8)

where nαis the normal vector of the cell surface, φis the mass flux, Sαis the cell surface area,

Cαis an adjustable coefficient on which the level of compression depends. The maximum

of Uris bounded to the maximum face velocity in the flow field and the direction is aligned

with nα(Cifani et al., 2016, Hoang et al., 2013).

Equations 3.1 and 3.2 are solved by using the PIMPLE algorithm. As outlined in Bayón

et al. (2015), pressure-velocity coupling is done by combining the inner corrector loops of

the PISO (Pressure Implicit with Splitting of Operators) algorithm (Issa, 1986) with outer

corrector loops of the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algo-

rithm (Patankar and Spalding, 1983). Detailed information on the algorithms can be found

in Jasak (1996). As advection schemes, a total variation diminishing scheme (interGamma)

(Jasak, 1996) combined with a flux corrected transport approach MULES (Multidimensional

Limiter for Explicit Solutions) (Damián, 2013) is used.

Simulations were carried out until a quasi-steady state was reached. This quasi-steady

state has been determined by sampling variables of interest such as the pressure and the

velocity at different points in the water and air phase.

In the first example, different RANS turbulence models, namely the Standard k-ε(Laun-

der and Sharma, 1974), k-ω(Wilcox, 1988) and k-ωShear Stress Transport (SST) (Menter,

1993, 1994) are compared. Wall functions have been used to describe the near-wall turbulent

flow. It was ensured that the y+-values for each simulation were in a range between 30 and

300 to guarantee the applicability of the wall functions.

3.2.3 Boundary conditions

Different boundary conditions were defined for single and two-phase simulations. The def-

inition of sidewalls was similar for single and two-phase simulations but varied depending

on the fact whether two or three-dimensional model setups were used. For two-dimensional

test cases, so-called empty boundary conditions, which are implemented in OpenFOAM to

describe sidewalls of two-dimensional geometries, were used. They declare that the side-

walls do not constitute solution directions within the defined domain. The sidewalls of

three-dimensional testcases were defined by no-slip conditions. The detailed set of bound-

ary conditions is explained in the chapters containing the respective test cases. However,

a short description about the definition of outlet boundary conditions in free surface CFD

models shall be given in the following.

In general, the definition of the outlet boundary condition in OpenFOAM using the in-

terFoam solver differs from the definition in shallow-water models. In these models, often

a water level is fixed at the outlet of the domain. In OpenFOAM, such a definition of a

constant water level is not possible, therefore alternatives have to be looked for. Possibilities

described in recent publications are to fix the pressure and the phase fraction value at the

outlet (Thorenz and Strybny, 2012) or to specify a velocity profile at the outlet. A third

option which circumvents the direct definition of a boundary condition in order to ensure a

certain water level is to set an atmospheric pressure condition at the outlet and integrate a

Chapter 3: Validation of hydrodynamic two-phase simulations

weir in the geometry close to the outlet (Bayón and Lopez-Jimenez, 2015).

In this paper, a certain water level h at the outlet has been obtained in two different ways

which will be referred to in the respective sections:

• Pressure boundary condition: prgh is defined as a stepwise function at the outlet. The

built-in function setFields is used to overwrite uniform values at the outlet which is

initially defined using a calculated boundary condition:

prgh =(0 if z <h

max(ρgh)if z ≥h

After the writing process, the type of boundary condition is changed to fixedValue. The

velocity at the outlet is defined as inletOutlet condition which applies a null Neumann

boundary condition in case of positive flux out of the domain and a user-specified

fixed velocity in case of negative flux into the domain. In case of negative flux, the

velocity is set to U = (0 0 0). The phase fraction value αis defined as null Neumann

boundary condition.

• Integration of weir: a weir structure in close proximity to the outlet is used as a means

to maintain the water level in the domain. At the outlet, a constant pressure is de-

fined using a Dirichlet boundary condition, the velocity is defined using the inletOut-

let condition. The phase fraction value αis specified using a null Neumann boundary

condition.

In order to avoid stability problems, pressure boundary conditions need an accurate defini-

tion of initial conditions, i.e. the initial water level, flow velocity and pressure. Therefore,

for simulations using this set of boundary conditions, suitable initial conditions have to be

specified. They will be further outlined in the subsequent sections. Depending on the RANS

turbulence model, an additional definition of boundary conditions for turbulent properties

such as the turbulent kinetic energy k and the turbulent dissipation εfor k-εor the turbulent

kinetic energy k and the specific dissipation ωfor k-ωmodels is needed depending on the

chosen model. For the boundary condition, an initial guess of the turbulent properties has

to be defined. During the simulation, these values are rewritten according to the simulation

results.

3.3 Applications

3.3.1 Single phase water flow over a ground sill

In the first example, single-phase water flow over a two-dimensional model hill in a duct has

been simulated and compared with experimental results obtained by Almeida et al. (1993).

Experiments were conducted in a two-dimensional duct with a height of hmax =0.17m

bounded by an upper and lower wall with a polynomial-shaped obstacle on its bottom (see

Figure 3.1). The mean centreline velocity at the inlet amounted to U0=2.147m/scaus-

ing pressurized flow throughout the domain. Velocity profiles of the simulations in x and

y-direction at four (velocities in x-direction) respectively two (velocities in y-direction) dif-

ferent locations were compared to the measurements: x−03 = -0.30 m(in front of the hill),

x00 = 0.00 m(top of the hill), x01 = 0.03 m(end of the hill) and x02 = 0.05 m(middle of recir-

culation zone). Figure 3.1 shows the model domain including the mesh that has been used

Chapter 3: Validation of hydrodynamic two-phase simulations

Figure 3.1: Domain of single-phase, two-dimensional hill flow including measurement loca-

tions, three selected indicator locations for grid convergence study and computational mesh

(Davroux et al., 1995).

and the different measurement stations. Simulations have been carried out using different

RANS turbulence models (Standard k-ε, k-ω, and k−ωSST).

A grid convergence study following Celik et al. (2008) has been conducted in order to

ensure mesh independence. Five different meshes with different cell sizes were tested using

the flow velocity in x-direction Uxin 15 locations in the nearfield of the hill structure. The

cell sizes and total number of cells of the analysed meshes were 0.0293 m (69,326 cells),

0.0346 m (31,102 cells), 0.0416 m (15,319 cells), 0.0495 m (7,338 cells), and 0.0581 m (3,855

cells), leading to a global refinement ratio between meshes of 2, which is above the minimum

recommended value of 1.3 (Celik et al., 2008). Figure 3.2 a) shows the resulting velocities

of three indicator locations φ4,φ7and φ10. The position of these three points is displayed

in Figure 3.1. Convergence was reached between the second finest and finest mesh. The

computed RMSE for the different refinement steps (i.e. step 1: between coarsest and sec-

ond coarsest mesh) is displayed in Figure 3.2 b). It shows a significant decline between the

third and fourth refinement step to 5.5e-5 m/s. These results lead to the conclusion that

the second finest mesh analysed has reached grid independence, so all subsequent analysis

is conducted on the mesh of 0.0346 m element size and a total number of 31,102 cells. For

the selected mesh, the average apparent order as defined by Celik et al. (2008) is 1.9, very

close to the model formal order. A variable time step in dependence of the Courant number

has been used, which converged to approximately ∆t = 0.0007 s. The model convergence

to a quasi-steady state can be observed after a few seconds. All simulations have been car-

ried out with a simulation time of 10 seconds. For all simulations, Intel Xeon IvyBridge

E5-2695v2 cores have been selected. Using parallel computations on 16 cores led to com-

putation times of approximately 34 minutes. In order to reach a good agreement with the

experimental setup, a velocity profile using the experimental data has been imposed as an

Chapter 3: Validation of hydrodynamic two-phase simulations

inlet boundary condition. The upper and lower walls were specified with no-slip conditions.

At the sidewalls, empty boundary conditions have been applied. At the outlet, the water

level was defined using a constant prgh-value defined by the maximum hydrostatic pressure:

prgh =max(ρgh). As the initial condition, the domain is filled with water but no velocities

are predefined. The velocity profile therefore develops during the simulation.

Simulated and measured velocity profiles are compared in Figures 3.3 (Ux) and 3.4 (Uy).

On the horizontal axis, the velocities are displayed in relation to the average flow velocity

U0. The height on the vertical axis is displayed in relation to the channel height hmax. The

imposed velocity profile at the inlet of the domain differs slightly from the experimental

results of the freestream profile (Figure 3.3(a)). Figures 3.3(b) to 3.3(d) show the resulting

velocity profiles in the reach of the hill structure and enable a qualitative analysis of the

accuracy of the different turbulence models. Small deviations can be found between simu-

lations and experimental results. The resulting root mean square error (RMSE) normalized

by U0for the different locations is listed in Table 3.1. The error has been calculated using

linear interpolation between the cell values of the simulation results. Overall, the Standard

k-eturbulence model leads to slightly better results in all locations than the other RANS

models analysed. The velocities in y-direction (Uy) show similar results (Figure 3.4). Table

3.1 shows that the error between the experimental results and the simulated cases is smaller

than the error of the Ux-velocities. In comparison to the velocities in x-direction, the veloc-

ities in y-direction are very small, being more challenging for numerical simulations and

measuring.

It can be concluded, that the chosen model set-up is able to describe a relatively complex

hydraulic test case for a single-phase flow problem appropriately and all turbulence models

show a good accuracy. Due to the good performance in this case as well as in the cases

referred to in Section 3.1, the Standard k-εturbulence model will be used in the following

cases.

Table 3.1: Water phase validation - RMSE/U0[-] between experimental results and simula-

tion

Observation point k-ek-ωk-ωSST

x−03 0.0515 0.0515 0.0515

x00 0.0232 0.0307 0.0237

x01 0.0595 0.1077 0.0681

x02 0.0361 0.0637 0.0361

x00 0.0116 0.0131 0.0130

x01 0.0156 0.0309 0.0193

Chapter 3: Validation of hydrodynamic two-phase simulations

0.85

0.9

0.95

1.05

1.1

1.15

10000

20000

30000

40000

50000

60000

70000

U/U0 [-]

No. of cells [-]

a) Velocities

φ4

φ7

φ10

0.001

0.002

0.003

0.004

RMSE/U0 [-]

Refinement step [-]

b) RMSE

Figure 3.2: Mesh sensitivity analysis, (a) Comparison of simulated indicator values φ(flow

velocities) for different grid sizes in three points within the domain; (b) RMSE between

different refinement steps over all 15 indicator values

Chapter 3: Validation of hydrodynamic two-phase simulations

0 0.2 0.4 0.6 0.8 1

y/hmax [-]

Ux/U0 [-]

k-epsilon

k-omega

k-omega SST

experimental data

(a)Uxvelocities at x−03 =−0.30m

0 0.2 0.4 0.6 0.8 1 1.2 1.4

y/hmax [-]

Ux/U0 [-]

k-epsilon

k-omega

k-omega SST

experimental data

(b)Uxvelocities at x00 =0.00m

-0.2 0 0.2 0.4 0.6 0.8 1 1.2

y/hmax [-]

Ux/U0 [-]

k-epsilon

k-omega

k-omega SST

experimental data

(c)Uxvelocities at x01 =0.03m

-0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2

y/hmax [-]

Ux/U0 [-]

k-epsilon

k-omega

k-omega SST

experimental data

(d)Uxvelocities at x02 =0.05m

Figure 3.3: Comparison of simulated velocities Uxusing different turbulence models

with experimental data for flow over a two-dimensional hill, (a) Uxvelocities at x−03 =

-0.30 m, (b) Uxvelocities at x00 = 0.00 m, (c) Uxvelocities at x01 = 0.03 m (d) Uxvelocities

at x02 = 0.05 m

Chapter 3: Validation of hydrodynamic two-phase simulations

0 0.02 0.04 0.06 0.08 0.1 0.12

y/hmax [-]

Uy/U0 [-]

k-epsilon

k-omega

k-omega SST

experimental data

(a)Uyvelocities at x00 =0.00m

-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08

y/hmax [-]

Uy/U0 [-]

k-epsilon

k-omega

k-omega SST

experimental data

(b)Uyvelocities x01 =0.03m

Figure 3.4: Comparison of simulated velocities Uyusing different turbulence models

with experimental data for flow over a two-dimensional hill, (a) Uyvelocities at x00 =

0.00 m, (b) Uyvelocities at x03 = 0.03 m

Chapter 3: Validation of hydrodynamic two-phase simulations

3.3.2 Free surface flow over a hill

After having analysed the water phase behaviour, a two-phase flow has been simulated in

order to investigate the water surface behaviour. The model setup was slightly changed.

A free surface flow was modelled in a two-dimensional model domain. Again, the system

was bounded by upper and lower walls. Under subcritical flow conditions the water level

drawdown was analysed and compared to the analytically computed water level deviation

by using continuity and Bernoulli’s equation. The computation time was approximately 1

h for parallel computation on 16 cores. In the following, three cases are investigated (Table

3.2). The two-dimensional model domains of cases 1 and 2 consisted of 68,542 cells with a

minimum cell length of 0.0024 m and a maximum cell length of 0.15 m. The domain of case

3 was discretized into 175,762 cells. The cell length ranged between 0.0035 m and 0.144 m.

A time step of ∆t=0.001 s was chosen.

The boundary conditions were defined as follows: The inlet was divided in two parts, con-

taining an inlet for the air phase and an inlet for the water phase. The air-inlet was specified

with a fixed total pressure, the water-inlet was specified with a fixed discharge. Since the

analysed test cases are two-dimensional, sidewalls were defined as empty boundaries. The

upper and lower walls were specified using no-slip conditions. The outlet was defined using

the pressure boundary condition mentioned in Section 3.2.3. As initial conditions for the

pressure and velocity, the freestream velocities (v1) and water levels (h1) as listed in Table

3.2 have been defined.

In the following, the results of different two-phase flow simulations for three different

model setups are presented. Variations were made concerning the discharge and water

level. A deeper analysis of factors such as the structure of the sill and the flow regime can

be found in Teuber et al. (2016).

In the three cases presented here, the water depth and flow velocity have been varied us-

ing the two-dimensional test case with a 30◦angular structure of the ground sill (Figure 3.5),

the flow regime was kept subcritical. Information about the setup including the analytically

calculated water level drawdown and the numerically simulated water level drawdown are

listed in Table 3.2, where v1is the approaching velocity and h1the water depth in front of

the sill. The domain of case 1 is displayed in Figure 3.5. The results show a similar trend

between the simulated results and the analytical solution. The reason for the higher draw-

down obtained by the numerical solution are additional energy losses due to the structure

of the sill that can be accounted for by the CFD simulation whereas the analytical solution

using continuity and Bernoulli’s equation leads to a one-dimensional solution that neglects

single energy losses caused by the structure of the sill. In addition, Teuber et al. (2016)

simulated a strictly supercritical setup leading to a rise of the water level above the sill.

In agreement with the previous observations, the water level deviation coincided with the

analytical solution.

It can be concluded that all the cases investigated show that the numerically calculated

deviation of the water level reasonably coincides with the analytical solution.

Chapter 3: Validation of hydrodynamic two-phase simulations

Figure 3.5: Free surface flow over a hill - Model domain of case 1

Table 3.2: Water-air-interface - Properties of the different test cases

Class Case 1 Case 2 Case 3

Length of domain 25 m 25m 35 m

Height of domain 2 m 2m 6 m

h11.0 m 1.0 m 3.0 m

v11.0 m/s 1.25 m/s 3.0 m/s

∆z0.2 m 0.2 m 0.2 m

∆h, analytical 0.036 m 0.070 m 0.110 m

∆h, numerical 0.042 m 0.090 m 0.140 m

Chapter 3: Validation of hydrodynamic two-phase simulations

3.3.3 Complex free surface flow in a sewer model

As a third step, the overall performance of the chosen model setup was tested using the

complex sewer geometry of Bayón et al. (2015). Data from an existing CFD model as well

as measured data gained from a 1:20 scale model has been compared to simulations using

a model setup describing a closed duct bounded by upper and lower walls as well as side-

walls. The simulation results reported in Bayón et al. (2015) were obtained using a similar

configuration as presented in this paper, also implemented in OpenFOAM with the differ-

ence in the definition of a simplified top and outlet boundary condition. The Standard k-e

model was used as closure to the turbulent stresses in both models.

The model domain describes a planned sewer stretch which has been developed due to a

necessary diversion of an existing sewer stretch in València, Spain (Figure 3.6). The chosen

stretch has a length of approx. 95 m and consists of an initial stretch, which was modelled

in order to achieve fully developed flow conditions. After the initial stretch, a 45◦constant

radius curve and a straight stabilisation stretch with a length of 10 m follow. Having passed

the curve, the water then reaches a spillway followed by a stilling basin which is designed in

order to control the occurrence of critical flow conditions, obtain a smooth change of slope

between the existing and new sewer stretch and to expand the channel width from 6.0 m to

7.5 m. An adjacent transition zone leads the flow into the ovoid channel with a length of 20

m which is the last part of the section. In order to force the flow back to a subcritical state

before reentering into the ovoid section, macro-roughness elements are placed in the stilling

basin and a hydraulic jump is forced. Due to its different features as well as the different

flow regimes, the chosen geometry leads to a model domain that is highly complex in terms

of hydraulic behaviour.

The innovation of the simulations presented in this paper compared to the existing CFD

simulations carried out by Bayón et al. (2015) is that here the whole domain was considered

as a closed duct leading to a realistic top boundary condition and another outflow bound-

ary condition to ensure stability. As mentioned before, the same case has been previously

simulated by Bayón et al. (2015) with a simplification for the closed conduit, i.e. an atmo-

spheric top boundary in order to avoid stability problems. In this paper, the top boundary

was not described with an atmospheric condition but as a wall with no-slip condition. The

lower walls as well as the sidewalls were defined using no-slip conditions in both cases. The

remaining boundary conditions were defined as follows by Bayón et al. (2015): The outlet

water level was set to 5.01 m and the discharge was defined as 100 m3/s. Accordingly, the

pressure and velocity outlet boundary condition were defined. At the outlet, the velocity

was defined to a constant value which forced the water level to a certain height. Within the

domain, an initial water level and flow velocity were defined according to these values.

In this paper, a weir structure has been added in close proximity to the outlet in order

to obtain the desired water level. This setup is very robust concerning initial conditions and

does not necessarily need a correct initialisation of water level and flow velocity. Stable sim-

ulations can even be obtained under dry initial conditions. In the following, the simulation

results using the weir outlet boundary condition will be compared to experimental mea-

surements as well as to existing CFD simulations as they have been carried out by Bayón

et al. (2015).

In Bayón et al. (2015) the accuracy of different mesh sizes ranging from 0.0861 m to 0.1236

m were compared. Mesh convergence has been reached for a mesh size of 0.103 m. In this

paper, simulations were performed on a grid with the same resolution, leading to an overall

number of 3,029,223 cells. The simulations demanded a computation time of 336 h on 96

Chapter 3: Validation of hydrodynamic two-phase simulations

cores. Simulations were performed with an adjustable time step in relation to the Courant

number which converged against ∆t = 0.0002 s. The simulations run stable on different

grid sizes, even if the domain is dry in the beginning of the simulations.

In order to compare the accuracy of the results for this complex test case, the water

level distribution at four different cross sections has been compared. Figure 3.7 shows the

location of the different monitoring points (1 - beginning of curve, 2 - end of the curve, 3 -

stabilisation stretch, 4 - spillway crest). The resulting water surface profiles in the different

cross sections of the two CFD models and the physical model are displayed in Figure 3.8.

The results obtained in this paper are called “new CFD model" in the following, while the

results of Bayón et al. (2015) are called “existing CFD model". A qualitative comparison

shows an overall good agreement of the new CFD model with the existing model. In the

beginning of the curve, the water level is closer to the physical model than the existing CFD

model. In the remaining monitoring points, the results of the new CFD model are similar to

the results of the existing model. This leads to the conclusion that the new CFD model with

a realistic top boundary condition is as accurate as the existing model using a simplified

boundary condition (atmospheric top boundary).

The results show that the closed model setup consisting of a weir structure as outlet

boundary condition, a high-resolution mesh and the Standard k-εturbulence model is stable

concerning initially dry conditions, flow transitions from super to subcritical flow as well as

high filling ratios and leads to reliable results.

Chapter 3: Validation of hydrodynamic two-phase simulations

Figure 3.6: Complex sewer geometry - Geometry of the sewer stretch (Bayón et al., 2015)

Chapter 3: Validation of hydrodynamic two-phase simulations

Figure 3.7: Complex sewer geometry - Overview of monitoring points for water surface

profiles

Chapter 3: Validation of hydrodynamic two-phase simulations

Height z [m]

1: Curve beginning 2: Curve end

-3 -2 -1 0 1 2 3

Height z [m]

Width w [m]

3: Stabilisation stretch

-3 -2 -1 0 1 2 3

Width w [m]

4: Spillway crest

Physical model

CFD-model (existing)

CFD-model with weir (new)

Figure 3.8: Complex sewer geometry - Comparison of water surface profile at different

locations (existing CFD model - results obtained by Bayón et al. (2015); new CFD model -

results obtained in this paper)

Chapter 3: Validation of hydrodynamic two-phase simulations

3.4 Conclusions

The aim of this study was to investigate OpenFOAM’s two-phase solver interFoam con-

cerning its ability to accurately describe complex hydraulic test cases in closed ducts when

a certain outlet water level is given. This paper presents three different application cases

where the interFoam solver has been used. Another aim was to evaluate the accuracy and

suitability of different RANS turbulence models.

In a first step, a single-phase water flow over a hill has been simulated and different

RANS turbulence models have been validated concerning their accuracy in describing the

eddy structure behind the hill. The flow velocities measured experimentally were success-

fully reproduced by the different models. In the following, the Standard k-εturbulence

model has been used for test cases describing the water surface behaviour under quasi-

steady state conditions.

As second test case, free surface flow over a hill has been simulated. Different pa-

rameters such as the discharge and the water level have been varied and the numerically

computed results were compared to an analytical solution obtained by using continuity and

Bernoulli’s equation. The simulations were able to predict a similar trend of the water level

drawdown compared to the analytical solution. Due to additional energy losses caused by

the specific structure of the sill, the analytically calculated drawdown is smaller than the

CFD simulations but reasonably coincides for different simulations. These results can be

seen to support the adequacy of the applied numerical modelling.

In a third step, complex free surface flow in a sewer model has been analysed in order

to investigate the stability of the simulations using the setup of a complex geometry. A

comparison with results of an existing CFD model, which used a simplified top boundary

condition (atmospheric instead of closed) and another outflow boundary condition, as well

as measured results from a 1:20 scale model showed a good agreement of this setup with

the existing results.

Summing up, the VOF approach implemented in OpenFOAM is capable of describ-

ing complex two-phase flows in closed ducts. Small differences in the accuracy have been

observed depending on the chosen turbulence models. The possibility of simulating even

complex closed systems opens up the chance to describe complex phenomena such as odour

and corrosion caused by hydrogen sulphide in concrete sewers. Future research will focus

on the detailed behaviour of the air phase.

Acknowledgement

Parts of the simulations were computed on the supercomputers of Norddeutscher Verbund

für Hoch- und Höchstleistungsrechnen in Berlin.

Funding

The funding provided by the German Research Foundation (DFG) within the Research

Training Group “Urban Water Interfaces" (GRK 2032) is gratefully acknowledged.

Chapter 4

Parameter study on hydrodynamic

two-phase simulations of flow over a

ground sill

This study was published as:

Teuber, K., Broecker, T., Barjenbruch, M. & Hinkelmann, R.: High-resolution numerical

analysis of flow over a ground sill using OpenFOAM, Proceedings of the 12th International

Conference on Hydroscience & Engineering (ICHE), Tainan, Taiwan, 2016.

This is the postprint version of the article.

The test cases’ setups are listed in Appendix B (single-phase flow over ground sill: B.1,

free-surface flow over ground hill B.2).

Abstract

In the field of hydraulic engineering, attention towards Computational Fluid Dynamics

(CFD) has increased within the last years. In this study, flow over a two-dimensional ground

sill is simulated and analyzed using the open source model OpenFOAM. Single-phase flow

simulations are compared to experimental results obtained by Almeida et al. (1993) and

two-phase flow simulations are compared to analytical solutions by using Bernoulli’s and

continuity equation. The results show that the model is capable of simulating such hydraulic

testcases. Different RANS and LES simulations were found to reproduce the analyzed flow

behavior well.

4.1 Introduction

Within the last years, Computational Fluid Dynamics (CFD) has gained importance in the

field of hydraulic engineering. Several publications, such as Bayón et al. (2015), Bayón

and Lopez-Jimenez (2015), Schulze and Thorenz (2014), Thorenz and Strybny (2012) have

investigated complex hydraulic testcases such as hydraulic jumps, and filling and emptying

of locks using the open source model OpenFOAM. In these publications the Volume of

Fluid (VoF) approach for two-phase flows has been used in order to describe free surface

flows. But the model also offers the possibility to simulate cases where two phases are of

importance. One application is the simulation of in-sewer processes (Edwini-Bonsu and

Chapter 4: Parameter study on hydrodynamic two-phase simulations

Steffler, 2004, Gessner et al., 2014, Hvitved-Jacobsen et al., 2013). In this work, a first step

of the validation process regarding two-phase flows in closed ducts such as sewer pipes is

made. The model is used to simulate flow over a two-dimensional ground sill. Validation

is first performed by comparing the results of single-phase flow with measurements by

Almeida et al. (1993). Two-phase flow simulations are performed for different two- and

three-dimensional model setups, variations are made concerning the structure of the sill,

discharge, water level and the flow regime.

4.2 Computational framework

4.2.1 Numerical model

Surface water flow is calculated by using the two-phase flow solver interFoam based on

a volume of fluid (VoF) approach for one- and two-phase flows. Both phases are consid-

ered as one fluid with rapidly changing fluid properties, therefore one set of Navier-Stokes

equations is solved. The phases are distinguished by an additional transport equation for

the volume fraction which is used as a marker to describe the distribution of the phases

throughout the domain. The equations can be formulated as follows (Rusche, 2003):

Mass conservation equation:

∇ · −→

U=0 (4.1)

Momentum conservation equation:

∂ρ−→

∂t+∇ · ρ−→

U−→

U=−∇prgh +∇ · µ∆−→

U+∇−→

U∇µ−−→

g·−→

x∇ρ(4.2)

Where prgh is the static pressure minus hydrostatic pressure:

prgh =p−ρ·g·h(4.3)

Volume of Fluid equation:

∂α

∂t+∇ · α−→

U+∇ · ((1−α)Urα)=0 (4.4)

with the following parameters:

ρ=αρw+ρa(1−α)(4.5)

µ=αµw+µa(1−α)(4.6)

(4.7)

where −→

Uis the velocity field [m/s]; ρis the density [kg/m3]; t is time [s]; p is pressure

[Pa]; µis dynamic viscosity [Ns/m2]; −→

qis acceleration vector due to gravity [m/s2]; −→

xis a

spatial position vector [m] ; αis a volume fraction or indicator function [-]; Ufis the relative

velocity between the phases [m/s]; the subscripts a and w denote the fluids air and water.

The indicator function αis defined as:

α=









1 fluid w

0<α<1 transitional region

0 fluid a

For single-phase simulations the volume fraction αis 1 and constant over the whole domain

and during the simulation time.

Chapter 4: Parameter study on hydrodynamic two-phase simulations

4.2.2 Turbulence modelling

Turbulence effects and their impact on the flow have been simulated using different Reynolds

averaged (RANS) turbulence models and Large Eddy Simulations (LES). From the wide

range of RANS models the Standard k-e(Launder and Sharma, 1974), k-ω(Wilcox, 1988)

and k-ωShear Stress Transport (SST) model (Menter, 1993, 1994) were used. As subgrid

scale model for the LES simulations the Smagorinsky model (Smagorinsky, 1963) was cho-

sen.

4.2.3 Boundary conditions

For all cases presented in this paper, a similar set of boundary conditions has been used.

The inlet of the domain has been subdivided in two components, an inlet for the air phase

and an inlet for the water phase. The height of the water inlet depends on the desired water

level and the flow is prescribed by using a fixed flow velocity or a discharge. The air inlet

is specified by using a total pressure boundary condition. The upper and lower walls of

the domain are defined using a no-slip condition. The outlet of the domain is specified by

using a pressure boundary condition. The water level in the domain is fixed by defining a

weir in close proximity to the outlet as outlined in Bayón et al. (2015). For three-dimensional

testcases, the sidewalls are determined using a no-slip condition. Two-dimensional testcases

are simulated with so-called empty boundary conditions which are implemented in Open-

FOAM to describe sidewalls of a two- dimensional geometry.

4.2.4 Geometry and mesh

Unstructured meshes with local refinements at the walls were set up using the open source

mesh generation tool gmsh. The single-phase flow cases (Figure 4.1) consist of 12,106 cells

for the RANS turbulence models and 89,284 cells for the LES simulations, leading to a cell

length between 0.0042 m and 0.056 m for the RANS simulations and between 0.00035 m

and 0.041 m for the LES simulations. The two-dimensional model domains of cases 1 and 2

consist of 68,542 cells with a minimum cell length of 0.0024 m and a maximum cell length of

0.15 m. The three-dimensional setup is based on the two- dimensional domain but extended

in z-direction by ten layers. The model consequently consists of 685,420 cells. The domain

of case 3 consists of 175,762 cells. The cell length ranges between 0.0035 m and 0.144 m.

Figure 4.1: Model setup and measurement locations of single-phase case.

Chapter 4: Parameter study on hydrodynamic two-phase simulations

4.3 Single-phase flow

In order to analyze the accuracy of the interFoam solver regarding flow behavior behind a

two-dimensional ground sill, a single-phase testcase has been implemented using different

turbulence models. Experimental data for this case has been obtained by Almeida et al.

(1993) and is available in the database of the European Research Community on Flow, Tur-

bulence and Combustion (ERCOFTAC) (Davroux et al., 1995). The domain consists of a

two-dimensional duct bounded by an upper and lower wall with a polynomial shaped ob-

stacle on its bottom. The mean centerline velocity at the inlet amounts to 2.147 m/s. At

four, respectively two different locations the velocities in x- and y-direction were compared

to the experimental results: x−03 = -0.30 m (in front of the sill, inlet profile), x00 = 0.00 m

(top of the sill), x01 = 0.03 m (end of the sill), x02 =0.05 m (recirculation zone). Turbulence

models used were the RANS and LES models previously outlined.

The results show that the chosen RANS models as well as the LES simulations lead

to a good approximation of the experimental results (Figure 4.2 and Figure 4.3), however,

the LES simulations are able to capture fluctuations as well which can be of interest when

analyzing eddy structures behind ground sills. One disadvantage of LES is the higher reso-

lution of the mesh that is needed in order to display large scale eddies which lead to much

higher computation times.

Chapter 4: Parameter study on hydrodynamic two-phase simulations

Figure 4.2: Velocity profiles at x−03 =

−0.30m(inlet profile).

Figure 4.3: Velocity profiles at x01 =

0.03m(end of sill).

Chapter 4: Parameter study on hydrodynamic two-phase simulations

4.4 Two-phase flow

When analyzing two-phase flow over a two-dimensional ground sill, two main aspects are

of interest: eddy structures behind the sill and the water level drawdown. Since the accu-

racy concerning eddy structures has already been analyzed in the previous case, the focus

of the two- phase flow cases will lie on the water level drawdown. Since the water level

drawdown is constant as soon as a quasi steady-state is reached, LES simulations are not

necessary. Therefore, the Standard k-e-turbulence model has been chosen in this part in

order to save computation time (approximately 1 h instead of 4 h for parallel computation

on 16 processors).

As a first step, the effect of the sill structure on the water level drawdown has been

analyzed. A two-dimensional model setup similar to case 1 as listed in Table 4.2 has been

used. Angular shaped structures with three different angles for the ground sill as well as

a round structure have been investigated (see Table 4.1). The maximum height of the sill is

∆z = 0.20m and similar for all cases.

Table 4.1: Overview over different sill structures.

round angular, 30◦

angular, 45◦angular, 90◦

Figure 4.5 shows the resulting water level drawdowns for the different structures. The

maximum water level drawdown is the smallest (∆h=0.03m) for the round structure. The

drawdown for the 30◦and 45◦angular structure are in a comparable order of magnitude

(∆h=0.04m). A significantly higher drawdown is caused by the ground sill with an angle

of 90◦(∆h=0.09m). With an analytically calculated drawdown of ∆h=0.036m(Table 4.2),

the sill structures with small angles or round structures show the highest accuracy. The ef-

fect of the ground sill is higher for steeper hill structures which can be explained by higher

single losses for steeper structures.

In a next step, changes have been made concerning the water depth and the flow velocity

using the 30◦angular structure of the ground sill. All cases have subcritical flow conditions

and do not show a flow transition over the ground sill.

The setups including the analytically calculated water level drawdown and the simu-

lated water level drawdown are listed in Table 4.2 (case 1 to 3), where v1is the approaching

velocity and h1the water depth in front of the sill. The results show a good agreement of the

Chapter 4: Parameter study on hydrodynamic two-phase simulations

Table 4.2: Subcritical two-phase flow: Properties of different testcases.

Case 1 Case 2 Case 3 Case 4

Length of domain 25 m 25 m 35 m 25 m

Height of domain 2 m 2 m 6 m 2 m

h11.0 m 1.0 m 3.0 m 0.3 m

v11.00 m/s 1.25 m/s 3.00 m/s 3.00 m/s

∆z 0.2 m 0.2 m 0.2 m 0.1 m

∆h, analytical 0.036 m 0.070 m 0.110 m 0.100 m

∆h, numerical 0.042 m 0.090 m 0.140 m 0.055 m

simulated results with the analytical solution obtained by using Bernoulli’s and continuity

equation. A reason for the slightly higher drawdown obtained by the numerical solution

are additional energy losses which have already been shown when analyzing the influence

of the structure due to the structure of the sill.

Case 1 has also been extended to a three-dimensional geometry with a width of 1 m and

sidewalls with no-slip condition.

The resulting water level drawdown is compared to the drawdown resulting from the

two-dimensional simulation in Figure 4.6. The Figure shows that the three-dimensionality

of the testcase causes a smaller drawdown and a shorter length of the water level drawdown.

The reason for this is the influence of the sidewalls. The no-slip condition at the sidewalls

causes additional continuous energy losses that were neglected in the two-dimensional ge-

ometry.

As a next step, a strictly supercritical setup (case 4) has been analyzed. The supercrit-

ical flow case has been simulated with a height of the ground sill of ∆z= 0.10 m and an

approaching flow velocity v1=3m/s(see Table 4.2). Compared to the analytical rise of the

water level (∆h=0.10m), the maximum increase achieved in the simulations (∆h=0.055m)

is considerably smaller. The high difference between analytical and numerical solution can

be considered reasonable since the high flow velocity of this testcase leads to higher single

losses at the sill structure and therefore a higher deviation between analytical and numerical

solution.

In the last case, the stability of the simulations has been analyzed under initially dry

conditions. A water level of h1=1mand a flow velocity of v1=1m/shas been chosen

at the inlet (similar to case 1) and at the outlet a free outflow without weir. Figure 4.4

shows the behavior of the two phases in the domain at different time steps. After a sim-

ulation time of 20 seconds a quasi-steady state is reached and subcritical conditions in the

upstream part of the domain and supercritical conditions in the downstream part of the

domain are reached. A flow transition occurs over the ground sill. Close to the inlet, a

disturbance of the water surface can be found. This disturbance develops due to an eddy in

the air phase that is caused by the increase of the water level from the ground sill upstream

to the inlet. Due to the inlet boundary patch the eddy is trapped in this point. This effect

could be moved further upstream by choosing an inlet in a higher distance to the ground sill.

Chapter 4: Parameter study on hydrodynamic two-phase simulations

Figure 4.4: Filling of the domain, segment of the computational domain for different time

steps.

Figure 4.5: Comparison of water level drawdown for different sill structures

Figure 4.6: Comparison of water level drawdown for two- and three-dimensional model

setup

Chapter 4: Parameter study on hydrodynamic two-phase simulations

4.5 Conclusions

In this study, flow over a two-dimensional ground sill has been analyzed using OpenFOAM.

First, a single-phase flow has been simulated and different turbulence models have been val-

idated concerning their accuracy in describing the eddy structure behind the sill. All models

analyzed led to a good accuracy, however, the LES turbulence model was capable to account

for fluctuations as well. One disadvantage of the model is the smaller necessary grid size

which leads to much higher computation times.

In a second step, two-phase flow has been simulated using the k-eturbulence model and

different parameters such as the structure of the ground sill, discharge, water level and flow

regime have been evaluated concerning their influence on the water level drawdown.

Due to additional energy losses caused by the angular structure of the sill as well as

losses due to sidewalls when a three-dimensional model is chosen, the analytically calcu-

lated drawdown is smaller than the two- dimensional simulations but reasonably coincides

for different simulations. Changes of the water level drawdown due to the structure of the

ground sill and sidewalls when three-dimensional testcases are computed are plausible. The

simulations are also stable for a filling case with initially dry conditions.

Summing up, the VoF approach implemented in OpenFOAM is capable of describing

flow over a ground sill and similar hydraulic cases. Future research aims to closer look at

the behavior of the air phase.

Acknowledgements

This work is carried out within the DFG Research Training Group Urban Water Interfaces

which is a joint initiative of TU Berlin and Leibniz-Institute of Freshwater Ecology and In-

land Fisheries (IGB) Berlin. The authors acknowledge the funding.

Parts of the simulations were computed on the supercomputers of Norddeutscher Ver-

bund für Hoch- und Höchstleistungsrechnen in Berlin.

Chapter 5

Single-phase transport simulations

This study was published as:

Teuber, K., Broecker, T., Jaydev, S.D., Goitom, G.M., Sielaff (née Grüneberger), M.,

Despot, D., Stephan, D., Barjenbruch, M. & Hinkelmann, R.: Multiphase CFD-Simulation

of Transport Phenomena in Sewer Systems, in: Mannina G. (eds) New Trends in Urban

Drainage Modelling. UDM 2018. Green Energy and Technology, 584-853, Springer, Cham.,

2019; 10.1007/978-3-319-99867-1_146.

This is the postprint version of the article.

The test cases’ setups are listed in Appendix B (transport in rectangular pipe: B.4, transport

around concrete probes B.5, complex sewer: B.3).

Abstract

This paper presents different computational fluid dynamics applications using the multi-

phase solver interFoam which is implemented in the open source software OpenFOAM.

The solver uses the volume of fluid approach. When modelling tracer transport in the prox-

imity of the interface between two phases, the problem of non-physical tracer spreading

across the interface has to be overcome. In this paper, two ways are presented to model

such systems successfully. First, tracer transport around concrete probes in the headspace

of a sewer pilot plant is considered. In this case a two-phase (water-air) system is assumed

by describing an idealized water surface as a boundary condition and a passive tracer is

applied. Second, flow in a complex sewer stretch containing a hydraulic jump is simulated

and a tracer is applied in the water phase. A multiphase transport approach based on the

Henry coefficient is used in this case and plausible results are obtained.

5.1 Introduction

In the field of hydraulic engineering, shallow water flow models using the depth-averaged

Navier-Stokes equations are often suitable to analyse flow phenomena. However, in certain

application areas where a hydrostatic pressure distribution within the flow field is not given,

these models reach their limits and alternatives have to be looked for. In the field of urban

drainage such structures can occur in complex sewer systems containing hydraulic jumps

due to changes of the flow regime. A thorough description of transport in the different

phases and mass transfer across the phase interface can be used to describe complex species

transport and biochemical processes, for example the development of odour and corrosion

Chapter 5: Single-phase transport simulations

due to hydrogen sulphide in sewer systems. Therefore, this study focuses on transport phe-

nomena in multiphase (water-air) systems.

When a passive advection-diffusion tracer is applied in the nearfield of the interface,

a non-physical spreading across the interface between two phases can occur. Figure 5.1

illustrates the problem using the example of a rectangular pipe which is partially filled

with water. After a few seconds the tracer spreads across the water surface into the air

due to diffusion. Since the interFoam solver only distinguishes between phases by using

an additional transport equation that modifies phase specific variables, a passive tracer can

spread into the air phase. This paper presents two ways to overcome this issue by using

different approaches.

Figure 5.1: Tracer transport in water phase of rectangular pipe with non-physical spreading

across water surface.

In the first part the headspace of a sewer pilot plant is being investigated with a simpli-

fied single-phase approach, assuming an idealized water surface. The second study presents

tracer transport in a sewer stretch containing complex hydraulic flow phenomena. In this

case the tracer is present in the water phase. A special formulation that regards the mass

transfer depending on the Henry coefficient has been introduced by Haroun et al. (2010a)

and has been applied in this study.

5.2 Material and methods

5.2.1 Hydrodynamic simulations

Surface water and air flow is calculated by using the two-phase flow solver interFoam based

on a volume of fluid approach for one- and two-phase flows as it is implemented in the open

source model OpenFOAM. Both phases are considered as one fluid with rapidly changing

fluid properties, therefore one set of Navier-Stokes equations is solved. The phases are

distinguished by an additional transport equation for the volume fraction which is used as

Chapter 5: Single-phase transport simulations

a marker to describe the distribution of the phases throughout the domain. The governing

equations can be found in Rusche (2003).

5.2.2 Transport simulations

For the first study, the transport of a passive tracer with a concentration C is examined

with an advection-diffusion equation that was implemented into the interFoam solver (see

Equation 5.1). The physical diffusivity Dphys as well as the turbulent Schmidt number Scturb,

which defines the turbulent diffusivity coefficient Dturb, have to be defined by the user

(Equation 5.2):

∂C

∂t+∇ · (uC) = (Dphys +Dturb)∂2C

∂x2(5.1)

with

Dturb =µturb/ρ

Scturb

(5.2)

5.2.3 Mass transfer

For the second study, the approach introduced by Haroun et al. (2010a) as it has been

implemented by Nieves-Remacha et al. (2015) has been used. The approach is based on the

interFoam solver and adds a transport equation as outlined in Equations 5.1 and 5.2. At the

interface the following boundary conditions have to be satisfied:

He =CL

(5.3)

i∇CL

i=DG

i∇CG

i(5.4)

Where Diis the diffusivity and Ciis the concentration in the respective phase. He denotes

the dimensionless Henry coefficient and superscripts L and G denote the liquid and gas

phase. To avoid spreading of the tracer across the water surface in the case presented in

section 5.3.2, the Henry coefficient has been set to 10−6which prevents mass exchange

across the water-air interface.

5.3 Results and discussion

5.3.1 Spreading of hydrogen sulphide in headspace of sewer pilot plant (air

phase)

In the first part the headspace of a sewer pilot plant is being investigated. In the pilot plant,

probes made of various types of concrete are installed to investigate their resistance against

sulphuric acid corrosion. The diameter of the pipe is 400mm. The original setup of the pipe

has a length of 25m, for the computational domain a length of 2m has been chosen.

Simulations are performed to detect whether the concentration around the concrete

probes is distributed homogeneously or whether peaks occur that might influence the cor-

rosion rate locally. The simulations are being performed with a simplified single-phase

approach, assuming an idealized water surface. Since the water surface is being described

using a slip boundary condition, only the upper half of the pipe has to be discretized with

Chapter 5: Single-phase transport simulations

the computational mesh. The computational domain is shown in Figure 5.2 and the result-

ing tracer concentration is illustrated in Figure 5.3. The results show an accumulation of the

tracer around the probes in the rear part of the pipe.

Figure 5.2: Top view of the model domain (pipeline containing the concrete probes).

Chapter 5: Single-phase transport simulations

Figure 5.3: Resulting tracer distribution.

Chapter 5: Single-phase transport simulations

5.3.2 Tracer transport in sewer pipe (water phase)

In the second case, a complex sewer stretch containing a hydraulic jump has been analysed.

The overall geometry has a length of roughly 90m and the sewer pipe consists of different

cross sections with widths up to 7.50m and a maximum height of 5.63m The detailed ge-

ometry has been outlined in Bayón et al. (2015). In this case the mass transfer formulation

by Haroun et al. (2010a) has been applied to keep the tracer in the water phase. Figure 5.4

shows that the tracer does not cross the water surface in the case of the simple rectangular

geometry introduced in Section 5.3.1. Figure 5.5 illustrates the movement of the tracer along

the complex sewer stretch and through the hydraulic jump. The results confirm that also in

cases with high turbulence the tracer remains in the water phase.

Figure 5.4: Movement and spreading of tracer along sewer.

Chapter 5: Single-phase transport simulations

Figure 5.5: Movement and spreading of tracer along sewer.

Chapter 5: Single-phase transport simulations

5.4 Conclusions

Within this study three-dimensional two-phase (water-air) flow and transport simulations

have been carried out in two different examples. These systems are difficult to model in

a way that a passive advection-diffusion tracer can cross the barrier of the water surface

without any physical constraints due to the formulation of the interFoam solver. In this

paper, two ways are presented to overcome this obstacle. First, to model the air phase as a

single phase system which can be valid for stratified flows in sewers and second, to use a

mass transfer approach based on the Henry coefficient and modify the Henry coefficient in

order to avoid mass transfer across the water surface.

Chapter 6

Mass transfer simulations under

equilibrium conditions: Validation

and solver extension

This study was published in Water Science and Technology (IWA Publishing) as:

Teuber, K., Broecker, T., Bentzen, T.R., Stephan, D., Nützmann, G. & Hinkelmann, R.:

Using computational fluid dynamics to describe H2S mass transfer across the water-air inter-

face in sewers, Water Science and Technology, 79 (10), 1934-1946, 2019; doi:10.2166/wst.2019.193.

©IWA Publishing 2019. The definitive peer-reviewed and edited version of this article is

published in doi:10.2166/wst.2019.193 and available at www.iwapublishing.com.

Reproduced from Teuber et al. 2019 Water Science & Technology 79(10) 1934-1946, with

permission from the copyright holders, IWA Publishing.

A detailled explanation of the solver extensions addressed in this paper is given in Ap-

pendix A. The test cases’ setups are listed in Appendix B (quasi one-dimensional tank: B.6,

rectangular pipe B.7, complex sewer: B.3).

6.1 Abstract

For the past 70 years, researchers have dealt with the investigation of odour in sewer systems

caused by hydrogen sulphide formations and the development of approaches to describe it.

The state-of-the-art models are one-dimensional. At the same time, flow and transport

phenomena in sewers can be three-dimensional, for example the air flow velocities in cir-

cular pipes or flow velocities of water and air in the reach of drop structures. Within the

past years, increasing computational capabilities enabled the development of more complex

models. This paper uses a three-dimensional two-phase Computational Fluid Dynamics

model to describe mass transfer phenomena between the two phases: water and air. The

solver has been extended to be capable to account for temperature dependency, the influence

of pH value and a conversion to describe simulated air phase concentrations as partial pres-

sure. Its capabilities are being explored in different application examples and its advantages

compared to existing models are demonstrated in a highly complex three-dimensional test

case. The resulting interH2SFoam solver is a significant step in the direction of describing

and analysing H2S emissions in sewers.

Chapter 6: Mass transfer: Solver extension and validation

6.2 Introduction

Wastewater in sewers undergoes a lot of physical and biochemical processes. One important

factor is the formation of hydrogen sulphide (H2S), which can cause health risks for sewer

workers. The tendency of more complex and longer sewer networks can lead to longer re-

tention times, which enhance the emission of H2S. Climate change at the same time causes

higher temperatures in the wastewater, which increases emission rates.

In the past 70 years extensive research has been performed to increase the knowledge

on H2S formations and to develop approaches, which describe the development of odour

in sewers (e.g. Gilchrist (1953), Thistlethwayte (1972)). The state-of-the-art models, which

have been developed within the last 20 years, are horizontal one-dimensional. These are the

SeweX model from Australia (Rootsey and Yuan, 2010, Rootsey et al., 2012) and the WATS

model from Denmark (Hvitved-Jacobsen et al., 2013). Both are not public domain.

An overview of existing model approaches has been given in Carrera et al. (2016) and

the need for further research has been highlighted.

To begin with, the mass transfer approach of the existing models is based on the so-called

two-film theory, which uses different assumptions. The WATS model additionally uses dif-

ferent approaches to account for turbulent H2S transfer rates across the water surface in

various applications. These different approaches are empirical or theoretical connections

between oxygen and H2S transfer on the one hand and empirical models linking H2S emis-

sions to flow properties in the pipe on the other hand (Carrera et al., 2016).

Wang et al. (2018) highlight the shortcomings of the two-film theory. It cannot account

for local changes of the flow regime or variations of fluid properties. Furthermore, the the-

ory is based on a constant liquid film that can change in real-life conditions due to flow

instabilities. The most limiting factor however is assumed to be the one-dimensionality of

the approach. More advanced approaches, the penetration theory and the surface renewal

theory, can account for the variability of the flux over time but do not account for local varia-

tions, the change of fluid properties or flow regimes (Wang et al., 2018). This has already led

to wide applications of CFD models for mass transfer applications in the chemical industry

(Wang et al., 2018).

Carrera et al. (2016) identified the models’ lack to describe mass transfer processes across

the water surface, the current approaches of which were considered to be simplified, espe-

cially when considering hydraulic structures such as gravity sewers, junctions and water

falls. Recent research on water falls or drop structures in sewer systems led to improved

formulations to account for the effect of local turbulence (Matias et al., 2017), but these ap-

proaches are still empirical equations which are fed into the model.

This short overview leads to the question whether a three-dimensional simulation model

could help in increasing the process understanding, especially when analysing complex and

turbulent flows in a sewer. Another benefit could be the in-depth analysis and design opti-

mization in hotspots of H2S emissions.

In order to address this question, a volume of fluid (VOF) approach as it is implemented

in OpenFOAM’s solver interFoam has been chosen to describe the two-phase flow of water

and air. This solver has already been used for a number of demanding hydraulic applica-

Chapter 6: Mass transfer: Solver extension and validation

tions (e.g. Thorenz and Strybny (2012), Bayón et al. (2015)) and enables a stable, robust and

accurate description of complex flow phenomena.

The VOF method is often used to describe mass transfer processes in CFD applications

(Wang et al., 2018). Therefore, Haroun et al. (2010a,b) have developed an approach to de-

scribe mass transfer processes across the interface between two fluids using the Henry co-

efficient for the VOF method. This approach has been implemented in OpenFOAM’s solver

interFoam by Nieves-Remacha et al. (2015), Yang et al. (2017) and Severin (2017), resulting

in a solver which will be called interHarounFoam in the following.

A short outline of the driving processes leading to H2S formation shall be given to

highlight important factors. When anaerobic conditions occur in the wastewater, sulphate-

reducing bacteria, which reside in the biofilms of sewer walls can reduce sulphate to sul-

phide (Sharma et al., 2008a). From the biofilm, sulphide is then diffused into the wastewater

as H2S. In the water, equilibrium conditions depending on the pH value and temperature

determine which amounts of sulphide are present as H2S and as bisulphide ion (HS-). The

air-water equilibrium, which can be described by the Henry coefficient for a volatile com-

pound such as H2S, can cause emissions of H2S from the water into the air phase. The

rate of the transfer process is influenced by factors such as the flow velocities within the

different phases, the pH value, temperature and the concentration of oxygen and nitrate.

The Henry coefficient describes the relative amount of a volatile compound in the gas phase

as a function of its relative occurrence in the water phase under equilibrium conditions and

at constant temperature. The temperature dependency of Henry’s law can be described by

different equations for example, the van’t Hoff equation. The concentration of H2S in the

air phase defines the intensity of odour (Hvitved-Jacobsen et al., 2013).

As this overview of the relevant processes shows, a sole consideration of the Henry coef-

ficient when describing H2S emissions is not sufficient. Therefore, relevant extensions have

been made to the solver, resulting in a new, specialized solver, interH2SFoam. This solver is

able to account for the temperature dependency of the Henry coefficient. Further extensions

enable the user to describe the equilibrium between HS-and pH value in the water phase

and to compute the partial pressure of H2Sgin the air phase in ppm in order to gain a better

comparability between simulations and measured values. The assessment of turbulent flow

effects on mass transfer will be subject to future research.

In the following, after an introduction on the methods used, the capabilities of the

interH2SFoam solver are explored in three simple application examples of vertical one-

dimensional flow. Then, mass transfer in a rectangular pipe is simulated. In a final example,

the new solver is applied to a highly complex sewer geometry.

6.3 Methods

6.3.1 Numerical model

OpenFOAM version 2.4.0 has been used for the work presented in this paper. Additionally,

a supplementary library called swak4Foam has been used to generate customized function

objects to calculate the equilibrium conditions between H2S and HS-as well as the partial

pressure in the air phase. This approach makes the use of this function optional for the

user. Depending on the framework of the model, the user can then decide whether these

Chapter 6: Mass transfer: Solver extension and validation

functions are needed or not. The temperature dependency on the Henry coefficient of H2S

has been directly implemented in the solver and makes a definition of the temperature as

an input parameter mandatory.

Hydrodynamic simulations

The mass transfer solvers are based on the two-phase flow solver interFoam which is based

on a VOF approach that considers both phases as one fluid with changing fluid properties.

One set of Navier-Stokes equations is solved. The volume fraction of a phase is stored as an

additional variable and the phases are distinguished by an additional transport equation.

The equations are defined as follows Rusche (2003):

Mass conservation equation:

∇ · −→

U=0 (6.1)

Momentum conservation equation:

∂ρ−→

∂t+∇ · ρ−→

U−→

U=−∇prgh +∇ · µ∆−→

U+∇−→

U∇µ−−→

g·−→

x∇ρ(6.2)

Where prgh is the static pressure minus hydrostatic pressure:

prgh =p−ρ·g·h(6.3)

Volume of Fluid equation:

∂α

∂t+∇ · α−→

U+∇ · ((1−α)Urα)=0 (6.4)

with the following parameters:

ρ=αρaq +ρg(1−α)(6.5)

µ=αµaq +µg(1−α)(6.6)

µi=µi,phys +µi,turb with i=aq,g (6.7)

where −→

Uis the velocity field [m/s]; ρis the density [kg/m3]; t is time [s]; p is the pressure

[Pa]; µis the dynamic viscosity [Ns/m2]; −→

qis the acceleration vector due to gravity [m/s2];

−→

xis a spatial position vector [m]; αis a volume fraction or indicator function [-]; Ufis the

relative velocity between the phases [m/s]; the subscripts aq and g denote the fluids water

(aq - aqueous) and air (g - gas). For the dynamic viscosity µ, the physical viscosity and the

turbulent viscosity are considered (see Equation 6.7).

The indicator function αis defined as:

α=









1 fluid aq

0<α<1 transitional region

0 fluid g

(6.8)

The water surface is defined as the area where α=0.5.

A turbulence model based on the Reynolds averaged Navier-Stokes equations (Standard

k-e) is applied to consider the turbulent part and the near-wall turbulence is modelled by

so-called wall functions. More advanced turbulence models, such as Large Eddy Simu-

lations (LES) or Direct Numerical Simulations (DNS), would offer the opportunity of re-

solving small-scale velocity variations but would come with the price of a highly increased

Chapter 6: Mass transfer: Solver extension and validation

computation time. As we expect that the application of more advanced turbulence effects

would not change the insights on the equilibrium conditions of the mass transfer simula-

tions addressed in this publication, a RANS turbulence model has been considered to be

sufficient as well as the best way to save computational resources. Even with the Standard

k-eturbulence model, the computation time of 10 seconds simulation for the complex sewer

geometry amounted to 12 hours on 80 parallel processors using the high performance com-

puting (HPC) clusters of TU Berlin.

The accuracy of the hydrodynamic simulations has been assessed in Teuber et al. (2019a).

Transport simulations

In general, the transport of a passive tracer with a concentration c is examined with an

advection-diffusion equation that can be implemented into the interFoam solver (see Equa-

tion 6.9). The physical diffusivity Dphys as well as the turbulent Schmidt number Scturb,

which defines the turbulent diffusivity coefficient Dturb, then have to be defined by the user

(Equation 6.10).

Advection-diffusion equation:

∂c

∂t+∇ · (−→

Uc) = (Dphys +Dturb)∇c(6.9)

with

Dturb =µturb/ρ

Scturb

(6.10)

Mass transfer

Mass transfer has been simulated using the approach defined by Haroun et al. (2010a,b) as

it has been implemented by Nieves-Remacha et al. (2015) and Severin (2017). The approach

is based on the interFoam solver and considers one additional transport equation for both

phases as outlined in Equations 6.9 and 6.10.

∂c

∂t+∇ · (−→

Uc) = ((Dphys +Dturb)∇c+φ)(6.11)

A concentration flux expression at the interface results in the following:

φ=−(Dphys +Dturb)c(1−He)

α+He(1−α)∇α(6.12)

In order to distinguish the species transport between the two phases, Henry’s law must be

fulfilled and the concentration flux must be consistent:

He =caq

(6.13)

(Dphys,aq +Dturb,aq)∇caq = (Dphys,g+Dturb,g)∇cg(6.14)

The concentrations and diffusion coefficients are considered as single-phase properties de-

pending on the phase fraction value α:

Chapter 6: Mass transfer: Solver extension and validation

c=αcaq +cg(1−α)(6.15)

Dphys =Dphys,aq ·Dphys,g

αDphys,aq + (1−α)Dphys,g(6.16)

The diffusion coefficients for Dphys,aq and Dphys,gare defined by the user. Note that these

coefficients are temperature dependent which has to be taken into account when defining

the values.

6.3.2 Henry coefficient

The Henry coefficient, also known as Henry constant, is a temperature dependent variable

which is reported in many different forms in literature and is often presented in different

units.

In this paper, three different definitions of the Henry coefficient are relevant for deriva-

tion and comparison with analytical solutions. Sander (2015) lists values of Henry coeffi-

cients in the unit [mol/(m3Pa)] and defines this Henry coefficient as Hcp. For the implemen-

tation in the interHarounFoam solver, the dimensionless Henry coefficient Hcc is relevant

(see Equation 6.13):

Hcc =He =caq

(6.17)

Where the Henry coefficient is expressed as the ratio between the concentration in the aque-

ous phase caq and the concentration in the gas phase cg.

Hcp can be converted to Hcc using the ideal gas law:

Hcc =Hcp ·R·T(6.18)

Where R is the universal gas constant 8.314 kgm2

s2molK and T is the temperature [K].

For H2S, the Henry coefficient at standard temperature (25◦C) results in:

Hcc

H2S=10−3mols2

m2kg ·8.314 kgm2

s2molK ·298.15K=2.479 (6.19)

6.3.3 Extensions

Temperature dependency of Henry coefficient

The Henry coefficient depends on the overall temperature in the domain. Therefore, the

temperature dependency has been added in a way that the solver takes one global temper-

ature value as an input parameter.

The temperature dependent Henry coefficient is computed using the van’t Hoff equation

following Sander (2015):

Hcp(T) = Hcpexp C1

T−1

TΦ (6.20)

Chapter 6: Mass transfer: Solver extension and validation

Where Cis a temperature coefficient which is depending on the enthalpy of dissolution and

defined as 2100 K (Sander, 2015), TΦis the standard temperature 298.15 K corresponding to

25◦C.

Equilibrium conditions

The equilibrium conditions are implemented following Hvitved-Jacobsen et al. (2013). The

aim is to describe the water-phase concentration of H2S depending on the amount of total

dissolved sulphide and the pH value since those values are usually measured in field inves-

tigations.

The dissociation of H2S is generally expressed by the following equilibrium:

H2Sg↔H2Saq ↔HS−+H+↔S2−+2H+(6.21)

The equilibrium between H2S in the gas phase (H2Sg) and H2S in the water phase (H2Saq)

is described by the Henry coefficient. In the water phase, an equilibrium between hydrogen

sulphide H2Saq and bisulphide ion (HS-) exists, where the total amount of both is described

as total dissolved sulphide.

Only H2Saq can be transferred across the air-water interface, not the ionized form HS-,

however, usually the concentration of total dissolved sulphide and the pH value are mea-

sured. Therefore it is useful to derive a way to calculate the concentration of cH2Saq when cS

and pH are given.

The equilibrium depends on the equilibrium constant Ka1(also known as acid dissociation

constant):

Ka1=cH+cHS−

cH2Saq

(6.22)

The dissociation can also be described using the negative logarithm of Ka1,pKa1(pKa1=

−logKa1), resulting in the Henderson-Hasselbalch equation:

log10

cH2Saq

cHS−=pKa1−pH (6.23)

For a temperature of 20◦C the equilibrium constant is pKa1=7.0. Between the ionized form

HS-and the sulphide ion another equilibrium exists in the water phase:

HS−↔S2−+H+(6.24)

Ka2=cH+cS−2

cHS−(6.25)

The value of pKa2=14.0 indicates that measurable amounts of the sulphide ion S2- only

exist at a value above a pH of about 12. Therefore, only the equilibrium value of Ka1is

important for wastewater and has been implemented for the interHarounFoam solver in

OpenFOAM using the utilities funkySetFields (for initial conditions) and funkySetBound-

aryFields (as boundary conditions) from swak4Foam.

Chapter 6: Mass transfer: Solver extension and validation

The step-by-step reformulation based on Hvitved-Jacobsen et al. (2013) of the equations

which results in the equation implemented in OpenFOAM is shown in the following, be-

ginning with the Henderson-Hasselbalch equation:

log10

cH2Saq

cHS−=pKa1−pH (6.26)

Solving the log-function and using the expression cS=cHS−+cH2Saq for the total dissolved

sulphide:

10pKa1−pH =cH2Saq

cHS−=cH2Saq

cS−cH2Saq

(6.27)

Rearranging leads to the mass concentration γH2Saq in [kg/m3]:

γH2Saq =cS·10pKa1−pH

1+10pKa1−pH (6.28)

This can be converted into a molar concentration [mol/m3] by dividing through the atomic

weight MS(0.032 kg/mol) of sulphur (S):

cH2Saq =γH2Saq

=γH2Saq

0.032kg/mol (6.29)

Thus, the resulting equation implemented in OpenFOAM is:

cH2Saq =

cS·10pKa−pH

1+10pKa−pH

32 (6.30)

Note that the equilibrium constants Ka1and Ka2are temperature dependent (Yongsiri et al.,

2004), which is not considered in the current version of the code. The value of Ka1is a

user-defined variable and the temperature dependency of the equilibrium constant has to

be accounted for by defining the corresponding Ka1value for the temperature analysed.

Calculation of partial pressure of H2Sgin ppm

The partial pressure of H2Sgis being computed using a function object in swak4Foam.

The input value in OpenFOAM is a tracer cH2Saq in [mol/m3], requiring a unit conversion:

c[mol

l] = 1mol

l=cH2Saq

1000 =1000mol/m3

1000 (6.31)

The conversion from ppm to atm is:

10−6atm =1ppm (6.32)

According to Hvitved-Jacobsen et al. (2013), the partial pressure of a trace quantity in the

air phase can be expressed by multiplying the molar concentration with the molar volume:

pH2Sg[atm] = cH2Saq [mol/l]·22.4l/mol (6.33)

Together, this leads to the following equation for conversion:

pH2Sg[ppm] = 106·cH2Saq [mol/l]

1000 ·22.4l/mol (6.34)

This conversion is only valid for the gas phase concentration, therefore the expression is

multiplied with (1-α) in order to keep the values constrained to the air phase.

Chapter 6: Mass transfer: Solver extension and validation

6.3.4 Case studies

In the following, three different cases will be used to explore the possibilities of the existing

solver, to validate the new features added to the solver and to show the importance of the

model compared to the existing model approaches. For all simulations, at a temperature of

25◦C, physical diffusivities for H2S in water of 2.2 ·10−9m2/sand in air of 1.74 ·10−5m2/s

are chosen.

The first setup is a quasi-one-dimensional cubic tank with the measures 1m x 1m x 1m

bounded by upper, lower and sidewalls with no-slip conditions. The tank is partially filled

with water (water depth d = 0.5 m). Both fluids water and air are at rest. As an initial condi-

tion, an H2S concentration of cH2Saq =1mol/m3in the water phase is given, the concentration

in the air phase is cH2Sg=0mol/m3. The domain is discretized with 100 cells in y-direction,

which is the vertical dimension of the domain, and 10 cells in x- and z-direction. At the

bottom wall, a concentration source is assumed, using a fixed value boundary condition of

1 mol/m3. The top wall as well as the sidewalls are defined with zeroGradient conditions.

This setup is used to illustrate the solver’s capabilities in a simple setup. In a first example,

mass transfer, as it can be described with the existing interHarounFoam solver, is shown in

a vertical one-dimensional case. Then, the extensions leading to the interH2SFoam solver

are demonstrated in different examples using this first setup.

In a second setup, mass transfer in a rectangular duct is analysed using two well-

documented examples of water-air pipe flow as they have been described by Bentzen et al.

(2016) (test cases no. 7 and 21). The investigated pipe has a length of 15.0m, a height of

0.26m and a width of 0.3m with two different water depths and slopes. The air phase is only

accelerated by the movement of the water surface. Bentzen et al. (2016) measured resulting

velocity profiles in detail using Laser Doppler Anemometry (LDA) velocity measurements.

The flow characteristics of the two test cases analysed are listed in Table 6.1. The setup

is a relatively simple three-dimensional setup of a pipe. It illustrates the applicability of

the model to regular pipes. The computational domain consists of 307,970 cells. The inlet

has been divided in two parts: one for the water phase and one for the air phase. For the

water phase, a fixed discharge has been defined, and the phase fraction value αhas been

defined to be α= 1. The pressure boundary condition has been defined as null Neumann

condition. For the air phase, a fixed pressure has been defined and the phase fraction value

has been set to α= 0. The velocity has been defined using a null Neumann condition. At

the outlet, a free outflow has been assumed. A fixed pressure has been defined and the

remaining boundary conditions were defined as null Neumann conditions. At the walls,

no-slip conditions were applied. Hydrodynamic simulations (without mass transfer) were

run for 200s, until quasi steady state conditions were reached, afterwards a concentration

cH2Saq =1mol/m3has been defined for the water at the inlet, assuming contaminated water

flowing into the domain. The upper fluid then has a concentration of cH2Sg=0mol/m3at

the inlet, all remaining boundaries were defined with null Neumann conditions.

The third setup describes a complex sewer geometry with an overall length of 93.3m, a

width ranging from 6.0m to 7.5m and a sewer height between 4.3m and 5.3m. The setup is

shown in Figure 6.6. This geometry has been simulated in OpenFOAM and compared to

experimental results from a 1:20 scale model by Bayón et al. (2015) and Teuber et al. (2019a).

The setup describes a highly three-dimensional pipe diversion, including bends and geome-

try changes as well as a hydraulic jump. The computational mesh consists of 3,029,223 cells.

The setup of boundary conditions is similar to the rectangular pipe of Bentzen et al. (2016).

Chapter 6: Mass transfer: Solver extension and validation

Table 6.1: Mass transfer in rectangular channel: flow properties of analysed test cases.

Test no. Duct slope

[%]

Water

depth[cm]

Uaq [m/s] Ug[m/s] Reynolds

number

Reaq

Reynolds

number

Reg

7 0.57 3.15 0.77 0.226 72,300 5,400

21 1.34 4.00 1.37 0.336 175,300 7,900

The hydrodynamic model has been simulated for 200s, until steady-state conditions were

reached. Then, a concentration cH2Saq =1mol/m3has been defined for the water phase and

the simulations using the interH2SFoam solver have been carried out for a simulation time

of 10s. A temperature of 25◦C is assumed.

The results of the numerical simulations are presented in the following Section.

6.4 Results and discussion

6.4.1 Saturation of H2S in a tank

Mass transfer modelling

In our first test case, we present the application of the model to a vertical one-dimensional

problem. It illustrates the advantage of the new model in describing vertical concentration

profiles in contrast to the existing horizontal one-dimensional approaches. The simplicity of

the test case enables a first illustration of the model’s capabilities. The simulation has been

carried out assuming normal temperature (25◦C).

Figure 6.1 shows the presence of the two phases within the domain (α= 1: water, α=

0: air) and the development of the concentration profile within the domain over time. After

t = 50s, a decrease of the overall concentration in the water phase can be observed. This is

due to the concentration jump at the interface, which has to be fulfilled by the solver. This

concentration jump occurs in the first second due to a direct flux of concentration across

the interface. After several seconds, the concentration in the water phase is re-established

by the source term at the bottom and after t = 1000s, a steady-state has developed and a

constant concentration profile is achieved. The concentration in the water phase is then

equal to the source term concentration and the air phase concentration is defined by the

Henry coefficient. A detailed validation of the flux under transient conditions has been

performed by Haroun et al. (2010a).

The concentration profile illustrates that the resulting air phase concentration is cH2Sg=

0.4034mol/m3, which is the expected concentration in the air phase when applying Henry’s

law for H2S:

cH2Sg=cH2Saq

Hcc

H2S

=1mol

2.479 =0.4034mol

m3(6.35)

Temperature dependency

In this test case, we will analytically analyse the temperature dependency of the Henry co-

efficient, which has been implemented. The application example is based on example 4.2

in Hvitved-Jacobsen et al. (2013). The Henry coefficient at a temperature of 15◦C is being

Chapter 6: Mass transfer: Solver extension and validation

Figure 6.1: H2S saturation in a tank (left: phase fraction value, right: concentration profiles

along the vertical axis over time).

calculated.

In this case, a temperature of 288.15 K has been chosen. The resulting Henry coefficient

can be determined as follows (following Equations 6.18 - 6.20):

Hcc(T) = Hcpexp C1

T−1

TΦRT (6.36)

Hcc(288.15) = 0.001 ·exp 2, 200 1

288.15 −1

298.15·8.314 ·288.15 (6.37)

H288.15 =3.083 (6.38)

Resulting in the following expected gas-phase concentration cH2Sg=0.324mol/m3:

cH2Sg(288.15) = 1

3.083 =0.324mol/m3(6.39)

Figure 6.2 shows the resulting concentration in the domain after t = 1000s. The result

agrees well with the expected concentration. The implemented temperature dependency

can therefore be considered as accurate.

Equilibrium conditions and unit conversion

In order to validate the solver extensions regarding the equilibrium conditions and the par-

tial pressure in the air phase, example 4.3 by Hvitved-Jacobsen et al. (2013) is simulated.

The resulting H2Saq and H2Sgconcentrations for a measured concentration of dissolved sul-

phide cS-and pH value have been simulated. Again, the basic setup of the case is the same

as for the first application example. A temperature of 15◦C, pKa1 = 7.0, pH = 7.0 and a

Chapter 6: Mass transfer: Solver extension and validation

Figure 6.2: H2S saturation in a tank for different temperatures (298.15K (cp. Figure 6.1, right)

and 288.15K) (left: phase fraction value, right: concentration profiles along the vertical axis).

dissolved sulphide concentration of 0.001 kg/m3are given.

In Hvitved-Jacobsen et al. (2013), the Henry coefficient for the given temperature of 15◦C

is assumed to be the same as the Henry coefficient from a previous calculation for a temper-

ature of 20◦C, i.e. 433 atm. For comparing the analytical solution with the simulated values,

the exact Henry coefficient for 15◦C has been calculated. Using this value and performing

the same calculation steps with the corrected Henry coefficient, the analytical solution leads

to a water phase H2S concentration of cH2Saq = 0.0075 mol/m3, a gas phase concentration of

cH2Sg= 0.0027mol/m3and a corresponding partial pressure of pH2Sg= 66 ppm.

Figure 6.3 shows the results of the numerical simulations. In the water phase, the con-

centration of H2S is cH2Saq = 0.0075mol/m3, in the air phase the concentration reaches a value

of cH2Sg= 0.0027mol/m3and a corresponding partial pressure (in the gas phase) of pH2Sg

= 66 ppm and thus agrees well with the analytical solution. The implemented approach

predicts the resulting concentrations accurately.

Chapter 6: Mass transfer: Solver extension and validation

Figure 6.3: Application example for HS-and H2S equilibrium and partial pressure of air

phase concentration (left: Phase fraction value profile over domain height, middle: concen-

tration profile, right: partial pressure of gas-phase concentration).

Chapter 6: Mass transfer: Solver extension and validation

6.4.2 Mass transfer in a rectangular channel

In this test case, we present mass transfer simulations in a rectangular pipe. Figures 6.4 and

6.5 present the resulting phase fraction, velocity and concentration profiles along the height

of the domain in the middle of the pipe. The simulated velocity profiles indicate a good

agreement with the measured values by Bentzen et al. (2016). For a pipe with a length of

15 m and the analysed flow velocities, the concentration profiles show that almost no mass

transfer across the water surface into the air phase can be observed. This can be explained by

small velocities in directions other than the main flow direction (i.e. in the yz-plane) which

cause advective transport to occur mostly in the main flow direction (x-direction). Further-

more, the small diffusion coefficients cause mainly advective transport. This example opens

the question how simulated mass transfer is influenced under highly turbulent conditions

or in cases with higher velocities in the yz-plane, which will be analysed in the next example.

Figure 6.4: Mass transfer in rectangular channel for test 7 (see Table 6.1) (left: phase fraction

value, middle: velocity, right: concentration profiles).

Chapter 6: Mass transfer: Solver extension and validation

Figure 6.5: Mass transfer in rectangular channel for test 21(see Table 6.1) (left: phase fraction

value, middle: velocity, right: concentration profiles).

Chapter 6: Mass transfer: Solver extension and validation

6.4.3 Complex sewer geometry

In a final example, the advantage of the new model are demonstrated by applying the

solver to a complex and highly three-dimensional sewer geometry. The existing models are

not public domain, therefore a direct comparison to the one-dimensional models is not pos-

sible, but the results of the CFD model will be used to highlight the advantages compared

to the concepts of the existing approaches.

The results of the simulations at t = 10s are displayed in Figure 6.6. Figure 6.6 a) gives

an overview of the computational domain and the water phase behaviour. The location of

highest turbulence occurs in the hydraulic jump, which is displayed in Figures 6.6 b) and 6.6

c). The velocity vectors in Figure 6.6 b) indicate the highly three-dimensional flow behaviour

in this location and show the complex water surface movement. In Figures 6.6 c) and 6.6

d), the isosurfaces of the resulting H2S concentration in the domain are displayed. The con-

centration range between 0 mol/m3and 1mol/m3 has been divided into 10 surfaces. The

value range in between is not displayed, leaving white spaces for better illustration of the

surfaces. The contour plots show, that a more diverse and highly three-dimensional concen-

tration profile develops at the reach of the hydraulic jump. This indicates the increased mass

transfer (i.e. higher distance of concentration isolines to the water surface) in the location of

the hydraulic jump.

Because the existing model approaches are not public domain, a direct comparison to

simulation results is not possible, however, the advantages of the new CFD based mass

transfer approach are the following:

1. In respect to the hydrodynamic behaviour, the new model can describe the three-

dimensional flow velocities in the air and water phase. The sewer geometry analysed

consists of a bent pipe structure with varying shapes and a hydraulic jump. A hydro-

dynamic one-dimensional approach would describe this geometry as one connection

pipe between beginning and end point. The flow velocity would be calculated as a

uniform value without accounting for the complexity of the geometry. The existing

model would not account for the highly complex interaction of water and air phase in

the hydraulic jump.

2. Regarding the mass transfer, the model would then account for advection and molec-

ular diffusion and for turbulence in the free-stream flow areas as well as in drop

structures in a very simplified way. This would lead to a simplified assumption of the

actual mass transfer occurring in the pipe, since the effect of turbulence on the mass

transfer is substantial. A validation of the actual mass transfer rate due to turbulence

effects is performed in Teuber et al. (in preparation).

Most sewer stretches in urban areas are not as complex as the previously shown example

and wide networks without high levels of turbulence justify the use of one-dimensional

models. However, locations of high turbulence can enhance H2S emissions and the three-

dimensional approach presented in this paper can help analyse the effect of local design

aspects on the resulting H2S emissions and improve the sewer network design.

Chapter 6: Mass transfer: Solver extension and validation

(a)

(b)

(c)

(d)

Figure 6.6: Mass transfer simulations in complex sewer geometry at t=10s (a) overview of

the domain filled with water under steady-state conditions, b) flow velocities and water

surface behaviour in hydraulic jump, c) tracer distribution in hydraulic jump, d) top view

on tracer distribution).

Chapter 6: Mass transfer: Solver extension and validation

6.5 Conclusions

H2S emissions and their consequences are an important topic when considering urban

drainage and the design of sewer networks. In the past, different model approaches, from

empirical to conceptual, have been developed in order to describe and predict H2S emissions

and resulting odour. These models are horizontal one-dimensional approaches, therefore

neglecting the occurrence of three-dimensional effects.

In this publication, a model approach has been introduced that can describe H2S emis-

sions across the water surface using a mass transfer approach based on the Henry coefficient,

which is implemented in the open source software OpenFOAM. Two-phase flow has been

simulated using a VOF method. The solver has been extended by different key features

that are crucial when describing H2S emissions. The temperature dependency of the Henry

coefficient has been taken into account. Equilibrium conditions between HS-and H2S are

described and enable the usage of the measured value for total dissolved sulphide and the

pH value as input parameters. The solver also computes the partial pressure of H2S in the

gas phase based on the simulated concentration of H2Sg.

The new solver has been applied to different simple test cases and the results have been

compared to analytical solutions. Furthermore, it has been applied to a highly complex

three-dimensional test case to highlight the advantages of the new model approach. Com-

pared to one-dimensional formulations, it can account for highly complex flow effects in a

sewer stretch and describe mass transfer in such environments. The analysis of the results

showed an increased mass transfer in the location of highest turbulence, which agrees with

existing observations. The exact quantification of local mass transfer rates has been vali-

dated in Teuber et al. (in preparation) and has led to a good agreement with experimental

results.

Overall, the new solver enables an analysis of mass transfer in complex three-dimensional

test cases, the description of which has so far only been possible with major simplifications.

Future research will deal with further extensions of the solver to account for temperature

effects in the fluids and reactive transport modelling.

6.6 Acknowledgements

The complex sewer geometry has been computed on the supercomputers of Norddeutscher

Verbund für Hoch- und Höchstleistungsrechnen in Berlin.

The funding provided by the German Research Foundation (DFG) within the Research

Training Group ‘Urban Water Interfaces’ (GRK 2032-1) is gratefully acknowledged.

We thank Prof. Arnau Bayón for providing the experimental data and the mesh for the

complex sewer geometry.

Chapter 7

Mass transfer under turbulent

conditions

This study will be submitted to Water Science and Technology (IWA Publishing) as:

Teuber, K., Broecker, T., Nützmann, G. & Hinkelmann, R.: CFD simulation of H2S mass

transfer under turbulent conditions in a stirring tank, in preparation for submission to Water

Science and Technology (IWA Publishing).

This is the preprint version of the article.

A detailled explanation of the solver extensions addressed in this paper is given in Appendix

A. The setup of the stirring tank test case is listed in Appendix B.8.

7.1 Abstract

The present paper deals with a study of H2S emissions under turbulent conditions in a stir-

ring tank. A new three-dimensional two-phase (water / air) solver which is able to describe

H2S emissions across the water surface is connected to a dynamic meshing functionality and

extended to compute the occurring mass transfer coefficient. H2S emissions are important

to be investigated because of the possible health risk for sewer workers and high costs for

sewer maintenance. The influence of turbulence on the mass transfer coefficient has been in-

vestigated in recent years but up to now there has not been any three-dimensional model to

simulate the direct exchange of H2S between water and air phase. Being able to describe the

mass transfer in a three- dimensional model can help to improve existing model concepts as

well as to better examine hotspots of H2S emissions in order to give recommendations for

design improvement.

7.2 Introduction

Cities in the modern world face problems caused by more and more centralized sewer sys-

tems. One of these problems is the occurrence of odour and corrosion due to hydrogen

sulphide emissions across the wastewater-air interface, which cause a health risk for sewer

workers as well as high costs for sewer maintenance.

In the past years, significant progress has been made in the field of understanding hy-

drogen sulphide (H2S) emissions in sewers. On the one hand, conceptual model approaches

have been developed to describe the occurrence of odour and corrosion (Hvitved-Jacobsen

Chapter 7: Mass transfer under turbulent conditions

et al., 2013). On the other hand, numerous experiments have been conducted to better un-

derstand driving factors such as the influence of turbulence on H2S mass transfer across the

water surface (Carrera et al. (2017), Matias et al. (2017), Wu (1995)).

Being able to model and directly quantify these emissions cannot only improve existing

model approaches by the possibility of investigating a wider range of cases but also to anal-

yse details of sewer construction such as the design of H2S hotspots and to find ways to

improve the design.

State of the art modelling approaches to describe these in-sewer processes are one-

dimensional approaches. However, previous research has shown that the flow processes

occurring are three-dimensional, for example air phase velocities in the headspace of a cir-

cular sewer (Edwini-Bonsu and Steffler, 2006), raising the question of how much more ac-

curate a three-dimensional Computational Fluid Dynamics (CFD) model could predict H2S

emissions. Furthermore, turbulent flow phenomena are ubiquitous in a sewer, starting in a

gravity pipe due to the flow velocities present, but also higher levels of turbulence can be

found in drop structures due to connections between pipe and manhole or due to obstacles.

A correct quantification of the effects of turbulence, which are highly three-dimensional, on

the H2S mass transfer is therefore crucial for the development of reliable models.

One experimental setup to analyse mass transfer under turbulent conditions that has

been investigated in previous publications are stirring tanks. A systematic sketch including

measures which are relevant for this publication can be found in Figure 7.1.

Figure 7.1: System sketch of stirring tank including relevant variables (Wu (1995), modified).

Wu (1995) analysed the influence of different factors such as the stirrer inundation (h/R

– ratio between water level above stirrer hand stirrer diameter R) and the stirring rate on O2

Chapter 7: Mass transfer under turbulent conditions

mass transfer. Carrera et al. (2017) investigated the influence of stirring rates on H2S mass

transfer and used CFD results from FLUENTTM 1.4 to extract hydrodynamic parameters

such as flow velocities. Being able to describe H2S emissions in a stirring tank with a CFD

model cannot only enable us to investigate more complex phenomena on a larger scale but

could also help avoiding lab experiments that come with a health risk and decrease the

number of experiments necessary.

In Teuber et al. (2019b), a new three-dimensional two-phase solver for the prediction of

H2S mass transfer called interH2SFoam has been introduced and tested on different applica-

tion examples. In this publication, we will extend a dynamic meshing solver in OpenFOAM

(interDyMFoam) by the functionalities of interH2SFoam as well as a functionality to compute

the mass transfer coefficient. This new solver called interDyMH2SFoam enables us to anal-

yse the effect of different stirring rates in a stirring tank on H2S mass transfer using a CFD

model.

In this paper, first an overview over the materials and methods is given. Then, the results

are evaluated and discussed. First, a graphic analysis is performed, then a sensitivity study

on parameters such as the grid size, the influence of the turbulent diffusivity, the choice of

turbulence model and the inundation of the stirrer is performed. In a quantitative analysis,

the simulated results are then compared to experimental results by Carrera et al. (2017) and

Wu (1995).

7.3 Methods

7.3.1 Geometry and mesh

In our model, the geometry analysed consists of a circular stirring tank. In the tank, a

Rushton turbine is placed with a defined stirring rate. The mesh has been discretized using

the OpenFOAM utility snappyHexMesh. The OpenFOAM solver interDyMFoam, which has

been extended here, is capable of handling mesh motion using both a static and a rotating

mesh domain. Therefore, the mesh motion has to be defined separately by specifying the

rotation speed.

7.3.2 Numerical model

OpenFOAM version 2.4.0 has been used for the work presented in this paper. A new solver

for two-phase simulations and mass transfer, which accounts for the temperature depen-

dency of the Henry coefficient, has been introduced in Teuber et al. (2019b) and is applied

to the turbulent case in this publication.

Hydrodynamic simulations

Hydrodynamic simulations are based on the two-phase flow solver interFoam which uses a

volume of fluid (VOF) approach that considers both phases as one fluid with changing fluid

properties. One set of Navier-Stokes equations is solved. The volume fraction of a phase is

stored as an additional variable and the phases are distinguished by an additional transport

equation. The equations are defined as follows (Rusche, 2003):

Mass conservation equation:

∇ · −→

U=0 (7.1)

Chapter 7: Mass transfer under turbulent conditions

Momentum conservation equation:

∂ρ−→

∂t+∇ · ρ−→

U−→

U=−∇prgh +∇ · µ∆−→

U+∇−→

U∇µ−−→

g·−→

x∇ρ(7.2)

Where prgh is the static pressure minus hydrostatic pressure:

prgh =p−ρ·g·h(7.3)

Volume of Fluid equation:

∂α

∂t+∇ · α−→

U+∇ · (1−α)−→

Urα=0 (7.4)

with the following parameters:

ρ=αρL+ρG(1−α)(7.5)

µ=αµL+µG(1−α)(7.6)

µ=µphys +µturb (7.7)

where −→

Uis the velocity field [m/s]; ρis the density [kg/m3]; t is time [s]; p is the pressure

[Pa]; µis the dynamic viscosity [Ns/m2]; −→

qis the acceleration vector due to gravity [m/s2];

−→

xis a spatial position vector [m]; αis a volume fraction or indicator function [-]; Ufis the

relative velocity between the phases [m/s]; the subscripts L and G denote the fluids water

(L - liquid) and air (G - gas). For the dynamic viscosity µ, the physical viscosity and the

turbulent viscosity are considered (see Equation 7.7).

The indicator function αis defined as:

α=









1 fluid L

0<α<1 transitional region

0 fluid G

(7.8)

The water surface is defined as the area where α= 0.5.

The accuracy of the hydrodynamic simulations has been assessed in (Teuber et al., 2019a).

Mass transfer simulations

Transport is considered by defining the H2S concentration as a passive tracer with an

advection-diffusion equation (Equation 7.9). The physical diffusivity Dphys as well as the

turbulent Schmidt number Scturb, which defines the turbulent diffusivity coefficient Dturb,

are defined by the user (Equation 7.10).

∂C

∂t+∇ · (−→

UC) = ∇ · (Dphys +Dturb)∇C(7.9)

with

Dturb =µturb/ρ

Scturb

(7.10)

Mass transfer at the wastewater-air interface is simulated using the approach defined by

Haroun et al. (2010a,b) as it has been implemented by Nieves-Remacha et al. (2015) and

Chapter 7: Mass transfer under turbulent conditions

Severin (2017). The approach is based on the interFoam solver and considers one additional

transport equation for both phases as outlined in Equations 7.9 and 7.10:

∂C

∂t+∇ · (−→

UC) = ∇ · ((Dphys +Dturb)∇C+φ)(7.11)

φ=−(Dphys +Dturb)C(1−He)

α+He(1−α)∇α(7.12)

It contains a modification of the original solver in its consideration of diffusivity. The origi-

nal approach by Haroun et al. (2010a) only considered molecular diffusivity Dphys, whereas

in this publication Dturb is included into the formulation. The difference can be explained

by the turbulence models used. Haroun et al. (2010a) performed simulations using Direct

Numerical Simulation (DNS), which is able to describe small scale turbulence effect without

any additional subgrid scale model, therefore no turbulent viscosity µturb is calculated. A

sensitivity study carried out in this paper will highlight the importance of including Dturb

in the test cases analysed in this publication.

In order to distinguish the species transport between the two phases, the following con-

ditions must be fulfilled at the interface:

He =CL

(7.13)

(Dphys,L+Dturb,L)∇CL= (Dphys,G+Dturb,G)∇CG(7.14)

Like density and viscosity, the concentrations and diffusion coefficients are considered as

single-phase properties depending on the phase fraction value α:

C=αCL+CG(1−α)(7.15)

The physical diffusivity Dphys is calculated by using a harmonic average:

Dphys =Dphys,L·Dphys,G

αDphys,L+ (1−α)Dphys,G(7.16)

The diffusion coefficients for Dphys,L and Dphys,G are defined by the user. Note that the tem-

perature dependency of these variables has to be accounted for.

The Henry coefficient depends on the temperature in the whole domain. Therefore, the

temperature dependency has been added in a way that the solver takes one global temper-

ature value as an input parameter.

The temperature dependent Henry coefficient is computed using the van’t Hoff equation

following Sander (2015):

He(T) = Hcpexp E1

T−1

TΦ (7.17)

He(T) = 10−3mols2

m2kg exp E1

T−1

TΦ·8.314 kgm2

s2molK ·298.15K(7.18)

Chapter 7: Mass transfer under turbulent conditions

Where Hecp is the Henry coefficient reported by Sander (2015) in the unit [mol/(m3Pa)], Eis

a temperature coefficient which is depending on the enthalpy of dissolution and defined as

2100 K (Sander, 2015), TΦis the standard temperature 298.15 K corresponding to 25◦C, Ris

the universal gas constant 8.314 kgm2

s2molK and Tis the temperature [K].

In order to combine dynamic meshing with mass transfer simulations, the respective

functionalities for mass transfer had to be transferred or connected to the interDyMFoam

solver, resulting in a new specialized solver interDyMH2SFoam.

Calculation of mass transfer coefficient and volumetric mass transfer coefficient

In the literature, two important physical values can be found, when analysing the mass

transfer coefficient. First of all, the mass transfer coefficient KLand second the volumetric

mass transfer coefficient KLa, which is a multiplication of KLawith a[m-1], the ratio be-

tween interfacial area and water volume. Many experiments report only KLabecause of

the difficulty of measuring the exact interfacial area in an experiment. In this respect, CFD

simulations offer the additional benefit of enabling the calculation of a.

The local mass transfer coefficient at the interface KL,local [m/(sm2)] is calculated accord-

ing to Haroun et al. (2010a) as

KL,local =−((Dphys +Dturb)∇C+φ)·nL

∆CL,local

(7.19)

The mass transfer coefficient for the area of the interface, KL, can then be computed by

integrating KL,local over the interfacial area:

KL=1

λZλ

0KL,localdξ(7.20)

where λis the interface length and ξthe curvilinear coordinate associated to the interface.

Furthermore, Carrera et al. (2017) propose a calculation of the mass transfer coefficient

based on the degassing technique, which has been used in this publication to cross validate

the simulated mass transfer coefficients. In the degassing technique, the measured decrease

of the H2S concentration in the water phase has been fitted to Equation 7.21:

CL,H2S=CL,0e−KL,H2Sa(t−t0)(7.21)

Where CL,0 is the initial concentration at t=t0. In this publication, the bulk H2S concentration

in the water phase as well as ahave been calculated and KL,H2Swas adjusted. To obtain

the mass transfer coefficient, an automated curve fitting using the programming language

Python has been carried out.

7.3.3 Turbulence

Turbulence and the impact of turbulence on flow can be simulated using different ap-

proaches. While DNS leads to most accurate results, the computational effort is highest

since the turbulence is discretized by using small enough cell sizes and time steps to dis-

play all vortices (Maric et al., 2014), leading to a high computational effort, that is in most

Chapter 7: Mass transfer under turbulent conditions

cases still not practically applicable today.

Using Reynolds averaging, the mean velocity is divided into an average velocity com-

ponent and a fluctuating velocity component resulting in a Reynolds stress tensor in the

Navier-Stokes equations (RANS). This tensor describes an additional eddy viscosity and

leads to new unknowns in the equations. These unknowns can be solved by using a two-

equation model, that adds two coupled transport equations in order to describe convection

and diffusion of turbulent energy. Depending on the model, these variables differ: k is the

turbulent kinetic energy, eis the turbulent dissipation, ωis the specific dissipation. For the

k-eturbulence models, three different formulations exist: the Standard (STD) k-emodel the

Realizable (Shih et al., 1995) and the RNG k-emodel using Re-Normalisation Group (RNG)

methods by Yakhot et al. (1992). RANS models used in this publication are the Standard k-e

(Launder and Sharma, 1974), the RNG k-e, as well as the k-ωShear Stress Transport (SST)

(Menter, 1993, 1994).

In turbulent flows, the viscous sublayer of a boundary layer close to the wall gets very

thin. In order to not discretize this sublayer by fine meshes, wall functions can be applied

when RANS turbulence models are used. These wall functions account for the fact that the

velocity profile normal to the wall follows a logarithmic profile, the log-law. The applicabil-

ity of these wall functions can be described using the dimensionless wall distance y+. It is

defined as:

y+=y

νυτ

(7.22)

where yis the absolute distance from the wall, υτdenotes the friction velocity and νdenotes

the kinematic viscosity. Depending on the chosen RANS turbulence model and the respec-

tive wall function, different values of y+and therefore different cell sizes in the nearfield

of the wall are necessary (Maric et al., 2014). Due to the log-law, RANS turbulence models

have smaller requirements in regard to cell size and time step than DNS and LES and are

therefore prone to lead to the smallest computational effort.

Large Eddy Simulation (LES) stands between DNS and RANS models in terms of com-

plexity. Large eddies are resolved spatially by choosing small enough cell sizes and small

eddies are modelled using a subgrid scale model which leads to smaller necessary cell sizes

than for RANS models (Maric et al., 2014).

7.3.4 Case study

The simulations are inspired by the experiments and simulations carried out by Carrera

et al. (2017). The domain is discretized as a cylindrical tank. A Rushton turbine is placed

above the bottom of the cylinder and two different levels of inundation are analysed. Fluid

properties can be found in Table 7.1. A Schmidt number of 1 is defined. An initial concen-

tration of 0.073 mol/m3, based on results by Carrera et al. (2017), is imposed.

Table 7.1: Fluid properties as defined in the CFD simulations.

Water phase Air phase

Density ρ[kg/m3] 1,000 1.20

Kinematic viscosity ν[m2/s] 1.0 10-6 1.48 10-5

Diffusion coefficient Dphys [m2/s] 1.8 10-9 1.74 10-5

Chapter 7: Mass transfer under turbulent conditions

In order to validate the amount of H2S emissions under steady state conditions for the

respective stirring speed, the first 25s of the simulations are carried out without any H2S

concentration present in the domain and by only performing hydrodynamic simulations

with the interDyMFoam solver. For this, an initial water level is defined depending on the

desired inundation depth using the setFields function in OpenFOAM.

After 25s, the OpenFOAM extension funkySetFields is used in order to apply a certain

H2S concentration Cin the areas were the fluid water is present (i.e. α>0.5). After the

H2S has been applied in the computational domain, another 10s of simulations are being

performed, this time by using the extended solver interDyMH2SFoam.

7.3.5 Initial and boundary conditions

The boundary conditions imposed to describe a stirring tank in the cylindrical geometry

consist of cylinder walls at the sides and the bottom of the cylinder, an upper open air

patch (top) as well as internal walls described by the shaft and stirrer. The usage of a

dynamic meshing algorithm also requires the definition of an internal boundary surface

between rotating and static mesh. The boundary surface between these regions is defined

as Arbitrary Mesh Interface (AMI). The boundary regions are displayed in Figure 7.2.

Figure 7.2: Computational grid used for the simulations including boundaries (red: stirrer,

grey: shaft, blue (outer cylinder): cylinder walls, blue (inner cylinder): AMI, blue (upper

plane): top).

At the walls, a no-slip condition is imposed. At the top, an atmospheric boundary

condition is defined to allow fluids to enter and leave the domain. This is achieved by

imposing a zero gradient (Neumann) condition to all variables except the pressure (prgh)

boundary condition, which is set to a fixed value of zero (Dirichlet) (atmospheric pressure).

For the variables of the turbulence model, i.e. k, eand ω, an arbitrary low value is defined

in the beginning of the simulations. During runtime, these values are corrected by the solver

using calculated values.

Chapter 7: Mass transfer under turbulent conditions

7.4 Results and discussion

7.4.1 Graphic analysis

A visual analysis of the model results led to the conclusion that a plausible behaviour of

the fluids in the stirring tank is reached. Figure 7.3 illustrates that the water surface forms

different vortex shapes depending on the stirring velocity.

Figure 7.3: Vortex shapes for different stirring rates (left: N=0.8333 1/s, middle: N=1.667

1/s, right: N=2.333 1/s).

Chapter 7: Mass transfer under turbulent conditions

7.4.2 Sensitivity analysis

A mesh and turbulence model sensitivity analysis has been conducted in order to determine

the best combination to achieve accurate results.

Regarding the grid convergence, the simulations reached grid convergence at a grid

size of 214,100 cells. The results of three different simulations performed on domains with

65,933, 214,100 and 374,127 cells were compared. The simulated Ux-velocities were com-

pared in six different locations and the cumulated Root Mean Square Error (RMSE) has

been calculated for each refinement step. The results yielded a RMSE=0.00306 m/s between

the coarse and medium sized mesh as well as an RMSE=0.00313 m/s between medium and

fine mesh, indicating a difference of 2.5% between the RMSE values of the different refine-

ment steps. Adding further refinement to the mesh would lead to a decrease of the y+-range

out of the acceptable range between 30 and 300 for the application of wall functions, which

is why the medium sized mesh has been chosen for further analysis. The final mesh is

presented in Figure 7.2.

The different turbulence models significantly influence the calculated mass transfer co-

efficient. Figure 7.4 shows the influence of the choice of turbulence model on the resulting

mass transfer coefficients. Both, RNG k-eas well as the SST k-ωturbulence model show

a significant increase of KLas well as KLa. A considerable difference can be found in the

behaviour of the STD k-emodel. The slope of the KLand KLacurves show a decreasing

trend for the observed stirring rates. The strong deviation can be explained by the differ-

ences in computed velocities and turbulent kinetic energy k, which influence the transport

and transfer of H2S in the model.

In a second step, the importance of including the turbulent diffusivity into the mass

transfer formulation has been analysed. Therefore, the simulation with a stirring rate of N

= 2.33 1/s and the RNG k-emodel has been performed with both, the interDyMH2SFoam

solver as well as a variation that does not include Dturb into the solution. For the first case,

the interDyMH2SFoam solver results in a mass transfer coefficient KL= 0.0040 m/s. Without

consideration of Dturb, a resulting KL= 0.0020 m/s is achieved, which is close to the value

computed for a stirring rate of N = 0.833 m/s, which amounts to KL= 0.0016 m/s. When

these two mass transfer coefficients are compared to the results and observations by Carrera

et al. (2017) as well as Wu (1995), it becomes clear, that a significant increase of the mass

transfer coefficient with increasing stirring rate is to be expected. All further discussion on

the results and model validation is therefore exclusively conducted by including the turbu-

lent diffusivity Dturb into the mass transfer approach and by using the results of the RNG

k-emodel as well as the k-ωSST model. Both models show the most similar results in the

computation of the mass transfer coefficient and coincide with the observations of Carrera

et al. (2017) as well as Wu (1995) and will be compared in more detail in a later step. In

both publications, at least a linear trend of increase has been observed with increasing stir-

ring rates, depending on the inundation depth and test case, even an exponential trend was

observed in some test cases.

Figure 7.5 shows the influence of inundation ratio h/R on the mass transfer coefficient

as well as the volumetric mass transfer coefficient. It can be seen that the value of the mass

transfer coefficient is not depending on the inundation ratio but the volumetric mass trans-

fer coefficient is. This is reasonable due to its dependency on the area to volume ratio a,

which changes due to the different volume, when a different water depth is being applied.

Chapter 7: Mass transfer under turbulent conditions

Figure 7.4: Influence of different turbulence models on the resulting mass transfer coefficient

(inundation ratio of h/R = 0.7).

The results furthermore agree with the observations by Wu (1995).

Overall, it has to be mentioned, that the sensitivity study is being complicated by the

fact that no experimental values were available for the exact chosen setup. The choice of

mesh and turbulence model are therefore made based on the experience of the authors and

literature from comparable cases.

Chapter 7: Mass transfer under turbulent conditions

Figure 7.5: Influence of stirrer inundation on the mass transfer coefficient.

Chapter 7: Mass transfer under turbulent conditions

Table 7.2: Setup of stirring tank in CFD simulations and known metrics of experimental

setup by Carrera et al. (2017) and Wu (1995).

CFD simulations Carrera et al. (2017) Wu (1995)

Turbine Rushton turbine

(6 blades)

Rushton turbine Rushton turbine

(6 blades, 4 baffles)

Tank diameter T [m] 0.250 approximately 0.123

(estimated from Reynolds

numbers,

a and figures in publication)

0.202

Stirrer diameter B [m] 0.165 0.100 (calculated

from Reynolds number)

0.101

Blade length L [m] 0.050 unknown 0.0253

Blade height W [m] 0.03 unknown 0.0202

Stirring rates N [1/s] 0.833, 1.333, 1.667, 2.000, 2.333 0.000, 0.833, 1.333, 1.667,

2.000, 2.333

0 - 900 (unknown intervals)

Area to volume ratio a [m-1] h/R = 0.2: 4.33 - 5.12

h/R = 0.7: 9.4 – 11.73

4.73 - 5.11 approximated

h/R = 0.10: 50.00

h/R = 0.62: 7.07

(assuming no surface

deformation due to stirring)

Reynolds number range [-] 22,687 - 63,525 8,333 - 23,333 0 - 153,015

Temperature [◦C] 20 20 25

7.4.3 Quantitative analysis

The H2S emissions in the stirring tank have been analysed for six different stirring rates with

two different inundation depths of the stirrer. The results were compared to experimental

investigations by Carrera et al. (2017), who calculated the KLaby using the degassing tech-

nique as outlined above, and to Wu (1995), who measured KLaof O2with the steady-state

sulphite feeding method in a stirring tank. The setups of the different domains are listed in

Table 7.2. The resulting Reynolds numbers were computed by

Re =ρLNB2

µL

(7.23)

Where ρLis the density of the water phase [kg/m3], µLis the dynamic viscosity of water

[kg

ms ], Bis the stirrer diameter [m] and Nis the stirring rate [1/s].

Figure 7.6 and Figure 7.7 illustrate the published results (as extracted from the publi-

cations). Due to the different range of Reynolds numbers, the authors see difficulties in

directly comparing the results to the observations by Carrera et al. (2017) and Wu (1995).

Therefore, the different results analysed are compared by plotting the Reynolds number

in a range between 20,000 and 65,000 against the local mass transfer coefficient KL,local in

Figure 7.8. To do so, the range of the results of Carrera et al. (2017) have been divided by

an estimated surface area derived from the given value for a and the figures shown in the

publication. Furthermore, the range of the Reynolds numbers given in the publication has

been used to derive the stirrer diameter which has been calculated to amount B=0.100 m.

The experimental results by Wu (1995) have been converted to KL,local of H2S by applying

the conversion factor KL,H2S

KL,O2=0.64 +/−0.24 that has been observed by Carrera et al. (2017),

resulting in the range of values displayed in Figure 7.8. Furthermore, the values have been

divided by a constant surface area calculated using the given tank diameter.

The results displayed in Figure 7.8 allow only a limited comparison of the results be-

cause of multiple factors of uncertainty such as the unknown exact geometry analysed by

Carrera et al. (2017) and the uncertainty when applying a conversion factor from O2to H2S

Chapter 7: Mass transfer under turbulent conditions

Figure 7.6: Influence of stirring rate and inundation ratio on mass transfer (Wu, 1995).

emissions for the results of Wu (1995) as well as the assumption of a constant water surface

area for different stirring speeds.

From the results by Carrera et al. (2017), the two highest stirring rates result in Reynolds

numbers in the area of interest (N=2 1/s leading to Re=20,000 and N=2.33 1/s leading to

Re=23,333). In general, the k-ωSST turbulence model computes slightly higher mass trans-

fer rates than the RNG k-turbulence model. But the simulated KL,local values agree well

with the results by Carrera et al. (2017) and Wu (1995).

100

Chapter 7: Mass transfer under turbulent conditions

Figure 7.7: Influence of stirring rate on mass transfer (Carrera et al., 2017).

101

Chapter 7: Mass transfer under turbulent conditions

Figure 7.8: Comparison of relation between Reynolds number and volumetric mass transfer

coefficient (Wu, 1995, Carrera et al., 2017).

102

Chapter 7: Mass transfer under turbulent conditions

7.5 Conclusions

The influence of turbulence on H2S mass transfer has been subject to intensive research

in the past years. So far, a general understanding has been gained from experimental in-

vestigations but the influence has not been simulated using a numerical three-dimensional

model.

In this publication, a CFD model has been developed to describe H2S mass transfer in

a stirring tank using a dynamic meshing functionality in order to predict the mass trans-

fer coefficient KLas well as the volumetric mass transfer coefficient KLaunder turbulent

conditions. After analysing the sensitivity of different parameters, namely, the choice of

turbulence model, the turbulent diffusivity, the mesh size and the stirrer inundation, the

resulting mass transfer coefficients in dependency of the Reynolds number were compared

to experimental findings by Carrera et al. (2017) and (Wu, 1995). The k-ωSST turbulence

model generally calculated the highest mass transfer rates, whereas the RNG k-emodel

yielded comparable results. The STD k-emodel was not able to account for the increase of

mass transfer with increasing stirring rates. The mass transfer coefficient KLwas found to

be independent of the stirrer inundation, while the volumetric mass transfer KLaincreased

with decreasing stirrer inundation. Furthermore, the importance of including the turbulent

diffusivity Dturb into the solver for this application has been shown.

A strong influence of the Reynolds number as well as the depth of stirrer inundation has

been observed in alignment with existing publications. When the local mass transfer coeffi-

cient KL,local of the simulated results is compared to the results from previous publications,

a good agreement can be found.

For future work, the model is currently being further validated using experimental re-

sults from a stirring tank at TU Berlin and will be used afterwards to describe a turbulent

H2S reactor. The simulation results will be used to gain a better understanding on the effect

of the hydraulic design on the measured H2S concentrations.

Nevertheless, these results demonstrate that the new model for three-dimensional sim-

ulation of mass transfer in dynamic meshing environments is able to accurately describe

the H2S mass transfer including its influencing parameters (stirrer inundation h/R, area to

volume ratio a, Reynolds number Re). This can help avoiding risky lab experiments for H2S

emissions in the future and in a later step help to improve the design of H2S hotspots in

sewer systems.

Acknowledgements

The funding provided by the German Research Foundation (DFG) within the Research

Training Group “Urban Water Interfaces" (GRK 2032-1) is gratefully acknowledged.

103

Chapter 8

Supplementary contributions

8.1 Advantages of three-dimensional flow simulations

This study was published as:

Teuber, K., Broecker, T., Elsesser, W. & Hinkelmann, R.: Beyond shallow water flow:

Navier-Stokes simulations with OpenFOAM, BIMoS Day: Shallow Water Flow Simulations,

Berlin International Graduate School in Model and Simulation based Research, TU Berlin,

Berlin, Germany, 2017.

This is the postprint version of the abstract submitted to the event.

Abstract

In the field of hydraulic engineering, shallow water flow models using the depth-averaged

Navier-Stokes equations are often suitable to analyze flow phenomena. However, in cer-

tain application areas where a hydrostatic pressure distribution within the flow field is not

given, these models reach their limits and alternatives have to be looked for. The model

used to describe free surface flow is the interFoam solver as it is implemented in the open

source software OpenFOAM. A collection of examples for these special applications using

single and two-phase flow approaches will be presented in this talk. Beginning with a sim-

ple example of flow over a ground sill, the talk moves on to complex flow and transport

simulations over ripples on a river bed. Apart from natural systems, the model can also

be used to simulate flow in technical systems such as a sewer stretch containing different

hydraulic structures such as weirs, flow transitions and quickly-varying shapes.

104

Chapter 8: Supplementary

8.2 Possibilities extending the interH2SFoam solver to multicom-

ponent reactive transport

This study was submitted as:

Dixit, A., Teuber, K., Barjenbruch, M., Stephan, D. & Hinkelmann, R.: Extension of a 3D

two-phase flow model to multicomponent reactive transport for odour and corrosion control

in sewer systems, Workshop on Applications of Multi-scale Approaches in Environmental

Chemistry (AMARE), Rennes, France, 2019.

This is the preprint version of the abstract submitted to the conference.

Abstract

Biological corrosion of sewers and sewage treatment plants constitutes a serious problem

and its effects result in the loss of billions of dollars every year. Changing demography and

more efficient use of water resources will lead to the reduction of the average volume of

waste water and leads to higher residence times in the sewer canals. Due to climate change,

i.e. warmer temperatures, the waste water in the canal will become more anaerobic. There-

fore, sewer networks with a concrete construction are subjected to various mechanisms that

subject it to rapid degradation. Due to the anaerobic conditions in sewage, sulfate present

in the waste water can be reduced to sulfide by sulfate-reducing bacteria residing in the

biofilms on the walls of the pipelines. For more than 70 years, researchers have been com-

mitted not only to study the processes for odour and corrosion but also creating empirical

and conceptual models for explanations. However, within the last 20 years a deeper under-

standing has been gained thanks to the efforts of research groups in Denmark and Australia.

Nearly all current models are confined to a one-dimensional approach which is very suitable

but is unable to sufficiently capture the turbulent effects. However, for processes which are

affected by the concentration profiles (e.g. H2S formation, mass transfer) and the air flow

as well as for the surroundings of drops, steps and hydraulic jumps, a three- dimensional

approach should be preferred accounting for water and gas phase. For this purpose, a high-

resolution three-dimensional model in OpenFOAM for water and air flow, multi-component

reactive transport and mass transfer between the water and air phase must be developed.

Teuber et al. (2019b) have developed a 3D two-phase (water, air) flow and transport model

that can account for temperature and pH dependency of the mass transfer of H2S. This

poster concentrates on the conceptual and computational extension of the model of Teuber

et al. (2019b) for reactive transport using and probably extending OpenFOAM libraries and

the model validation.

105

Chapter 8: Supplementary

8.3 Research on further multiphase CFD applications carried out

at the Chair of Water Resources Management and Modeling of

Hydrosystems, TU Berlin

This study was published as:

Broecker, T., Teuber, K., Elsesser, W. & Hinkelmann, R.: Multiphase Modeling of Hy-

drosystems Using OpenFOAM, in: Gourbesville P., Cunge J., Caignaert G. (eds) Advances

in Hydroinformatics, Springer Water, 1013-1029, Springer, Singapore, 2018.

This is the abstract of the book chapter (postprint).

Abstract

This paper presents three computational fluid dynamics applications regarding multiphase

modeling of hydro systems with the open source software OpenFOAM. The first model

investigates flow processes of groundwater and surface water using an integral approach

which solves the three-dimensional Navier–Stokes equations, extended by the consideration

of porosities. For the validation, seepages through homogeneous dams with impervious

foundations were compared with analytical and numerical solutions. A further application

examines the water–air interface in sewer systems.The focus of the model lies on the de-

scription of in-sewer water–air flow and transformation processes, reaeration and hydrogen

sulfide emission which highly depend on the three-dimensionality of the hydraulic behavior

in the closed duct. A test case analyzing the hydraulic behavior in a sewer stretch showed

a good agreement of the numerical results with measured water levels. In the third model,

fluid–structure interaction is investigated applying FOAM Extend Project. Calculations of

the fluid phase are linked with the solid phase via a coupling algorithm to achieve an equi-

librium state. To describe the time-varying position of the fluid boundary, caused by the

structural response, dynamic meshes are considered. A technical case, consisting of the air

flow around a thin tower as well as a natural case, describing the water flow around aquatic

vegetation and its response, were examined.

106

Chapter 8: Supplementary

8.4 Application of free-surface flow and transport simulations to

natural river beds

This study was published as:

Broecker, T., Elsesser, W., Teuber, K., Özgen, I., Nützmann, G. & Hinkelmann, R.: High-

resolution simulation of free-surface flow and tracer transport over streambeds with ripples,

Limnologica, 68: 46-58, 2018.

This is the abstract of the article (postprint).

Abstract

This study presents a novel high-resolution simulation of free-surface flow and tracer reten-

tion over a streambed with ripples based on varying ripple morphologies, surface hydraulics

and the transport of a tracer pulse from surface water to surface dead zone. For the simu-

lations, the computational fluid dynamics (CFD) model OpenFOAM was used to solve the

three-dimensional Navier-Stokes equations in combination with an implemented transport

equation. Pressure gradients at the streambed were used to account for hyporheic exchange,

assuming water flow from high pressure zones to low pressure zones. Flow velocities, rip-

ple sizes and spacing showed to significantly affect these pressure gradients, but also the

transport of a passive tracer at the streambed, which was not investigated so far. Due to

the velocity field, large parts of the tracer mass were transported alongside the main stream

above the ripples. Tracer mass reaching the space between the ripples was temporarily re-

tained due to low velocities and recirculations. It was shown that the retention is depending

on the ripple size and space between the ripples as well as on the flow velocity. Decreasing

ripple sizes and higher flow velocities lead to a smaller tracer retention. Furthermore we

showed that the ripple length to height ratio controls the generation of recirculation zones

which affect the residence time of the tracer significantly. Ripple spacing leads to temporar-

ily higher tracer concentration at the streambed, but smaller tracer retention. We conclude

that the impact of the streambed morphology on the hydraulics in combination with tracer

retention should be addressed for a comprehensive understanding of compound movement,

exchange and transformation within the hyporheic zone.

107

Chapter 9

Synthesis

This Chapter synthesizes the specific outcomes of Chapters 1 to 8.

9.1 Conclusions

Sewer networks are complex systems containing locations of highly complex turbulent mul-

tiphase flows and complex biochemical processes, for example the emission of H2S. Existing

model approaches tend to describe these processes in one-dimensional models, leading to

simplified estimations of the overall emissions, which are valid for long sewer networks.

They are not suitable for predicting the small-scale effects of certain geometrical features in

the domain, especially in highly turbulent locations. A three-dimensional two-phase model

can be used to investigate the local effects of the sewer geometry on a smaller scale. In-

cluding mass transfer into the model enables a detailled analysis of locations with high

emissions.

The main outcome of this work are two OpenFOAM solvers that have been developed by

extending existing solvers, namely the H2S mass transfer solver for (waste)water-air sim-

ulations with static meshes (interH2SFoam) (Chapter 5), and the interDyMH2SFoam solver

for dynamic meshes (Chapter 6). The two-phase simulation approach is based on the VOF

method and mass transfer is described based on the work of Haroun et al. (2010a). After

performing a thorough analysis and validation on the applicability of the hydrodynamic be-

haviour of the VOF method in closed conduits (Chapters 3 and 4), the first solver has been

studied in several test cases based on the equilibrium conditions of concentrations in the

two phases, ranging from benchmark tests with analytical solutions to a real world appli-

cation of a complex sewer geometry. The second solver has been used to validate turbulent

mass transfer in a laboratory scale experiment.

The simulation results are promising, the equilibrium conditions as well as the mass transfer

rates could be reproduced by the solvers.

A short overview of supplementary work related to the research topic was given in Chapter

9.1.1 General outcomes

Before more specific conclusions of the different parts, hydrodynamics, transport and mass

transfer are drawn, some general conclusions can be drawn:

108

Chapter 9: Synthesis

• It is possible to describe water-air mass transfer in sewer systems or - more general, in

closed conduits - by using a three-dimensional CFD model.

• The three-dimensional two-phase model is based on the VOF method and accounts

for general mass transfer with an approach by Haroun et al. (2010a).

• Using a two-phase CFD model requires detailled knowledge about the geometry and

flow conditions at the model boundaries in order to obtain stable simulations. It fur-

thermore needs a high-quality and high-resolution three-dimensional computational

mesh requiring more preprocessing than a one-dimensional model approach. Fur-

thermore, the mesh quality significantly influences the quality of the results and grid

convergence needs to be assessed. The choice of turbulence model is crucial for the

model results and needs to be carefully evaluated for each test case.

• The three-dimensional simulations lead to a considerably higher computation time

than one-dimensional approaches due to their high complexity, making the use of

high performance computers necessary. The simulations offer the benefit of detailled

insight into the local flow and mass transfer processes.

• Overall, the solvers yield accurate results for the flow properties, mass transfer rates

and equilibrium conditions.

• Within the DFG Research Training Group UWI, different interfaces in sewer systems

have been identified when analysing H2S emissions. These interfaces are namely

the biofilm-(waste)water interface, the (waste)water-air interface as well as the air-

(biofilm)-(concrete) wall interface. Three different projects focusing on different inter-

faces were established. Although the three projects focus on different interfaces, they

can profit from each other’s work in various ways. Within the work of this thesis,

a major step towards the collaboration between the different projects has been made

by solver customizations that can lead to a direct exchange of results between exper-

iments and simulations. In Chapter 6, the solver interH2SFoam was customized to

the direct demands of the experimental site. Therefore, the pH value and the total

dissolved sulphide have been defined as input parameters for the model and the H2S

concentration in the air phase has been calculated as partial pressure in the unit ppm

to enable direct comparison with measurements. Furthermore, the single-phase trans-

port of a tracer in the sewer pilot plant, as analysed in Chapter 5, led to valuable

insights of concentration distributions between the concrete probes placed in the pilot

plant. These insights can help interpret future results when corroded concrete probes

are analysed and help design experiments.

9.1.2 Outcomes of hydrodynamic modelling

In Chapters 3 and 4, hydrodynamic three-dimensional two-phase simulations in closed con-

duits are performed using the VOF method as it is implemented in OpenFOAM’s solver

interFoam. The results are compared to experiments, analytical solutions and existing CFD

models in order to assess the applicability of the solver for closed conduits and complex

geometries as they can be found in sewer systems. The following conclusions are drawn for

the hydrodynamic behaviour of the two-phase CFD simulations in closed conduits:

• Three different test cases were analysed. First, a two-dimensional single-phase flow

of water has been analysed. The domain consisted of a two-dimensional channel

bounded by an upper and a lower wall with a hill structure on its bottom. Results of

109

Chapter 9: Synthesis

simulations with different RANS turbulence models were compared to experiments

by Almeida et al. (1993). Overall, the Standard k-eturbulence model yielded best

agreement and was chosen for the remaining simulations.

• For the free-surface flow over a ground sill, a comparison to analytical solutions using

the continuity and Bernoulli equation led to a good agreement.

• The example of free-surface flow over a ground sill has been deeper analysed by

varying the geometric shape of the sill and the flow conditions. Overall the simulations

led to plausible results.

• A complex sewer geometry as it has been described by Bayón et al. (2015) has been

simulated and compared to experiments from a 1:20 scale model as well as to results

from an existing CFD model. The difference of this model to the existing CFD model is

that the model yielded stable results although it has been simulated as a closed duct.

The previous CFD model had to be simulated with an open atmospheric boundary

because of stability problems.

• For the complex sewer geometry, the problem of instabilities for the setup of a closed

duct had to be overcome. Therefore, the water level at the outlet boundary has not

been defined by a boundary condition but by adding a weir structure in the close

proximity of the outlet. Finding a way to describe complex flows in closed conduits

has been a crucial part for this thesis due to the overall aim of investigating H2S

emissions in sewers. Being capable of describing such closed systems enables the

future analysis of H2S spreading in a sewer or in connecting points.

9.1.3 Outcomes of single-phase transport modelling

In Chapter 5, the outcomes of single-phase transport simulations are presented. Three

different test cases are analysed. In a first test case, the limitations of adding a standard

advection-diffusion equation to simulate single-phase transport phenomena in a two-phase

system are shown. Then, two different ways of overcoming this problem are presented.

The first solution is to describe the domain as a single-phase system and describing the

interface as a boundary condition. The second approach uses the mass transfer approach

developed by Haroun et al. (2010a) with a very low Henry coefficient. This approach enables

a simulation of phase-constrained transport in a two-phase system. For the single-phase

transport, the following conclusions are drawn:

• In a first step, a rectangular pipe with a length of 15m, a height of 1m and a width

of 1m has been simulated as a free-surface flow with the interFoam solver and an

additional standard advection-diffusion transport equation. The aim was to simulate

tracer spreading in the water phase of the pipe. The simulation results showed a

spreading of the tracer into the air phase when the tracer is placed in closed proximity

to the water surface. This spreading can be considered unphysical when thinking of

a dye tracer. Numerically it can be explained by the governing equations of the VOF

method: Only one set of Navier-Stokes equations is solved and the phases are only

accounted for by an additional transport equation. The solver is therefore not able

to account for the change of phases accurately when an advection-diffusion transport

equation is applied.

• As an alternative solution, the interHarounFoam solver which accounts for mass trans-

fer by using a single-phase formulation and which depends on the Henry coefficient

110

Chapter 9: Synthesis

has been applied to the same setup of a rectangular pipe. The Henry coefficient was

specified with a very low value (10−6) and the tracer remained in the water phase.

• This solution has been furthermore applied to the complex sewer geometry. This

setup has been chosen due to the high levels of turbulence present in the domain. The

results were able to prove that the tracer travelled through locations of high turbu-

lence without spreading into the air phase. Therefore, the conclusion can be drawn

that applying the interHarounFoam solver to two-phase simulations with single-phase

transport phenomena is one solution to retain the tracer in one phase.

• As an alternative solution, a single-phase simulation approach has been presented.

This approach can be of interest when the behaviour of only one phase is of special

interest. As an example, tracer spreading in the air phase of the pilot plant has been

simulated by only considering the air phase. The water surface behaviour has been

assumed to be stratified and the resulting air phase velocities at the interface and in

the headspace were assumed to be following observations from existing publications.

The velocities at the interface were described by applying a slip-condition.

• The results showed a tracer accumulation behind the concrete probes. This observation

can be used for interpreting the results of concrete corrosion.

• The simulations furthermore showed that the singe-phase approach is a way to de-

scribe transport in one phase in a computationally efficient way because in this case

only one half of the pipe needs to be discretized. It comes with the limitation of being

only applicable for uniform water levels.

9.1.4 Outcomes of mass transfer modelling

In Chapters 6 and 7, the capabilities of the mass transfer solver introduced by Haroun et al.

(2010a) are explored and novel extensions are presented. This is done by perfoming mass

transfer simulations for equilibrium conditions and comparing the results to analytical solu-

tions. The advantages of a three-dimensional formulation are shown by applying the solver

to the complex sewer geometry. The functionalities of the extended solver are then applied

to a dynamic meshing solver which is used to investigate the solver’s capabilities of describ-

ing turbulent mass transfer. This is done by simulating H2S mass transfer in a stirring tank

and comparing the mass transfer rates to results of two laboratory scale experiments. The

detailled conclusions are listed in the following:

Extension of interHarounFoam solver

• In a first step, several extensions to the interHarounFoam solver were developed in

order to create a solver that is better suited to describe the specifics of H2S mass

transfer across the (waste)water-air interface.

• The temperature dependency of the Henry coefficient has been included by imple-

menting the van’t Hoff equation. The Henry coefficient then depends on a temper-

ature for the whole computational domain which has to be defined by the user as a

part of the preprocessing.

111

Chapter 9: Synthesis

• The equilibrium between H2S and HS-in the wastewater has been implemented fol-

lowing Hvitved-Jacobsen et al. (2013). The model takes the values of the total dis-

solved sulphide and the pH and computed the resulting H2S concentration in the

water phase. This function is and has been integrated by using the OpenFOAM utility

swak4Foam, the user can therefore decide whether the H2S concentration should be

computed or manually defined.

• Furthermore, the partial pressure of H2S in the air phase is being computed by an-

other extension. Because H2S concentrations in the air phase are measured as a partial

pressure and are then interpreted as odour intensity, this functionality offers the user

the possibility of directly comparing simulation results to measured values. This ex-

tension therefore results in the description of odour in the air phase of a sewer.

• In a final step, the functionalities of the interH2SFoam solver have been connected to

a dynamic mesh solver (interDyMFoam) in OpenFOAM in order to describe moving

geometries resulting in a new solver interH2SDyMFoam, which has been used for the

validation of turbulent mass transfer simulations.

Validation and application of the interH2SFoam solver

• The solver was tested under equilibrium conditions in a tank partially filled with wa-

ter. As transported pollutant, H2S was chosen and the respective Henry coefficient

as well as diffusivities were applied. The development of equilibrium conditions af-

ter a certain amount of time has been observed and the ratio between water phase

and air phase concentration could be related to the Henry coefficient applied. The

extensions regarding temperature dependency, equilibrium conditions and calcula-

tion of the partial pressure in the air phase were compared to analytical solutions by

Hvitved-Jacobsen et al. (2013).

• The solver was then applied to a rectangular pipe with a length of 15 m and two dif-

ferent flow velocities. The resulting air flow velocities were compared to experimental

results by Bentzen et al. (2016). The mass transfer simulations showed a minor impor-

tance of the mass transfer simulations in this case due to the strong advective transport

in the domain.

• Therefore, the solver has been applied to the complex sewer geometry, which contains

high levels of turbulence and three-dimensional flow velocities. Especially in the point

of highest turbulence, the hydraulic jump, a higher mass transfer could be observed

by analysing contour plots of the tracer distribution within the domain.

• The analysis in the rectangular pipe and the complex sewer geometry were able to

show, that especially in complex geometries the application of the three-dimensional

model approach can lead to valuable insights and could be used for the analysis of

H2S hotspots or for design optimizations.

Validation under turbulent conditions

• The interH2SDyMFoam solver has been used to analyse turbulent mass transfer in a

stirring tank and the results have been compared to the findings by Wu (1995) and

Carrera et al. (2017).

112

Chapter 9: Synthesis

• A circular tank with a Rushton turbine placed above the tank bottom has been simu-

lated with different stirring rates and stirrer inundations, where the stirrer inundation

is described by the ratio between water level above the Rushton turbine and stirrer

diameter.

• A sensitivity analysis has been perfomed regarding the choice of turbulence model,

the turbulent diffusivity and the stirrer inundation. The mass transfer coefficient KL

as well as the volumetric mass transfer coefficient KLa were compared. A sensitivity

analysis of the turbulence models STD k-e, RNG k-eand k-ωSST showed a good

performance of the RNG k-eand k-ωSST turbulence models. The stirrer inundation

influenced the volumetric mass transfer coefficient KLa, whereby for smaller inunda-

tions a higher KLa has been observed. The results of the mass transfer coefficient KL

were independent of the stirrer inundation. This observation of the behaviour of KLa

agreed with the observations by Wu (1995) and is reasonable when the fact is con-

sidered that by decreasing the stirrer inundation the water volume and therefore the

interfacial area to volume ratio is increased.

• The simulated mass transfer rates were furthermore directly compared to the results

by Wu (1995) and Carrera et al. (2017). A comparison on the basis of the Reynolds

numbers showed a good agreement when the local mass transfer coefficient KL,local is

compared. KL,local is independent of the interfacial area and the water volume in the

domain. The Reynolds number in a stirring tank further accounts for the diameter

of the stirrer of the turbine. These factors appear to be the most influencing factors

for mass transfer in a stirring tank and therefore led to a good comparability of the

different setups analysed.

• Overall, the work enabled a direct comparison between experimental and simulated

results of turbulent mass transfer in a stirring tank. The simulated results agreed well

with the experimental results and can therefore be considered as plausible, enabling a

usage of the solver for more complex application cases.

9.1.5 Final notes

In general, two new CFD solvers for H2S mass transfer in (waste)water-air systems have been

developed and validated under turbulent conditions. The solvers are capable to describe

transport and mass transfer of H2S in sewers. By modifying the Henry coefficient and

removing the functionalities for equilibrium conditions between H2S and HS-as well as the

calculation of the partial pressure, it could be furthermore applied to different species such

as oxygen or methane.

Additionally, the interFoam solver has been analysed regarding its ability to describe the

hydrodynamics of water-air flow in closed conduits and the possibility to describe single-

phase transport.

9.2 Outlook

A number of open issues are identified and could be addressed in future research in order

to enhance the model’s applicability to real-life conditions:

• Further validation is needed in order to test the approaches under field conditions:

First, they can be validated using the setup of a stirring tank as presented in Chapter

7, which is currently being prepared at TU Berlin. The advantage is that in this case

113

Chapter 9: Synthesis

the exact specifics of the geometry are known. Furthermore, a more advanced experi-

mental setup using an H2S mass transfer reactor will be set up. The reactor uses real

wastewater taking one step further into the direction of a model that can be applied

under real life conditions. Furthermore, data from the BWB pilot plant will be used to

validate the model.

• Regarding the temperature dependency, the model can be further extended to account

for the temperature dependency of fluid properties such as the density, viscosity and

H2S diffusivity which currently have to be defined by the user.

• When thinking about a network-wide application of the model in terms of a decision-

support tool, the topics upscaling and parametrization should be addressed: Three-

dimensional simulations lead to very high computation times compared to one-dimen-

sional simulations. Looking at sewer systems that consist of networks of multiple

kilometers of pipes and corrosion processes that take months to years, the simulations

are not feasible. Extracting important parameters from the CFD model and feeding

them into a one-dimensional approach, could improve the existing approaches. In the

WATS model, this is partially already done by accounting for drop structures in the

sewer network. It could be done for other points of high turbulence as well as, for

example, hydraulic jumps in connection points between manholes and pipes.

• In order to include additional interfaces of the sewer, the biofilm can be integrated

into the model: This can be done in a simple way by setting a boundary condition

which is already possible with the current version of the model. A more complex

consideration would demand the model to account for reactive transport in order to

include the complex processes at the biofilm-(waste)water interface.

• Reactive transport modelling can furthermore help including the effect of counter-

measures: The interHarounFoam solver implemented by Nieves-Remacha et al. (2015)

offers the possibility of adding reactive transport by including the production term

into the mass transport equation. This possibility should be first explored. If the ap-

proach is too simplified, OpenFOAM offers a solution for more complex reactions in

the chemFoam solver. However, this solver is only designed for single cell, i.e. zero-

dimensional, cases and a way to apply it to multi cell problems with multiple phases

would be subject to future research.

The future work regarding further extension of the model is illustrated in Figure 9.1.

114

Chapter 9: Synthesis

Figure 9.1: Sewer processes (following Hvitved-Jacobsen et al. (2013))

115

Appendix A

Flow charts of solver extensions

116

Appendix A

A.1 Integration into the solution procedure

A.1.1 interH2SFoam solver

Figure A.1 illustrates the solution procedure carried out within the interFoam solver.

Figure A.1: Flow chart of interFoam solver (following Devolder et al. (2015), Lopes et al.

(2017))

In Figure A.2, the solution procedure of the interH2SFoam solver is shown. The exten-

sions made to the interFoam solver (see Figure A.1) are highlighted in green.

117

Appendix A

Figure A.2: Flow chart of interH2SFoam solver (following Devolder et al. (2015), Lopes et al.

(2017)), solver extensions compared to interFoam solver in green

118

Appendix A

A.1.2 interDyMH2SFoam solver

Figure A.3 illustrates the solution procedure carried out within the interDyMFoam solver,

which is based on the interFoam (see Figure A.1) solver but includes a solution step for

mesh motion.

Figure A.3: Flow chart of interDyMFoam solver (following Devolder et al. (2015), Lopes

et al. (2017))

In Figure A.4, the solution procedure of the interDyMH2SFoam solver is shown. The

extensions made to the interDyMFoam solver (see Figure A.3) are highlighted in green.

119

Appendix A

Figure A.4: Flow chart of interDyMH2SFoam solver (following Devolder et al. (2015), Lopes

et al. (2017)), solver extensions compared to interFoam solver in green

120

Appendix A

A.2 Code extensions and integration into the object-oriented frame-

work

A.2.1 File structure

In OpenFOAM, solver modifications should be stored in a specific user directory. In contrast

to the official OpenFOAM solvers of the current release, which are stored in the directory

$WM_PROJECT_DIR (full path of the directory illustrated in Figure A.5), custom solvers

can be called via $WM_PROJECT_USER_DIR (full path of the directory illustrated in Fig-

ure A.6) (Nilsson, 2019).

In order to create the interH2SFoam solver, the folder containing the interFoam solver

should be copied into a subdirectory of $WM_PROJECT_USER_DIR with the same file

structure as the original interFoam solver (/applications/solvers/multiphase).

Then, modifications are to be made to files createFields.H and interFoam.C. The file CEqn.H

will be newly created. In subfolders Make and interDyMFoam, further modifications will be

made. The other subfolders of the original interFoam solver refer to other derived solvers

and can be deleted. This refers to the folders interMixingFoam,LTSInterFoam and porousIn-

terFoam. The remaining files are listed in Figure A.6. The files that need to be modified in

a next step are highlighted in green. The file structure of the interDyMH2SFoam solver is

illustrated in Figure A.7. As shown in Figures A.6 and A.7, the interDyMH2SFoam solver

is placed as a subdirectory within the interH2SFoam solver directory.

In this thesis, transport simulations with an advection-diffusion transport equation were

simulated using a modified interFoam solver called passiveScalarInterFoam. The file struc-

ture is comparable to the interH2SFoam solver and is illustrated in Figure A.8.

In the solver directory, the file interFoam.C is then renamed to interH2SFoam.C or pas-

siveScalarInterFoam.C. The same applies for the interDyMFoam solver: The file interDyM-

Foam.C is renamed to interDyMH2SFoam.C. The modifications within the files will be out-

lined in the following Sections.

121

Appendix A

Figure A.5: File structure of interFoam solver

122

Appendix A

Figure A.6: File structure of interH2SFoam solver

123

Appendix A

Figure A.7: File structure of interDyMH2SFoam solver

124

Appendix A

Figure A.8: File structure of passiveScalarInterFoam solver

125

Appendix A

A.2.2 Make/files

The first step is to modify the files within the Make directory in order to ensure a correct

initialization of the solver. In Listing A.1, the content of the file files in subfolder Make are

illustrated for the interH2SFoam solver (Nilsson, 2019). For the interDyMH2SFoam and

passiveScalarInterFoam solver, the modifications are made accordingly.

interH2SFoam.C

EXE = $( FOAM_USER_APPBIN )/ interH2SFoam

Listing A.1: Modified version of Make/files

126

Appendix A

A.2.3 createFields.H

The relevant field variables and variables to extend the interFoam solver are defined within

the file createFields.H (Listing A.2). In this file, the tracer is being defined as a volScalarField

C. The temperature-dependent Henry coefficient is being computed by looking up the user-

defined temperature from the case directory. Furthermore, the diffusivity is being read from

the case directory.

Info << " Reading field C" << endl ;

volScalarField C

(

IOobject

(

"C",

runTime . timeName () ,

mesh,

IOobject :: MUST_READ ,

IOobject :: AUTO_WRITE

mesh

);

Info << " Reading transportProperties \n" << endl ;

IOdictionary transportProperties

(

IOobject

(

"transportProperties",

runTime . constant () ,

mesh,

IOobject :: MUST_READ_IF_MODIFIED ,

IOobject::NO_WRITE

)

);

// water phase ( phase 1)

dimensionedScalar DT1

(

transportProperties . lookup (" DT1 ")

);

// air phase

dimensionedScalar DT2

(

transportProperties . lookup (" DT2 ")

);

dimensionedScalar Schmidtnumber

(

transportProperties . lookup (" Schmidtnumber ")

);

volScalarField DT

(

"DT",

DT1 * DT2 / ( alpha1 * DT1 + (1 - alpha1 ) * DT2)

127

Appendix A

);

// Constant temperature in domain

dimensionedScalar Temp

(

transportProperties . lookup (" Temp ")

);

// Temperature dependent Henry coefficient for H2S , normal temperature of

25.15 degree

dimensionedScalar He

(

"He",

1/(0.001 * Foam :: exp (2200.0*(1/ Temp -1/298.15) ) * 8.314 * Temp )

);

Listing A.2: Initialization of relevant variables in createFields.H

128

Appendix A

A.2.4 CEqn.H

Transport and mass transfer are described within the file CEqn.H which is displayed in

Listing A.3.

{

surfaceScalarField phiCi =

(

fvc :: interpolate (( DT +( nut/ Schmidtnumber ))) *

(1-He)

/ ( fvc :: interpolate ( alpha1 )+(1 - fvc :: interpolate ( alpha1 ))* He)

)

* fvc :: snGrad ( alpha1 )

) * mesh . magSf () ;

solve

(

fvm :: ddt(C)

+ fvm :: div (phi , C, "div (phi ,C)")

- fvm :: laplacian ( fvc :: interpolate ( DT ), C , " laplacian (C)")

- fvm :: laplacian (( nut / Schmidtnumber ) , C, " laplacian (C)")

+ fvm :: div (phiCi , C , "div (phi ,C)")

mesh . solver ("C")

);

}

Listing A.3: Mass transfer equation in CEqn.H

129

Appendix A

A.2.5 interH2SFoam.C

For the mass transfer equation to be solved, it has to be included into the file interH2SFoam.C.

This is done by calling the include command. Furthermore, the turbulent viscosity is being

stored in a volScalarField so that it can be called when solving the mass transfer equation.

In Listing A.4, the content of interH2SFoam.C is displayed and the changes highlighted in

red.

---------------------------------------------------------------------------*\

========= |

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

\\ / O peration |

\\/ M anipulation |

-------------------------------------------------------------------------------

License

This file is part of OpenFOAM .

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

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

the Free Software Foundation , either version 3 of the License , or

(at your option ) any later version .

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

ANY WARRANTY ; without even the implied warranty of MERCHANTABILITY or

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

for more details .

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

along with OpenFOAM . If not , see < http :// www . gnu .org / licenses / >.

Application

interFoam

Description

Solver for 2 incompressible , isothermal immiscible fluids using a VOF

( volume of fluid ) phase - fraction based interface capturing approach .

The momentum and other fluid properties are of the " mixture " and a single

momentum equation is solved .

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

For a two - fluid approach see twoPhaseEulerFoam .

\*---------------------------------------------------------------------------

# include " fvCFD .H"

# include " CMULES .H"

130

Appendix A

# include " EulerDdtScheme .H"

# include " localEulerDdtScheme .H"

# include " CrankNicolsonDdtScheme .H"

# include " subCycle .H"

# include " immiscibleIncompressibleTwoPhaseMixture .H"

# include " turbulenceModel .H"

# include " pimpleControl .H"

# include " fvIOoptionList .H"

#include "fixedFluxPressureFvPatchScalarField.H"

# include " convectionScheme .H"

//*************************************

int main(int argc , char * argv [])

{

# include " setRootCase .H"

# include " createTime .H"

# include " createMesh .H"

pimpleControl pimple ( mesh );

# include " initContinuityErrs .H"

# include " createFields .H"

# include " readTimeControls .H"

# include " createPrghCorrTypes .H"

# include " correctPhi .H"

# include " CourantNo .H"

# include " setInitialDeltaT .H"

//***********************************

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

while ( runTime . run ())

{

# include " readTimeControls .H"

# include " CourantNo .H"

# include " alphaCourantNo .H"

# include " setDeltaT .H"

runTime ++;

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

// --- Pressure - velocity PIMPLE corrector loop

while ( pimple . loop ())

{

# include " alphaControls .H"

# include " alphaEqnSubCycle .H"

131

Appendix A

mixture.correct();

# include " UEqn .H"

// --- Pressure corrector loop

while ( pimple . correct ())

{

# include " pEqn .H"

}

if ( pimple . turbCorr ())

{

turbulence -> correct ();

}

// interH2SFoam : read nut for mass transfer equation and solve mass

transfer equation by including "CEqn .H"

volScalarField nut(" nut ", turbulence -> nut () );

# include " CEqn .H"

runTime . write () ;

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

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

<< nl << endl ;

}

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

return 0;

}

// *************************************************************************

Listing A.4: Extension of interFoam.C to interH2SFoam.C, see highlighted in red

132

Appendix A

A.2.6 Integration of interDyMH2SFoam

Similarly to the above mentioned steps, the interDyMFoam folder has to be modified to

account for mass transfer and be renamed to interDyMH2SFoam. As for the modification

of the interFoam solver, the directory has to be renamed from interDyMFoam to inter-

DyMH2SFoam. After resetting the solver using wclean, the file Make/files has to be modified

similarly to Listing A.1 but with the new solver name interDyMH2SFoam. Then, the inter-

DyMFoam.C has to be renamed to interH2SDyMFoam.C and Listing A.5 has to be written

into interH2SDyMFoam.C. Finally, the solver has to be recompiled using wmake to account

for the changes within the code.

volScalarField nut(" nut ", turbulence -> nut () );

# include " CEqn .H"

Listing A.5: Extension of interDyMH2SFoam.C

133

Appendix A

A.2.7 Passive scalar transport

For the passiveScalarInterFoam solver, the file CEqn.H contains a standard advection-diffusion

equation as displayed in Listing A.6. The other changes within the file structure are to be

made in a similar way as for the interH2SFoam solver.

fvScalarMatrix CEqn

(

fvm :: ddt(C)

+ fvm :: div (phi , C)

- fvm :: laplacian (D , C)

- fvm :: laplacian (( nut / Schmidtnumber ), C)

);

CEqn . solve ();

Listing A.6: CEqn.H for passiveScalarInterFoam

134

Appendix A

A.3 Function objects and boundary conditions

A.3.1 Equilibrium conditions

Equilibrium conditions between H2Sand HS−are computed using a custom boundary con-

dition. The boundary condition is implemented using the OpenFOAM extension swak4Foam

and is displayed in Listing A.7. It is defined as a boundary condition for the tracer concen-

tration in the file Cwithin the subfolder 0of the case directory.

bottom

{

type groovyBC ;

variables "pH =7.5; CHS =1.0; pKa =7.0; ";

valueExpression "( CHS * pow (10 , pKa - pH ) /(1+ pow (10 , pKa - pH )) )/32";

// placeholder

value 0.001;

}

Listing A.7: Boundary conditions for equilibrium conditions

135

Appendix A

A.3.2 Calculation of partial pressure in ppm

The extension is added to the model by adding a function object using the OpenFOAM

extension swak4Foam. The code is added to the controlDict file in the system folder of the

case directory. The relevant libraries have to be imported and the function object has to

be defined. The code displayed in Listing A.8 needs to be added to the bottom of the

controlDict file.

libs (

" libOpenFOAM .so"

" libswakFunctionObjects .so"

"libsimpleFunctionObjects.so"

);

functions

{

conc_ppm

{

type expressionField ;

fieldName ppm ;

valueType internalField ;

aliases {

alphawater alpha . water ;

}

expression "(1 - alphawater ) * 1 e6 * C /1000 * 27 * 0.008206 ";

autowrite true;

verbose true ;

}

Listing A.8: Function object for calculation of partial pressure in ppm in controlDict

136

Appendix B

Test case overview

137

Appendix B

B.1 Single-phase flow over a ground sill

Table B.1: Model setup of ERCOFTAC test case

Title: Single-phase flow over a ground sill

General

Referred to in chapters 2.4.2, 3.3.1, 4.3

References Almeida et al. (1993), Davroux et al. (1995)

Published in Teuber et al. (2017), Teuber et al. (2016), Teuber et al. (2019a)

Domain discretization

Dimensions two-dimensional (length: 0.85m, height (max): 0.17m)

Mesh generator GMSH

Number of cells 31,102

Turbulence models Standard k-e, Standard k-ω, k-ωSST, LES

Hydrodynamic simulations

Solver interFoam

Time step variable, converged against 0.005 s RANS simulations and against 0.001 s for LES

Simulation time 10s

Single-phase transport simulations

None performed

Mass transfer simulations

None performed

138

Appendix B

Table B.2: Boundary conditions for hydrodynamic simulations of single-phase flow over a ground sill (Turbulence properties calculated

using https://www.cfd-online.com/Tools/turbulence.php)

alpha.water [-] prgh [ kg

m·s2] U [m/s] k [m2

s2]e[m2

s3] / ω[1

inlet inletOutlet

inletValue: uniform 1

value: uniform 1

zeroGradient fixedValue

value: parabolic

function following

experimental results,

applied using

funkySetBoundary

turbulentIntensity-

KineticEnergyInlet

intensity: 0.03

value: uniform 0.00622

turbulentMixingLength-

DissipationRateInlet

mixingLength: 0.0119

value: uniform

0.00371

ω:

fixedValue

value: $internalField

outlet zeroGradient fixedValue

value: uniform 1667.7

zeroGradient zeroGradient zeroGradient

upper and

lower wall

zeroGradient zeroGradient fixedValue

value: uniform (0 0 0)

kqRWallFunction

value: uniform 0

epsilonWallFunction

value: uniform 0

ω:

omegaWallFunction

value: $internalField

sidewalls empty empty empty empty empty

initial condi-

tions

uniform 1 uniform 0 uniform (0 0 0) uniform 0.00622 e:

uniform 0.00371

ω:

uniform 6.62906

139

Appendix B

B.2 Two-phase flow over a two-dimensional hill

Table B.3: Model setup, two-phase flow over a two-dimensional hill

Title: Two-phase flow over a two-dimensional hill

General

Referred to in chapters 2.4.2, 3.3.2, 4.4

References -

Published in Teuber et al. (2017), Teuber et al. (2016), Teuber et al. (2019a)

Domain discretization

Dimensions two-dimensional:

cases 1,2: length: 25m, height (max): 2m

case 3: length: 35m, height (max): 6m

three-dimensional:

extension of case 1: length: 25m, height (max): 2m, width: 1m

Mesh generator GMSH

Number of cells 68,542 (cases 1 and 2, comparable cell numbers for ground sill variations referred to in

chapter 4.4), 175,762 (case 3), 685,420 (3D case based on case 1)

Turbulence models Standard k-e

Hydrodynamic simulations

Solver interFoam

Time step variable, converged against 0.001 s

Simulation time 100s

Single-phase transport simulations

None performed

Mass transfer simulations

None performed

140

Appendix B

Table B.4: Boundary conditions for hydrodynamic simulations of two-phase flow over a hill (case 1, 2D setup, turbulence properties were

calculated using https://www.cfd-online.com/Tools/turbulence.php)

alpha.water [-] prgh [ kg

m·s2] U [m/s] k [m2

s2]e[m2

s3]

inlet_air inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

U: U

phi: phi

rho: rho

psi: none

gamma: 1

value: uniform 0

zeroGradient fixedValue

value: $internalField

fixedValue

value: $internalField

inlet_water inletOutlet

inletValue: uniform 1

value: uniform 1

zeroGradient fixedValue

value: uniform (1 0 0)

(for other cases

respective velocities)

fixedValue

value: $internalField

fixedValue

value: $internalField

outlet zeroGradient fixedValue

value: uniform 0

inletOutlet

inletValue:

uniform (0 0 0)

value: uniform (0 0 0)

zeroGradient zeroGradient

upper and

lower wall

zeroGradient zeroGradient fixedValue

value: uniform (0 0 0)

zeroGradient zeroGradient

sidewalls empty empty empty empty empty

initial condi-

tions

uniform 0

define initial water

level of 1m using

setFields

(volScalarFieldValue

alpha.water 1)

uniform 0 uniform (0 0 0)

define initial flow

velocity of 1 using

setFields

(volVectorFieldValue

U (1 0 0))

uniform 0.00375 uniform 0.00054

141

Appendix B

B.3 Complex sewer geometry

Table B.5: Model setup of complex sewer geometry

Title: Complex sewer geometry

General

Referred to in chapters hydrodynamic simulations: 2.4.2, 3.3.3

single-phase transport: 5.3.2

mass transfer: 6.3.4, 6.4.3

References Bayón et al. (2015)

Published in hydrodynamic simulations: Teuber et al. (2017), Teuber et al. (2016), Teuber et al. (2019a)

single-phase transport: Teuber et al. (2019c)

mass transfer: Teuber et al. (2019b)

Domain discretization

Dimensions three-dimensional (length: approx. 95 m, height (max): 7.50m)

Mesh generator snappyHexMesh

Number of cells 3,029,223

Turbulence models Standard k-e

Hydrodynamic simulations

Solver interFoam

Time step variable, converged against 0.00019 s

Simulation time 200s

Single-phase transport simulations

Solvers interH2SFoam with modified Henry coefficient

Simulation time 10s, starting from quasi steady-state reached by hydrodynamic simulations

Tracer diffusivity water: 10−5m2/sair: 10−5m2/s

Schmidt number 1.0

Henry coefficient 10−6

Mass transfer simulations

Solver interH2SFoam

Simulation time 10s, starting from quasi steady-state reached by hydrodynamic simulations

Tracer diffusivity water: 10−3m2/sair: 10−3m2/s

Schmidt number 1.0

Temperature 298 K

142

Appendix B

Table B.6: Boundary conditions for simulations of complex sewer geometry (turbulence properties were calculated using https://www.cfd-

online.com/Tools/turbulence.php)

alpha.water [-] prgh [ kg

m·s2] U [m/s] k [m2

s2]e[m2

s3] C [mol

m3] (only applies

to single-phase transport

and mass transfer cases)

inlet_air inletOutlet

inletValue:

uniform 0

value: uniform 0

totalPressure

p0: uniform 0

U: U

phi: phi

rho: rho

psi: none

gamma: 1

value: uniform 0

pressureInletOutlet-

Velocity

phi: phi

tangentialVelocity:

uniform (0 0 0)

value: uniform (0 0 0)

zeroGradient zeroGradient zeroGradient

inlet_water inletOutlet

inletValue:

uniform 1

value: uniform 1

fixedFluxPressure flowRateInletVelocity

volumetricFlowRate:

100

fixedValue

value: $internalField

turbulentMixingLength-

DissipationRateInlet

mixingLength: 0.005

value: uniform 0.0375

zeroGradient

outlet zeroGradient fixedValue

value: uniform 0

inletOutlet

inletValue:

uniform (0 0 0)

value: uniform (0 0 0)

inletOutlet

inletValue:

$internalField

value: $internalField

inletOutlet

inletValue:

$internalField

value: $internalField

zeroGradient

bottom zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0 0)

kqRWallFunction

value: uniform 0.063;

epsilonWallFunction

value: uniform 0.0375

single-phase transport:

zeroGradient

mass transfer:

fixedValue

value: 1

top wall and

side walls

zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0 0)

kqRWallFunction

value: uniform 0.063;

epsilonWallFunction

value: uniform 0.0375

zeroGradient

initial condi-

tions

uniform 0 uniform 0 uniform (0 0 0) uniform 0.063 uniform 0.0375 uniform 0, and for

- single-phase transport:

use setFields to define

box where C=1

- mass transfer: use

expression in funkySet-

Fields to define C=1

where alpha.water >0.5

143

Appendix B

B.4 Single-phase transport in rectangular pipe

Table B.7: Model setup of single-phase transport in rectangular pipe test case

Title: Single-phase transport in rectangular pipe

General

Case description Tracer transport along rectangular channel

case 1: interFoam extended by standard advection-diffusion equation

case 2: interH2SFoam with modified Henry coefficient

Referred to in chapters 5.1, 5.3.2

References -

Published in Teuber et al. (2019c)

Domain discretization

Dimensions three-dimensional (length: 15m, height (max): 1m, depth: 1m)

Mesh generator blockMesh

Number of cells 45,000

Turbulence model Standard k-e

Hydrodynamic simulations

None performed

Single-phase transport simulations

Solvers case 1: interFoam extended by advection-diffusion equation

case 2: interH2SFoam with modified Henry coefficient

Simulation time 150s

Tracer diffusivity water: 10−3m2/sair: 10−3m2/s(high diffusivity in order to force fast spreading during

transport along the channel)

Schmidt number 1.0

Mass transfer simulations

None performed

144

Appendix B

Table B.8: Boundary conditions for simulations of single-phase transport in rectangular pipe (turbulence properties were calculated using

https://www.cfd-online.com/Tools/turbulence.php)

alpha.water [-] prgh [ kg

m·s2] U [m/s] k [m2

s2]e[m2

s3] C [mol

m3]

inlet_air inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

U: U

phi: phi

rho: rho

psi: none

gamma: 1

value: uniform 0

zeroGradient zeroGradient zeroGradient zeroGradient

inlet_water inletOutlet

inletValue: uniform 1

value: uniform 1

fixedFluxPressure flowRateInletVelocity

volumetricFlowRate:

0.05

turbulentIntensity-

KineticEnergyInlet

intensity: 0.05

value: uniform

0.00004

turbulentMixing-

LengthDissipation-

RateInlet

mixingLength: 0.038

value: uniform 5e-7

zeroGradient

outlet zeroGradient zeroGradient outletPhaseMean-

Velocity

Umean: 0.1

alpha: alpha.water

value: uniform (0.1 0 0)

zeroGradient inletOutlet

inletValue:

$internalField

value: $internalField

zeroGradient

bottom zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0 0)

kqRWallFunction

value: uniform

0.00004;

epsilonWallFunction

value: uniform 5e-7

zeroGradient

top wall and

side walls

zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0 0)

kqRWallFunction

value: uniform

0.00004;

epsilonWallFunction

value: uniform 5e-7

zeroGradient

initial condi-

tions

uniform 0

use setFields to define

water phase

(alpha.water =1)

within bounding box

between (0 0 -1) and

(15 0.5 2)

uniform 0 uniform (0.1 0 0) uniform 0.00004 uniform 5e-7 uniform 0

use setFields to de-

fine initial tracer con-

centration (C=1)

within bounding box

between (1 0.1 0.4) (1.5

0.4 0.6)

145

Appendix B

B.5 Concrete probes in sewer pilot plant

Table B.9: Concrete probes in sewer pilot plant

Title: Concrete probes in sewer pilot plant

General

Case description Single-phase tracer transport around concrete probes in sewer pilot plant

Referred to in chapters 5.3.1

References -

Published in Teuber et al. (2019c)

Domain discretization

Dimensions three-dimensional (length: 2m, height (max): 0.3m, depth: 0.5m)

Mesh generator SnappyHexMesh

Number of cells 1,041,772

Turbulence model laminar

Hydrodynamic simulations

Solver interFoam

Simulation time 50s

Time step variable, converged against 0.024 s

Single-phase transport simulations

Solvers interFoam extended by advection-diffusion equation

Simulation time 65s, starting from quasi steady-state reached by hydrodynamic simulations

Tracer diffusivity air: 10−9m2/s

Schmidt number 1.0

Mass transfer simulations

None performed

146

Appendix B

Table B.10: Boundary conditions for transport simulations around concrete probes (turbulence properties were calculated using

https://www.cfd-online.com/Tools/turbulence.php)

alpha.water [-] prgh [ kg

m·s2] U [m/s] C [mol

m3]

inlet inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

U: U

phi: phi

rho: rho

psi: none

gamma: 1

value: uniform 0

pressureInlet-

OutletVelocity

phi: phi

tangentialVelocity:

uniform (0 0 0)

value: uniform (0 0 0)

zeroGradient

outlet zeroGradient fixedValue

value: uniform -0.02

inletOutlet

inletValue: uniform (0

0 0)

value: uniform (0 0 0)

zeroGradient

watersurface zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0.5

zeroGradient

pipe zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0 0)

zeroGradient

concrete

probes

zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0 0)

zeroGradient

initial condi-

tions

uniform 0 uniform 0 uniform (0 0 0) uniform 0

use setFields to de-

fine initial tracer con-

centration (C=1)

within bounding box

between (0 1.379 0)

(100 2.12 100)

147

Appendix B

B.6 Quasi one-dimensional cubic tank

Table B.11: Model setup of quasi one-dimensional cubic tank

Title: Quasi one-dimensional cubic tank

General

Case description Validation of mass transfer solver extensions using a cubic tank setup and analysing the

equilibrium conditions at steady state

case 1: analysis of concentration profile at normal temperature (298.15 K)

case 2: analysis of concentration profile at different temperature (293.15 K)

case 3: analysis of computed equilibrium conditions between HS−and H2Susing new

boundary condition and analysis of calculated partial pressure of H2Sin air phase

Referred to in chapters 6.3.4, 6.4.1

References -

Published in Teuber et al. (2019b)

Domain discretization

Dimensions quasi one-dimensional (height: 1m, width: 1m, depth: 1m)

Mesh generator blockMesh

Number of cells 10,000

Turbulence models laminar

Hydrodynamic simulations

None performed (no fluid flow in domain)

Single-phase transport simulations

None performed

Mass transfer simulations

Solver interH2SFoam

Simulation time 1000s

Time step 0.1 s

Tracer diffusivity water: 10−2m2/sair: 10−2m2/s(a high diffusivity was chosen because only the equilib-

rium conditions established were of interest)

Temperature 298.15 K and 293.15 K

148

Appendix B

Table B.12: Boundary conditions for simulations of quasi one-dimensional cubic tank

alpha.water [-] prgh [ kg

m·s2] U [m/s] C [mol

m3]

bottom zeroGradient zeroGradient fixedValue

value: uniform (0 0 0)

- cases 1 and 2:

fixedValue

value: uniform 1

- case 3:

groovyBC according

to Listing A.7

top wall and

side walls

zeroGradient zeroGradient fixedValue

value: uniform (0 0 0)

zeroGradient

initial condi-

tions

use setFields to define

a water level of 0.5 m

(alpha.water =1 for

h<0.5m,

alpha.water =0

elsewhere)

uniform 0 uniform (0 0 0) use setFields to define

initial concentration

in water phase (C=1

for h<0.5m,C=0

elsewhere)

149

Appendix B

B.7 Mass transfer in rectangular duct

Table B.13: Model setup of mass transfer in rectangular duct

Title: Mass transfer in rectangular duct

General

Case description Mass transfer simulations according to test cases 7 (here: case 1) and 21 (here: case

2) as carried out by Bentzen et al. (2016). Hydrodynamic simulations were carried out

until a quasi steady-state has been reached, then the mass transfer simulations have been

performed.

Referred to in chapters 6.3.4, 6.4.2

References Bentzen et al. (2016)

Published in Teuber et al. (2019b)

Domain discretization

Dimensions three-dimensional (length: 15m, height: 0.26m, width: 0.3m)

Mesh generator GMSH

Number of cells 307,970

Turbulence model laminar

Hydrodynamic simulations

Solver interFoam

Time step variable, converged against 0.007 s (case 1), 0.019 (case 2)

Simulation time 200s

Single-phase transport simulations

None performed

Mass transfer simulations

Solver interH2SFoam

Simulation time 50s, starting from quasi steady-state reached by hydrodynamic simulations

Tracer diffusivity water: 2.2 ·10−9m2/sair: 1.74 ·10−5m2/s

Schmidt number 1.0

Temperature 298 K

150

Appendix B

Table B.14: Boundary conditions for mass transfer in rectangular duct

alpha.water [-] prgh [ kg

m·s2] U [m/s] C [mol

m3] (only applies

to mass transfer cases)

inlet_air inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

U: U

phi: phi

rho: rho

psi: none

gamma: 1

value: uniform 0

pressureInletOutlet-

Velocity

phi: phi

tangentialVelocity:

uniform (0 0 0)

value: uniform (0 0 0)

zeroGradient

inlet_water inletOutlet

inletValue: uniform 1

value: uniform 1

fixedFluxPressure flowRateInletVelocity

volumetricFlowRate:

0.0072 (case 1, 0.0164

for case 2)

fixedValue

value: uniform 1

outlet zeroGradient fixedValue

value: uniform -0.025

inletOutlet

inletValue:

uniform (0 0 0)

value: uniform (0 0 0)

inletOutlet

inletValue:

uniform (0 0 0)

value: uniform (0 0 0)

bottom zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0 0)

zeroGradient

top wall and

side walls

zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0 0)

zeroGradient

initial condi-

tions

uniform 0 uniform 0 uniform (0 0 0) uniform 0, use expres-

sion in funkySetFields

to define C=1 where

alpha.water >0.5

151

Appendix B

B.8 Stirring tank

Table B.15: Model setup of stirring tank

Title: Mass transfer in stirring tank

General

Case description Mass transfer simulations according to experiments as carried out by Carrera et al. (2017).

Hydrodynamic simulations were carried out with different stirring rates (case 1: 50 rpm,

case 2: 80 rpm, case 3: 100 rpm, case 4: 120 rpm, case 5: 140 rpm) until a quasi steady-state

has been reached, then the mass transfer simulations have been performed.

Referred to in chapter 7

References Carrera et al. (2017)

Published in Teuber et al. (in preparation)

Domain discretization

Dimensions three-dimensional (tank diameter: 0.25m, blade length: 0.05m, blade height: 0.03m)

Mesh generator SnappyHexMesh

Number of cells 214,100

Turbulence model RNG k-e

Hydrodynamic simulations

Solver interFoam

Time step variable, converged against 0.007 s (case 1), 0.019 (case 2)

Simulation time 25s

Single-phase transport simulations

None performed

Mass transfer simulations

Solver interH2SFoam

Simulation time 10s, starting from quasi steady-state reached by hydrodynamic simulations

Tracer diffusivity water: 1.4 ·10−9m2/sair: 1.5 ·10−5m2/s

Schmidt number 1.0

Temperature 298 K

152

Appendix B

Table B.16: Boundary conditions for stirring tank

alpha.water [-] prgh [ kg

m·s2] U [m/s] k [m2

s2]e[m2

s3] C [mol

m3] (only applies

to mass transfer cases)

top inletOutlet

inletValue: uniform 0

value: uniform 0

totalPressure

p0: uniform 0

U: U

phi: phi

rho: rho

psi: none

gamma: 1

value: uniform 0

pressureInletOutlet-

Velocity

phi: phi

tangentialVelocity:

uniform (0 0 0)

value: uniform (0 0 0)

inletOutlet

inletValue: uniform

1e-5

value: uniform 1e-5

inletOutlet

inletValue:

$internalField

value: $internalField

inletOutlet

inletValue: 0

value: 0

static walls zeroGradient fixedFluxPressure fixedValue

value: uniform (0 0 0)

kqRWallFunction

value: uniform 1e-5

epsilonWallFunction

value: uniform 1e-5

zeroGradient

rotating wall

elements

zeroGradient fixedFluxPressure movingWallVelocity

value: $internalField;

kqRWallFunction

value: uniform 1e-5

epsilonWallFunction

value: uniform 1e-5

zeroGradient

AMI cyclicAMI cyclicAMI cyclicAMI cyclicAMI cyclicAMI cyclicAMI

initial condi-

tions

use setFields to define

a water level of 0.5 m

(alpha.water =1 for

h<0.5m,

alpha.water =0

elsewhere)

uniform 0 uniform (0 0 0) uniform 1e-5 uniform 1e-5 uniform 0, use expres-

sion in funkySetFields

to define C=0.073

where alpha.water >

0.5

153

Bibliography

Almeida, G., Durao, D., and Heitor, M.: Wake flows behind two-dimensional model hills,

Experimental Thermal and Fluid Science, 7(1), 87–101, 1993.

Ashley, R.: Solids in Sewers: Characteristics, Effects and Control of Sewer Solids and Asso-

ciated Pollutants (Scientific and Technical Report), IWA Publishing, 2004.

Auguet, O., Pijuan, M., Guasch-Balcells, H., Borrego, C., and Gutierrez, O.: Implications of

downstream nitrate dosage in anaerobic sewers to control sulfide and methane emissions,

Water research, 68, 522–532, 2015.

Barjenbruch, M.: Prevention of odour emergence in sewage networks, Water science and

technology, 47(7-8), 357–363, 2003.

Barjenbruch, M., Hinkelmann, R., Huhnt, W., Krämer, T., Nehrig, M., and Rühmland, S.:

An Online-Monitoring and Operating System to Prevent Odour and Corrosion in Sewer

Networks, Feasibility Study, 2008.

Barlag, C., Hinkelmann, R., Helmig, R., and Zielke, W.: Adaptive methods for modelling

transport processes in fractured subsurface systems, in: 3rd International Conference on

Hydroscience and Engineering, Cottbus, Center of Computational Hydroscience and En-

gineering, The University of Mississippi, vol. 284, 1998.

Bates, P., Lane, S., and Ferguson, R.: Computational Fluid Dynamics: Applications in Envi-

ronmental Hydraulics, John Wiley & Sons, Ltd., Chichester, UK, 2005.

Bayón, A. and Lopez-Jimenez, P.: Numerical analysis of hydraulic jumps using OpenFOAM,

Journal of Hydroinformatics, 17(4), 662–678, 2015.

Bayón, A., Vallés-Morán, F., and López-Jiménez, P.: Numerical analysis and validation of

South Valencia sewage collection system diversion, in: 36th IAHR World Congress, The

Hague, The Netherlands, 2015.

Bentzen, T., Østertoft, K., Vollertsen, J., Fuglsang, E., and Nielsen, A.: Airflow in Grav-

ity Sewers–Determination of Wastewater Drag Coefficient, Water Environment Research,

88(3), 239–256, 2016.

Berberovi´c, E., van Hinsberg, N., Jakirli´c, S., Roisman, I., and Tropea, C.: Drop impact onto

a liquid layer of finite thickness: Dynamics of the cavity evolution, Physical Review E,

79(3), 036306, 2009.

Bjerre, H. L., Hvitved-Jacobsen, T., Teichgräber, B., and Schlegel, S.: Modeling of aerobic

wastewater transformations under sewer conditions in the Emscher River, Germany, Wa-

ter Environment Research, 70(6), 1151–1160, 1998.

154

Bibliography

Boon, A. and Lister, A.: Formation of sulphide in rising main sewers and its prevention by

injection of oxygen, Prog. Wat. Tech., 1975.

Caretto, L., Gosman, A., Patankar, S., and Spalding, D.: Two calculation procedures for

steady, three-dimensional flows with recirculation, in: Proceedings of the third interna-

tional conference on numerical methods in fluid mechanics, pp. 60–68, Springer, 1973.

Carrera, L., Springer, F., Lipeme-Kouyi, G., and Buffiere, P.: A review of sulfide emissions in

sewer networks: overall approach and systemic modelling, Water Science and Technology,

73(6), 1231–1242, 2016.

Carrera, L., Springer, F., Lipeme-Kouyi, G., and Buffiere, P.: Sulfide emissions in sewer

networks: focus on liquid to gas mass transfer coefficient, Water Science and Technology,

75(8), 1899–1908, 2017.

Celik, I., Ghia, U., and Roache, P.: Procedure for estimation and reporting of uncertainty

due to discretization in CFD applications, Journal of Fluids, 130(7), 078 001, 2008.

Cifani, P., Michalek, W., Priems, G., Kuerten, J. G., van der Geld, C., and Geurts, B. J.:

A comparison between the surface compression method and an interface reconstruction

method for the VOF approach, Computers & Fluids, 136, 421–435, 2016.

Damián, S.: Description and utilization of interFoam multiphase solver, URL:

http://infofich. unl. edu. ar/upload/3be0e16065026527477b4b948c4caa7523c8ea52. pdf,

2012.

Damián, S.: An extended mixture model for the simultaneous treatment of short and long

scale interfaces, Ph.D. thesis, Universidad Nacional Del Litoral, Argentina, 2013.

Davroux, A., Hoa, C., and Laurence, D.: Flow over a 2D hill—reference solutions for k–and

second moment closure turbulence models, in: ERCOFTAC Workshop, Karlsruhe, 1995.

De Belie, N., Monteny, J., and Taerwe, L.: Apparatus for accelerated degradation testing of

concrete specimens, Materials and Structures, 35(7), 427–433, 2002.

Devolder, B., Schmitt, P., Rauwoens, P., Elsäßer, B., and Troch, P.: A Review of the Implicit

Motion Solver Algorithm in OpenFOAM to Simulate a Heaving Buoy, 2015.

DIN 1045-2:2008-08: Tragwerke aus Beton, Stahlbeton und Spannbeton–Teil 2: Beton–

Festlegung, Eigenschaften, Herstellung und Konformität–Anwendungsregeln zu DIN EN

206-1, Beuth Verlag Berlin, 2008.

DIN E 197-1: 2011-11: Zement – Teil 1: Zusammensetzung, Anforderungen und Kon-

formitätskriterien von Normalzement; Deutsche Fassung EN 197-1:2011., Beuth Verlag

Berlin, 2011.

DIN EN 196-1:2005-05: Prüfverfahren für Zement – Teil 1: Bestimmung der Festigkeit;

Deutsche Fassung EN 196-1:2005, Beuth Verlag Berlin, 2005.

E DIN 19573:2013-06: Mörtel für Neubau und Sanierung von Entwässerungssystemen

außerhalb von Gebäuden, Beuth Verlag Berlin, 2013.

Edwini-Bonsu, S. and Steffler, P.: Air flow in sanitary sewer conduits due to wastewater

drag: a computational fluid dynamics approach, Journal of Environmental Engineering

and Science, 3(5), 331–342, 2004.

155

Bibliography Appendix B

Edwini-Bonsu, S. and Steffler, P.: Dynamics of air flow in sewer conduit headspace, Journal

of Hydraulic Engineering, 132(8), 791–799, 2006.

Fourie, C. and Alexander, M.: A dynamic hydrochloric acid test for evaluating acid resis-

tance of sewer pipe concrete, in: Workshop on Performance of Cement-based Materials in

Aggressive Aqueous Environments-Characterisation, Modelling, Test Methods and Engi-

neering Aspects, pp. 104–113, RILEM Publications SARL, 2008.

Freudenthal, K., Koglatis, J., Otterpohl, R., and Behrendt, J.: Prediction of sulfide formation

in sewer pressure mains based on the IWA Anaerobic Digestion Model No. 1 (ADM1),

Water Science and Technology, 52(10-11), 13–22, 2005.

Ganigué, R., Jiang, G., Sharma, K., Chen, J., Vuong, L., and Yuan, Z.: Online Control of

magnesium hydroxide dosing for sulfide mitigation in sewers: algorithm development,

simulation analysis, and field validation, Journal of Environmental Engineering, 142(12),

04016069, 2016.

Gessner, M., Hinkelmann, R., Nützmann, G., Jekel, M., Singer, G., Lewandowski, J., Nehls,

T., and Barjenbruch, M.: Urban water interfaces, Journal of Hydrology, 514, 226–232, 2014.

Gilchrist, F.: Microbiological studies of the corrosion of concrete sewers by sulphuric acid

producing bacteria, South African Industrial Chemist, 7, 214–215, 1953.

Greenshields, C. J.: OpenFOAM user guide, OpenFOAM Foundation Ltd, version 2.4.0, 3(1),

47, 2015.

Grengg, C., Mittermayr, F., Baldermann, A., Böttcher, M., Leis, A., Koraimann, G., Grunert,

P., and Dietzel, M.: Microbiologically induced concrete corrosion: A case study from a

combined sewer network, Cement and Concrete Research, 77, 16–25, 2015.

Haroun, Y., Legendre, D., and Raynal, L.: Volume of fluid method for interfacial reactive

mass transfer: application to stable liquid film, Chemical Engineering Science, 65(10),

2896–2909, 2010a.

Haroun, Y., Legendre, D., and Raynal, L.: Direct numerical simulation of reactive absorption

in gas–liquid flow on structured packing using interface capturing method, Chemical

Engineering Science, 65(1), 351–356, 2010b.

Higuera, P., Lara, J., and Losada, I.: Simulating coastal engineering processes with Open-

FOAM®, Coastal Engineering, 71, 119–134, 2013.

Hinkelmann, R.: Efficient Numerical Methods and Information-Processing Techniques for

Modeling Hydro- and Environmental Systems, Lecture Notes in Applied and Computa-

tional Mechanics, Springer Berlin Heidelberg, 2005.

Hoang, D., van Steijn, V., Portela, L., Kreutzer, M., and Kleijn, C.: Benchmark numerical

simulations of segmented two-phase flows in microchannels using the Volume of Fluid

method, Computers & Fluids, 86, 28–36, 2013.

Hüttl, R.: Flugaschebetone mit hohem Säurewiderstand im Kühlturm-und Wasserbau–

Historie, Stand der Technik und Bewertungskriterien, VGB Power-Tech, 4, 2008.

HV-Consult: Description of the WATS model, URL http://hvcon.com/Download/

Description%20of%20the%20WATS%20model.pdf, 2019.

156

Bibliography

Hvitved-Jacobsen, T., Jütte, B., Nielsen, P., and Jensen, N.: Hydrogen sulphide control in

municipal sewers, in: Pretreatment in chemical water and wastewater treatment, pp. 239–

247, Springer Berlin Heidelberg, 1988.

Hvitved-Jacobsen, T., Raunkjær, K., and Nielsen, P. H.: Volatile fatty acids and sulfide in

pressure mains, Water Science and Technology, 31(7), 169–179, 1995.

Hvitved-Jacobsen, T., Vollertsen, J., and Nielsen, A.: Sewer processes: microbial and chemi-

cal process engineering of sewer networks, CRC press, Boca Raton, USA, 2013.

Issa, R.: Solution of the implicitly discretised fluid flow equations by operator-splitting,

Journal of computational physics, 62(1), 40–65, 1986.

Jasak, H.: Error analysis and estimation for the finite volume method with applications to

fluid flows., Ph.D. thesis, Imperial College London (University of London), 1996.

Ji, Z.: Water Quality and Eutrophication. Hydrodynamic and Water Quality: Modeling

Rivers, Lakes, and Estuaries, John Wiley & Sons, Inc., Hoboken, New Jersey, 2008.

Jiang, G., Sun, J., Sharma, K., and Yuan, Z.: Corrosion and odor management in sewer

systems, Current Opinion in Biotechnology, 33, 192–197, 2015.

Kaempfer, W. and Berndt, M.: Estimation of service life of concrete pipes in sewer networks,

Durability of Building Materials and Components, 8, 36–45, 1999.

Keough, S.: Optimising the Parallelisation of OpenFOAM Simulations, Tech. rep., Defence

Science and Technology Organisation Fishermans Bend Australia Maritime Div., 2014.

Kinyua, M., Zhang, J., Camacho-Céspedes, F., Tejada-Martinez, A., and Ergas, S.: Use of

physical and biological process models to understand the performance of tubular anaer-

obic digesters, Biochemical Engineering Journal, pp. 35–44, 2016.

Kortelainen, J.: Meshing tools for open source CFD-a practical point of view, VTT, Espoo,

Finland, Tech. Rep, 2009.

Launder, B. and Sharma, B.: Application of the energy-dissipation model of turbulence

to the calculation of flow near a spinning disc, Letters in Heat and Mass Transfer, 1(2),

131–137, 1974.

Leigh, N. G. and Lee, H.: Sustainable and Resilient Urban Water Systems: The Role of

Decentralization and Planning, Sustainability, 11(3), 918, 2019.

Lopes, P., Leandro, J., and Carvalho, R. F.: Self-Aeration Modelling Using a Sub-Grid

Volume-Of-Fluid Model, International Journal of Nonlinear Sciences and Numerical Sim-

ulation, 18, 559–574, 2017.

Maric, T., Hopken, J., and Mooney, K.: The OpenFOAM technology primer, Sourceflux,

2014.

Matias, N., Nielsen, A., Vollertsen, J., Ferreira, F., and Matos, J.: Reaeration and hydrogen

sulfide release at drop structures, in: 8th International Conference on Sewer Processes

and Networks, Rotterdam, The Netherlands, 2016.

Matias, N., Nielsen, A., Vollertsen, J., Ferreira, F., and Matos, J.: Erratum: Water Science

and Technology 75 (10), 2257–2267: Liquid-gas mass transfer at drop structures, Natércia

Matias, Asbjørn Haaning Nielsen, Jes Vollertsen, Filipa Ferreira and José Saldanha Matos,

Water Science and Technology, 76(6), 1584–1594, 2017.

157

Bibliography Appendix B

Matsché, N., Saracevic, E., Bertrán de Lis, F., and Brooks, L.: Korrosions-und Geruchsprob-

leme in Abwasserdruckleitungen, Institut für Wassergüte, Ressourcenmanagement und

Abfallwirtschaft der Technischen Universität Wien, 2005.

Matta, E.: Multi-dimensional flow and transport modeling of a surface water body in a

semi-arid area, Ph.D. thesis, Technische Universität Berlin, Germany, 2018.

Menter, F.: Zonal two equation k-ωturbulence models for aerodynamic flows, in: 23rd Fluid

Dynamics, Plasmadynamics, and Lasers Conference, p. 2906, 1993.

Menter, F.: Two-equation eddy-viscosity turbulence models for engineering applications,

AIAA Journal, 32(8), 1598–1605, 1994.

Morgan, G.: Application of the interFoam VOF code to coastal wave/structure interaction,

Ph.D. thesis, University of Bath, 2013.

Nielsen, P., Raunkjær, K., and Hvitved-Jacobsen, T.: Sulfide production and wastewater

quality in pressure mains, Water Science and Technology, 37(1), 97–104, 1998.

Nieves-Remacha, M., Yang, L., and Jensen, K.: OpenFOAM computational fluid dynamic

simulations of two-phase flow and mass transfer in an advanced-flow reactor, Industrial

& Engineering Chemistry Research, 54(26), 6649–6659, 2015.

Nilsson, H.: How to implement your own application, URL http://www.tfd.chalmers.se/

~hani/kurser/OS_CFD_2011/implementApplication.pdf, 2019.

ÖNORM B 4710-1: Beton–Teil 1: Festlegung, Herstellung, Verwendung und Konformität-

snachweis (Regeln zur Umsetzung der ÖNORM EN 206-1 für Normal-und Schwerbeton),

Wien: Österreichisches Normungsinstitut, 2007.

Parker, C.: The corrosion of concrete: 1. The isolation of a species of bacterium associated

with the corrosion of concrete exposed to atmospheres containing hydrogen sulphide.,

Australian Journal of Experimental Biology and Medical Science, 23(2), 14–17, 1945.

Patankar, S. and Spalding, D.: A calculation procedure for heat, mass and momentum

transfer in three-dimensional parabolic flows, in: Numerical Prediction of Flow, Heat

Transfer, Turbulence and Combustion, pp. 54–73, Elsevier, 1983.

Petersen, L. and Lohaus, L.: Entwicklung eines Prüfstandes für Parameterstudien

zur Säurebeständigkeit von Hochleistungsbetonen, Proceedings of 16th International

Baustofftagung (Ibausil), pp. 0637–0644, 2006.

Pomeroy, R. and Parkhurst, J.: The forecasting of sulfide build-up rates in sewers, in: Eighth

International Conference on Water Pollution Research, pp. 621–628, Elsevier, 1978.

Rootsey, R. and Yuan, Z.: New insights into sewer odour and corrosion, in: 6th International

Conference on Sewer Processes and Networks, Queensland, Australia, 2010.

Rootsey, R., Melchers, R., Stuetz, R., Keller, J., and Yuan, Z.: Taking control of odours

and corrosion in sewers, in: Proceedings of Australia’s National Water Conference and

Exhibition (OzWater 2012), Sydney, Australia, pp. 8–10, Citeseer, 2012.

Rusche, H.: Computational fluid dynamics of dispersed two-phase flows at high phase

fractions, Ph.D. thesis, Imperial College London (University of London), 2003.

158

Bibliography

Rutz, M.: Laboratory investigation of physical disintegration of gross solids in the sewer,

MSc. Thesis. ETH Zürich, 2016.

Sander, R.: Compilation of Henry’s law constants (version 4.0) for water as solvent., Atmo-

spheric Chemistry & Physics, 15(8), 2015.

Schankat, M.: DiaTrans: A Multi-component Model for Density-driven Flow, Transport and

Biogeochemical Reaction Processes in the Subsurface, Citeseer, 2010.

Schulze, L. and Thorenz, C.: The multiphase capabilities of the CFD toolbox OpenFOAM for

hydraulic engineering applications, in: ICHE 2014. Proceedings of the 11th International

Conference on Hydroscience & Engineering, pp. 1007–1016, 2014.

Schulze, L. and Thorenz, C.: Mehrphasenmodellierung im Wasserbau, in: Wasserbauwerke-

Vom hydraulischen Entwurf bis zum Betrieb, pp. 53–58, Bundesanstalt für Wasserbau,

2015.

Severin, T.: Computational Fluid Dynamics Assisted Design of Thin-Layer Cascade Photo-

bioreactor Components, Ph.D. thesis, Technische Universität München, 2017.

Sharma, K., De Haas, D., Corrie, S., O’Halloran, K., Keller, J., and Yuan, Z.: Predicting

hydrogen sulfide formation in sewers: a new model, Journal of the Australian Water

Association, 35(2), 132–137, 2008a.

Sharma, K., Yuan, Z., de Haas, D., Hamilton, G., Corrie, S., and Keller, J.: Dynamics and

dynamic modelling of H2S production in sewer systems, Water Research, 42(10-11), 2527–

2538, 2008b.

Shih, T., Liou, W. W., Shabbir, A., Yang, Z., and Zhu, J.: A new k-eeddy viscosity model for

high reynolds number turbulent flows, Computers & Fluids, 24(3), 227–238, 1995.

Shuard, A., Mahmud, H., and King, A.: Comparison of two-phase pipe flow in OpenFOAM

with a mechanistic model, IOP Conference Series: Materials Science and Engineering,

121(1), 2016.

Smagorinsky, J.: General circulation experiments with the primitive equations: I. The basic

experiment, Monthly Weather Review, 91(3), 99–164, 1963.

Teuber, K., Broecker, T., Barjenbruch, M., and Hinkelmann, R.: High-Resolution Numerical

Analysis of Flow Over a Ground Sill Using OpenFOAM, in: Proceedings of the 12th

International Conference on Hydroscience & Engineering (ICHE), Tainan, Taiwan, 2016.

Teuber, K., Sielaff, M., Despot, D., Stephan, D., Barjenbruch, M., and Hinkelmann, R.: Mod-

eling and Measuring of Interfaces in Sewer Systems, in: Proceedings of the 37th IAHR

(International Association for Hydro-Environment Engineering and Research) World

Congress, Kuala Lumpur, Malaysia, 2017.

Teuber, K., Broecker, T., Bayón, A., Nützmann, G., and Hinkelmann, R.: CFD-modelling of

free-surface flows in closed conduits, Progress in Computational Fluid Dynamics, 19(6),

2019a.

Teuber, K., Broecker, T., Bentzen, T., Nützmann, G., and Hinkelmann, R.: Using compu-

tational fluid dynamics to describe H2S mass transfer across the water-air interface in

sewers, Water Science and Technology, 79(10), 1934–1946, 2019b.

159

Bibliography Appendix B

Teuber, K., Broecker, T., Jaydev, S. D., Goitom, G. M., Sielaff, M., Despot, D., Stephan, D.,

Barjenbruch, M., and Hinkelmann, R.: Multiphase CFD-Simulation of Transport Phenom-

ena in Sewer Systems, in: In: Mannina G. (eds) New Trends in Urban Drainage Modelling.

UDM 2018. Green Energy and Technology., Springer Cham., 2019c.

Teuber, K., Broecker, T., Nützmann, G., and Hinkelmann, R.: CFD simulation of H2S mass

transfer under turbulent conditions in a stirring tank, Water Science and Technology, in

preparation.

Thaker, J. and Banerjee, J.: CFD simulation of two-phase flow phenomena in horizontal

pipelines using openfoam, in: Proceedings of the 40th National Conference on Fluid

Mechanics & Fluid Power, 2013.

Thistlethwayte, D.: Control of sulphides in sewerage systems, Ann Arbor Science, 1972.

Thorenz, C. and Strybny, J.: On the numerical modelling of filling-emptying systems for

locks, in: 10th International Conference on Hydroinformatics, 2012.

Ubbink, O.: Numerical prediction of two fluid systems with sharp interfaces, Ph.D. thesis,

University of London London, UK, 1997.

Vollertsen, J., Nielsen, A., Yang, W., and Hvitved-Jacobsen, T.: Effects of in-sewer processes:

a stochastic model approach, Water science and technology, 52(3), 171–180, 2005.

Vollertsen, J., Nielsen, A., Jensen, H., Wium-Andersen, T., and Hvitved-Jacobsen, T.: Corro-

sion of concrete sewers—the kinetics of hydrogen sulfide oxidation, Science of the total

environment, 394(1), 162–170, 2008.

Wang, C., Xu, Z., Lai, C., and Sun, X.: Beyond the standard two-film theory: Computational

fluid dynamics simulations for carbon dioxide capture in a wetted wall column, Chemical

Engineering Science, 184, 103–110, 2018.

Ward, M., Hamer, G., McDonald, A., Witherspoon, J., Loh, E., and Parker, W.: A sewer

ventilation model applying conservation of momentum, Water Science and Technology,

64(6), 1374–1382, 2011.

Wilcox, D.: Reassessment of the scale-determining equation for advanced turbulence mod-

els, AIAA Journal, 26(11), 1299–1310, 1988.

Wu, H.: An issue on applications of a disk turbine for gas-liquid mass transfer, Chemical

Engineering Science, 50(17), 2801–2811, 1995.

Yakhot, V., Orszag, S., Thangam, S., Gatski, T., and Speziale, C.: Development of turbu-

lence models for shear flows by a double expansion technique, Physics of Fluids A: Fluid

Dynamics, 4(7), 1510–1520, 1992.

Yang, L., Nieves-Remacha, M. J., and Jensen, K. F.: Simulations and analysis of multiphase

transport and reaction in segmented flow microreactors, Chemical Engineering Science,

169, 106–116, 2017.

Yongsiri, C., Vollertsen, J., Rasmussen, M., and Hvitved-Jacobsen, T.: Air–water transfer of

hydrogen sulfide: an approach for application in sewer networks, Water Environment

Research, 76(1), 81–88, 2004.

160