Numerical simulation of single rising bubbles influenced by soluble surfactant in the spherical and ellipsoidal regime Numerische Simulation des Einflusses grenzflächenaktiver Substanzen auf aufsteigende Einzelblasen im sphärischen und ellipsoiden Regime Master-Thesis von Matthias Steinhausen Tag der Einreichung: 1. Gutachten: Prof. Dr. Dieter Bothe 2. Gutachten: Prof. Dr.-Ing. Peter Stephan Departement of Mathematics Mathematical Modelling and Analysis Numerical simulation of single rising bubbles influenced by soluble surfactant in the spherical and ellipsoidal regime Numerische Simulation des Einflusses grenzflächenaktiver Substanzen auf aufsteigende Einzelblasen im sphärischen und ellipsoiden Regime Vorgelegte Master-Thesis von Matthias Steinhausen 1. Gutachten: Prof. Dr. Dieter Bothe 2. Gutachten: Prof. Dr.-Ing. Peter Stephan Tag der Einreichung: Bitte zitieren Sie dieses Dokument als: URN: urn:nbn:de:tuda-tuprints-82963 URL: http://tuprints.ulb.tu-darmstadt.de/8296 Dieses Dokument wird bereitgestellt von tuprints, E-Publishing-Service der TU Darmstadt http://tuprints.ulb.tu-darmstadt.de tuprints@ulb.tu-darmstadt.de Die Veröffentlichung steht unter folgender Creative Commons Lizenz: Namensnennung 4.0 International https://creativecommons.org/licenses/by/4.0/ – Abstract – Surfactants are surface active agents that accumulate at fluid interfaces and influence interfacial properties, e.g. the surface tension. For single rising bubbles, even a small amount of surfactant causes Marangoni forces that influence the bubble rise significantly. In this work, Direct Numerical Simulations (DNS) with an Arbitrary Lagrangian-Eulerian (ALE) Interface-Tracking method are performed. The use of a subgrid-scale model enables the simulation of realistic time and length scales and the comparison with experiments. The resolution require- ments close to the interface are examined using 2D simulations to reduce the computational costs further. Then, 3D simulations of single rising bubbles under the influence of Triton-X100 are carried out, investigating different bubble diameters and initial surfactant bulk concentrations. The 3D simulations provide new insights into the transition from a helical motion into a zig-zag motion, which can only be observed in the presence of a surfactant. Additionally, the reciprocal influence of the local surfactant distribution on the interface and the vortex structures for path-unstable bubbles are analysed. Finally, the local surfactant distribution on the interface is modelled using a data-driven approach. The model is based on the DNS data obtained from the 3D simulations and is in good agreement with the validation data. In future work, the derived model can be used to improve existing simplified models for the simulation of bubbly flows under the influence of surfactant. Tenside sind oberflächenaktive Substanzen, die sich an fluiden Grenzflächen anlagern und die Eigenschaften der Grenzfläche, wie beispielsweise die Oberflächenspannung, beeinflussen. Bei aufsteigenden Einzelblasen führt bereits eine geringe Menge an Tensiden zu Marangoni-Kräften, die den Blasenaufstieg signifikant beeinflussen. In dieser Arbeit werden Direkte Numerische Simulationen (DNS) mit einer Arbitrary Lagrangian-Eulerian (ALE) Interface-Tracking Methode durchgeführt. Die Verwendung eines Subgridskalen-Modells ermöglicht die Simula- tion von realistischen Zeit- und Längenskalen und den Vergleich mit experimentellen Daten. Die Auflösungsan- forderungen in der Nähe der Grenzfläche werden mithilfe einer 2D Studie untersucht, um den Rechenaufwand weiter zu reduzieren. Daraufhin, werden 3D-Simulation von Einzelblasen unter dem Einfluss von Triton-X100 durchgeführt und verschiedene Blasendurchmesser und Tensidkonzentrationen in der Flüssigphase untersucht. Die Simulationen geben neue Einblicke in die Blasendynamik. Besonders interessant ist dabei der Übergang von einem helikalen Aufstiegspfad zu einem Zig-Zag Pfad, der nur in Gegenwart von Tensiden beobachtet werden kann. Zusätzlich, wird die wechselseitige Beeinfussung der lokalen Tensidverteilung auf der Grenzfläche und den Wirbelstrukturen von pfadinstabilen Blasen analysiert. Zuletzt wird die lokale Tensidverteilung auf der Gren- zfläche mithilfe eines daten-basierten Ansatzes modelliert. Das hergeleitete Modell basiert auf den DNS Daten der durchgeführten 3D Simulationen und ist in guter Übereinstimmung mit den Validierungsdaten. Das hergeleitete Modell kann im nächsten Schritt dazu verwendet werden, um bestehende Modelle für den Einfluss von Tensiden auf Blasenströmungen zu verbessern. i Table of Contents Abstract i Table of Contents iii List of Figures v List of Tables vii Glossary ix 1 Introduction 1 2 Governing equations and algorithms 3 2.1 Mathematical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.1 Hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.2 Surfactant transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.3 Surfactant influence on the surface tension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1.4 Forces acting on the interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.5 Dimensionless quantities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.1 Spatial discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.2 Temporal discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.3 OpenFOAM-specific boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.4 Pressure-velocity coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Machine learning algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.1 Machine learning domains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.2 Machine learning workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.3 Feature selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.4 Multilayer perceptron . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3 Study on the mesh requirements 15 3.1 Mesh study for a clean 2D bubble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.1 Radial grid resolution: liquid phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.2 Radial grid resolution: gaseous phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1.3 Tangential grid resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 3D meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.3 Mesh requirements: conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4 Parallelisation study 23 4.1 Decomposition techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.2 Comparison of the decomposition techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.3 Domain dependency study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.3.1 Manual decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.3.2 Scotch decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.4 Parallelisation study: conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 5 Simulation results and discussion 29 5.1 Experimental studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.2 Mesh sensitivity study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.2.1 Mesh sensitivity for dB = 0.8 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.2.2 Mesh sensitivity for dB = 1.3 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.2.3 Mesh sensitivity for dB = 2.0 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.3 Bubble path and terminal velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.3.1 Bubble path and terminal velocity for dB = 0.8 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.3.2 Bubble path and terminal velocity for dB = 1.3 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.3.3 Bubble path and terminal velocity for dB = 2.0 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 iii 5.4 Forces acting on the interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.4.1 Forces acting on the interface for dB = 0.8 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.4.2 Forces acting on the interface for dB = 1.3 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.4.3 Forces acting on the interface for dB = 2.0 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.5 Local bulk velocity and surface fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.5.1 Local bulk velocity and surface fields for dB = 0.8 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.5.2 Local bulk velocity and surface fields for dB = 1.3 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.5.3 Local bulk velocity and surface fields for dB = 2.0 mm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.6 Surfactant distribution on the interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.6.1 Surfactant distribution on the interface for a rectilinear rise . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.6.2 Surfactant distribution on the interface for a zig-zag rise . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.6.3 Surfactant distribution on the interface for a helical rise . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 6 Modeling of the local surfactant distribution on the interface 53 6.1 Data pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.1.1 Data filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.1.2 Data averaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.2 Feature selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.2.1 Linear feature correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.2.2 Feature importance: random forest regressor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 6.2.3 Sequential backward selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.2.4 Feature selection: conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.3 Model architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.4 Model training . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 6.5 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 6.6 Influence of training data distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.7 Model of the local surfactant distribution: conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 7 Summary and Outlook 67 Bibliography 69 Acknowledgements xii Thesis Statement xiv iv List of Figures 1.1 Applications and effects of surfactants in bubbly flows. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Bubble path under the influence of surfactant. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2.1 Sketch of the computational domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Overview of the numerical solution procedure of the solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Pressure residual compared with terminal velocity vy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4 Fluctuation in the drag force, nmin = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.5 Structure of a multilayer perceptron. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.1 Influence of different mesh parameters on the grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Total pressure force acting on the interface for different outer radial resolutions, δh/∆Rout,1. . . . . . . . . . . 17 3.3 Total pressure force at the interface for different ratios ∆Rin,1/∆Rout,1. . . . . . . . . . . . . . . . . . . . . . . . . 18 3.4 Force balance normal to the interface with different radial resolutions, δh/∆Rout,1. . . . . . . . . . . . . . . . . 19 3.5 Bubble deformation for different radial resolutions, δh/∆Rout,1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.6 Force balance normal to the interface with different tangential resolutions, dB/∆Tan. . . . . . . . . . . . . . . 20 3.7 Computational mesh for the 3D simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.8 Interface mesh of a 2.0 mm bubble. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.1 Manual decomposition on ten processors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.2 Single processor domain of the scotch decomposition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.3 Scotch decomposition on ten processors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.4 Runtime comparison between the manual and scotch decomposition. . . . . . . . . . . . . . . . . . . . . . . . . 25 4.5 Runtime comparison for different cell ratios χProc for the manual decomposition. . . . . . . . . . . . . . . . . . 26 4.6 Runtime comparison for different cell ratios χProc with scotch decomposition and varying processor weights. 27 4.7 Runtime comparison for different cell ratios χProc with scotch decomposition and varying interface domain. 28 5.1 Trajectories of rising bubbles under the influence of Triton X-100. . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.2 Mesh comparison of the bubble terminal velocity, dB = 0.8 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.3 Mesh comparison of the bubble terminal velocity, dB = 1.3 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5.4 Mesh comparison of the bubble path, dB = 1.3 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5.5 Mesh comparison of the bubble terminal velocity, dB = 2.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.6 Mesh comparison of the bubble path, dB = 2.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 5.7 Region of the mesh showing fluctuations at the interface before the crash. . . . . . . . . . . . . . . . . . . . . . 33 5.8 Terminal velocity under the influence of Triton X-100. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.9 Rise path under the influence of Triton X-100, dB = 1.3 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.10 Rise path under the influence of Triton X-100, dB = 2.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.11 Sketch of the direction of the drag and lift force. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.12 Drag forces acting on the bubble, dB = 0.8 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.13 Drag forces acting on the bubble, dB = 1.3 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.14 Lift forces acting on the bubble, dB = 1.3 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.15 Drag forces acting on the bubble, dB = 2.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.16 Lift forces acting on the bubble, dB = 2.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.17 Lift force direction along the bubble path for dB = 2.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.18 Lift force contributions in the top view of the helical path , dB = 2.0 mm. . . . . . . . . . . . . . . . . . . . . . . 43 5.19 Local pressure field and surfactant distribution on the interface, dB = 0.8 mm. . . . . . . . . . . . . . . . . . . 44 5.20 Local velocity field around the bubble, dB = 1.3 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.21 Local velocity field around the bubble, dB = 2.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.22 Local surfactant distribution, zig-zag rise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.23 Local surfactant distribution, helical rise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.24 Interface velocity influenced by the bulk velocity field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.25 New coordinate system to evaluate the surfactant distribution on the interface. . . . . . . . . . . . . . . . . . . 49 5.26 Surfactant distribution for a straight rise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.27 Surfactant distribution for a zig-zag rise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.28 Surfactant distribution for a helical rise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 6.1 Local distributions of pressure, velocity and surfactant on the interface, t = 0.01 s and dB = 0.8 mm. . . . . 54 6.2 Filtered local surfactant concentration on the interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.3 Local surfactant concentration on the interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.4 Scatterplot cΣ, c̄Σ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.5 Correlation matrix with cΣ, dB = 0.8 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 6.6 Scatterplot c̃Σ, ϕ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 6.7 Correlation matrix with c̃Σ, dB = 0.8 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 v 6.8 Feature importance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 6.9 Sequential backward selection: r2-score. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.10 Terminal velocity, dB = 0.8 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 6.11 Local surfactant concentration on the interface for low bulk concentrations, test data. . . . . . . . . . . . . . . 60 6.12 Local surfactant concentration on the interface for high bulk concentrations, test data. . . . . . . . . . . . . . 61 6.13 Local surfactant concentration on the interface for low bulk concentrations, validation data. . . . . . . . . . . 63 6.14 Local surfactant concentration on the interface for high bulk concentrations, validation data. . . . . . . . . . 64 6.15 Validation scores for different combinations of input data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.16 Influence of the initial surfactant concentration on cΣ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 vi List of Tables 2.1 Numerical schemes used for the simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Boundary conditions of the numerical domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.1 Fluid properties of the gaseous and liquid phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Default setup for the 2D mesh study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.3 Parameters for the outer radial grid resolution study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4 Parameters for the inner radial grid resolution study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5 Parameters for the tangential grid resolution study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.6 Grid resolution and mesh size for dB = 0.8 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.7 Grid resolution and mesh size for dB = 1.3 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.8 Grid resolution and mesh size for dB = 2.0 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.1 Number of cells on the decomposed processors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.2 Parameters for the scotch decomposition domain study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 5.1 Surfactant properties (Triton-X100). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 6.1 Setup parameters for the MLP in tensorflow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 vii Glossary Greek symbols Symbol Unit Description β m s−1 Mass transfer coefficient δh m Hydrodynamic boundary layer thickness ∆ m Cell thickness ε - Threshold κ m−1 Curvature λ - Learning rate µ kg m−1 s−1 Dynamic viscosity ν m2 s−1 Kinematic viscosity ρ kgm−3 Density σ kgm s−2 Surface tension coefficient σ - Standard deviation Σ(t) - Deformable interface ϕ - Polar angle χProc - Ratio between the number of cells on the processor holding the interface and the mean number of cells on an outer processor ψ - Azimuthal angle Ω - Liquid domain ∇ - Divergence Latin symbols Symbol Unit Description a mol m−2 Langmuir constant A m2 Face area c mol m−3 Surfactant concentration in the bulk cΣ mol m−2 Surfactant concentration on the interface cΣeq mol m−2 Equilibrium concentration cΣ∞ mol m−2 Saturated surfactant concentration d m Diameter D m2 s−1 Molecular diffusivity E - Loss function f kg s−2 m−1 Area specific force vector F kgm s−2 Integral force vector acting on the bubble g ms−2 Gravitation constant I - Identity matrix j molm−3 s−1 Diffusive flux lPN m Distance between two face centers L m Characteristic length n - Number of processors n - Normal vector N - Number of cells p kgm−1 s−2 Pressure pdyn kgm−1 s−2 Dynamic pressure ptot kgm−1 s−2 Total pressure r m Radius R kgm2s−2mol−1K−1 Universal gas constant Re - Reynolds number R - Correlation coefficient s molm−2 s−1 Sorption source term Svisc kgm−1 s−2 Viscous stress tensor ix Symbol Unit Description Sh - Sherwood number t s Time T K Temperature U ms−1 Characteristic velocity v ms−1 Barycentric velocity V m3 Volume w - Weights of the MLP Indicees Index Description B Bubble Bot Bottom/Wake region ca Capillary pressure Drag Drag contribution Eq Equator region G Gaseous phase loc Local value L Liquid phase Lift Lift contribution ma Marangoni max Maximum value min Minimum value model Model estimate num DNS value Proc Processor Rin Component close to the interface, in the gaseous phase Rout Component close to the interface, in the liquid phase Σ Interfacial value Tan Tangential cells on the interface Top Top region Total Total visc Viscous ∥ Acting tangential to the interface ⊥ Acting normal to the interface |Σ Defined in the bulk phase at the interface ·̃ Non-dimensional quantity ·̄ Mean quantity x Abbreviations Abbreviation Description ALE Arbitrary Lagrangian-Eulerian CFD Computational Fluid Dynamics CFL Courant-Friedrichs-Lewy DNS Direct Numerical Simulation OpenFOAM Open Field Operation And Manipulation FV Finite Volume FA Finite Area MLP Multilayer Perceptron PBiCG Preconditioned biconjugate Gradient PCG Preconditioned Conjugate Gradient PISO Pressure Implicit with Splitting of Operators SGS Subgrid-scale VOF Volume-Of-Fluid xi 1 Introduction Bubbly flow plays an important role in a variety of technical applications such as bubble column reactors, flotation processes or waste-water treatment. Within this applications, surface active agents, so-called surfactants, are present. The surfactants are either a side effect of contamination or added on purpose to influence the interaction between the phases. A well-known example of bubbly flow under the influence of surfactant is froth flotation, which was first developed in mineral processing to separate hydrophobic and hydrophilic materials. The raw ores are crushed and dissolved in water to form a slurry. Then, a surfactant is used to render the desired minerals hydrophobic. The slurry is transferred into a bubble column, the so-called flotation cell, that is aerated. The hydrophobic particles attach to the bubbles that rise to the surface of the flotation cell where they form a froth that can be extracted. Nowadays froth flotation is also used in paper recycling and waste-water treatment. The Direct Numerical Simulation (DNS) of a flotation cell requires to resolve a wide range of length scales. These scales range from the column size down to scales smaller than a single bubble. The high computational costs of such simulations make them impracticable for real applications. Therefore, scale-reduced approaches have been developed that model the phenomena occurring at small scales like the influence of surfactant on a single rising bubble. Figure 1.1 shows the schematic of a flotation cell and the influence of surfactants on a bubble swarm. (a) Flotation cell reproduced from [17]. (b) Bubble swarm under the influence of methyl isobutyl carbinol (MIBC) and sodium chlorid (NaCl) [9]. Figure 1.1: Applications and effects of surfactants in bubbly flows. One of the most fascinating effects involving surfactant is that even small traces of it can lead to microscopical forces that cause a tremendous change in the macroscopic flow pattern. Figure 1.1b shows an experiment performed by Finch et al. [9] displaying the influence of surfactant on a bubble swarm. At the left-hand side, a bubble swarm under the influence of methyl isobutyl carbinol (MIBC) is shown. An concentration of 10 ppm MIBC already has a significant effect on the flow pattern, the void fraction, and the bubble size distribution. This effect, however, complicates the modelling of contaminated bubbles. While for a uncontaminated bubble rising in water a trailing vortex in the wake occurs at Reynolds numbers above 600, for contaminated bubbles this effect is already observed at Reynolds numbers as low as 20 [19, 38]. The vortex influences the bubble rise velocity and path, leading to terminal velocities for contaminated bubbles that are up to two times smaller than for uncontaminated ones [7]. For instance, figure 1.2 shows the influence of surfactant on the bubble path from experiments performed by Tagawa et al. [31]. While the uncontaminated bubble (left) shows a helical rise with a helix-width of approximately one bubble diameter, small traces of surfactant may increase (b) or decrease (c) the lateral motion, or even change the motion from helical to zig-zag (d). 400 300 200 100 –2 0 2 –2 0 2 500 0 400 300 200 100 –2 0 2 –2 0 2 500 0 400 300 200 100 –2 0 2 –2 0 2 500 0 400 300 200 100 –2 0 2 –2 0 2 500 0 (a) (b) (c) (d) Figure 1.2: Bubble path under the influence of surfactant [31]. (a) 0 ppm; (b) 25 ppm; (c) 50 ppm; (d) 75 ppm. 1 Bubbles rising in impure liquids show terminal velocities that vary between the terminal velocity of particles with a fully mobile and an immobile (rigid) interface. In the absence of impurities, a rising bubble is characterised by a mobile interface, i.e., the fluid elements at the interface are able to move and can be exchanged or displaced. This leads to smaller velocity gradients in the surrounding liquid compared to liquid motion around a solid particle. Consequently, less energy is dissipated, which leads to higher rise velocities of uncontaminated bubbles. In the presence of contamination, however, surfactant adsorbs onto the interface and is accumulated in the rear part of the bubble, reducing the mobility of the interface [23, 18]. Within this area, the interface acts almost like a solid surface. This phenomenon paired with experimental observations led to the idea of the so-called stagnant cap model, presented by Davis and Acrivos [8], which is described in the following. The interface is divided into two segments symmetric around the rise velocity vector: one fully covered with surfactant, acting like a solid particle, and one entirely clean, with a mobile interface. While the rear part of the bubble, the stagnant cap, shows a zero velocity in a coordinate system moving with the bubble centre, the bubble front is characterised by zero shear stresses. The stagnant cap angle denotes the polar angle where the interface changes from fluid-like to solid-like behaviour. This strict separation between a stagnant cap and the bubble front is a strong idealisation and is only valid for cases where the transition zone from high to low surface contamination is small compared to the bubble size. Furthermore, in most applications, the condition of symmetry around the rise velocity vector is rather an exception than the rule [25]. An uneven surfactant distribution on the interface leads to surface tension gradients that cause so-called Marangoni forces acting from areas with low to areas with high surface tension. The Marangoni forces acting on the interface are balanced by shear forces in the liquid phase. Due to the high contamination at the rear part and the low one in the front of the bubble, the main contribution of these forces acts opposed to the bubble rise direction and increases the overall drag force. In the case of path instability, however, the surfactant distribution on the interface is non-symmetric along the direction azimuthal to the bubble movement direction. This asymmetry leads to additional lift forces acting on the bubble. This effect has been recently encountered by Pesci et al. [25] and is not yet completely understood. Advances in computational fluid dynamics enabled the simulation of two-phase flows with deforming interfaces. The current simulation techniques are still under development and lead to high computational costs [11]. The simulation of complex processes, e.g. flotation cells, requires simplified simulation methods. To capture the effect of contamination, the surfactant distribution on the interface needs to be modelled. A recent study by Fleckenstein and Bothe [11] models the Marangoni stresses at the interface based on the assumptions of stationary bubble shapes and surfactant distribu- tions. A Navier-type jump condition at the interface is derived that accounts for the Marangoni stresses based on the tangential interface velocity in a coordinate system co-moving with the bubble’s barycenter and a friction parameter. The friction parameter is based on the diffusion constant of the interfacial surfactant transport and the Gibbs elasticity, which is related to the surfactant distribution on the interface. In the model, the Gibbs elasticity is considered constant for an individual bubble which implies a linear dependency of the change of surface tension and the surfactant concen- tration on the interface. This assumption is only valid for very small concentration differences. To extend the model’s range of validity for unsteady surfactant distributions on the interface and higher surfactant concentrations, the angular dependency of the surfactant concentration at the interface can be accounted for in the Gibbs elasticity. This, however, requires knowledge about the local and temporal surfactant distribution on the interface. Within this work, the local, time-dependent concentration profiles are reconstructed based on a statistical learning algorithm that leverages DNS data of single rising bubbles under the influence of surfactant. In the first part of this work, the influence of surfactants on single rising bubbles in the spherical and ellipsoidal regime is examined, extending the investigations done by Pesci et al. [24, 25]. An Arbitrary Lagrangian-Eulerian (ALE) interface tracking approach [22, 14, 33, 34] is used combined with a recently introduced subgrid-scale (SGS) model for the surfactant adsorption onto the interface [37]. This model reduces the resolution requirements close to the interface significantly, and thereby enables the simulation of realistic systems and the comparison of the numerical results with experimental data. Despite the application of the SGS model, the simulations have high computational costs due to resolution requirements at the interface and the simulation of long physical times. Therefore, a 2D setup is used to examine the tangential and radial mesh resolution requirements of the 3D simulations. Additionally, different parallelisation techniques are investigated, comparing a manual decomposition technique [36] with a decomposition based on the scotch algorithm [6]. The obtained simulation results are compared with an experimental study performed by Tagawa et al. [31], and the bubble terminal velocities, bubble paths and the local surfactant distribution on the interface are evaluated. In the second part of this work, the simulation results are assessed to model the surfactant distribution on the interface using a data-driven approach. 2 2 Governing equations and algorithms 2.1 Mathematical model In this work, the two-phase flow problem is modelled as follows: the fluid domain is separated into a liquid and a gaseous phase by the moving interface Σ(t) with unknown time-dependent shape and location. The interface is a surface of zero thickness and, therefore, the model is usually referred to as sharp interface representation. Furthermore, the assumptions of incompressible Newtonian fluids, isothermal conditions and the absence of phase change and chemical reactions are made. The governing equations are based on the conservation of mass, momentum, and surfactant concentration. Additionally, only the surfactant transport in the liquid phase and on the interface is taken into account. 2.1.1 Hydrodynamics The two-phase flow problem is described by the two-phase Navier-Stokes equation for incompressible Newtonian fluids. The continuity and momentum equation for the gaseous and liquid phase are given by equations (2.1) and (2.2): ∇ · v= 0 , (2.1) ∂t (ρv) +∇ · (ρv⊗ v) = −∇p+∇ · Svisc +ρg , (2.2) where v is the barycentric velocity, p the pressure, ρ the density of the gas or the liquid and g the gravitational constant. The viscous stress tensor is given by Svisc = µ � ∇v+ (∇v)T � , with µ being the dynamic viscosity of the respective fluid. Additionally, the two phases are coupled via jump conditions at the interface: JvK= 0 , (2.3) v · nΣ = vΣ · nΣ , (2.4) JpI− SviscK · nΣ = σκnΣ +∇Σσ , (2.5) where vΣ is the interface velocity and κ = −∇Σ · nΣ the surface curvature, with ∇Σ being the surface gradient. For contaminated systems, the surface tension coefficient σ is a function of the surfactant concentration σ = σ � cΣ � . A jump of a physical quantity across the interface is denoted by J·K and is defined as: JφK(t,x) = lim h→0+ (φ (t,x+ hnΣ)−φ (t,x− hnΣ)) , x ∈ Σ(t) . (2.6) The boundary condition at the outer domain is given by a zero velocity Dirichlet boundary condition and reads: v= 0 on ∂Ω . (2.7) 2.1.2 Surfactant transport The model of the surfactant transport on the interface is essential to simulate the influence of contamination on the bubble motion. In this work, only the surfactant transport in the liquid domain Ω and on the interface Σ(t) is taken into account. The concentration in the gaseous phase is assumed to be zero. The local surfactant transport equations, derived for example in [2, 30], are: ∂t c +∇ · (cv+ j) = 0 in Ω \ Σ(t) , (2.8) ∂ Σt cΣ +∇Σ · cΣvΣ + jΣ = sΣ on Σ(t) , (2.9) with c being the molar concentration of surfactant in the bulk phase, cΣ being the molar surface concentration of surfac- tant on the interface, and j and jΣ denoting the diffusive fluxes in the liquid domain and on the interface, respectively. The sorption term sΣ reads: sΣ + Jj · nΣK= 0 on Σ(t) . (2.10) Equation (2.9) is a dynamic boundary condition for equation (2.8). The initial conditions are: c (0,x) = c0 (x) , x ∈ Ω(0) , (2.11) cΣ (0,x) = cΣ0 (x) , x ∈ Σ(0) . (2.12) 3 The system of equations (2.8) to (2.10) is not closed. Additional relations are required that define the diffusive fluxes and source terms as functions of the primitive variables. The surfactant can be considered diluted thus the diffusive fluxes j, jΣ can be modelled by Fick’s law: j= −D∇c in Ω(t) , (2.13) jΣ = −DΣ∇ΣcΣ in Σ(t) , (2.14) where D and DΣ represent the diffusion constant in the respective fluid and on the interface, respectively. At the outer boundary of the liquid phase, a zero flux boundary condition is imposed: j · n= 0 on ∂Ω(t) . (2.15) Additionally, the sorption term sΣ in equation (2.9) needs to be modelled. Two limiting sorption processes can be distinguished: diffusion-controlled and kinetically-controlled sorption [20]. In this work, a diffusion-controlled sorption model has been adopted [16]. In this case, the ad- and desorption rates are locally in equilibrium: sads � c|Σ, cΣ � = sdes � cΣ � . (2.16) This local relationship between the surfactant concentration on the interface cΣ and the surfactant bulk concentration close to the interface c|Σ is called adsorption isotherm and needs to be accounted for in the numerical solution1. For the Langmuir adsorption model, the adsorption isotherm reads: cΣ = cΣ∞ c/a 1+ c/a , (2.17) where c is the surfactant bulk concentration close to the interface, a the Langmuir equilibrium constant in mol/m3 and cΣ∞ is the saturated surfactant concentration, which is the maximum number of molecules per area on the interface and a surfactant specific constant. Additionally, the equilibrium concentration cΣeq can be calculated using the maximum bulk concentration in equa- tion (2.17), which is the initial surfactant bulk concentration c0: cΣeq = cΣ∞ c0/a 1+ c0/a . (2.18) 2.1.3 Surfactant influence on the surface tension The presence of surfactant on the interface changes the interfacial surface tension and thereby influences the forces acting on the interface. For fast and slow sorption mechanisms, the effect of surfactant on the surface tension is given by the equation of state, and reads: σ−σ0 = Π � cΣ � . (2.19) The function Π � cΣ � is defined with respect to the sorption model employed [24]. In the current work, the Langmuir model is used and, therefore, the surface tension equation of state is given by: σ = σ0 + RT cΣ∞ ln � 1− cΣ cΣ∞ � , (2.20) where R is the universal gas constant and T is the absolute system temperature in Kelvin. 1 The notation ·|Σ denotes the trace of a quantity defined in Ω at the interface. 4 2.1.4 Forces acting on the interface Equation (2.5) can be used to derive the forces acting on the interface. The normal force balance at the interface is given by equation (2.21), while the tangential force balance leads to equation (2.22): JptotKnΣ + 2JµK (∇Σ · v)nΣ = σκnΣ , (2.21) −Jµn · ∇vK− Jµ (∇Σv) · nK− JµKnΣ (∇Σ · v) =∇Σσ , (2.22) where L denotes the liquid and G the gaseous phase. The area-specific forces f j i can be divided into: • the Marangoni force, which is the result of a non-uniform surfactant distribution on the interface fma =∇Σσi . (2.23) • the capillary pressure force, which increases either with the surface tension σ or with the curvature κ of the interface. The capillary pressure force is large for small, uncontaminated bubbles and for deformed bubbles with interfacial regions of high curvature fca = σikinΣ,i . (2.24) • the dynamic pressure, which is the kinetic energy per volume of fluid and thereby increases with the velocity of fluid elements fpdyn = � pdyn,Gi − pdyn,Li � nΣ . (2.25) • the total pressure force, which is the sum of the dynamic, ambient and hydrostatic pressure force. For a freely rising bubble, the ambient pressure difference is zero. fptot = � ptot,Gi − ptot,Li � nΣ . (2.26) • the viscous force, which acts between the two phases (liquid and gaseous). It is determined by the velocity gradient at the interface and can be separated into parts normal and tangential to the interface: – The normal viscous force is defined as: fvisc ⊥ = 2µG � ∇Σ · v � inΣ,i − 2µL � ∇Σ · v � i nΣ,i . (2.27) – While the tangential viscous force reads: fvisc ∥,Li = µL � (n · ∇v)Li + (∇Σv)Li · n+ nΣ (∇Σ · v) � , (2.28) fvisc ∥,Gi = µG � (n · ∇v)Gi + (∇Σv)Gi · n− nΣ (∇Σ · v) � , (2.29) fvisc ∥,i = fvisc ∥,Li + fvisc ∥,Gi . (2.30) By integrating the area specific forces f j i over the bubble surface, the global forces acting on the interface are obtained. Due to the discrete representation the integration transforms into a sum over all faces N f : F j = N f ∑ i=1 f j i Ai . (2.31) The jump condition for the global forces acting on the interface reads: Fptot + Fvisc = Fca + Fma . (2.32) Note that the jump condition (2.32) is only fulfilled approximately in the numerical procedure. Thus, the equality of terms on the left and right side of equation (2.32) may be used as a quality measure for the coupling of the two phases. 5 2.1.5 Dimensionless quantities To describe the physical properties and behaviour of the simulated bubbles in a case-independent manner, it is useful to introduce dimensionless quantities. The dimensionless quantities are described in the following: • Reynolds number: The Reynolds number describes the ratio between inertia and viscous forces and is an indicator for the appearance of distinct flow patterns. It is defined as: Re = U · L ν = v · dB ν , (2.33) with U being the characteristic velocity, L the characteristic length of the flow problem and ν the kinematic viscosity of the fluid. Within this work, the characteristic velocity is given by the barycentric velocity v and the characteristic length by the bubble diameter dB. Additionally to the global Reynolds number describing the motion of the bubble center, a local Reynolds number is defined, capturing the motion on the interface. The local Reynolds number is defined as: ReΣ = vΣ · dB ν , (2.34) with vΣ being the relative interface velocity of the bubble in a coordinate system moving with the bubble centre. • Sherwood number: The Sherwood number is a measure for the species transfer rate and represents the ratio of convective and diffusive mass transfer. It is defined as: Sh= β L D , (2.35) with β being the mass transfer coefficient, L the characteristic length of the transfer problem and D the diffusion coefficient. Within this work, the Sherwood number is used to characterise the surfactant transfer at the interface. In order to extract the Sherwood number from the simulation results, the local Sherwood number is defined as: Shloc = (n · ∇c)|Σ · dB c0 . (2.36) with (n · ∇c)|Σ being the gradient in normal direction to the interface, c the bulk concentration close to the interface and c0 the concentration far away from the interface. The global Sherwood number is determined using the area-weighted integral of the local Sherwood number [10]: Sh= 1 ∂ V ∫ ∂ V ShlocdA = ∑N i=0 Shloc,i · Ai ∑N i=0 Ai , (2.37) where Shloc,i is the cell specific local Sherwood number, Ai the cell area and N the number of cells at the interface. • Dimensionless time: In order to define a dimensionless time, the gravitational time tgr = p dB/g is used [4]. The dimensionless time is therefore given by: t̃ = t tgr = t p dB/g , (2.38) with dB being the bubble diameter and g the gravitational constant. • Dimensionless acceleration: For freely rising bubbles in water, the driving force acting on the bubble results from the buoyancy. Thus, the acceleration is normalized using the gravitational constant g. The dimensionless acceleration is defined as: ã = a g . (2.39) • Dimensionless area: Due to forces acting on the interface the bubble is deformed, resulting in a change of the bubble area. In order to compare the area deviation of the different cases, the bubble area is normalised using the area of a sphere with equivalent volume Asphere: Ã= A Asphere = A 1/6πd3 B . (2.40) 6 2.2 Numerical setup The equations presented in section 2.1 are solved numerically using the open-source platform OpenFOAM1 (Open Field Operation And Manipulation). OpenFOAM is a free to use toolbox for the development of numerical solvers for the solution of continuum mechanical problems. The main focus lies on computational fluid dynamics (CFD). To solve the surfactant transport equations in the bulk and on the interface, the arbitrary Lagrangian-Eulerian (ALE) Interface Tracking method is used, which was originally presented by Hirt et al. [12] and further developed by Muzaferija and Perić [22] and Tuković and Jasak [34]. The interface is represented by a surface mesh that divides the computa- tional domain into two subdomains: the liquid and the gaseous phase. The two subdomains are coupled via boundary conditions for pressure and velocity derived from the jump conditions at the moving interface Σ(t), as described in the previous chapter. A detailed description of the solver can be found in [24, 25]. In the following, the basic ideas are outlined. To reduce the required number of cells, a moving reference frame fixed to the bubble centre is used, i.e., the bubble remains in the centre, while the fluid flows around it. Figure 2.1 shows a sketch of the computational domain. inflowing liquid rB gas liquid space interface Figure 2.1: Sketch of the computational domain. In the simulation the interface is duplicated, consisting of two surface meshes: interface and interfaceShadow being the boundaries of the liquid and gas phase, respectively. 2.2.1 Spatial discretisation The equations in the bulk phases and on the interface are solved by applying a Finite Volume (FV) and Finite Area (FA) method, respectively. Table 2.1 reports the schemes used to discretise each term of the transport equations. Table 2.1: Numerical schemes used for the simulation. (a) FV schemes. Scheme Setting gradScheme Gauss linear divScheme default Gauss linear divScheme div(phi,U) Gauss GammaVDC 0.5 divScheme div(phi,C) Gauss limitedLinear 1.0 lapacianScheme Gauss linear corrected interpolationScheme default linear interpolationScheme C limitedLinear phi 1.0 snGradScheme corrected (b) FA schemes. Scheme Setting gradScheme Gauss leastSquares2 divScheme div(Us) Gauss leastSquares2 div(phi,Cs) Gauss upwind lapacianScheme Gauss linear corrected interpolationScheme default leastSquares2 1 Link: https://www.openfoam.com/ 2 The leastSquares edge interpolation scheme is part of an in-house development and not yet published. 7 https://www.openfoam.com/ 2.2.2 Temporal discretisation For the temporal discretisation, the second order backward-differencing time scheme, known in OpenFOAM as backward, with a fixed time step ∆t is used. The simulation time step ∆t is chosen with regard to the numerical stability of the interface [34], which is a more restrictive condition for surface tension driven flows than the CFL (Courant-Friedrichs- Lewy) condition. As a stability criterion the time step should satisfy: ∆t <  √ρmin (lPN) 3 2πσ , (2.41) with min (lPN) being the minimum distance between two face centers on the interface. 2.2.3 OpenFOAM -specific boundary conditions In addition to the numerical schemes, boundary conditions for the numerical domain need to be defined. The computa- tional domain contains three boundaries: the outer boundary of the liquid domain, here called space, the interface and the interfaceSchadow, see figure 2.1. Table 2.2 lists the OpenFOAM-specific boundary conditions used in the simulations. Table 2.2: Boundary conditions of the numerical domain. Boundary U motion U p c interfaceSchadow fixedValue fixedValue fixedGradient zeroGradient interface fixedGradient fixedValue fixedValue fixedValue space inletOutlet fixedValue zeroGradient inletOutlet The initial surfactant concentration on the interface is zero, while in the bulk phase, the surfactant is uniformly dis- tributed over the liquid domain with a given concentration c0. 2.2.4 Pressure-velocity coupling For the pressure-velocity coupling a modified version of the iterative pressure implicit with splitting of operators (PISO) algorithm is used [13]. The linear solver tolerance for the pressure (PCG) and velocity (PBiCG) solvers is set to 10−6, while for the concentration (PBiCG) a much lower bound of 10−12 is deployed. Figure 2.2 shows a schematic overview of the numerical solution procedure of the subgrid-scale (SGS) solver used in this work and first presented in [25]. In the SGS solver implemented by Pesci et al. [25], the pressure-velocity coupling with a modified PISO algorithm is embedded in an outer loop. Besides the PISO algorithm, the interface displacement, mesh and fluxes, as well as the interface tangential and normal jump conditions, and surfactant transport at the interface are updated, see figure 2.2. Thereby, the number of these outer iterations is considered constant. Thus a compromise between a steady simulation and runtime has to be made. Figure 2.3 shows the terminal velocity vy and the initial pressure residual for the 2D simulation of an uncontaminated bubble over time. In the acceleration phase of the bubble the residual rises. Once the bubble decelerates, the residual drops and approaches zero. The residual trend is an indication for the numerical effort needed to solve the problem. While in the first acceleration phase, the simulation is demanding, requiring a large num- ber of outer iterations, for later times fewer outer iterations are sufficient to ensure a satisfying pressure-velocity coupling. To stabilise the simulation in the acceleration phase and to reduce the computational effort in the quasi-steady state, a residual control is implemented varying the number of outer iterations. Therefore, the initial residuals for pressure and velocity are stored and compared to user-defined tolerances. Once the initial residuals fall below the tolerance, the outer loop is terminated and the simulation continues with the next time step. The break condition reads: � Resp ≤ Tolp ∧Resv ≤ Tolv � ∨ (n≥ nmax) , (2.42) with Res denoting the residual, Tol the given tolerance, n being the number of iterations and nmax the maximum number of outer iterations to avoid an infinite loop. In a first attempt, a pressure and velocity tolerance of 10−6 with nmax = 30 is used. The residual control was tested for a 3D bubble with a bubble diameter of dB = 1.3 mm. 8 Field initialization Initialize/Update c, cΣ, p, U from the previous time level a. Compute interface displacement & mesh and fluxes update b. Update tangential component of momentum jump condition at Σ c. Surfactant surface transport d. Update normal component of momentum jump condition at Σ e. PISO Loop (modified Rhie-Chow interpolation) f. Mesh update Further outer iteration? Compute SGS model parameter Surfactant bulk transport Fast sorption Fast sorption Next Time Step { sΣ fast = −DSGS(n · ∇c)Σ { c|Σ = f−1(cΣ) Time loop Outer loop Yes No Figure 2.2: Overview of the numerical solution procedure of the solver [25]. With this setup fluctuations in pressure and velocity occur once the number of outer iterations falls below three, see figure 2.4a. The reason for the fluctuation is unclear and needs further investigation. A possible explanation might be the influence of additional steps performed inside the outer loop, which are not captured by the pressure and velocity residuals, e.g. the mesh motion. To avoid the fluctuation, within this work, a minimal number of outer iterations nmin is specified, changing the break condition of the loop to: �� Resp ≤ Tolp ∧Resv ≤ Tolv � ∧ n≥ nmin � ∨ (n≥ nmax) . (2.43) Figure 2.4b shows the fluctuation in the total drag force of a 1.3 mm bubble under the influence of surfactant for the residual control with nmin > 0 and nmin = 0. With more than three minimal outer iterations, the fluctuations are gone. 9 0.00 0.05 0.10 0.15 0.20 0.25 0.30 t in s 0.00 0.05 0.10 0.15 0.20 0.25 0.30 v y in m /s 0.00 0.05 0.10 0.15 0.20 0.25 0.30 t in s 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 In ita l P re ss ur e Re sid ua l Figure 2.3: Pressure residual development over time compared with the terminal velocity vy of a 2D bubble with dB = 1.3 mm. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 t in s 0.90 0.95 1.00 1.05 1.10 ||F dr ag ||/ (V b g) nmin = 0 (a) No minimum number of outer iterations nmin = 0. 0.200 0.205 0.210 0.215 0.220 0.225 0.230 t in s 0.96 0.97 0.98 0.99 1.00 1.01 ||F dr ag ||/ (V b g) nmin = 0 nmin = 3 nmin = 4 nmin = 5 (b) Detailed view with varied number of minimum outer iterations, nmin. Figure 2.4: Fluctuation in the drag force for dB = 1.3 mm, c0 = 0.008 mol/m2. 10 2.3 Machine learning algorithms Within this work, machine learning techniques are used to derive a model for the local surfactant distribution on the interface of a single rising bubble. In the following, the basic machine learning algorithms and workflows are introduced. 2.3.1 Machine learning domains Machine learning can be separated into three main fields: • Supervised learning: In supervised learning the data has defined input and output parameters, the so-called features and labels. The goal of supervised learning is to train a model based on feature-label-pairs to predict the labels for unseen features. Supervised learning can be further subdivided into classification, where the output label is defined by a discrete value, e.g. true or false, and regression, where the output is continuous. • Unsupervised learning: In unsupervised learning the given data structure is unknown. The main objective is to find a structure within the data. A popular example is clustering. • Reinforcement learning: The goal of reinforcement learning is to develop a system that improves its performance based on interactions with the environment. It is commonly used, for example, in self-driving cars or humanoid robots. The modelling of the surfactant concentration on the interface is a supervised learning problem, more accurately a regression problem. Therefore, the following chapters will refer to regression problems. Additional information regarding the other machine learning fields can be found, for instance, in [26]. 2.3.2 Machine learning workflow In the following, the main steps when building and training a machine learning model are discussed. For further infor- mation, the reader is referred to [26]. The first step to generate a machine learning algorithm is the gathering of data. Within this work, the raw data is a collection of results from several DNS of rising bubbles under the influence of surfactant. From the simulations the raw data is extracted and pre-processed in several steps: • Data filtering: Data filtering is a suitable approach to reduce outliers in the data. Usually, DNS are not susceptible to outliers. The performed simulations, however, show outliers in the first time steps due to the discretisation of the surfactant transport on the interface. Therefore, a data filtering of the first time steps is performed, as described in section 6.1.1. • Feature selection: To model the surfactant concentration on the interface cΣ, several input features can be ex- tracted from the simulation, e.g. the rise velocity v or the mean surfactant concentration on the interface c̄Σ. In feature selection, the most relevant input parameters are selected to reduce the model’s complexity. • Feature scaling: The input parameters of the model can have scales ranging over orders of magnitude. While the Sherwood number Sh reaches values up to 800, the mean surfactant concentration c̄Σ is of order O (10−6). Therefore, the parameters are scaled to a range between 0 and 1 to avoid very small or very big input parameters. This scaling leads to a better-conditioned optimisation problem for the model training and thereby speeds up the training process and increases the model’s accuracy. Furthermore, feature scaling reduces rounding errors that are caused by the limited accuracy of the numerical data. • In the last step, the processed data is split into a training and a validation data set. The training data set is used to train the model, while the validation data set is unknown to the model and only used to evaluate the model performance after training. Within this work, entire simulations of different initial surfactant bulk concentrations c0 are used for validation. After data extraction and pre-processing, the machine learning model is trained. The training of the model is an iterative procedure, in which a defined criterion, the so-called loss, is minimised by adjusting the model’s weights. Finally, the model performance is assessed using the validation data. Once satisfactory results are obtained the required model can be used to make predictions on unknown input data. 11 2.3.3 Feature selection Too many input features make the model prone to overfitting, due to the curse of dimensionality [26], while too few input features instead may reduce the potential model performance. In the following, the basic concepts of the algorithms that are used to select significant input parameters are introduced. Further details about the presented techniques can be found in [26]. • Decision trees and random forest regression A random forest regressor is a machine learning algorithm that is based on multiple decision trees. A decision tree defines a series of conditions to make a prediction. Therefore, it uses the input features and creates data-based decisions to explain the output labels of the data. The feature importance of a random forest regressor is a measure for feature relevance based on the average impurity reduction computed from all decision trees within the forest. The impurity of a node is the criterion based on which the (locally) optimal condition for a split of the tree is chosen. Further information can be found in the ski-kit learn documentation1. In feature selection, decision trees are a fast approach to get a first impression of the data. The algorithm is simple to implement and has the feature importance as an internal property. The feature importance, however, is dependent on the data structure. For multiple linear correlated features, the algorithm produces poor results. Therefore, the sequential backward selection (SBS) algorithm is used for these cases. • Sequential backward selection The SBS algorithm reduces the dimensionality of an initial feature subspace by evaluation the model performance iteratively for different feature combinations. The algorithm can be divided into four steps: 1. The algorithm is initialised with the full available feature space. 2. The algorithm determines the feature that leads to the lowest performance loss of the model if removed from the feature space. 3. The feature with the lowest performance loss is removed. 4. The algorithm restarts from 2. until only one feature or a predefined minimum number of features is left. In SBS, a machine learning model is required to evaluate the performance of the predictions with the given feature subset. In general, any machine learning model may be used. Here, the K-nearest neighbour algorithm (KNN) is employed. The KNN algorithm is based on saving the entire training data set and using it to make predictions based on interpolation. The basic steps of the algorithm are: 1. Define a number of neighbouring points k and a distance metric for the interpolation. 2. Find the k nearest neighbours for the given feature vector. 3. Interpolate between the k neighbouring points to compute the label. A nice visualisation of the KNN algorithm can be found, for instance, on YouTube2. When performing SBS with the KNN algorithm, it is essential to either split the data into a training and a validation data set and use the validation data to assess the performance of the model or to use a simple average for the interpolation. Otherwise, the prediction performance can be misleading, because perfect scores are reached independently of the selected features. The KNN algorithm can be understood as a database of values that are assessed to make predictions. Thereby, the data is stored and used for estimation, instead of learning a function. When one of the input features has a value corresponding to exactly one output label, the assignment of this feature to the label is unique. This means that for every feature value only one output is defined. If the KNN algorithm is used combined with a distance metric scaling with the inverse distance as weight, the algorithm returns a perfect score, because the algorithm is not predicting, but just returning the stored values from the database. The perfect score, however, does not depict an actual functional correlation between the feature and the label but is a result of overfitting. A simple fix is the use of a validation data set to calculate the feature performance. 1 http://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html, visited 10/2018 2 https://www.youtube.com/watch?v=3lp5CmSwrHI, visited 10/2018 12 http://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html https://www.youtube.com/watch?v=3lp5CmSwrHI Finally, scoring and error metrics are required to validate the quality of the model. Popular metrics for regression are the mean squared or the mean absolute error and the r2-score1. For the SBS the r2-score is used, which is based on the mean squared error: r2-score(y, ŷ) = 1− ∑nsamples−1 i=0 (yi − ŷi)2 ∑nsamples−1 i=0 (yi − ȳ)2 , (2.44) with yi being the true label of the i-th sample, ŷi being the predicted value of the i-th sample and ȳ being the mean of all labels. The best possible r2-score is 1.0, while a value of 0.0 corresponds to a model that always returns the mean. 2.3.4 Multilayer perceptron A Multilayer perceptron (MLP) is a popular algorithm to learn non-linear relationships between input and output data [28]. A MLP is characterised by an input and an output layer combined with multiple hidden layers. The hid- den layers are densely connected and typically employ a sigmoid activation function2 [1]. The importance of the connection between two neurons is presented by weights, which are adjusted during the training process. Figure 2.5 shows the basic structure of a MLP. Figure 2.5: Structure of a multilayer perceptron [21]. Formally, a MLP is defined as a function with a weight matrix Wk and a bias bk for each layer k. The forward pass through a MLP with one hidden layer and an input vector x ∈ Rn and the output y ∈ R is computed as: y(x) =W2 · f (W1 · x+ b1) + b2 , (2.45) where W1 ∈ Rm and W2,b1,b2 ∈ R are the model parameters. The weights of the input and hidden layer are represented by W1 and W2, respectively. The biases b1 and b2 are added onto the hidden layer and the output layer, respectively. The function f (·) : R→ R is the, so-called, activation function of the neuron. The backpropagation algorithm is used to adjust the network’s weights and biases for optimal representation of the training data. Optimal network parameters are those which lead to a minimum loss, defined, for example, as the mean squared error: E = 1 N N ∑ i=1 (yi − ŷi) 2 , (2.46) where yi is the actual label, ŷi is the network estimate and N the number of data points fed to the network. 1 http://scikit-learn.org/stable/modules/model_evaluation.html, visited 10/2018 2 A sigmoid function is a bounded, differentiable, real function that is defined for all real input values and has a non-negative derivative at each point. 13 http://scikit-learn.org/stable/modules/model_evaluation.html The incremental update of model parameters via backpropagation includes the following three steps: 1. The network is fed with data, and the output of every neuron in each layer is computed. 2. The network output is compared with the desired labels and the difference is considered as the network’s error. 3. The error is ’backpropagated’ to the input layer, starting with the output layer. Thereby, the weights of the network are adjusted according to their impact on the error. The update of a weight is calculated by: w(n)ri j = w(n−1) ri j −λ · ∂ E ∂ w(n−1) ri j , (2.47) where wri j is the weight from the input neuron j to the hidden neuron i and r indicating the layer, E being the loss of the prediction, n the number of iterations and λ the step width or learning rate. In chapter 6 a MLP is implemented to model the local surfactant distribution on the interface. For the implementation and training, the open-source software library tensorflow1 is used, which was originally developed by Google and enables high-performance numerical computations of neural networks and other machine learning algorithms. 1 https://www.tensorflow.org/, visited 10/2018 14 https://www.tensorflow.org/ 3 Study on the mesh requirements The simulation of rising bubbles under the influence of surfactants has high demands regarding the mesh resolution at and close to the interface. In a Finite Volume/Finite Area framework methodology, the computational domain is divided into control volumes in the bulk and control areas on the interface. In order to determine the resolution requirements close to the interface in the radial and tangential direction, a 2D mesh study is carried out. The fluid properties for these simulations are chosen according to the experimental conditions reported in the study by Tagawa et al. [31]. The liquid phase consists of water, while the gaseous phase is nitrogen at T = 295.65 K. Table 3.1: Fluid properties of the gaseous and liquid phase [35]. ρL kg/m3 ρG kg/m3 µL kg/(ms) µG kg/(ms) σ0 N/m 997.7 1.140 0.943 · 10−3 1.77 · 10−5 0.07235 3.1 Mesh study for a clean 2D bubble A 2D study of a clean nitrogen bubble rising in water with diameters of 0.8, 1.3 and 2.0 mm is carried out to determine the mesh resolution requirements. First, the grid convergence in the radial direction of the outer and inner mesh is investigated using the total pressure force acting on the interface as a metric. The total pressure can be divided into the dynamic, ambient and hydrostatic pressure, see section 2.1.4. For a freely rising bubble, the ambient pressure difference is zero, and the hydrostatic pressure force is approximately constant over time. Therefore, the change in total pressure force is dominated by the dynamic pressure, which is the kinetic energy per volume of fluid and thereby directly related to the rise velocity of the bubble. Compared to the terminal velocity, however, it shows pronounced deviations. Secondly, the tangential mesh quality is evaluated using the jump conditions at the interface decomposed in a normal and a tangential contribution. For an uncontaminated bubble the surface tension σ at the interface is constant and the jump condition in equation (2.5) simplifies to: JptotI− SviscK · nΣ = σκnΣ , (3.1) For a closed surface, the right hand-side of equation (3.1) cancels out: Fca = ∫ A σκnΣdA= 0 , (3.2) leading to the normal jump condition at the interface: JFptot − Fvisc ⊥ K= 0 . (3.3) An insufficient tangential resolution results in numerical errors in the computation of the surface integrals of equa- tion (3.1). Table 3.2 shows the default settings for the meshes. If not stated otherwise, these are used in all simulations. Figure 3.1 displays the effect of the parameters on the mesh. The meshes are created using the open-source software salome-meca1. Table 3.2: Default setup for the 2D mesh study. NTan NRin NRout 180 36 28 1 https://www.code-aster.org/spip.php?article303, visited 10/2018 15 https://www.code-aster.org/spip.php?article303 NRin NRout NTan Figure 3.1: Exemplary computational grid of a 2D bubble showing the different mesh parameters. NTan defines the tan- gential grid resolution, while NRin and NRout define the radial resolution in the gaseous and in the liquid do- main, respectively. 3.1.1 Radial grid resolution: liquid phase To resolve the hydrodynamic boundary layer, the size of the first cells close to the interface is significant. Therefore, different cell sizes at the interface have been compared, varying only the number of cells in radial direction NRout, stretching the cell size at the interface by a factor of 1.5. The case setups are shown in table 3.3. The radial cell width at the interface ∆Rout,1 is normalised using the hydrodynamic boundary layer thickness δh. The correlation presented by Tomiyama et al. [32] provides an estimate of the Reynolds number for a clean bubble. The boundary layer thickness is then given by: δh ≈ dBp 2Re . (3.4) Figure 3.2 shows the total pressure force acting on the interface normalised with the buoyancy force VB∆ρg. Other parameters, e.g. the terminal velocity vy , show a comparable behavior. For dB = 2.0 mm the 2D simulation gets unstable after t ≈ 0.1 s. Therefore, the simulation results for t < 0.1 s are used to evaluate the radial grid resolution. With increasing mesh resolution the drag force acting on the bubble is computed more accurately, leading to an higher terminal velocity of the bubble. Mesh independent results are obtained for δh/∆Rout,1 > 5.5 and δh/∆Rout,1 > 3 for the small bubble with dB = 0.8 mm and for the bigger bubbles with dB = 1.3 mm and dB = 2.0 mm, respectively. The corresponding number of radial cells NRout and the cell width at the interface ∆Rout,1 are displayed in table 3.3. These values provide a reference for the minimum radial grid resolution in 3D. Table 3.3: Parameters for the outer radial grid resolution study. (a) dB = 0.8 mm; δh = 45 µm. NRout ∆Rout,1 µm δh/∆Rout,1 8 19.7 2.3 12 13.5 3.3 20 8.2 5.5 28 5.9 7.6 (b) dB = 1.3 mm; δh = 42 µm. NRout ∆Rout,1 µm δh/∆Rout,1 8 31.1 1.4 12 21.8 1.9 20 13.4 3.1 28 9.6 4.4 (c) dB = 2.0 mm; δh = 57 µm. NRout ∆Rout,1 µm δh/∆Rout,1 8 49.4 1.2 12 33.6 1.7 20 20 2.9 28 14.7 3.9 16 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 t in s 0.02 0.04 0.06 0.08 ||F p t ot ||/ (V b g) h/ Rout, 1 = 2.3 h/ Rout, 1 = 3.3 h/ Rout, 1 = 5.5 h/ Rout, 1 = 7.6 (a) dB = 0.8 mm. 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 t in s 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 ||F p t ot ||/ (V b g) h/ Rout, 1 = 1.4 h/ Rout, 1 = 1.9 h/ Rout, 1 = 3.1 h/ Rout, 1 = 4.4 (b) dB = 1.3 mm. 0.00 0.02 0.04 0.06 0.08 0.10 t in s 0.000 0.002 0.004 0.006 0.008 0.010 ||F p t ot ||/ (V b g) h/ Rout, 1 = 1.2 h/ Rout, 1 = 1.7 h/ Rout, 1 = 2.9 h/ Rout, 1 = 3.9 (c) dB = 2.0 mm. Figure 3.2: Total pressure force acting on the interface for different outer radial resolutions, δh/∆Rout,1. 17 3.1.2 Radial grid resolution: gaseous phase By varying the number of radial cells in the inner domain NRin, the resolution of the inner mesh (gas phase) is investigated. A bubble with a diameter of dB = 1.3 mm is considered. The simulations for the other bubble diameter show similar results. Table 3.4 lists the parameters used in this study. Table 3.4: Parameters for the inner radial grid resolution study, dB = 1.3 mm. NRin ∆Rin,1 µm ∆Rout,1 µm ∆Rin,1/∆Rout,1 12 37.4 9.6 3.9 20 23.1 9.6 2.4 28 16.6 9.6 1.7 36 13.0 9.6 1.4 48 9.8 9.6 1.02 The influence of the inner mesh resolution is small compared to the outer one. The deviation in the total pressure force amplitude between the default setting with ∆Rin,1/∆Rout,1 = 1.7 and the finest mesh ∆Rin,1/∆Rout,1 = 1 is smaller than 0.2%. Nevertheless, a too coarse inner mesh in relation to the outer one leads to deviations of about 1%, see ∆Rin,1/∆Rout,1 = 3.9 in figure 3.3. For the 3D simulation a ratio ∆Rin,1/∆Rout,1 < 2 is recommended. 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14 0.15 t in s 0.0320 0.0325 0.0330 0.0335 0.0340 0.0345 0.0350 0.0355 0.0360 ||F p t ot ||/ (V b g) Rin, 1/ Rout, 1 = 3.9 Rin, 1/ Rout, 1 = 2.4 Rin, 1/ Rout, 1 = 1.7 Rin, 1/ Rout, 1 = 1.4 Rin, 1/ Rout, 1 = 1 Figure 3.3: Total pressure force at the interface, dB = 1.3 mm, for different ratios∆Rin,1/∆Rout,1. 18 3.1.3 Tangential grid resolution In the following, the fulfilment of the normal force balance at the interface is examined. Therefore, the normal contribu- tion of the global forces acting on the interface in equation (2.32), derived from the normal jump condition, is analysed, in other words, the difference between the total pressure and the normal viscous force normalised with the total pres- sure force. For dB = 0.8 mm, this difference is approximately 10−10 and thereby smaller than the solver tolerance. For dB = 1.3 mm the difference between the total pressure and normal viscous force increases with the mesh resolution. It reaches a maximum of about 1% at t = 0.15 s. For dB = 2.0 mm the force difference at t = 0.1 s is greater than 10% of the total pressure force, see figure 3.4. In the following, only the results for dB = 2.0 mm are further examined. 0.00 0.02 0.04 0.06 0.08 0.10 t in s 0.00 0.02 0.04 0.06 0.08 0.10 0.12 ||F p t ot Fvi sc ||/ ||F p t ot || h/ Rout, 1 = 1.2 h/ Rout, 1 = 1.7 h/ Rout, 1 = 2.9 h/ Rout, 1 = 3.9 Figure 3.4: Force balance normal to the interface, dB = 2.0 mm, with different radial resolutions, δh/∆Rout,1. The increasing deviation with mesh resolution can be explained by the greater bubble deformation with higher terminal velocities and bubble sizes. Due to the bubble deformation, the interfacial mesh is deformed. As a result, the tangential resolution is not sufficient in the region of highest curvature, leading to different bubble shapes with respect to the various δh/∆Rout,1, see figure 3.5. δh/∆Rout,1 = 1.2 δh/∆Rout,1 = 1.7 δh/∆Rout,1 = 2.9 δh/∆Rout,1 = 3.9 δh/∆Rout,1 ↑ Figure 3.5: Bubble deformation for different radial resolutions δh/∆Rout,1 at t = 0.1 s, dB = 2.0 mm. The increased deformation results in an increased error in the normal jump condition. It is likely that with increasing deformation the regions of higher curvature are not resolved. Therefore, the principle that the total surface tension force on a closed surface cancels out is not fulfiled, and neither is the normal jump condition, see equation (3.2). To determine the required tangential resolution the number of cells in the tangential direction NTan is varied. Table 3.5 and figure 3.6 show the comparison of different tangential resolutions. As suspected, an increase of the tangential resolution decreases the error for the normal jump condition at the interface. For dB/∆Tan = 131 the deviation is less than 1% of the total pressure force after t = 0.06 s. For a 3D mesh, this resolution leads to a large number of faces at the interface. Therefore, a tangential grading is suggested, increasing the number of cells close to the region of maximum curvature for bubbles with a high deformation. 19 Table 3.5: Parameters for the tangential grid resolution study, dB = 2.0 mm. NTan ∆Tan µm dB/∆Tan 140 45.5 44 180 35.4 56 220 29.0 69 280 21.3 94 380 16.8 119 420 15.3 131 0.00 0.02 0.04 0.06 0.08 0.10 t in s 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 ||F p t ot Fvi sc ||/ ||F p t ot || dB/ Tan = 44 dB/ Tan = 56 dB/ Tan=69 dB/ Tan = 94 dB/ Tan = 119 dB/ Tan = 131 Figure 3.6: Force balance normal to the interface, dB = 2.0 mm and δh/∆Rout,1 = 3.9, with different tangential resolutions, dB/∆Tan. 20 3.2 3D meshes The initial shape of the bubble is a sphere positioned in the centre of a spherical domain. The outer domain radius is twenty times larger than the radius of the bubble. The computational grid is divided into two meshes, one for the gaseous phase consisting of polyhedral cells and one for the liquid phase consisting of prismatic cells with a polygonal base, see figure 3.7. In the liquid domain, the outer boundary mesh (space) is created by projecting the interface grid onto the outer surface. (a) Full domain. (b) Enlarged view of the bubble region. Figure 3.7: Computational mesh for the 3D simulation. The inner and outer domain are displayed, and a part of the interface mesh is shown in a darker color. The radial cell thickness of the first cell at the interface ∆Rout,1 is chosen according to the outcome of the parameter study for the radial grid resolution, see section 3.1. The criteria for the radial mesh resolution in the liquid domain are ∆Rout,1 < 9 µm, ∆Rout,1 < 14 µm and ∆Rout,1 < 20 µm for the 0.8, 1.3 and 2.0 mm bubble, respectively. The radial cell thickness in the gaseous phase is chosen with respect to the one in the outer domain and satisfies ∆Rin,1/∆Rout,1 < 2. Regarding the tangential resolution, a uniform mesh fulfiling the requirements set by the tangential grid resolution study results in a large number of faces on the interface. Such a grid would imply high computational costs. Thus, the solution adopted here consists of a division of the interface in three regions: top, equator and bottom, see figure 3.8. For these regions, different tangential resolutions are employed: (i) In the inflow region (top domain) the resolution requirements are lower than for the other two regions. Therefore, a coarser grid is used to lower the computational costs. (ii) In the equator region, high local curvatures can occur for bubbles with high deformation, resulting in increased resolution requirements, see chapter 3.1.3. Thus, the resolution is adapted according to the expected bubble deformation. (iii) The bubble wake is located in the bottom region of the bubble. The tangential grading at the interface is reflected in the liquid domain, thus, a finer grid in this region enables a higher resolution of the wake compared to the inflow region. equator bottom top Figure 3.8: Interface mesh of a 2.0 mm bubble undeformed (left) and deformed at t = 0.08 s (right). 21 The average edge length of the faces in the top, equator and bottom regions are chosen such that: (i) the total number of faces is limited to a prescribed number NΣ, (ii) there is a smooth transition between the face sizes in the different areas and (iii) the tangential grid resolution requirement is satisfied as well as possible. Figure 3.8 shows the mesh of a bubble with dB = 2.0 mm and tangential grading in its undeformed and deformed state after t = 0.08 s, before the bubble starts oscillating. For dB = 0.8 mm, only little deformation is expected. Therefore, a uniform mesh is used. The simulations are carried out using a coarse and a fine mesh. Table 3.6 lists the approximate resulting face and cell sizes of the meshes. Table 3.6: Grid resolution and mesh size for dB = 0.8 mm. Mesh ∆Tan µm ∆Rout µm NΣ NTotal coarse 30 8.9 2.526 193.878 fine 25 6.7 3.604 350.328 For the 1.3 mm and 2.0 mm bubble, a grading between the top, equator and bottom region is used. As for the smaller bubble, the simulations are executed on a coarse and a fine mesh. Additionally, for dB = 1.3 mm two contaminated cases are considered. For these cases, the flow field is strongly affected by the presence of surfactants with vortices developing in the bubble wake. Thus, a mesh with a refined bottom region and initial surfactant bulk concentrations of c0 = 4.4 · 10−4 mol/m3 and c0 = 0.002 mol/m3 is examined in order to asses mesh independence, see section 5.2. Tables 3.7 and 3.8 list the resulting mesh statistics for dB = 1.3 and 2.0 mm. Table 3.7: Grid resolution and mesh size for dB = 1.3 mm. Mesh ∆Tan,Eq µm ∆Tan,Bot µm ∆Tan,Top µm ∆Rout µm NΣ NTotal coarse 30 45 60 13.9 2.563 197.678 medium 30 40 50 9.5 3.631 352.394 fine 30 25 50 7.6 4.839 517.607 Table 3.8: Grid resolution and mesh size for dB = 2.0 mm. Mesh ∆Tan,Eq µm ∆Tan,Bot µm ∆Tan,Top µm ∆Rout µm NΣ NTotal coarse 35 60 90 13 3.620 352.584 fine 30 40 80 13 6.092 591.870 3.3 Mesh requirements: conclusion In this chapter, the mesh resolution requirements at the interface in radial and tangential direction have been assessed. These results serve as reference values for the creation of the 3D meshes used for the following simulations. 22 4 Parallelisation study The preformed 3D simulations require high computational efforts due to the resolution requirements at the interface, see chapter 3. Additionally, in most of the contaminated cases relatively long physical times of approximately 1 s have to be simulated to reach a quasi-steady state. Due to the time step criterion, see equation (2.41), a total number of about 5 · 105 time steps is necessary to reach such a physical time of 1 s. Similar simulations that have been carried out in earlier studies required runtimes of 30 to 60 days while distributed on six processors [25]. To reduce the runtime, two different decomposition methods are examined. Additionally, various numbers of processors are considered to choose the most efficient way to distribute the computational effort. For the parallelisation study, an uncontaminated 3D bubble with a diameter of dB = 1.4 mm and a uniform mesh with NΣ ≈ 4500 faces on the interface and NTotal ≈ 350000 cells is considered. The simulation is executed for 1000 time steps with a fixed number of five outer iterations. Due to the different performances of the assigned processors on the Lichtenberg cluster, the simulation time can differ. Nevertheless, the study should be able to provide trends for the runtime. Within the interface tracking methodology, the interface (and its counterpart bounding the gas phase) cannot be distributed over multiple processors. Thus, the decomposition method has to account for this restriction. The decomposition methods examined are • a manual decomposition, with processor regions defined by the user [36]. • a scotch decomposition, where the processor regions are chosen automatically by the scotch library [6]. For the scotch decomposition, a speed up in runtime is desired since no restriction on the number of subdomains is prescribed. Additionally, the scotch decomposition requires less user intervention and is better suited for non-uniform meshes, see section 4.2. 4.1 Decomposition techniques The manual decomposition divides the computational domain in an inner part containing the interface and some cells around it, the so-called outer domain. The inner part can be further decomposed into a region containing the interface (like a hollow sphere) and the remaining part. The outer domain can be split into two, four or eight subdomains, each one cutting the volume in direction of the xy-, xz- and yz-planes. Figure 4.1 shows a cut through the domain with the regions coloured by the processor number after a manual decomposition on ten processors. Figure 4.1: Cut through the yz-plane of a manual decomposition on ten processors. The processor domain containing the interface is marked in red. 23 In the case of the scotch decomposition, a restriction needs to be specified to ensure that the interface is not shared by multiple processors. Therefore, a small hollow sphere containing the interface is defined to stay on a single processor, see figure 4.2. The hollow sphere defines the smallest domain that needs to be included in the single processor region containing the interface. Nevertheless the scotch decomposition algorithm can include additional cells to that region. The rest of the domain is decomposed automatically to optimise the number of cells per processor and the number of faces per processor boundary according to the scotch algorithm [6]. Figure 4.3 shows a cut through the yz-plane for the scotch decomposition on 10 processors. The processor containing the interface is coloured in red. interface rRin rRout Figure 4.2: Schematic drawing of the restriction applied to the scotch decomposition for the processor domain holding the interface. Figure 4.3: Cut through the yz-plane of the scotch decomposition on ten processors. The processor holding the interface domain is marked in red. 24 4.2 Comparison of the decomposition techniques To compare the decomposition techniques, the scotch decomposition is used with a hollow sphere containing the interface with radii of rRout = 0.72 mm and rRin = 0.69 mm and no further restrictions. For the manual decomposition the number of cells in the processor domains are evenly distributed. Table 4.1 shows the cell distribution for varying numbers of processors nProc. Table 4.1: Number of cells on the processor containing the interface NProc,Σ compared with the mean number of cells on the other processor regions NProc,/∈Σ for a varying number of processors, nProc. (a) Manual decomposition Method nProc NProc,Σ NProc,/∈Σ χProc manual 2 164000 177000 1.08 manual 5 66000 70000 1.06 manual 6 62000 65000 1.05 manual 10 31000 37000 1.2 (b) Scotch decomposition Method nProc NProc,Σ NProc,/∈Σ χProc scotch 2 170000 171000 1.01 scotch 5 86000 59000 0.7 scotch 6 66000 57000 0.86 scotch 10 43000 34000 0.79 scotch 12 19000 28500 1.5 scotch 14 18000 24000 1.3 Figure 4.4 compares the runtime of the two decomposition methods on a different number of processors scaled by the runtime on a single core. As a reference, a half linear scaling of y = 1/2 · x is shown in orange. Both decomposition techniques perform slightly worse than the given reference and do not differ significantly for more than five processors. For more than 10 processors, the gain in runtime reduces even further and the runtime stays approximately constant. 2 4 6 8 10 12 14 nProc 1 2 3 4 5 6 7 t ru n, n P ro c= 1/t ru n y = 1/2 x Scotch Manual Figure 4.4: Runtime comparison between the manual and scotch decomposition. For uniform meshes the manual decomposition gives overall better results compared to the scotch decomposition. Nevertheless, the processor domains need to be specified by the user and are restricted to a decomposition in xy-, xz-, and yz-direction, see section 4.1. In case of a non-uniform mesh with top, equator and bottom grading, as used for the 1.3 mm and 2.0 mm bubble, this decomposition leads to an uneven cell distribution on the processor domains not containing the interface. The number of cells in these domains differs up to a factor of 8. Therefore, the performance of the manual decomposition decreases drastically. The scotch decomposition, on the other hand, shows comparable runtime results due to an even cell distribution. Furthermore, the decomposition can be done on any number of processors giving more flexibility to the user. 25 4.3 Domain dependency study The area near the interface region has the highest requirements regarding the computational effort. Therefore, the computational time on the processor containing the interface is expected to be greater than on the other ones with the same amount of cells. A further study is performed to gain an insight into the influence of the cell ratio on the runtime: χProc = NProc,/∈Σ NProc,Σ , (4.1) with NProc,/∈Σ being the mean number of cells in the processor domains not containing the interface and NProc,Σ being the number of cells on the processor containing the interface. Cell ratios χProc of 1, 1.5 and 2 are examined for the manual and scotch decomposition on 2, 5, 6 and 10 processors. For the scotch decomposition, two different approaches to control the number of cells on the interface are presented. 4.3.1 Manual decomposition For the manual decomposition, the outer radius of the processor domain containing the interface is chosen in order to fulfil the prescribed χProc. Figure 4.5 displays the obtained results. For less than five processors, the decomposition with χProc = 2 shows the best results. The cases of nProc = 5 and 6, χProc = 1 and 1.5 do not differ significantly from one another. In case of nProc = 10, the performance does not decrease continuously with χProc. In fact, χProc = 1.5 shows the worst performance, indicating counteracting effects on the runtime: • A higher χProc increases the decomposition efficiency. This effect can be explained due to the increased complexity of the physical processes on and close to the interface and thereby increased computational effort to solve the numerical problem. An indication for this hypothesis is the overall better performance for χProc = 2 on nProc > 5. • Cutting the processor domain too close to the interface leads to badly conditioned matrices, slowing down the iterative solver and thereby increasing the time to solve the numerical problem on the processor containing the interface. This hypothesis is based on the performance increase for χProc = 1 compared to χProc > 1 for a decom- position on 10 processors. The higher ratios lead to a processor domain close to the interface consisting of only five or three radial cell layers outside the interface in the processor domain for χProc = 1.5 and 2, respectively. In order to verify these hypotheses and to obtain more meaningful results, further decomposition studies are needed, which would exceed the scope of this work. 2 3 4 5 6 7 8 9 10 nProc 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 t ru n, n P ro c= 1/t ru n y = 1/2 x Manual: Proc = 1 Manual: Proc = 1.5 Manual: Proc = 2 Figure 4.5: Runtime comparison for different cell ratios χProc for the manual decomposition. 26 4.3.2 Scotch decomposition To control the processor region containing the interface in the scotch decomposition, processor weights wProc can be defined that control the cell ratio χProc. Inside the scotch algorithm higher processor weights lead to an increasing number of cells on these processors. For the scotch decomposition, two different approaches are compared. 1. The radii of the hollow sphere containing the interface is kept constant while the processors outside the interface are weighted according to χProc. 2. The processor containing the interface is forced to contain only the hollow sphere by applying high processor weights for the domains not containing the interface, while the outer and inner radius of the hollow sphere rRin and rRout are varied. Table 4.2 shows the different parameters for the scotch decomposition while figures 4.6 and 4.7 show the results in terms of runtime for the different domain decomposition strategies. The first decomposition technique does not show a recognisable trend according to χProc and the performance is comparable to the unweighted decomposition. The second strategy shows slightly better results with a maximum performance for χProc = 1 and 2. On 10 processors a maximum speed up factor of 4.4 is reached. Table 4.2: Parameters used for the scotch decomposition domain study. A value of 0 for the ratio between the processor weights corresponds to a processor containing the interface consisting only of the volume inside the hollow sphere. (a) Varying the processor weight. nProc χProc wProc,Σ wProc,/∈Σ rRin mm rRout mm 2 1 6/7 0.72 0.69 5 1 3/5 0.72 0.69 6 1 4/7 0.72 0.69 10 1 1/3 0.72 0.69 2 1.5 6/9 0.72 0.69 5 1.5 4/11 0.72 0.69 6 1.5 2/4 0.72 0.69 10 1.5 1/15 0.72 0.69 2 2 2/5 0.72 0.69 5 2 1/5 0.72 0.69 6 2 1/7 0.72 0.69 10 2 0 0.72 0.69 (b) Varying the hollow sphere containing the interface. nProc χProc wProc,Σ wProc,/∈Σ rRin mm rRout mm 2 1 0 2.0 0 5 1 0 0.85 0 6 1 0 0.88 0.69 10 1 0 0.77 0.69 2 1.5 0 1.4 0 5 1.5 0 0.75 0 6 1.5 0 0.8 0.69 10 1.5 0 0.73 0.69 2 2 0 1.2 0 5 2 0 0.714 0 6 2 0 0.75 0.69 10 2 0 0.714 0.69 2 3 4 5 6 7 8 9 10 nProc 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 t ru n, n P ro c= 1/t ru n y = 1/2 x No Weights Weight: Proc = 1 Weight: Proc = 1.5 Weight: Proc = 2 Figure 4.6: Runtime comparison for different cell ratios χProc with scotch decomposition and varying processor weights. 27 2 3 4 5 6 7 8 9 10 nProc 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 t ru n, n P ro c= 1/t ru n y=1/2*x No Weights Domain: Proc = 1 Domain: Proc = 1.5 Domain: Proc = 2 Figure 4.7: Runtime comparison for different cell ratios χProc with scotch decomposition and varying interface domain. 4.4 Parallelisation study: conclusion The parallelisation technique using the OpenFOAM scotch algorithm is not able to improve the runtime for uniform meshes. On the contrary, the manual decomposition shows slightly better performance. Even though the desired per- formance boost on uniform meshes was not reached, the less restrictive scotch decomposition enables a satisfactory decomposition of non-uniform meshes, which was not possible using the manual decomposition. Furthermore, the scotch decomposition provides higher flexibility and requires less user intervention. In the current work, a combination of both techniques is used. While for uniform meshes the manual decomposition on 10 processors with a ratio of χProc ≈ 2 is used, for non-uniform meshes the scotch decomposition on 10 processors with a sharp processor domain holding the interface (2nd decomposition method presented in section 4.3.2) and a ratio of χProc ≈ 2 is favoured. Due to the restriction that the interface needs to stay on one processor, a noticeable performance increase for fur- ther decomposition strategies is unlikely. This restriction limits the capability of the decomposition techniques included in OpenFOAM and also the maximum number of processors that can be used. The possibility to cut the interface in multiple parts allows new decomposition approaches and promises performance improvements. Nevertheless, such an implementation exceeds the scope of this work. 28 5 Simulation results and discussion The objective of this work is to simulate single rising nitrogen bubbles in clean and contaminated water. The surfactant employed is Triton-X100, which is known to follow a fast sorption mechanism [5]. The surfactant is dissolved in the liquid phase and can ad- and desorb onto the interface via the sorption mechanisms. To reduce the resolution requirements close to the interface, a SGS model for the surfactant transport in the bulk phase is used [37, 25]. The sorption process is modelled using the Langmuir model with the respective adsorption isotherm and surface tension equation of state. Four different initial surfactant bulk concentrations are considered: a relatively small one c0 = 4.4·10−4 mol/m3 corresponding to an experimental study by Tagawa et al. [31], two intermediate ones, c0 = 0.002 and 0.008 mol/m3, and a high one c0 = 0.05 mol/m3. Table 5.1 lists the surfactant specific properties. Table 5.1: Surfactant properties (Triton-X100). cΣ∞ mol/m2 a mol/m3 D m2/s DΣ mol/m2 T K 2.9 · 10−6 6.6 · 10−4 2.6 · 10−10 2.6 · 10−7 315.65 The simulation setup is described in section 3.2, while the simulation results are discussed in the following. To assess their quality, experiments examining the path instability of rising bubbles under the influence of Triton-X100 performed by Tagawa et al. [31] are considered as a reference. 5.1 Experimental studies The experimental study performed by Tagawa et al. [31] examines the path of single nitrogen bubbles rising in water under the influence of Triton-X100. Thereby, bubble diameters of 1.3, 2.0 and 3.1 mm are studied with a surfactant bulk concentration of 27 ppm, which corresponds to a molar concentration of 4.4 · 10−4 mol/m3. Figure 5.1 shows the experimentally obtained trajectories. –5 –5 0 0 5 5 (a) 800 600 400 200 1000 0 z( m m ) y (mm) x (mm) –10 –10 –10–10 0 0 10 10 10 10 (b) y (mm) x (mm) 0 0 (c) 800 600 400 200 1000 0 x (mm) y (mm) 800 600 400 200 1000 0 Zigzag Helical Zigzag Helical Straight Zigzag Zigzag Helical Figure 5.1: Trajectories of rising bubbles under the influence of Triton X-100 [31]: (a) dB = 1.3 mm, (b) dB = 2.0 mm, (c) dB = 3.1 mm. The bubble with dB = 1.3 mm shows a straight–helical–zig-zag transition. The onset of helical motion is after approx- imately 200 mm rise height. The bubble follows a helical trajectory with a constant diameter before it transitions into a zig-zag motion. For dB = 2.0 mm, first a helical rise is observed with a decreasing helix diameter that becomes a zig-zag motion with a constant amplitude. The biggest bubble shows a zig-zag–helical–zig-zag transition. The first zig-zag state shows an increasing motion amplitude, followed by a helical motion with decreasing helix diameter. The zig-zag motion after approximately 400 mm rise height shows a constant amplitude. Due to the graphical representation, it is difficult to obtain exact values for the motion amplitudes. Nevertheless, the experimental data can be used for a qualitative comparison and discussion of the numerical results. In this work, the bubble diameter dB = 3.1 mm is not considered, because of the high computational effort. 29 5.2 Mesh sensitivity study In order to assess the mesh sensitivity, the bubble terminal velocity vy and the bubble path are examined on a coarse and a fine mesh. Tables 3.6 to 3.8 list the corresponding mesh setups. 5.2.1 Mesh sensitivity for dB = 0.8 mm Figure 5.2 shows the terminal velocity under the influence of different initial surfactant bulk concentrations for the smallest bubble investigated. The obtained results on both meshes do not differ significantly. The maximum difference between the coarse and fine mesh is encountered for the most contaminated bubble and is smaller than 0.1% of the terminal velocity. Due to