ABSTRACT. This work presents a python package built on top of FEniCS, FEniCS Mechanics, developed to carry out computational mechanics simulations. The goal of this package is to facilitate defining and numerically solving mechanics problems through creating a rapid prototyping environment for users of all programming levels and mechanics knowledge. Problems in continuum mechanics are generally described by the weak form for the balance of linear momentum, \begin{gather} \int_R \boldsymbol\xi\cdot\rho\mathbf{a}\;dv + \int_R \pd{\bs\xi}{\mathbf{x}}\cdot\mathbf{T}\;dv = \int_R \boldsymbol\xi\cdot\rho\mathbf{b}\;dv + \int_{\Gamma_q}\bs\xi\cdot\bar{\mathbf{t}}\;dv, \end{gather} where R denotes the current configuration of the body, $\rho$ the density of the material, $\mathbf{a}$ the acceleration, $\mathbf{T}$ the Cauchy stress tensor, $\mathbf{b}$ the body force applied, $\bar{\mathbf{t}}$ a prescribed Cauchy traction on $\Gamma_q\subset\partial R$, and $\boldsymbol\xi$ an arbitrary, kinematically admissible, vector field. Note that Equation (1) is applicable to any continua since a constitutive relation is yet to be specified. This facilitates the implementation of various materials within the same framework. The user can also turn incompressibility on, triggering a mixed function space formulation that results in a block system. Furthermore, inverse elastostatics has been implemented as described by Govindjee and Mihalic [1]. This method is used to determine the unloaded geometry when the loaded configuration under a given load is known. In addition to supporting fluid and solid mechanics problems, the simulations can be specified as quasi-static or time-dependent, where a generalized-$\alpha$ method is used for time integration.

The input provided by the user is a python dictionary defining the problem. This dictionary contains the subdictionaries material, mesh, and formulation. The constitutive equation and relevant material parameters are given in material. Mesh information, as well as function space discretization is given in mesh. Lastly, boundary conditions, time integrating parameters, domain formulation, and initial conditions are given in formulation. The user then creates a problem object by providing this dictionary. During instantiation of this object, the values in the configuration dictionary, and their combination, are checked to ensure they are consistent with physical processes. Then, a solver object can be created and used to compute the solution to Equation (1).

To illustrate a simple use of this package, consider the quasi-static loading of an incompressible unit cube modeled as a neo-Hookean material. The material parameters are $\kappa = 10$ GPa, and $\mu = 1.5$ MPa, where $\kappa$ and $\mu$ are the bulk and shear moduli of natural rubber, respectively. The cube is clipped at x = 0 m, while a traction of $\bar{\mathbf{t}} =$ (1 MPa) $\mathbf{e}_1$ is applied at x = 1 m. The code for this problem is given below: \begin{python} import dolfin as dlf import fenicsmechanics as fm

config = {'material': { 'const_eqn': 'neo_hookean', # constitutive equation 'type': 'elastic', # choose between elastic and viscous 'incompressible': True, # toggle incompressibility 'density': 1.0, 'kappa': 10e9, # bulk modulus [Pa] 'mu': 1.5e6}, # shear modulus [Pa] 'mesh': { 'mesh_file': 'mesh.h5', # mesh file in xml or hdf5 format 'mesh_function': 'mesh_function.h5', 'element': 'p2-p1'}, # finite element choices, # p1, p2, p2-p1, etc. 'formulation': { 'time': {'unsteady': False}, # quasistatic or unsteady using gen. alpha 'domain': 'eulerian', # lagrangian or eulerian 'inverse': False, # toggle inverse finite elastostatics 'bcs': { 'dirichlet': { 'displacement': [dlf.Constant([0.,0.,0.])], 'regions': [1]}, # surface tags 'neumann': { 'regions': [2], # surface tags 'types': ['cauchy'], # choose between piola, cauchy, # and pressure 'values': [dlf.Constant([1e6,0.,0.])]} # traction vector [Pa] }}}

problem = fm.SolidMechanicsProblem(config) # set up solid mechanics problem solver = fm.SolidMechanicsSolver(problem) # set up solid mechanics solver solver.full_solve(fname_disp='results/displacement.pvd') # Solve PDE and save \end{python} Running the above script in parallel yielded the results shown in Figure \ref{subfig:cube_elongation}, where the maximum displacement is 0.263 m. \begin{figure}[h] \begin{subfigure}{0.47\textwidth} \centering \includegraphics[width=0.8\textwidth]{figures/cube_elongation-neo_hooke.png} \caption{Elongation of a unit cube (1m$\times$1m$\times$1m) with the unloaded configuration shown in transparent gray.} \label{subfig:cube_elongation} \end{subfigure} \hfill \begin{subfigure}{0.47\textwidth} \centering \includegraphics[width=\textwidth]{figures/ventricle.png} \caption{The idealized left ventricle in its unloaded (left) and loaded (right) states.} \label{subfig:ventricle} \end{subfigure} \end{figure}

The code has been verified using several test cases and comparing simulation results against in-house computational tools, as well as against demos within FEniCS repositories. Furthermore, inflation of an idealized left ventricle was compared with results obtained from other groups [2], and is shown in Figure 1a. The simulation of the left ventricle was accomplished with minimal changes to the configuration dictionary shown above.

FEniCS Mechanics will be available through git repositories once further quantitative verification is performed for quality assurance. This code makes modeling of custom mechanics problems even more accessible to those with a limited programming background or knowledge of mechanics.

References [1] Govindjee, Sanjay, and Paul A. Mihalic. “Computational Methods for Inverse Finite Elastostatics.” Computer Methods in Applied Mechanics and Engineering 136.1-2 (1996): 47-57. [2] Land S et al. “Verification of Cardiac Mechanics Software: Benchmark Problems and Solutions for Testing Active and Passive Material Behaviour.” Proceedings of the Royal Society of London A 471 (2015): 20150641.

ABSTRACT. The configuration of a so-called intrinsic beam model with shear deformation is defined by the position vector \boldsymbol{x} and the associated orientation tensor \boldsymbol{\alpha}\in\mathrm{SO}(3). Let identify the position and orientation in the deformed configuration with an appended prime, \boldsymbol{x}' and \boldsymbol{\alpha}', with \boldsymbol{x} and \boldsymbol{\alpha} the corresponding position and orientation in the reference configuration. The corresponding linear form reads

\begin{split}\int_{L}\left(\delta\boldsymbol{\varepsilon}\boldsymbol{T}+\delta\boldsymbol{\beta}\boldsymbol{M}\right)\mathrm{d}s & =\int_{L}\left(\delta\boldsymbol{x}'\boldsymbol{t}+\boldsymbol{\varphi}_{\delta}\boldsymbol{m}\right)\mathrm{d}s+\sum_{i}\delta\boldsymbol{x}'_{i}\boldsymbol{F}_{i}+\boldsymbol{\varphi}_{\delta i}\boldsymbol{C}_{i}\end{split}

where s is the arc length and

\begin{split}\boldsymbol{\varepsilon} & =\boldsymbol{\alpha}'^{T}\boldsymbol{x}'_{/s}-\boldsymbol{\alpha}^{T}\boldsymbol{x}'_{/s}\\ \boldsymbol{\beta} & =\boldsymbol{\alpha}'^{T}\mathrm{ax}(\boldsymbol{\alpha}'_{/s}\boldsymbol{\alpha}'^{T})-\boldsymbol{\alpha}^{T}\mathrm{ax}(\boldsymbol{\alpha}_{/s}\boldsymbol{\alpha}^{T})\\ & =\boldsymbol{\alpha}^{T}\boldsymbol{\varPhi}^{T}\mathrm{ax}(\boldsymbol{\varPhi}_{/s}\boldsymbol{\varPhi}^{T}) \end{split}

are the generalized deformation measures energetically conjugated to the internal forces and moments T and M; the axial operator \boldsymbol{a}=\mathrm{ax}(\boldsymbol{A}) is the operator extracting the vector \boldsymbol{a} that characterize the skew symmetric part of tensor \boldsymbol{A}, i.e. \boldsymbol{a}\times=1/2(\boldsymbol{A}^{T}+\boldsymbol{A}) ; \boldsymbol{\varPhi}=\boldsymbol{\alpha}'\boldsymbol{\alpha}^{T} is the rotation tensor bringing \boldsymbol{\alpha} into \boldsymbol{\alpha}'. Being \boldsymbol{\alpha} and \TT{\varPhi} orthogonal, tensors \boldsymbol{\alpha}'_{/s}\boldsymbol{\alpha}'^{T}, \boldsymbol{\alpha}_{/s}\boldsymbol{\alpha}^{T} and \boldsymbol{\varPhi}_{/s}\boldsymbol{\varPhi}^{T} are skew-symmetric. Vectors t and m are the the external forces and moments per unit length, respectively; vectors F_{i} and C_{i} are the concentrated external forces and moments. Finally, the virtual rotation vector \T{\varphi}_{\delta} characterizes the virtual rotation of the orientation tensor in the deformed configuration, \boldsymbol{\varphi}_{\delta}=\mathrm{ax}(\delta\boldsymbol{\alpha}'\boldsymbol{\alpha}'^{T})=\mathrm{ax}(\delta\boldsymbol{\varPhi}\boldsymbol{\varPhi}^{T}) . The internal forces and moments, \boldsymbol{T} and \boldsymbol{M} , are functions of the generalized strains \boldsymbol{\varepsilon} and \boldsymbol{\beta}. A linear constitutive law can often be assumed.

The rotation tensor \TT{\varPhi} belongs to the Special Orthogonal Group \text{SO}(3). Tensors fields belonging to the intrinsically nonlinear group \text{SO}(3) needs to be treated differently than vectors/tensors that belongs to linear spaces, such as the position or displacement fields. Thus, the solution of linear forms involving finite rotations can be challenging when using finite element libraries that assume – from the ground up – that all unknown fields do belong to a linear space. In a nutshell, a parametrization of \text{SO}(3) is often chosen, and the parameters are interpolated using standard finite elements, that are linear with respect to the unknown field. The rotation vector \T{\varphi} is one of the most often chosen parametrization. The corresponding rotation tensor \TT{\alpha} can be computed using the exponential map,\TT{\varPhi}=\exp(\T{\varphi}\times)=\sum_{n=0}^{\infty}\frac{1}{n!}(\T{\varphi}\times)^{n}=I+\frac{\sin\varphi}{\varphi}\T{\varphi}\times+\frac{1-\cos\varphi}{\varphi^{2}}\T{\varphi}\times\T{\varphi}\times with \varphi the norm of \T{\varphi}; the virtual rotation vector \T{\varphi}_{\delta} can be computed from the variation of the rotation parameter \delta\T{\varphi} as \T{\varphi}_{\delta}=\TT{\varGamma}(\T{\varphi})\delta\T{\varphi} , with tensor \TT{\varGamma} defined by\TT{\varGamma}(\T{\varphi})=\text{dexp}(\T{\varphi}\times)=\sum_{n=0}^{\infty}\frac{1}{(n+1)!}(\T{\varphi}\times)^{n}=I+\frac{1-\cos\varphi}{\varphi^{2}}\T{\varphi}\times+\frac{\varphi-\sin(\varphi)}{\varphi^{3}}\T{\varphi}\times\T{\varphi}\times; Equation [WeakForm-PLV] and its linearization need to be written as a function of \T{\varphi}, with test and trail functions \delta\T{\varphi} and \partial\T{\varphi}, respectively; tensor \TT{\varGamma} enters the formulation because \T{\varphi}_{\delta}=\TT{\varGamma}\delta\T{\varphi} , \T{\varphi}_{\partial}=\TT{\varGamma}\partial\T{\varphi} and \mathrm{ax}(\boldsymbol{\varPhi}_{/s}\boldsymbol{\varPhi}^{T})=\TT{\varGamma}\T{\varphi}_{/s} . The parametrization is singular for 2\pi: given an orientation tensor \TT{\alpha} the corresponding rotation vector is defined up to a coaxial rotation of 2\pi; furthermore, and more importantly, tensor \TT{\varGamma} is singular for \varphi=2\pi. Note also that the coefficients \sin(\varphi)/\varphi, (1-\cos(\varphi))/\varphi^{2} and (\varphi-\sin(\varphi))/\varphi^{3} need to be evaluated near \varphi=0 using their Maclaurin series; the expansion must be truncated to an order that is high enough to guarantee the sought precision not only for the coefficients themselves, but also for their second derivative, that is required to compute the linearization \partial\delta\T{\beta} of \delta\T{\beta}=\boldsymbol{\alpha}^{T}\delta(\boldsymbol{\varPhi}^{T}\mathrm{ax}(\boldsymbol{\varPhi}_{/s}\boldsymbol{\varPhi}^{T}))=(\boldsymbol{\alpha}'^{T}\TT{\varGamma}\T{\varphi}_{/s})\times\boldsymbol{\alpha}'^{T}\TT{\varGamma}\delta\T{\varphi}+\TT{\alpha}'^{T}\otimes\T{\varphi}_{/s}:\delta\TT{\varGamma}+\TT{\alpha}'^{T}\TT{\varGamma}\delta\T{\varphi}_{/s} .

Summing up, and after writing a python library that allows to deal with finite rotations and the above sketched parametrization, the above formulation can be implemented in FEniCS provided that

• the expressions of \TT{\varPhi} and \TT{\varGamma} are computed correctly, with a truncated Macalurin series for \varphi\approx0 ;

• the magnitude of the rotation vector is kept lower that 2\pi, with ideally \varphi\in[-\pi,\pi].

The first requirement can be accomplished by resorting to Uflacs' excellent handling of conditional expressions. For example, coefficient \sin(\varphi)/\varphi can be computed asdef a(phi2): return conditional(le(phi2, ath), aexp(phi2), sin(sqrt(phi2))/sqrt(phi2)) where phi2 =\varphi^{2}, ath is the threshold for switching from the Maclaurin series to the trigonometric expansion, and aexp(phi2) computes the Maclaurin truncated series. Using these conditionals withing the library methods written to compute \TT{\varPhi} and \TT{\varGamma} allows to easily write the linear form Eq. [WeakForm-PLV] and have the corresponding bilinear form automatically generated.

Dealing with the second requirement is a bit more tricky. The chosen approach is suitable if the rotation field is approximated using Lagrange elements. In other words, it is doomed to fail for e.g. Nedelec elements. First of all, a C++ Expression computes either the rotation vector \T{\varphi} or, if \varphi>\pi, its complementary vector \T{\varphi}_{c}=\T{\varphi}(1-2\pi/\varphi). This Expression is used to replace the values of \T{\varphi} with its complementary rotation angle when \varphi>\pi. This well known trick [1] allows to keep the rotation vector magnitude under control, and works perfectly well if all the rotation vector unknowns of an element get changed into their complementary. Unfortunately one must take care of the elements for which only some of the nodal rotation vector unknowns have been changed into their complementary, see e.g. [2]. This would result in an incorrectly interpolated field. To fix this problem, a custom Assembler is required, whose assemble method checks, for each element, whether only some of the rotation vector unknowns have been changed. If this is the case then it changes back, in place, the data that gets passed to the element's tabulate_tensor. Then, after calling tabulate_tensor, is modifies the assembled vector/tensor in order to account that the locally changed rotation vector unknowns are actually a function of the complementary vector. This is accomplished by accounting for the derivative of the complementary vector\frac{\partial\boldsymbol{\varphi}_{C}}{\partial\boldsymbol{\varphi}}=\left(1+\frac{2\pi}{\varphi}\right)\boldsymbol{I}-\frac{2\pi}{\varphi^{3}}\boldsymbol{\varphi}\otimes\boldsymbol{\varphi}, either on the left for the test functions \delta\T{\varphi} or on the right for the trial functions \partial\T{\varphi}.

This rotation handling strategy is not free from problems, see e.g. [3] and references therein. Nonetheless, the resulting finite elements, whose complex Jacobian matrix is computed automatically, are able to withstand finite rotations of any magnitude and to correctly solve complex nonlinear problems. As an example, the very first deformed configurations of a classic benchmark test are plotted in Fig. [fig:Wrench-beam-testcase].

The library is general enough to be suitable for the implementation of intrinsic shell models with the so-called drilling degree of freedom. In this case, however, special care should be given to alleviate the shear locking problem, that is rather easily dealt with for beam finite elements.

[float Figure:

[Figure 1: Wrench beam testcase, first load steps. ] ]

References

[1] M. Ritto-Corrêa, D. Camotim: On the differentiation of the Rodrigues formula and its significance for the vector-like parameterization of Reissner-Simo beam theory, IJNME, 55 (2002), 1005-1032.

[2] J. Mäkinen: Total Lagrangian Reissner's geometrically exact beam element without singularities, IJNME, 70 (2007), 1009-1048.

[3] T. Merlini, M. Morandini: On successive differentiations of the rotation tensor: An application to nonlinear beam elements, JOMMS, 8 (2013), 305-340.

ABSTRACT. Mathematical modeling and numerical solution of the applied problems in highly heterogeneous media is a actual, necessary and interesting task. A classical numerical method for solving such problems should use a computational grid that resolve all small heterogeneities. For such problems, a numerical homogenization or multiscale methods are used. For the processes in a periodic heterogeneous media, where the length of the period is small in comparison with other parameters of the problem, the asymptotic homogenization methods is used to construct coarse-scale approximation and calculate effective properties that takes into account small scale heterogeneities. Nowadays, a multiscale methods are becoming popular, for example, heterogeneous multiscale method (HMM), multiscale finite element method (MsFEM), variational multiscale method (VMS) and others. Multiscale methods should combine the simplicity and efficiency of a coarse-scale models, and the accuracy of microscale approximations. In this work, we construct a reduced order model using the generalized multiscale finite element method (GMsFEM). This method involves two basic steps: (1) the construction of multiscale basic functions that take into account small scale heterogeneities in the local domains and (2) the construction of the coarse scale approximation. For the construction of the multiscale basic functions we solve spectral problems in local domains. Spectral problems help to identify the most important characteristics of the solution. In contrast to the available techniques, this method allows to avoid the limitations associated with idealization and limitations on the applicability of the method. This method also more general technology that takes into account the different scale processes. The construction of basic functions occurs independently for each local domain, doesn't require the exchange of information between processors, and has a high parallelization efficiency. Using constructed multiscale basic functions, we construct a mathematical model on a coarse grid that allows to significant reduce the solution time, the amount of used memory, and can be used to perform calculations for a given configuration of heterogeneous properties. We consider several applied problems with different types of heterogenuities. As first problem, we consider gas filtration in the fractured poroelastic medium and present multiscale method for solution of the problem on a coarse grid. Next, we consider wave propagation in the elastic fractured media, and finally present multiscale method for the solution of the pore-scale electrochemical processes in the litium-ion batteries.

ABSTRACT. We consider multi-level Monte Carlo methods (MLMC) [1] for accelerating the computation of stochas- tic eigenvalue problems arising from PDEs in the presence of uncertainty. Eigenvalue problems are impor- tant in a number of applications where randomness in the underlying system has a significant impact on the response, such as structural vibration. In these systems, designing to coincide with or avoid natural frequencies is important. Monte Carlo (MC) methods can quantify the impact of uncertainties on eigenvalue problems in a simple way, but are considered to be too expensive in many practical cases. We propose employing MLMC methods to accelerate the computation of the statistics of eigenvalue problems to make MC techniques viable. MLMC methods compute the expectation and variance of the quantity of interest, Q, by solving the governing equation on grids of varying fidelity. The expectation of Q at the finest level, L, is calculated using the following telescopic sum:

E[Q_L] = E[Q_0] + \sum_{l=0}^{L}E[Q_l − Q_{l−1}]

where increasing l represents levels of increasing mesh refinement. More expensive ‘finer grids’ are sampled less often than coarser grids because the variance of the ‘finer’ levels is smaller than the variance between the ‘coarser’ levels. This achieves equivalent accuracy as MC methods but at a much lower cost. The total variance of the estimator can also be calculated by summing the variance of the estimator at each level.

We consider uncertainty in our systems from random coefficients in the governing equations. We solve the eigenvalue problem using a discrete realisation of the random field and add point masses to the weak form of the equation. This is similar to the approach used in other random vibration research, e.g. [2]. We discuss the development of the FEniCS PointSource code, which now enables correct addition of the point sources to mass matrices. We use the SLEPc eigenvalue solver combined with the FEniCS libraries to compute the natural frequencies for each random realisation and schedule and run the MLMC algorithm using the podS library in parallel on HPC and Microsoft Cloud facilities. We argue that the demonstrated efficiency gains of the multi-level method make the approach compet- itive with alternatives to MC for eigenvalue problems. This enables us to compute quantities of interest when the underlying system has many degrees of freedom without having to make assumptions that reduce the accuracy of the result.

[1] M. Giles: Multilevel Monte Carlo methods, Acta Numerica, 24 (2015), 259-328. [2] R. S. Langley and A. W. M. Brown: The ensemble statistics of the energy of a random system subjected to harmonic excitation, Journal of Sound and Vibration, 275 (2004), 823-846.

ABSTRACT. The simulation of the propagation of elastic waves in unbounded domain is still an open problem. Different kinds of absorbing boundary conditions (ABC) has been introduced, such as Gaussian taper method, high-order ABC or paraxial approximations. Among them, one of the most successful approach is the perfectly matched layers (PML) method, introduced by Berenger in the 90s in the context of electrogamnetism, is now widely used in elastodynamics, acoustic, poroelasticity and other research fields. A PML is an absorbing layer model for linear wave equations that theoretically absorbs, quite efficiently, propagating waves of all non-tangential angles of incidence and of all non-zero frequencies.

This work is inspired by the article of Kucukcoban & Kallivokas (2013), in which they propose a mixed and hybrid formulation of the elastodynamics equations using the PML method. We show that this method can be easily implemented in Fenics. Some results are presented to demonstrate feasibility, stability and accuracy of this implementation.

ABSTRACT. An acoustic metamaterial (AMM) is an engineered composite material, often constructed using a periodic lattice, that behaves as an effective material with unconventional dynamic effective properties. AMMs can manipulate and control acoustic and/or elastic waves in ways that are not possible with conventional materials, and are therefore of interest for applications such as acoustic lenses that defy the diffraction limit, material slabs that prevent certain frequencies from propagating, and acoustic cloaks that redirect acoustic signals around an object. Reconfigurable AMMs allow even greater control over the direction of energy propagation by altering the structure through an external stimulus, such as mechanical deformation, which fundamentally changes the effective dynamic properties. In this work, a periodic arrangement of pre-curved beams in a honeycomb configuration, commonly referred to as negative stiffness honeycombs (NSH), is studied as a candidate reconfigurable AMM. This work seeks to explore the degree to which exotic dynamic properties and extreme changes in the modal dispersion can be elicited by uniaxial pre-strain of the NSH. This talk will highlight the use of FEniCS to compute the dispersion properties of NSH, as well as describe how this model can fit into a design framework for finding optimal NSH configurations for particular applications.

ABSTRACT. In this contribution we present a novel theoretical analysis of PCD (pressure-convection-diffusion) preconditioner for approximation of Navier-Stokes equations, present some important implementational details, and provide an open-source implementation suitable for deployment on HPC systems.

ABSTRACT. The wetting of liquid films on solid surfaces occurs in many industrial processes like coating or painting and apparatuses involving trickle films. At the intersection of the gas-liquid interface with the solid surface a moving contact line is formed. Applying the common no-slip boundary condition at the solid surface, a non-physical divergent stress at the contact line occurs. One possibility to circumvent this difficulty in the context of continuum mechanics is the coupling of the incompressible Navier-Stokes equations with the Cahn-Hilliard (CH) equation. The CH equation models the interface between the fluids with a diffuse interface of positive thickness and describes the distribution of the different fluids by a smooth indicator function. Especially, the CH equation allows the contact line to move naturally on the solid surface due to a diffusive flux across the interface, which is driven by the gradient of the chemical potential. The Cahn-Hilliard-Navier-Stokes (CHNS) model forms a tightly coupled and nonlinear system of four partial differential equations.

In this talk, we present a new FEniCS-based library for massive parallel simulations using CHNS models for problems involving moving contact lines. The library is programmed mainly in Python, whereas some time critical parts are in C++. Furthermore, Parallel HDF5 is applied for in- and output and PETSc's SNES with MUMPS is used for the solution of the nonlinear system. The library aims to be easily extendable and to serve as a platform for rapid benchmarking of different CHNS models and solution schemes. For efficient simulations of the thermodynamically consistent CHNS model by [1] for large density ratios, we employ a stable solution scheme as proposed in [2]. Besides the structure of the library, we discuss the compilation of the FEniCS toolchain and our own implementations on a Cray XC40/XC30 supercomputer system. Furthermore, we present some promising physical results on moving contact lines and discuss the good scalability on the supercomputer.

[1] Helmut Abels, Harald Garcke, and Günther Grün. "Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities." Mathematical Models and Methods in Applied Sciences 22.03 (2012): 1150013. [2] Harald Garcke, Michael Hinze, and Christian Kahle. "A stable and linear time discretization for a thermodynamically consistent model for two-phase incompressible flow." Applied Numerical Mathematics 99 (2016): 151-171.

ABSTRACT. We investigate an exactly mass conserving higher order interior penalty discontinuous Galerkin method for the incompressible and variable density Navier-Stokes equations. Our intended application is low viscosity flows with density fields containing large and sharp jumps such as simulation of multi phase flows with a water/air free surface embedded in the computational domain.

Sharp jumps in the solution gives rise to convective instabilities due to Gibbs phenomena. For our application we see these instabilities clearly in the momentum equation when used with a discontinuous density field. In most cases the simulations will almost immediately break down due to exponential growth of the velocity in localised regions near the discontinuity due to the factor 1000 jump in density.

In the context of DG methods there are existing methods known as slope limiters for eliminating these non-linear convective instabilities in scalar transport equations. We extend the existing work with an approximate slope limiter for divergence free, solenoidal, vector fields and compare it with an existing hierarchical Taylor based scalar limiter applied to each of the velocity components.

We investigate the combination of the solenoidal approximate limiter with the stable hierarchical Taylor limiter and establish a higher order method for two phase flows that is mass conserving and has no unphysical smearing of the fluid properties near the free surface as well as no Gibbs instabilities due to the jump in density.