Log In Sign Up

A stable and optimally convergent LaTIn-Cut Finite Element Method for multiple unilateral contact problems

In this paper, we propose a novel unfitted finite element method for the simulation of multiple body contact. The computational mesh is generated independently of the geometry of the interacting solids, which can be arbitrarily complex. The key novelty of the approach is the combination of elements of the CutFEM technology, namely the enrichment of the solution field via the definition of overlapping fictitious domains with a dedicated penalty-type regularisation of discrete operators, and the LaTIn hybrid-mixed formulation of complex interface conditions. Furthermore, the novel P1-P1 discretisation scheme that we propose for the unfitted LaTIn solver is shown to be stable, robust and optimally convergent with mesh refinement. Finally, the paper introduces a high-performance 3D level-set/CutFEM framework for the versatile and robust solution of contact problems involving multiple bodies of complex geometries, with more than two bodies interacting at a single point.


page 21

page 24

page 25

page 30

page 33

page 34

page 35

page 36


A high order unfitted finite element method for time-Harmonic Maxwell interface problems

In this paper, we propose a high order unfitted finite element method fo...

An extended mixed finite element method for elliptic interface problems

In this paper, we propose an extended mixed finite element method for el...

An Enriched Immersed Finite Element Method for Interface Problems with Nonhomogeneous Jump Conditions

This article presents and analyzes a p^th-degree immersed finite element...

A cut finite element method for a model of pressure in fractured media

We develop a robust cut finite element method for a model of diffusion i...

Immersed Virtual Element Methods for Elliptic Interface Problems

This article presents an immersed virtual element method for solving a c...

BFEMP: Interpenetration-Free MPM-FEM Coupling with Barrier Contact

This paper introduces BFEMP, a new approach for monolithically coupling ...

1 Introduction

Unfitted or non-conforming finite element methods uncouple the description of the geometry from the representation of the solution field itself. Typically, the geometry of the computational domain is projected over a regular background grid. In this setting, boundaries or interfaces between objects cut through elements of the corresponding mesh. Non-conforming methods are attractive for applications where coupling analysis codes and third party meshing libraries is either impractical, numerically expensive and/or prone to errors. Nonetheless, a number of specific challenges needs to be addressed for unfitted numerical solutions to be computable. Firstly, integrals need to be calculated over cut elements, which requires specialised numerical quadratures. Secondly, unfitted approaches require stabilisation. This is because combinations of degrees of freedom may be poorly controlled in regions where contributions from cut elements to integral forms are small. Finally, in the context of multiple interacting materials, enrichment of the finite element solution space is required in order to allow for numerical jumps or kinks to develop over embedded interfaces. Failure to do so may severely impair the convergence rate of the unfitted finite element solver.

Over the last two decades, several encompassing frameworks have been developed to provide guidance for the development of unfitted finite element solvers. The eXtended Finite Element Method (XFEM) [36, 18, 7, 21] relies on the Partition of Unity Method [34] to enrich the approximation space. Within this popular framework, non-conforming discontinuities can be introduced in an elegant and robust manner. Field discontinuities can also be enabled by making use of Embedded Discontinuity Methods [39, 37, 31], whereby discontinuous shape functions are directly defined at the element level. The CutFEM framework [22, 10], which is fundamental to the developments reported in this paper, uses an alternative concept to introduce non-conforming field discontinuities, and to perform computations involving non-conforming geometrical domains. In the context of multiple interacting materials, CutFEM solvers enrich the solution space by overlaying (i.e. doubling or tripling) cut elements. Subsequently, penalty techniques such as the Nitsche method are employed to disallow unwanted discontinuities, for instance to represent jumps in fluxes whilst preserving the continuity of the primal field [22]. CutFEM solvers are usually equipped with a ghost-penalty stabilisation technique [9], which results in optimal convergence rates whilst ensuring that the condition number of the system matrix degrades at the same rate as that of the conforming finite element method when the mesh is refined.

Extending unfitted finite element methods to nonlinear problems, and especially to problems involving complex interface behaviour, remains a challenge that needs to be addressed on a case-by-case basis. The present paper focuses on problems involving unilateral contact between elastic solids. This class of problems is important to engineers, and they are typically encountered when modelling joints between components in large-scale engineering assemblies. Numerically, unilateral contact is notoriously difficulty to treat, even in the context of fitted finite element approaches. References on the matter include [45, 30, 42, 1, 33]. Of particular interest to us are LaTIn-based solvers [13, 28, 8, 32, 29, 2], which treat unilateral contact through a two-field formulation of the interface equations, and resolve these fields iteratively using a dedicated two-stage solver. A number of researchers have demonstrated the versatility of LaTIn solvers, owing to the fact that the predictor stage (so-called “linear stage”), is weakly dependent on the type of interface conditions, and completely independent of the current interface state. As a result, switching between nonlinear interface conditions is relatively painless [24]. Moreover, dedicated space-time Model Order Reduction techniques may take advantage of the predictor invariance to reduce the cost of time-dependent or multi-query computations [29, 24, 19]. However, a noteworthy issue related to the hybrid-mixed LaTIn formulation is the difficulty of discretising displacement and traction interface fields in a consistent manner. Whilst a piecewise constant discretisation of interface fields seems to be physically sound when associated with P1 finite element spaces for the bulk equations, such a discretisation scheme allows for the onset of numerically uncontrolled deformation modes whose associated amplitude creep up with increasing LaTIn iteration count, and eventually alter the convergence rate of the finite element solver. A claimed remedy to this numerical plague is to introduce additional degrees of freedom for the bulk fields, by means of either h-refinement or p-refinement. An alternative to this numerically expensive procedure has recently been proposed in [20], as discussed later on. Although both approaches seem to regularise the contact problem efficiently, the effect on the convergence of the mixed finite element solver with mesh refinement has not been thoroughly studied.

Over the years, several attempts to develop stable unfitted finite element solvers for unilateral contact have been met with success. XFEM-based solutions have been proposed in [16, 25, 40, 17, 32, 20, 38]. Of particular importance to us are LaTIn-based unfitted solvers, as the LaTIn is also a fundamental technological element for the present contribution. The first of these solvers was proposed in [16] to solve crack propagation problems (i.e. the historically favoured playground for XFEM). In [40], the XFEM-LaTIn idea was extended in order to simulate fatigue crack propagation. In [32], the authors compared a LaTIn formulation to a penalty/Newton approach, and claimed the superiority of the latter in terms of convergence rate and stability. Finally, the authors of [20] proposed a regularised, LaTIn-type formulation of contact for XFEM. The main idea of the approach is to relax the kinematic continuity condition between the bulk primal field and its interface trace and to enforce it weakly via a regularised mortar method.

In contrast to XFEM-based developments, the CutFEM framework is relatively underdeveloped in the context of unfitted unilateral contact. We mention the recent work described in [14, 15], which extends the Nitsche-type treatment of interface conditions that is usually employed in CutFEM solvers to the case of unilateral contact. The resulting fitted finite element solver is proven to be optimally convergent, stable and numerically efficient. A very first nonconforming CutFEM strategy for unilateral contact is discussed in [11], whereby the contact conditions are treated by an augmented Lagrange multiplier formulation, and the discrete equations are stabilised by a ghost-penalty regularisation.

The present paper proposes a novel LaTIn-CutFEM solver for unilateral contact between multiple interacting solid bodies. The solid bodies will be described by linearised elasticity, which implies that the location of contact interfaces is known a priori. Our algorithm relies on the traditional CutFEM enrichment strategy [22]. However, contact conditions over nonconforming interfaces will be enforced via the LaTIn method. The resulting algorithm is a two-stage solver, whereby the linear stage (i.e. linear predictor) consists in solving a series of independent linear problems for each of the solid bodies, whilst the local stage (i.e. interface corrector) is nonlinear. The method naturally inherits the coarse-grain parallel characteristics of LaTIn-based domain decomposition methods [27, 28, 24], together with the previously mentioned versatility of LaTIn solvers with regards to the type of interface conditions being dealt with (i.e. unilateral contact with Coulomb friction [8, 12], frictional sliding [5, 6], cohesive fracture models with transition to unilateral frictionless contact [24]).

Two algorithmic elements of the proposed solver require particular care. Firstly, the LaTIn-CutFEM scheme needs to be stabilised. We regularise each of the subproblems of the linear stage using a ghost penalty technique [9]. We will show that, when using the appropriate scaling for the ghost penalty terms, the condition number of the system matrices “degrades optimally” with mesh refinement (i.e. it degrades at the same rate as that of the conforming finite element method). Secondly, we need to be careful when discretising the contact interface fields. As mentioned before, a naive discretisation of the LaTIn hybrid-mixed formulation leads to the onset of spurious interface solutions. We will show that the P1-P1 discretisation scheme that we propose (i.e. P1 for the bulk primal field, P1 for both its trace and the trace of its normal flux) yields optimal convergence rates with mesh refinement, and eliminates the onset of instabilities. Our scheme “nonlocalises” the local stage of the LaTIn, as it relies on a two-scale treatment of the interface fields and nonlinear equations. We will also show that the nonlinear solver of the LaTIn converges fast, in the sense that it quickly leads to iterates whereby the discretisation error is larger than the algorithmic error, typically within twenty iterations. To our best knowledge, these convergence results have never been reported in the context of LaTIn solvers, even in the context of conforming geometries.

Our formal developments are accompanied by a high-performance computer implementation. The core of our implementation is the finite element C++/Python library FEniCS [4], which, in particular, proposes a range of high-level tools to rapidly developed finite element solvers. The CutFEM C++ library, LibCutFEM, developed by S. Claus and A. Massing, and partially described in [10], defines additional tools that are specific to unfitted finite elements. This library forms the basis for the LaTIn-CutFEM code. We will present results in two and three space dimensions, which is made seamless by the FEniCS framework. Some of our examples exhibit a high level of geometric complexity, which is handled through a robust implementation of multiple level sets. In particular, we will study the convergence of contact problems exhibiting multiple interacting bodies (i.e. more than two) at a single point, (see other contributions on the topic, in a different context, in [5, 6]).

This paper is organised as follows. The LaTIn hybrid-mixed formulation of contact is presented in a continuous setting in Section 2, together with the associated two-stage nonlinear algorithm. The non-conforming regularised LaTIn-CutFEM formulation is presented in Section 3. Section 4 is dedicated to numerical investigations, where we demonstrate the optimal convergence and stability properties of the proposed framework. We also showcase the versatility of our high-performance computer implementation.

2 LaTIn method for unilateral contact: the continuous setting

2.1 Domain decomposition

We consider a problem of linear elasticity over a 2 or 3-dimensional spatial domain composed of non-overlapping subdomains , where . The domain decomposition satisfies . Let us denote the set of all subdomains by . Subdomains may be multiply-connected. The boundary of subdomain is denoted by , for any , whilst denotes the boundary of computational domain . Furthermore, we define interface domains for any such that the two indexes satisfy and the result of intersection is of non-zero measure. The union of these interfaces is denoted by . The set of pairs of indexes such that is defined as a valid interface is denoted by , and the corresponding set of interfaces is denotes by . We also denote the union of all the interfaces of subdomain by , for any subdomain index . We denote the outward normal pointing from to on by . We further note that we have , and we denote the outward normal to over by .

Figure 1: Schematic representation of a composite material made of three distinct phases.

As an illustration for this set of definitions, consider the case of the multiphase domain depicted in Figure 1. The cubic domain comprises three subdomains , and . There are three interfaces between the domains: , and . Therefore, the interface of individual subdomains reads as , and .

As further detailed in the next subsection, subdomain field variables will be governed by the equations of linear elasticity, while contact equations to be solved over interfaces will couple subdomain problems to one another. The LaTIn framework [26] is particularly well suited to the treatment of such problems where complex nonlinear equations need to be solved over embedded interfaces [8, 29, 3]. In this particular context, the LaTIn method relies on two main ingredients:

  • a hybrid description of bulk and interface fields, and a mixed representation of the interface unknowns,

  • a two-step iterative solution algorithm for the solution of the global nonlinear problem formulated in its hybrid-mixed form.

Both these steps are described below.

2.2 Hybrid-Mixed LaTIn formulation of unilateral contact problems

Hybrid fields

For every subdomain , where , we introduce a displacement field , where (2D problem) or (3D problem).
For every domain interface with , we define two field variables:

  • a force field that represents a density of tractions (i.e. flux) applied to subdomain by adjacent subdomains through interface .

  • a displacement field that represents the trace of over .

As detailed below, will be governed by linear elasticity, while over an interface between two domains and , where , interface models will provide a link between the restriction of fields , to and the restriction of fields , to .

Bulk equations and boundary data

For every subdomain , , the displacement field is required to be sufficiently regular (continuous in particular), and to satisfy the following equations:

  • Static equilibrium equation and constitutive relation:
    For all , the displacement field satisfies


    In equation (1), is a body force and is the outer normal of domain . The stress function is defined as , where the strain is the symmetric part of the displacement gradient, i.e.