A stabilized finite element method using a discontinuous level set approach for the computation of bubble dynamics

Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method

The aim of this work is to provide a solver for viscoelastic multi-phase flows within the Bounded Support Spectral Solver (BoSSS) currently under development at the Chair of Fluid Dynamics at the Technical University of Darmstadt. The discretisation in BoSSS consists of a high-order Discontinuous Galerkin (DG) method for single-phase flow and a high-order eXtended Discontinuous Galerkin (XDG) method for the multi-phase purpose. The solver shall be used to investigate numerically the behaviour of viscoelastic droplets. The macroscopic Oldroyd B model which is used in a wide range of applications is chosen as the constitutive model. A detailed derivation of the system of equations including the modeling principles for the Oldroyd B model is presented. A DG discretisation of the system of equations including the Local Discontinuous Galerkin (LDG) method is presented after introducing the field of the DG method. The derivation of appropriate flux functions for the constitutive equations and the extra stress tensor are one of the key derivations of this scientific work. Difficulties arising in the numerical solution of viscoelastic flow problems for higher Weissenberg numbers for different discretisation methods are due to the convection dominated, mixed hyperbolic-elliptic-parabolic nature of the system of equations. Several strategies are presented which overcome these problems and are known from the literature. A key achievement of this scientific work is the application of the LDG method, originally developed for a hyperbolic system of equations for a Newtonian fluid, on the viscoelastic system of equations which renders methods for preserving ellipticity unnecessary. Furthermore, various strategiesn to enhance and to support convergence of the solution of the DG discretised system are presented. These are the Newton method with different approaches determining the Jacobian of the system, a homotopy continuation method based on the Weissenberg number for a better initial guess for the Newton method, and a troubled cell indicator combined with an artificial diffusion approach or an adaptive mesh refinement strategy, respectively. For the completeness of this work the XDG method is presented using a sharp interface approach with a signed distance level-set function as it is implemented in BoSSS. The single-phase solver is combined with these methods and appropriate flux functions for the interface are implemented to enable multi-phase applications for viscoelastic fluid. Several numerical experiments are conducted to verify and to validate the viscoelastic singlephase solver and to show the capability of the viscoelastic multi-phase solver to simulate viscoelastic droplets. Advantages and disadvantages of the implementation and an outlook for future research can be found in the conclusion.

A VOSET method combined with IDEAL algorithm for 3D two-phase flows with large density and viscosity ratio

A simplified method is proposed to implement a wetting boundary condition on curved surfaces within the conservative phase-field lattice-Boltzmann (LB) simulation framework. It combines the idea of Huang et al. [“An alternative method to implement contact angle boundary condition and its application in hybrid lattice-Boltzmann finite-difference simulations of two-phase flows with immersed surfaces,” Eur. Phys. J. E 41, 17 (2018)] to find the order parameter on the other side of the wall with the conservative Allen–Cahn equation (CACE) for interface evolution solved by the LB equations. It inherits the advantage of the original method using the Cahn–Hilliard equation to avoid complicated interpolations under different geometries. By using the CACE, the boundary condition for the chemical potential is circumvented (making it more simplified), and the overshooting of the order parameter is also greatly suppressed, enabling it to simulate two-phase flows with solid objects of various shapes and wettabilities at large density and viscosity ratios. Several two-dimensional, axisymmetric, and three-dimensional problems, including some previously studied by experiments, were simulated and the capability of the proposed method is demonstrated through its good agreement with theoretical predictions and/or experimental observations.

A level set discontinuous Galerkin method for free surface flows

We present a new numerical method to simulate the structural patterns emerging from the long-term large-deformation tectonic flows in both two and three spatial dimensions.  The domains of different material properties are each represented by a level set function discretized on a Eulerian mesh with the discontinuous Galerkin method. The level sets are advected by a velocity field provided by a coupled Stokes flow solver. Our method accurately captures the material interface by the adaptive mesh refinement, reduces the computational expenses compared to the traditional particle-in-cell method and offers straightforward handling of geometric splitting and merging.  Under the unified finite element framework, our method promises the flexibility in the choice of mesh geometry as well as the potential for extending to complex rheology.  With passive tracers geat and around areas of interest, the finite strain of the flow field can be integrated through any time interval within the total simulation time.  The strain ellipsoids thus obtained offers the possibility for ground-truthing the simulated deformation patterns with the field structural analysis.  Our results demonstrate identical physical behaviour when compared with established structural geology and geodynamic benchmarks.The style of the crustal dynamics on the Archean Earth has been subject to controversy on whether a vertical tectonic style in the form of Rayleigh-Taylor instability, induced by an inverted density profile, prevails in the early history of the Earth and if so, how the transition to the present-day plate tectonics, characterized by dominantly horizontal movement, is manifested in the rock record.  Equipped with our modelling scheme, we construct numerical models to simulate the lithological distributions and deformation patterns resulted from a synchronous operation of vertical tectonism and horizontal shearing. The latter can be viewed as a possible result of some far-field tectonic boundary condition (e.g. oblique convergence).  Many aspects of the simulation in terms of the map pattern, foliation/lineation trend and strain distribution compare favorably with the field observations in Neoarchean granitoid-greenstone terranes in the Superior Province as well as worldwide.  Therefore, it is concluded that the vertical and horizontal tectonism are not mutually exclusive tectonic regimes  The symbiosis of both tectonic processes is a viable mechanism for establishing the crustal architecture and the deformation pattern we see today in many Neoarchean terranes and might represent a transition from the former to the latter in the Neoarchean.

Many multiphase flows relevant to natural phenomena and technical applications are characterized by large density ratios between the phases or fluids. Numerical simulations of such flows are especially challenging if the phase interface has a complex shape and is subject to large shear. This scenario is typical for atomization of liquids under ambient conditions. In this paper we discuss some of the reasons why one-fluid level-set based methods are prone to become unstable for high-density ratio/high shear atomizing flows and present a consistent rescaled momentum transport method that addresses the identified shortcomings. We present results obtained with the new method for a number of high-density ratio test cases, including the advection of a 1 000 000:1 density ratio impulsively accelerated drop, a 1000:1 density ratio damped surface wave, and the collapse of a water column in air under ambient conditions. The new method shows significantly improved results compared to standard level set based single-fluid methods.

In this paper, an integrated smoothed particle hydrodynamics (SPH) model for complex interfacial flows with large density ratios is developed. The discrete continuity equation and acceleration equation are obtained by considering the time derivative of the volume of particle and Eckart's continuum Lagrangian equation. A continuum surface force model is used to meet the fact that surface force may not be distributed uniformly on each side of the interface. An improved boundary condition is imposed to model wall free‐slip and no‐slip condition for interfacial flows with large density ratios. Particle shifting algorithm (PSA) is added for interfacial flows by imposing the normal correction near the interface, called as Interface‐PSA. Then four representative numerical examples, including droplet deformation, Rayleigh‐Taylor instability, dam breaking, and bubble rising, are presented and compared well with reference data. It is demonstrated that inherent interfacial flow physics can be well captured, including surface tension and the dynamic evolution of the complex interfaces.

Discontinuous Galerkin (DG) methods are a class of finite element methods that use discontinuous basis functions. This particular feature enables the use of non-conforming meshes and facilitates th...

With increasing computer capabilities, numerical modeling of two-phase flows has developed significantly over the last few years. Although there are two main categories, namely 'one' fluid and 'two' fluid methods, the 'one' fluid methods are more commonly used for tracking or capturing the interface between two fluids. Level set (LS), volume-of-fluid (VOF), front tracking, marker-and-cell (MAC) and lattice Boltzmann (LB) methods are all 'one' fluid methods. It is clear that there is no perfect method; each method has advantages and disadvantages which make it more appropriate for one kind of problem than for others. For instance, a LS method will accurately compute the curvature and the normal to the interface, but tends to loss mass which is physically incorrect. On the other hand, a VOF method will conserve mass up to machine precision, but the computation of the curvature and normal to the interface is not as accurate. In order to minimize the disadvantages of these methods, several authors have used two or more methods together to model two-phase flows. This is the case for the CLSVOF (Couple Level Set Volume Of Fluid) method, where LS and VOF are coupled together in order to better capture the interface. In CLSVOF, the level set function is used to compute the interface curvature and normal to the interface, while the volume of fluid function is used to capture the interface. For two-phase flows in microchannels, surface tension forces play an important role in determining the dynamics of bubbles whereas gravitational forces are generally negligible. Also it is very important to consider the interaction between the boundaries and the fluids by prescribing or computing the correct contact angle between them. The commercial CFD code FLUENT allows the use of static constant angles, or the use of User Defined Functions (UDF) to compute the dynamic contact angles. It is inappropriate to use a static contact angle to model cases involving moving contact lines. For such cases a dynamic contact angle scheme should be implemented. In this study, FLUENT was used to model adiabatic and diabatic, time dependent two-phase flows. Since FLUENT already contains a VOF method, a LS method was implemented and coupled with VOF into FLUENT via UDFs. Furthermore, since the LS function, used to compute the surface tension force, ceases to be a signed distance to the interface even after one time step, a re-initialization equation was solved after each time step. This involved using a fifth order WENO (Weighted Essentially Non Oscillatory) scheme to discretize the space derivatives (otherwise oscillations of the interface occurred), and a first order Euler method for the time integration. In another part of the study, a 3D dynamic contact angle model based on volume fraction, interface reconstruction, and experimentally available advancing and receding static contact angles was also developed and implemented into FLUENT via UDFs. Several validations for the developed CLSVOF method and dynamic contact angle model are presented in this thesis, these includes a static bubble, a bubble rising in a stagnant liquid for Morton numbers ranging from 102 to 10-11, droplet deformation due to a vortex flow field, droplets spreading over a wall under the gravity effect and droplets sliding over a wall due to gravity. These validations demonstrated the high accuracy and the stability of our methods for modeling these phenomena. A heat and mass transfer model was also implemented into the commercial CFD code FLUENT for simulating of boiling (and condensation) heat transfer. Several simulations were presented with water and R134a as working fluids. The influence of the contact angle and the wall superheat was also studied.

The computation of flows with large density contrasts is notoriously difficult. To alleviate the difficulty we consider a discretization of the Navier-Stokes equation that advects mass and momentum in a consistent manner. Incompressible flow with capillary forces is modeled and the discretization is performed on a staggered grid of Marker and Cell type. The Volume-of-Fluid method is used to track the interface and a Height-Function method is used to compute surface tension. The advection of the volume fraction is performed using either the Lagrangian-Explicit / CIAM (Calcul d’Interface Affine par Morceaux) method or the Weymouth and Yue (WY) Eulerian-Implicit method. The WY method conserves fluid mass to machine accuracy provided incompressibility is satisfied. To improve the stability of these methods momentum fluxes are advected in a manner “consistent” with the volume-fraction fluxes, that is a discontinuity of the momentum is advected at the same speed as a discontinuity of the density. To find the density on the staggered cells on which the velocity is centered, an auxiliary reconstruction of the density is performed. The method is tested for a droplet without surface tension in uniform flow, for a droplet suddenly accelerated in a carrying gas at rest at very large density ratio without viscosity or surface tension, for the Kelvin-Helmholtz instability, for a 3mm-diameter falling raindrop and for an atomizing flow in air-water conditions.

SummaryIn this paper, a numerical method, which is about the coupling of continuous and discontinuous Galerkin method based on the splitting scheme, is presented for the calculation of viscoelastic flows of the Oldroyd‐B fluid. The momentum equation is discretized in time by using the Adams‐Bashforth second‐order algorithm, and then decoupled via the splitting approach. Considering the Oldroyd‐B constitutive equation, the second‐order Runge‐Kutta approach is selected to complete the temporal discretization. As for the spatial discretizations, the fundamental purpose is to make the best of finite element method (FEM) and discontinuous Galerkin (DG) method to handle different types of equations. Specifically speaking, for the subequations, FEM is chosen to treat the Poisson and Helmholtz equations, and DG is employed to deal with the nonlinear convective term. In addition, because of the hyperbolic nature, DG is also utilized to discretize the Oldroyd‐B constitutive equation spatially. This coupled method avoids resorting to extra stabilization technique occurred in standard FEM framework even for moderately high values of Weissenberg number and also reduces the complexity compared with unified DG scheme. The Oldroyd‐B model is applied to investigate several typical and challenging benchmarks, such as the 4:1 planar contraction flow and the lid‐driven cavity flow, with a wide range of Weissenberg number to illustrate the feasibility, robustness, and validity of our coupled method.

The standard weakly compressible Smoothed Particle Hydrodynamics (WCSPH) is successfully applied to multi-phase problems involving fluids with similar densities, but when density ratio increases at some order of magnitude, serious instability phenomena occur at the interface. Several remedies have been proposed based on numerical correctives that deviate from standard formulation, increasing the algorithm complexity and, sometimes, the computational cost. In this study, the standard SPH has been adapted to treat free-surface multi-phase flows with a large density ratio through a modified form of the governing equations which is based on the specific volume (i.e. the inverse of particle volume) instead of density: the former is continuous across the fluid interface while the latter is not and generates numerical instability. Interface sharpness is assured without cohesion forces; kernel truncation at the interface is avoided. The model, relatively simple to implement, is tested by simulating two-phase dam breaking for two configurations: kinematic and dynamic features are compared with reference data showing good agreement despite the reduced computational effort.

The cross-sectional shape of the lift pipe of the bubble pump is considered to have a significant impact on the merging of rising bubbles. The effects of cross-sectional shapes of pipes and Reynolds number on the rise and coalescence of near-wall bubbles in lithium bromide solutions are investigated. A hybrid thermal multiphase model incorporating the D3Q15 pseudopotential multiphase lattice Boltzmann model with a large density ratio and the finite difference method is utilized. And the volume translated Peng-Robinson (VTPR) state equation is adopted. It is found that the merging time of rising bubbles near the walls of rectangular, circular, elliptical, and triangular cross-section pipes increases in order at the same Reynolds number. Under the same cross-sectional conditions of the pipes, a suitable increase in the Reynolds number will accelerate the merging of bubbles near the wall. This study may offer valuable insights for informing the design of pipe types for bubble pumps in absorption refrigeration systems.

The field of the mathematical and numerical analysis of systems of nonlinear pdes involving interfaces and free boundaries is a burgeoning area of research. Many such systems arise from mathematical models in ma- terial science and fluid dynamics such as phase separation in alloys, crystal growth, dynamics of multiphase fluids and epitaxial growth. In applications of these mathematical models, suitable performance indices and appropriate control actions have to be specified. Mathematically this leads to optimiza- tion problems with pde constraints including free boundaries. It is now timely to consider such control problems because of the maturity of the field of com- putational free boundary problems. The aim of the mini-workshop was to bring together leading experts and young researchers from the separate fields of numerical free boundary problems and optimal control in order to estab- lish links and to identify suitable model problems to serve as paradigms for progressing knowledge of optimal control of free boundaries.