The subject of this thesis is the derivation and analysis of numerical approximations of multi-component, multi-phase field systems. Recent approximations of solutions to such models are mostly based on explicit time stepping schemes and require the computation of many time steps. Implicit methods exhibit inherent numerical challenges, in particular due to the non-smoothness of the underlying energy functionals.

Our focus lies on the derivation of numerical approximations within the thermodynamically consistent context with high efficiency and robustness. We aim to exploit the special mathematical structure of the model and the underlying thermodynamics without introducing additional regularizations.

We introduce the thermodynamic and multi-phase setting in chapter 2 and continue by motivating and presenting a thermodynamically consistent multi-component, multi-phase field model in chapter 3. Based on Rothe’s method, we obtain a semi-discretization allowing for adaptive meshes in chapter 4 and the implicit problems are analyzed. In chapter 5, a full discretization with adaptive finite elements based on hierarchical a posteriori error estimation is set up. We transition to a purely algebraic formulation and present the iterative approximation of solutions with a nonsmooth Schur–Newton multigrid approach in chapter 6. Finally, in chapter 7, we perform numerical experiments to underline the thermodynamical consistency and numerical efficiency of our method.

Discrete-continuous hybrid models are a popular means for describing elastic membrane-mediated particle interactions in and on lipid bilayers. Here, the continuous part is usually given by an approximation of the lipid membrane by an infinitely thin and sufficiently smooth hypersurface, whose elastic energy is determined by a Canham-Helfrich type functional. The discrete component results from modeling non-membrane particles as rigid discrete entities, which, depending on their configuration, induce local constraints on the membrane along the membrane-particle interfaces. In this context, the interaction potential describes the optimal elastic energy of such hybrid systems with a fixed particle configuration. Correspondingly, the energy minimization principle yields that stationary particle configurations are given by the local minima of the interaction potential. The main goal of this work is the proof of differentiability of the interaction potential for a selected class of models. This is accomplished using a variational approach that is already established in the literature in order to develop and apply robust numerical optimization methods for computing stationary particle configurations. Correspondingly, an additional focus is the derivation of a numerically accessible representation of the gradient, including its discretization and relevant numerical analysis. The proof of differentiability is brought forward by an application of the implicit function theorem. The basis for this is so-called boundary preserving domain transformations, which are induced by suitable families of vector fields and which locally admit the reformulation of the minimization problem that is implicitly defined by the interaction potential with respect to a fixed particle configuration. This subsequently enables the representation of the gradient as a volume integral using matrix analysis methods. The discretization of the partial differential equations for describing optimal membrane shapes is done via finite element methods. For particle methods with so-called curve restrictions a fictitious domain stabilized Nitsche method is developed, and for models with point value restrictions a conforming Galerking discretization is made possible by local QR transformations of the nodal finite element basis. For both cases suitable a priori error estimates are proven, and in addition also error estimates for the volume representation of the gradient are shown within that context. These developed methods open up the domain of efficient simulation of macro structures by isotropic and anisotropic particles, which is illustrated with the aid of various example applications and by means of perturbed gradient methods.

Partial differential equations with random coefficients (random PDEs) is a very developed and popular field. The variety of applications, especially in biology, motivate us to consider the random PDEs on curved moving domains. We introduce and analyse the advection-diffusion equations with random coefficients on moving hypersurfaces. We consider both cases, uniform and log-normal distributions of coefficients. Furthermore, we will introduce and analyse a surface finite element discretisation of the equation. We show unique solvability of the resulting semi-discrete problem and prove optimal error bounds for the semi-discrete solution and Monte Carlo samplings of its expectation. Our theoretical findings are illustrated by numerical experiments. In the end we present an outlook for the case when the velocity of a hypersurface is an uniformly bounded random field and the domain is flat.

In this thesis we introduce a novel framework for uncertainty quantification in problems with random coefficients. The developed framework utilizes the ideas of multilevel Monte Carlo (MLMC) methods and allows for exploiting the advantages of adaptive finite element techniques. In contrast to the standard MLMC method, where levels are characterized by a hierarchy of uniform meshes, we associate the MLMC levels with a chosen sequence of tolerances. Each deterministic problem corresponding to a MC sample on a given level is then approximated up to the corresponding accuracy. This can be done, for example, using pathwise a posteriori error estimation and adaptive mesh refinement techniques. We further introduce an adaptive MLMC finite element method for random linear elliptic problems based on a residual-based a posteriori error estimation technique. We provide a careful analysis of the novel method based on a generalization of existing results, for deterministic residual-based error estimation, to the random setting. We complement our theoretical results by numerical simulations illustrating the advantages of our approach compared to the standard MLMC finite element method when applied to problems with random singularities.

Dynamic large deformation contact problems arise in many industrial applications like auto mobile engineering or biomechanics but only few methods exists for their numerical solution, all having their advantages and disadvantages. In this thesis the numerical solution of large deformation contact problems is tackled from an optimisation point of view and an application of this approach within a femoroacetabular impingement analysis is described. In this thesis we use a non-smooth Hamilton principle and Fréchet subdifferential calculus to derive a weak formulation of the problem. The resulting subdifferential inclusion is discretised in time by constructing a contact-stabilised midpoint rule. For the spatial discretisation the state-of- the-art dual mortar method is applied which results in non-convex constrained minimisation problems that have to be solved solved during each time step. For the solution of these problems an inexact filter trust-region method is derived which allows to use inexact linearisations of the non-penetration constraints. This method in combination with fast monotone multigrid method is then shown to be globally convergent.

This work is concerned with the proof of optimal error bounds for the discretization of $H^1$-elliptic minimization problem with solutions taking values in a Riemannian manifold. The discretization is done using Geodesic Finite Elements, a method of arbitrary order that is invariant under isometries. The discretization error is considered both intrinsically in a specially introduced Sobolev-distance as well as extrinsically. Optimal estimates of $H^1$- and $L^2$-type are shown, that have been observed experimentally in previous works of other authors. Using the Rothe method consisting of an implicit Euler method for the time discretization and Geodesic Finite Elements for the spatial discretization, error estimates for $L^2$-gradient flows of $H^1$-elliptic energies are derived as well. The core of the work is formed by the discretization error estimates for minimization problems in instrinsic $H^1$- and $L^2$-distances. To derive these, inverse estimates and interpolation errors for Geodesic Finite Elements and their discrete variations are shown. Using a nonlinear Cea's Lemma, this leads to the $H^1$-error estimate for minimizers of $H^1$-elliptic energies. A generalization of the Aubin-Nitsche-Lemma shows optimal $L^2$-error estimates for (essentially) semilinear energies, as long as the dimension of the domain of the minimizer is limited to $d<4$ for technical reasons. All results are illustrated using harmonic maps into a smooth Riemannian manifold satisfying certain curvature bounds as an example.

In this work, the model of rate-and-state friction, which can be viewed as central to the numerical simulation of earthquakes, is considered from a mathematical point of view. First, a framework is presented through which a general class of such friction laws can be understood and analysed. A prototypical viscoelastic problem of earthquake rupture is then formulated, both in strong and in variational form. Analysis of this problem is difficult, since the incorporation of rate-and-state friction leads to a coupling of variables. In a time-discrete setting, nonetheless, results on existence, uniqueness, and continuous parameter dependence of solutions can be obtained. The principal idea is to reformulate the variable interdependence as a fixed point problem and to prove convergence for a corresponding iteration. With that in mind, next, a numerical algorithm is proposed that resolves the coupling through a fixed point iteration. Since it puts a state-of-the-art solver and adaptive time stepping to use, it is not only stable but also fast. Its applicability to problems of interest is demonstrated in the penultimate chapter, which focuses on simulations of megathrust earthquakes that form at the base of a subduction zone. The main assumptions made throughout this work are summarised and discussed in the last chapter.

The ban of lead in electronics solder by EU directives results in the technological challenge to develop lead-free alternatives with comparable life span and processing properties. Numerical simulations of the microstructure evolution may contribute to identify promising candidates and thus focus the immense experimental effort. Aim of this work is on the one hand to develop a numerical framework for the efficient and robust simulation of the microstructure evolution in binary alloys combining adaptive finite element methods with fast solvers for the Cahn-Hilliard model. On the other hand we will extend the existing fast solvers for the discrete scalar Cahn-Hilliard equation to the vector-valued case. After some preliminary remarks on phase diagrams, phase separation, and phasefield models in Chapter 1 we will firstly discuss anisotropic Allen-Cahn equations in Chapter 2. Alle-Cahn-like problems arise as subproblems in the Nonsmooth Schur-Newton (NSNMG) method for Cahn- Hilliard equations in Chapters 3 and 4. Here we prove existence and uniqueness of solutions to the anisotropic Allen-Cahn equation with logarithmic potential using the theory of maximal monotone operators. For the numerical solution we introduce an adaptive spatial mesh refinement cycle for evolution problems and several variants of implicit Euler time discretization. We prove stability for the latter and numerical experiments conclude the chapter. Chapter 3 combines existing and newly developed numerical tools to a simulation software for microstructure evolution in binary alloys. Key ingredients are the adaptive mesh refinement cycle of Chapter 2, the NSNMG solver, a quantification algorithm for measuring "coarseness" of microstructures and a quotient space multigrid method for indefinite problems. An application of this software to simulate the microstructure evolution in a eutectic AgCu alloy shows only marginal impact of elastic stresses on coarsening in the setting considered; while the use of a smooth interpolant of the logarithmic potential affects the coarsening dynamics considerably. In the final chapter we consider the multicomponent Cahn-Hilliard equation and derive a unified formulation for the discrete problems which allows a direct application of the NSNMG method. Existence and uniqueness of discrete solution are proved and numerical examples illustrate the robustness of the scheme with respect to temperature, mesh size, and number of components.

Phase field models are a widely used approach to describe physical processes that are characterized by thin interfacial regions between large almost homogeneous domains. Important application areas of phase field models are transition processes of the state of matter and the separation of alloys. A fundamental property of these models is, that the transition and separation of phases is driven by a double-well potential with distinct minima for the different phases. Already the pioneering work of Cahn and Hilliard used a temperature dependent logarithmic potential that is differentiable with singular derivatives. If the temperature tends to zero it degenerates to the non-differentiable obstacle potential. The goal of this thesis is to develop methods for the efficient numerical solution of such equations that are also robust for nonsmooth potentials and anisotropic surface energies. These methods are derived for the Cahn-Hilliard equation that are prototypic for a multitude of such models. The main result of the thesis is the development of a fast iterative solver for nonlinear saddle point problems like the ones that arise from finite element discretization of Cahn-Hilliard equations. The solver relies on a reformulation of the problem as dual minimization problem whose energy functional is differentiable. The gradient of this functional turns out to be the nonlinear Schur complement of the saddle point problem. Generalized linearizations for the Schur complement are derived and used for a nonsmooth Newton method. Global convergence for this 'Schur Nonsmooth Newton' method and inexact versions is proved using the fact the equivalence to a descent method for the dual minimization problem. Each step of this method requires the solution of a nonlinear convex minimization problem. To tackle this problem the 'Truncated Nonsmooth Newton Multigrid' (TNNMG) method is developed. In contrast to other nonlinear multigrid methods the TNNMG method is significantly easier to implement and can also be applied to anisotropic problems while its convergence speed is in general comparable or sometimes even faster. Numerical examples show that the derived methods exhibit mesh independent convergence. Furthermore they turn out to be robust with respect to the temperature including the limiting case zero. The reason for this robustness is, that all methods do not rely on smoothness but on the inherent convex structure of the problems.

A frequent problem during numerical computations consists in the uncertainty of certain model parameters due to measuring errors or their high variability. In the last years, one could observe an increasing interest in the quantification of these uncertainties and their effects to the solution of numerical simulations; a powerful tool which has been proven to be an efficient approach in this context is the so-called polynomial chaos method which is based on a spectral decomposition of the covariance function of the uncertain parameters and a representation of the solution in a polynomial basis. The aim of this thesis is the application of this method to the Richards equation modeling groundwater flow in saturated and unsaturated porous media. The main difficulty consists in the saturation and the hydraulic conductivity appearing in the time derivative and in the spatial derivatives, since both depend nonlinearly on the pressure. Considering uncertain parameters like random initial and boundary conditions and, in particular, a random permeability leads to a stochastic variational inequality of second kind with obstacle conditions and a nonlinear convex functional as superposition operator. Considering variational inequalities in the context of uncertain parameters and the polynomial chaos method is new, and we start by deriving a weak formulation of the problem and approximating the parameters by a Karhunen-Loève expansion. The existence of a unique solution u in a tensor space can be proven for the time-discrete problem by reformulation as a convex minimization problem. We proceed by discretizing with finite elements and polynomial ansatz functions and by approximating the convex functional with Gaussian quadrature. The convergence of the solution of the discretized problem to the solution u is proved in a special case for a stochastic obstacle problem. Moreover, we perform numerical experiments to determine the discretization error. In the second part of this thesis, we develop an efficient numerical method to solve the discretized minimization problems. It is based on a global converging Block Gauß Seidel method and exploits a transformation which decouples the stochastic coefficients and connects the stochastic Galerkin with the stochastic collocation approach. This also allows us to establish a multigrid solver to accelerate the convergence. We conclude this thesis by demonstrating the power of our approach on a realistic example with lognormal permeability and exponential covariance.

The thesis presents a new model for the numerical simulation of the mechanics of the human knee. In this model bones are described using linear elasticity. Ligaments instead are modelled as one-dimensional Cosserat rods. The simulations give insight into the mechanical behavior of human joints. This can be helpful for a number of applications. For example, it is possible to estimate the long-term effect of certain surgical interventions. Also, the design of prosthetic devices can be improved. The main mathematical focus is on the correct formulation of the coupling conditions between one- and three- dimensional objects. Starting from the case of two three-dimensional objects, for which coupling conditions can be derived rigorously, conditions for the multidimensional case are formulated. A solution algorithm for this coupled problem is presented, and the existence of solutions is shown under certain symmetry assumptions. For the subproblems, large contact problems and minimization problems on Riemannian manifolds have to be solved. For both problems, robust and efficient numerical methods are introduced. Numerical experiments show the applicability for real-world problems.

The thesis presents a new method for the solution of saturated-unsaturated groundwater flow problems in heterogeneous porous media. Concretely, highly nonlinear degenerate elliptic problems arising from a certain time discretization of the Richards equation are the basis of this work. The problems are considered as homogeneous in subdomains where a single soil prevails and, therefore, the parameter functions do not depend on space. These nonlinearities, however, may jump across the interfaces between the subdomains and, thus, account for the heterogeneous setting of different soil types in different subdomains. As a consequence, non-overlapping domain decomposition problems in which subproblems are coupled via nonlinear transmission conditions are obtained. In this work these problems are solved without any linearization. By Kirchhoff transformation the homogeneous subproblems are transformed into convex minimization problems. Here, additional constraints like Signorini-type boundary conditions, which occur on seapage faces around lakes, can be taken into account. Finite elements are chosen for the space discretization, and convex analysis is applied as the solution theory. Finally, monotone multigrid methods provide efficient solvers which are robust with respect to degenerating soil parameters. In order to deal with the coupling of the homogeneous subproblems, nonlinear Dirichlet-Neumann and Robin methods are used. Here, the thesis provides new convergence results for these iterations applied to nonlinear elliptic problems in 1D as well as well- posedness results, which generalize existing linear theory. On the other hand, detailed numerical experiments demonstrate that the methods can also be applied successfully to problems in 2D. Finally, based on the artificial viscosity method, an upwind discretization with finite elements is developed in order to account for gravity. Hence, stability of the numerical solutions is obtained. In a closing numerical example the Richards equation is solved in 2D with four different soils and coupled to a surface water reservoir. The result demonstrates the applicability of the developed solution technique to a heterogeneous problem with realistic hydrological data.

The focus of this thesis ist the numerical computation of flow in special geometries dominated by jumps in the flow coefficients and large differences in the scales of the main flow pathes and the surrounding materials. These characteristics result in difficulties in the numerical computation of the modelling equations. By means of groundwaterflow in fractured porous media we present a hierarchical domain decomposition method for the numerical computation of flow. Under certain assumptions this new method converges independently of the fracture width, the refinement depth and the jump in the flow coefficient. The theoretical results are confirmed by practical computations of a model problem and a fracture network. Thus for a new class of completely overlapping domain decomposition methods multigrid efficiency is shown for a class of problems, for which so far no comparably theoretically validated method existed.

In this work, we consider the numerical simulation of contact problems. Since the numerical realization of contact problems is of high importance in many application areas, there is a strong demand for fast and reliable simulation method. We introduce and analyze a new nonlinear multigrid method for solving contact problems with and without friction. As it turns out, by means of our new method nonlinear contact problems can be solved with a computational amount comparable that of linear problems. In particular, in our numerical experiments we observe our method to be of optimal complexity. Moreover, since we do not use any regularization techniques, the computed discrete boundary stresses as well as the computed displacements turn out to be highly accurate. The new method is based on the succesive minimization of the associated energy functional in direction of properly choosen functions. We show the global convergence of our method and give several numerical examples in two and three space dimensions, illustrating the robustness and the performance of the method. In addition to the theoretical analysis, the method has been implemented in an object oriented way. We explain the concepts of our implementation and show the flexibility of our approach by deriving a nonlinear algebraic multigrid method. To include frictional effects, we use a discrete fixed point iteration. As a faster alternative, also a Gauss-Seidel like iteration scheme is proposed. Both methods are compared in numerical examples. The resulting nonlinear algorithm turns out to be fast and reliable. Finally, we consider the case of contact between elastic bodies. Here, the information transfer at the interface is realized by means of non conforming domain decomposition methods (mortar methods). This gives rise to a non-linear Dirichlet Neumann Algorithm.