A contact algorithm for voxel-based meshes using an implicit boundary representation

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 352 (2019) 276–299 www.elsevier.com/locate/cma

A contact algorithm for voxel-based meshes using an implicit boundary representation Alexander Leichnera ,∗, Heiko Andräa , Bernd Simeonb a

Fraunhofer Institute for Industrial Mathematics, Department of Flow and Material Simulation, Kaiserslautern, Germany b TU Kaiserslautern, Department of Mathematics, Germany Received 1 October 2018; received in revised form 2 April 2019; accepted 3 April 2019 Available online 24 April 2019

Highlights • • • • •

A numerical method on structured meshes for contact in 2D and 3D is presented The concept of intermediate surfaces as contact surfaces is considered. Implicit surface dynamics is used for contact detection and kinematics. A modified Nitsche’s approach enforces the contact constraints. Geometrical approximations are based on implicit surface representations.

Abstract Numerical methods operating on structured grids have become popular since they offer good run time performance and are able to process directly voxel-based digital data from image recordings. Hence, the general framework of these fast solvers presupposes an unfitted boundary approximation avoiding complicated meshing of bodies. This allows an efficient handling with geometrical issues. Nevertheless, contacts between deformable solids are hard to deal with in the presence of this boundary representation. For this difficulty we suggest the usage of an implicit boundary representation combined with a modified saddle point formulation, resembling Nitsche’s approach. Both ideas give an elegant approach for discretizing the contact terms and enable a simple contact detection. Moreover, we suggest an intermediate surface as new reference contact area which fits well into our proposed method. In the end we present numerical results and analyze the accuracy and convergence rate. Furthermore, we demonstrate the current application range of our approach for problems with multiple contacts. In this paper we focus only on frictionless contact with small deformations. c 2019 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ⃝ (http://creativecommons.org/licenses/by/4.0/). Keywords: Contact mechanics; Regular grid discretization; Implicit geometry; Level set method; Nitsche’s method; Intermediate contact surface

1. Motivation and introduction As mentioned, digital processing of image scans [1] has become an attractive alternative to experiments on real samples. In our case we focus on the simulation of structural analysis of solids. Appropriate scans are then obtained ∗

Corresponding author. E-mail address: [email protected] (A. Leichner).

https://doi.org/10.1016/j.cma.2019.04.008 c 2019 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons. 0045-7825/⃝ org/licenses/by/4.0/).

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

277

at micro scale and are used to determine macroscopic material behavior or quality. Materials of interest are mostly composites, rocks and metal alloys. Beyond that, characteristics of fabrics or textiles are depending on their structure at micro scale, too, and can be engineered to satisfy certain purposes. In contrast to composites, textiles, built by fiber networks, have void as matrix or background material and therefore allow many possibilities for contact. Hence, the occurrence of collisions between fibers is inevitable in stress analysis. Together with a voxel-based discretization, simulation of contact means then a challenging task. Certainly, the main problem is to find the actual collisions in a multibody setting. The search can be classified into two problems: First, we have to determine which bodies are in contact and second, where exactly on their surfaces the penetrations have to be prevented. Additionally, the reaction forces must be computed precisely in order to achieve realistic stress results. Today, the finite element method [2] can cope with such contact problems. However, appropriately fitted meshes as geometrical approximations are usually required. Standardly, the contact conditions are then enforced by means of Lagrange multipliers [3] or via the penalty method [4]. The resulting algorithms prerequire that the contact area is currently known. Hence, these methods apply a detection for contact in a superordinate loop which terminates when the correct area is found. Methods like the primal–dual mortar approach [5], the scalable FETI method [6] or the monotone multigrid method [7] demonstrated that they are powerful tools for contact problems in computational solid mechanics. Yet, the generation of fitted meshes for standard finite element methods can be very expensive in computational time. Therefore, the development of methods which use an unfitted mesh has been carried out in the past [8]. Various names have been given these methods, like extended finite element [9], fictitious domain [10] or immersed boundary [11] method. These novel approaches can be applied directly on voxel data from CT scans or digital images in general. Moreover, transformation techniques like the FFT [12] are then applicable to problems discretized on structured grids benefiting the solver’s runtime. Admittedly, the surface representation by uniform cells is poor and does not allow a continuous description of geometrical properties, like the outward normal. In order to remedy that, we suggest the usage of an implicit boundary approximation by level set functions [13,14]. These have two great advantages: First, appropriate methods are already available for the manipulation of those. And second, they are most efficiently used on structured grids, which suits our purposes best. Beyond that, we exploit non-standard methods for handling the contact terms [15,16]. The approach, we follow here, is inspired by Nitsche’s method [17] and is considered for the enforcement of the contact conditions. Together with the already mentioned level set description, this approach fits well into the framework of contact treatment on now implicitly defined geometries. In fact, we will see that our method resembles the standard penalty method [18]. Anyhow, the innovation is based on the new applications of level set methods in the context of contact mechanics. The combination of both ideas, the level set and Nitsche’s approach, results in an efficient procedure for contact problems placed in a structured grid. This is due to the fast and accurate contact detection by means of level sets. Moreover, we establish a new contact reference surface which describes contact more intuitively than usual ways and yields results which are competitive with standard procedures. The outline of this paper is now explained in the following: In Section 2 we introduce the governing equations for elasticity and the balance of momentum for of small deformations. Then we extend the model by certain constraints and derive the weak form for contact problems. After that, we introduce in Section 3 the level set function and explain briefly its background from a mathematical point of view. Then we turn our attention to the usage of level sets for contact detection, the description of gaps and the determination of the correct contact surface. Afterwards, Nitsche’s method is presented via a perturbed Lagrangian formulation and the connection to level set functions will be emphasized. In Section 4, the establishment of the total linear system will be demonstrated in the set-up of a regular grid discretization. In particular, we show how to reconstruct (continuous) surfaces and volumes within uniform cells. The last Section 5 provides numerical tests of our proposed approaches. There, we present solutions to standard problems from literature, like the Hertz problem and compare our results with the theory. Finally, we show more numerical experiments with complex geometries in three dimensions. 2. Problem setting We begin with some notational preliminaries: At first we define a set I := {1, . . . , NΩ } for indexing the solids Ω (i) which are finite, simply connected domains with Lipschitz boundaries in Rd where d ∈ {2, 3}. Independent

278

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299 (i)

of d, we call the boundary ∂Ω (i) ⊂ Ω also surfaces for all i ∈ I. By the definition of the solids, there exists a finite bounding box B containing all Ω (i) where i = 1, . . . , NΩ . Moreover, we assume that Ω (i) ∩ Ω ( j) = ∅ for all (i, j) ∈ I × I, such that we may write: ⋃ ˙ (i) Ω =: Ω ⊂ B ⊂ Rd . i∈I

2.1. Strong formulation In this contribution we are primarily interested in boundary value problems (BVPs) for quasi-static, small deformations of solids. We define the deformation d as a bijective map from the initial, undeformed state Ω to the deformed configuration Ω˘ , i.e. d : Ω → Ω˘ , d(x) = x + u(x),

(1)

where u(x) is the displacement at x ∈ Ω . Unless otherwise mentioned, index superscripts in brackets denote restrictions on the corresponding domain, e.g. u(i) := u|Ω (i) . We consider only linear elastic, isotropic materials for which we give the stress–strain relation: σ (u(i) ) = C(i) : ε(u(i) ), ) 1 ( (i) ε(u(i) ) := ∇u + (∇u(i) )T , 2 where ε is the linearized strain tensor and σ denotes the Cauchy stress tensor. The material tensor C(i) describes the homogeneous material characteristics of the body Ω (i) , and is defined via the Young’s modulus E (i) and Poisson’s ratio ν (i) . For each Ω (i) we introduce the BVP which consists of the balance of momentum as well as Dirichlet and Neumann boundary conditions. So the strong form of this problem reads: Find for all i ∈ I a displacement vector u(i) : Ω (i) → Rd such that −∇ · σ (u(i) ) = b(i) in Ω (i) , u(i) = uˆ (i) on Γd(i) , σ (u ) · n (i)

(i)

= ˆt

(i)

on

(2)

Γn(i) ,

where Γn(i) ⊂ ∂Ω (i) and Γd(i) ⊂ ∂Ω (i) denote the Neumann and nonempty Dirichlet boundary surface, respectively. ˙ Γd(i) = ∂Ω (i) , hence Γn(i) ∩ Γd(i) = ∅. The functions uˆ (i) and ˆt (i) are prescribed Moreover, we assume that Γn(i) ∪ along the corresponding boundaries and n(i) is the outward normal on ∂Ω (i) . The boundary conditions ensure that solving the partial differential equation is a well-posed problem and results in a unique solution [19]. In case contact occurs on ∂Ω (i) , the corresponding possible contact area is denoted with Γc(i) ⊂ ∂Ω (i) . In addition to the BVP (2), frictionless contact problems require the fulfillment of the Karush–Kuhn–Tucker (KKT) conditions on each contact surface: ϱn (d (i) ) ≥ 0, σn (u(i) ) ≤ 0, σn (u(i) ) ϱn (d (i) ) = 0 on Γc(i) .

(3)

These conditions are well known from problems of nonlinear optimization [20] where contact problems can also be placed. The terms σn := n · σ · n and ϱn (u) are called contact stress (or pressure) and global gap function, respectively. The latter is defined as continuous function in terms of the deformations of both bodies ϱn (d (i) (x)) := min gn (d (i) (x), d ( j) (P ( j) (x))). j∈I\{i}

(4)

The local gap function gn in (4) returns the signed distance between x ∈ ∂Ω (i) and P ( j) (x) ∈ ∂Ω ( j) under consideration of the corresponding deformations. For the undeformed setting it can be defined as gn : ∂Ω × ∂Ω → R, gn (x, y) = sign(n( y)T · (x − y)) ∥x − y∥.

(5)

The deformed case is treated in Section 3.4. In definition (4) P ( j) is then the closest point projection: P ( j) : Ω → ∂Ω ( j) , P ( j) (x) = arg min ∥x − y∥. y∈∂Ω ( j)

(6)

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

279

Fig. 1. Depiction of some contact problems.

For simplification, we also write ϱn (u, x) instead of ϱn (d, x) in this paper. In order to return to the KKT conditions (3), we briefly interpret them: The first inequality describes the non-penetration condition for a body with index i w.r.t. all j ∈ I \ {i}, i.e. gn (u(i) , u( j) ) ≥ 0. The second condition ensures that only inward contact forces act over the contact area. Finally, the last condition is also referred to as complementary slackness since the contact pressure is non-zero, if and only if, the global gap function vanishes. In summary, the strong form of mechanical contact problems is the BVP in (2) subject to the conditions in (3). 2.2. Weak formulation Now we reformulate the contact problem in strong form (2)–(3) as a saddle point problem. For this purpose, we define the following sets { } [ ]d V := u ∈ H 1 (Ω ) | u = uˆ on Γd , (7) Λ N := {λn ∈ L 2 (Γc ) | λn ≤ 0},

(8)

where V is the space of admissible displacements, and Λ N is the convex cone [4,21] for the Lagrange multipliers λn . For the sake of a simple problem formulation, we consider homogeneous Dirichlet conditions in this section, hence uˆ (i) = 0 for all i ∈ I. Note that V is decomposed body-wise and the cone Λ N contact-wise, i.e. for all Ω (i) ( j) and Ω ( j) in contact λn is acting on Γc(i) and Γc , see Fig. 1a. The actual determination of bodies in contact is then explained in Section 3.3. Similar to Ω , the surface Γc might be understood as a union of all Γc(i) . Finally, the weak formulation of the contact problem (2)–(3) reads in global sense, i.e. without index notation: Find (u, λn ) ∈ V × Λ N , such that a(u, v) + b(v, λn ) = f (v), b(u, µn ) = 0,

∀v ∈ V, ∀µn ∈ Λ N ,

where the bilinear forms a(·, ·) and b(·, ·) and linear operator f (·) are defined again in global sense: ∫ ∫ ∫ ˆt · v dγ , a(u, v) := σ (u) : ε(v) dx, f (v) := b · v dx + Ω Ω Γn ∫ b(v, λn ) := ϱn (v)λn dγ .

(9)

(10)

Γc

In general, the saddle point problem (9) can be written depending on the Lagrangian L which is defined as the sum of potential energy and the contact contribution πc : 1 L(u, λn ) := a(u, u) − f (u) + πc (u, λn ). (11) 2

280

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

Note that the weak formulation (9) originates from a optimization problem and is the result of a variational treatment of the Lagrangian L via Gˆateaux differentiation. However, we omit to give these details and refer to [4,19]. 3. Implicit geometry and Nitsche’s method for contact Since we want to avoid the generation of boundary-aligned meshes, we consider a different way of representing surfaces in space. In common Lagrangian methods, e.g the standard finite element method, positions of nodes which build surface elements are precisely known. Here, we use an implicit boundary representation such that boundaries are identified as surface contours within the bounding box B. The surface dynamics is then described from the Eulerian point of view. Furthermore, we present an approach by Nitsche for treating the contact conditions (3). Our proposed modification resembles the standard penalty method with an additional constant term. Nonetheless, this approach is still related to the Lagrangian functional in (11). Indeed, the direct involvement of the level set idea in our proposed Nitsche’s approach is clearly recognizable. 3.1. Level set function A level set function Φ [13,14] is a function characterized by the shape of a domain Ω where in this contribution it is confined to a bounding box B ⊂ Rd . In general, it is described as ⎧ ⎪ ⎨< 0, x ∈ Ω , (12) Φ : B → R, Φ(x) = 0, x ∈ ∂Ω , ⎪ ⎩ > 0, x ∈ B \ Ω . Accordingly, an isosurface or contour of value Ψ ∈ R is defined as the set L Ψ (Φ) := {x ∈ B | Φ(x) = Ψ }. The surface given by L Ψ (Φ) is referred to as surface contour or level set. Hence, the boundary of the domain Ω equals the zero level set ∂Ω = L 0 (Φ). Indeed, we use this fact in order to reconstruct boundaries of bodies in a structured grid. Furthermore, level set functions enable the computations of the outward normal n on ∂Ω : n(x) = ∇Φ(x)/∥∇Φ(x)∥,

x ∈ ∂Ω .

(13)

In the following sections we use an abbreviation for the set of points which satisfy certain conditions for level set functions. For locations inside Ω , we write: {Φ < 0} := {x ∈ B | Φ(x) < 0} = Ω .

(14)

Analogously, the complement can be identified B \ Ω as {Φ ≥ 0}. 3.2. Implicit surface dynamics Level sets are able to represent surfaces, even after deformations. This becomes very important for tracking contact and boundary surfaces. For this purpose, we focus on a certain position x ∈ B. By introducing the (fictitious) time variable t ∈ T, the level set through x is able to evolve. The finite interval T := [tα , tω ] covers the time frame from initial configuration Ω till deformed state Ω˘ in equilibrium. For the sake of simplicity, we assume that tα = 0 and tω = 1 because we deal with quasi-static processes. Therefore, we redefine d as a function of space and time, in particular d(x, t) = x + u(x)t. The level set function Φ becomes time-dependent, too, i.e. Φ = Φ(x, t). Clearly, the motion towards the deformed body is governed by the deformation field d from (2) such that it holds by means of (1) d(Ω , tω ) = Ω˘ = {Φ(·, tω ) < 0}. We consider the unsteady level set L Ψx (t) (Φ(·, t)) where Ψ x (t) := Φ(d(x, t), t) and further denote Φα := Φ(·, tα ) and Φω := Φ(·, tω ). The motion of the initial level set L Ψx (0) (Φα ) is given as an evolution in time: L Ψx (∆t) (Φ(·, ∆t)) → L Ψx (1) (Φω ), for ∆t → 1.

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

281

In general, level sets are not steady in time if u ̸= 0. The level set value Ψ x (t), however, is constant for all t ∈ T since it tracks the front L Ψx (Φ(·, t)). Hence, we obtain by differentiating Ψ x (t) w.r.t. the time t the partial differential equation (PDE) ∂t Φ + ∇Φ · u = 0,

(15)

which is called level set equation and is similar to the transport equation. Since Φ(·, tα ) = Φα (·), surface tracking can be regarded as an initial value problem. The established methods [13,14] for solving (15) use upwind schemes and achieve different order of accuracy. In general, these methods are called level set methods (LSM) and compute solutions to Hamilton–Jacobi equations. For our research the standard level set method is sufficient, although many improvements were presented in the last years [22–25]. The involvement of the displacements as motion field in (15) sums up a direct connection between the methods operating from the Lagrangian point of view and those from the Eulerian one. Various extensions and modifications of level set methods are nowadays established and well-known: For instance, Sethian presented in [26] an efficient algorithm, the fast marching method (FMM), for the stationary version of (15): ∥∇Φ∥F = 1

in Ω ,

(16)

where F is the velocity in normal direction (F = u · n). In our algorithm we apply the FMM mainly for the initialization [27] of gn (x, P(x)) where P is the closest point projection to ∂Ω and x ∈ B. The definitions (5) and (6) establish a direct relation between the signed distance and the gap function. Indeed, if F = 1, the FMM applied to any level set function of the domain Ω will return the signed distance function gn . In general, the LSM does not conserve the signed distance property, i.e. ∥∇Φα ∥ = 1 does not automatically apply to Φω . Hence, the FMM can be thought of as a reinitialization method [28] for Φω , too. Apart from deforming level sets, the PDE in (15) might be rewritten, such that it treats extrapolation problems, too. The corresponding initial value problem for extending an arbitrary quantity q off a zero isosurface L 0 (Φ) reads then { ∂q ∇Φ q0 (x), if x ∈ Ω , + · ∇q = 0, q(x, tα ) = (17) ∂t ∥∇Φ∥ 0, if x ∈ B \ Ω . For our purposes, we seek the steady-state solution of this unsteady problem [29]. Considering the time interval T, the solution qω at tω represents the extrapolated counterpart of q0 , i.e. in terms of projections qω (x) = q0 (P(x)). Note that we are only interested in a solution within the exterior of Ω , hence x ∈ B \ Ω . In particular, the quantity q0 is extrapolated constantly in rays given by the normal n = ∇Φ/∥∇Φ∥ along the surface (see Fig. 3). Similar to (15), the problem (17) can be solved by the LSM, too. An alternative method solving the stationary version of (17) can be found in [30]. The extrapolation approach is crucial for the continuation of our Eulerian approach: Note that u is defined only in Ω , but the level sets might evolve into B \ Ω due to the formulation in (15). For a graphical illustration, see Fig. 2b. Thus, extending u according to (17), the level set method applied to problem (15) can be continued in B \ Ω. As a final remark we want to note that evaluations and operations in the context of level set functions do not have to be carried out in the whole domain B. On the contrary, it is sufficient to restrict the application of all mentioned methods in the immediate vicinity of boundaries or narrow bands [31]: Nβ (Φ) = {x ∈ B | −β ≤ Φ(x) ≤ β}.

(18)

This will allow the algorithm to benefit from minimal memory consumption and optimal performance [26]. Besides, level set or distance data in the far interior or exterior of bodies are rarely interesting in practice. The width β is set to a few multiples of the grid size h. If extrapolation is considered, the bandwidth can be increased as large as needed. 3.3. Description of contact by means of level sets In this section, we show how implicit boundary descriptions can be advantageously used for both the contact detection and an alternative description of the gap function gn in (5). In principle, the following ideas are extendable

282

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

Fig. 2. Level set functions and applications in Cartesian grids.

Fig. 3. Extrapolation by LSM.

to problems with large deformations. For this purpose, we first explain how level set functions support the search for sufficiently close bodies which are likely to get into contact in a multibody setting. Second, a novel approach for determining the contact surface is introduced. And finally, we show how the gap function gn can be reinterpreted as a level set function and how to linearize it. In this section, we label the level set function for domain Ω (i) with Φ (i) and require that Φ (i) is a signed distance function, i.e. it additionally satisfies ∥∇Φ (i) ∥ = 1, for all i ∈ I. 3.3.1. Contact detection In order to deal with contact between Ω (i) and Ω ( j) , we aim to find an environment in close neighborhood to both objects, see Fig. 4a. Since in this setting Nβ (Φ (i) ) and Nβ (Φ ( j) ) contain each body’s surfaces, a simple approach would be to determine the set intersection Nβ (Φ (i) ) ∩ Nβ (Φ ( j) ). Then, it captures for an appropriate choice of β ( j) the right contact surfaces, Γc(i) and Γc , and the proximity zone between those bodies. This procedure is sufficient for detecting all possible contacts. A more accurate alternative, however, can be achieved by a dilation of level set functions. This approach is based on the union Ω (i) ∪ Ω ( j) , whose level set function is given by (i j)

Φmin = min(Φ (i) , Φ ( j) ). (i j)

A shift given by Φmin − ξ for ξ > 0, and a subsequent reinitialization [28] result in the desired dilated signed (i j) distance function Φ˜ min . This dilation supports a contact search by “closing” small gaps between bodies: Shifting

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

283

Fig. 4. Detecting close bodies.

(i j) (i j) back by ξ , i.e. Φ˜ min + ξ , yields the original level set function Φmin , apart from the region in the vicinity to both (i j) bodies. There, we will find that Φ˜ min + ξ has a negative sign which can be used as an indicator for contact, see Fig. 4b. In particular, for an appropriate ξ the following set (i j)

(i j)

Ωc(i j) := {Φ˜ min + ξ ≤ 0} ∩ {Φmin > 0}

(19) (i j)

is a better alternative to the intersection of narrow bands since Ωc captures the possible contact surfaces more efficiently. To recognize this, check the illustration in Fig. 4b against that in Fig. 4a where the narrow band intersection might be excessively large compared to the size of the contact surfaces. From the point of (i j) (i j) implementation, the above set can be easily determined by evaluating the level set functions Φ˜ min + ξ and Φmin . Focusing on narrow bands, this set determination can be narrowed down to operations close to surfaces. This should prevent wasting performance for the contact search in the deep interiors of the domains Ω (i) or outside of them. (i j) In regard to complexity, the computational effort for finding all Ωc depends linearly on the amount of bodies NΩ , assuming their narrow bands have similar volumes, see Algorithm 1. In comparison, the naive approach, checking each body against all remaining ones for collisions, has quadratic effort. 3.3.2. Determination of the contact surface The choice for a proper contact surface is commonly made by a master–slave relation [5] between contacting bodies: Once a side is chosen as master, the closest point projection is defined as a map from the slave side onto the former. The integration of the corresponding contact terms is usually performed on the slave’s surface. This is a biased formulation of contact, however, an unbiased extension would alternately consider both bodies as slave and master. The resulting terms for the contact contribution are then averaged in the final variational problem, see [32–34]. Instead of referring individual surfaces to each other, we introduce an artificial surface as a neutral alternative. This choice solves then the issue about the bias against slave or master side. With regard to later discretization, this surface has to be constructed fast and must fit in the previous presented contact detection approach. For this reason, we propose an intermediate hypersurface like in [35–37], fulfilling all mentioned criteria. First, we define an intermediate level set function between Φ (i) and Φ ( j) by Φ (i) − Φ ( j) . (20) 2 This level set function has a special property: From the above definition one can immediately conclude that points (i j) x ∈ B satisfying Φint (x) = 0 have equal distance to both ∂Ω (i) and ∂Ω ( j) . The corresponding zero level set (i j)

Φint :=

(i j)

(i j)

Γint := L 0 (Φint )

(21)

284

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

Algorithm 1 Initialization and Contact Detection 1: 2:

for all i ∈ I do Initialize Φ (i) (i)

▷ A possibility is to use analytical formulas (i)

3:

Φ

4:

Create narrow band Nβ (Φ (i) )

← FMM(Φ )

▷ Apply fast marching method from (16)

5:

end for

6:

Set IC ← ∅ and uh ← 0

7:

for all i ∈ I do

8:

▷ Bandwidth β is few grid sizes large

for all j ∈ I with Nβ (Φ (i) ) ∩ Nβ (Φ ( j) ) ̸= ∅ do

9:

if ( j, i) ∈ / IC then

10:

(i j) Ωc

Construct

▷ This means new contact is found according to (19)

▷

12:

Φ int(i j) ← FMM((Φ (i) − Φ ( j) )/2)|Ω (i j)

13:

for l = i, j do p ← ∇Φ

14:

u(l)

⊂ Nβ (Φ (i) ) ∩ Nβ (Φ ( j) ) ▷ See definition (21)

(i j)

h|Ωc

▷ For initialization of active set, see Alg. 2 (l)

(i j)

▷ Provoking penetration

· ∇Φ int

← d sign( p)∇Φ int(i j)

▷ Small translatory motion, e.g. d = nh, n ∈ N small

end for

16:

18:

(i j)

̸= ∅, since Ωc

c

15:

19:

(i j) Ωc

IC ← IC ∪ {(i, j)}

11:

17:

▷ This loop operates only in Nβ (Φ (i) )

end if end for end for

Fig. 5. Graphical interpretation of intermediate hypersurface as interface.

separates both bodies Ω (i) and Ω ( j) with Ω (i) ∩ Ω ( j) = ∅ into two domains like in Fig. 5b. Then, the following (i j) portion of intermediate surface Γint offers a further alternative besides the master and slave sides: (i j)

Γc(i j) := Γint ∩ Ωc(i j) .

(22)

This surface is the first guess for a contact surface between Ω (i) and Ω ( j) . In deformed state it contains both Γc(i) ( j) and Γc , if equilibrium is achieved. Because it is in general not aligned to any of the two surfaces ∂Ω (i) and ∂Ω ( j) , it can be considered as a fictitious interface.

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

285

(i j)

In undeformed state, the surface Γc acts as reference surface which supports the active set strategy. For this (i j) active set, we locate points on Γint which have to be considered in the KKT constraints (3). One possibility (i j) (i j) (i j) (i j) is to search on Γω,c (deformed equivalent to Γc ) for penetration, i.e. locate the set {Φω,min ≤ 0} ∩ Γω,c (i j) ( j) where Φω,min = min(Φω(i) , Φω ). Since the resulting system of equations is usually defined w.r.t. the undeformed (i j) configuration, one has to trace back Γω,c by d −1 . Eventually, this procedure gives an exact localization of an active set, which consists currently of points violating the non-penetration condition, i.e. gn < 0. In Section 3.5.2, however, we introduce a different criterion. 3.4. Gap function in context of level sets The gap between bodies is given by the distance between two closest points on each contacting surface, see the definition of the (global) gap function in (4). With regard to level set functions, the gap function between Ω (i) and Ω ( j) can be rewritten w.r.t. the deformed state: gn (d (i) (x, tω ), d ( j) (P ( j) (x), tω )) = Φω( j) (d (i) (x, tω )), x ∈ Γc(i) .

(23)

3.4.1. Reformulation of the closest point projection ( j) Projections between two given contact surfaces Γc(i) and Γc are necessary in order to couple both bodies in the mathematical model (2). In most cases the closest point projection is used since it is the obvious choice for contact problems with small deformations. Again, the approach by implicit boundaries is able to simplify the search for closest neighboring points. To this ( j) end, we seek a y ∈ Γc such that the distance to the isosurface L Ψ (Φ ( j) ) for x ∈ Γc(i) is minimized, where now Ψ = Φ ( j) (x). The closest point projection can then be explicitly expressed by P ( j) (x) = x − Φ ( j) (x)

∇Φ ( j) (x) , ∥∇Φ ( j) (x)∥

(24)

such that y = P ( j) (x). Since ∇Φ ( j) (x)/∥∇Φ ( j) (x)∥ is the outward normal on L Ψ (Φ ( j) ), Eq. (24) is an orthogonal projection from y onto the surface L Ψ (Φ ( j) ) and, hence, a closest point projection. As a final remark, we note that the projection in (24) can be regarded as Newton iteration for solving Φ ( j) ( y) = 0 with starting point x. Assuming ∥∇Φ ( j) (x)∥ = 1, the scheme reads x k+1 = x k − Φ ( j) (x k )∇Φ ( j) (x k ), for k = 1, 2, . . . and x 0 = x. Since x is close ( j) to Γc , the above scheme needs only few iterations to converge. 3.4.2. Linearization of the gap function In Section 3.3.2 we introduced the notion of an intermediate surface as contact reference surface. Now we focus on the usage of that surface for the saddle point problem (9). (i j) ( j) First, we state the fact that Γc is generally not identical to Γc(i) or to Γc in undeformed configurations, because (i j) of the definition in (22). In order to integrate the gap function gn (u(i) , u( j) ) on Γc , we extend its definition domain (i) by extrapolating each displacement field u : uˇ (i) : B → Rn , uˇ (i) = u(i) , |Ω (i) analogously for u( j) . For instance, the extrapolation approach in (17) can be exploited for that. Recall from (5) and (i j) (6) that the measure of the gap depends on the closest point projection P. With regard to Φint , we define the gap vector function as: g(u(i) (P (i) (x)), u( j) (P ( j) (x))) := P ( j) (x) − P (i) (x) + u( j) (P ( j) (x)) − u(i) (P (i) (x)), (i j)

(i j)

(25) (i j)

where x ∈ Γc . We assume that Φint is a signed distance function, therefore it holds ∥∇Φint ∥ = 1. Motivated by (23), we use the following ansatz for the right-hand side of (25): [ ] (i j) (i j) (i j) ˇ ˇ g(u(i) (P (i) (x)), u( j) (P ( j) (x))) ≈ Φω(i) (x + {u(x)} tω ) + Φω( j) (x + {u(x)} tω ) ∇Φint (x), (26) ( ) ( ) (i j) where {q}(i j) := α q (i) + (1 − α) q ( j) is defined as the weighted average along Γc for any quantity q and ( j) α ∈ [0, 1]. Often, one simply chooses α = 0.5. After linearizing Φω(i) and Φω on the right-hand side of (26) w.r.t.

286

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

Fig. 6. Zoom into narrowest place between Ω (i) and Ω ( j) in Fig. 5b.

time, we can reassign terms in (25) to their linearized counterpart: [ ] (i j) P ( j) (x) − P (i) (x) ≈ ∇Φint (x) Φ (i) (x) + Φ ( j) (x) ,

(27)

ˇ u( j) (P ( j) (x)) − u(i) (P (i) (x)) = Ju(x)K (i j) . Γ int

(i j)

The double brackets J·K denote the jump along the intermediate surface Γint : ( ) ( ) (i j) (i j) ( j) (i) ˇ ˇ ˇ Ju(x)K x + θ ∇Φ (x) − u x − θ ∇Φ (x) , (i j) := lim u int int Γ int

θ →0

(i j)

(i j)

ˇ := JuK ˇ Γ (i j) · ∇Φint , jint (u)

(28) (29)

int

(i j)

ˇ is the normal jump. It holds that uˇ (i) (x) = u(i) (P (i) (x)) for x ∈ B \ Ω (i) . So uˇ (i) implicitly depends where jint (u) on the closest point projection P (i) , which can be easily evaluated by (24). (i j) Now we redefine the local gap function as projection of the gap vector g from (26) onto ∇Φint , i.e. gn := (i j) g · ∇Φint . For our purposes, this definition makes sense since we aim to avoid a biased gap formulation. According ( j) (i j) to (27), the gap between the contacting surfaces Γc(i) and Γc at x ∈ Γc can be expressed by the normal jump (i j) ˇ (29) and signed distance functions Φ (i) and Φ ( j) : jint (u) (i j)

ˇ gn (u(i) (P (i) (x)), u( j) (P ( j) (x))) ≈ jint (u(x)) + 2Φ (i) (x),

(30)

(i j)

where we used Φ ( j) = Φ (i) on Γc . In contrast to the initial gap definition in (4), the linearized gap function in (30) (i j) measures the gap w.r.t. the intermediate surface Γc . It still involves applications of point projection twice, but now on a single surface. Because it has equal distance to both bodies’ surfaces, the projection’s length ∥x − P ( j) (x)∥ (i j) ( j) for x ∈ Γc is approximatively halved in comparison to x ∈ Γc(i) , equivalently for ∥x − P (i) (x)∥ and x ∈ Γc . Figs. 6a and 6b depict the difference between the gap vector (25) and the approximated versions (26) or (27). The arrows in Fig. 6a represent the gradients of the level set functions on either side of the intermediate surface. (i j) Hence, the direction of g is defined by the locations of both projections P (i) (x) and P ( j) (x) for x ∈ Γc . On the other hand, the principal direction of the gap function in Fig. 6b depends on the gradient of the intermediate level (i j) set function Φint . Moreover, the initial distance is estimated by evaluating the signed distance functions Φ (i) and Φ ( j) . Note that both figures illustrate the situation in undeformed state. 3.5. Nitsche’s approach for contact problems In the following, we continue to consider contact only between Ω (i) and Ω ( j) , since the subsequent ideas can be easily extended to a multibody scenario.

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

287

3.5.1. Modified perturbed Lagrangian We follow the approach from [38] and define a contact contribution πc (see the definition of the Lagrangian L in (11)) depending on a perturbation term and the intermediate contact surface: ∫ ∫ ϵn ˇ (i j) )2 dγ . (λn − {σn (u)} (31) πc (u, λn ) = gn (uˇ (i) , uˇ ( j) ) λn dγ − (i j) (i j) 2 Γc Γc Through differentiating (31) w.r.t. the dual variable, the solution λn of the minimization problem can be expressed explicitly via the contact pressure and gap function: [ ] 1 (i) ( j) (i j) ˇ λn = {σn (u)} + gn (uˇ , uˇ ) , (32) ϵn − (i j)

where [z]− := min(0, z), for z ∈ R. If equilibrium is achieved, the gap vanishes, i.e. gn = 0 on Γc . Hence, (i j) the Lagrange multiplier and the contact pressure on the contact surface coincide, i.e. λn = {σn (u)}(i j) on Γc . ( j) Moreover, balance of forces (or Newton’s third law actio = reactio) implies σn(i) = σn . It is easy to reason that the multiplier in (32) satisfies the KKT conditions (3) (instead of the contact pressure σn ). Then the resulting contact contribution by Nitsche reads ∫ ∫ )2 ( 1 ˇ (i j) dγ , πcN (u) = gn (uˇ (i) , uˇ ( j) ) {σn (u)} (33) gn (uˇ (i) , uˇ ( j) ) dγ + (i j) 2ϵn Γc(i j) Γc (i j)

which is the result of inserting (32) into (31), assuming Γc is the correct contact surface. Since the standard perturbation approach [4] contains the first integral in (33), it follows that the modified Lagrangian functional in (31) leads to a generalization of the penalty method. In the definition of the contact contribution (33) the contact pressure is involved. With regard to the final weak form, a variational analysis of the second integral results in a lengthy version of the original penalty method [21]. We propose to approximate the contact pressure with a constant, i.e. τn ≈ σn (u), known for instance from previous iterations or computations. If we consider the linearized gap function from (30), the problem from (9) becomes for ϵ˜n := ϵn /2: Find u ∈ V such that ⟩ ⟩ 1 ⟨ (i j) 1 ⟨ (i) (i j) ⟩ ⟨ (i j) (i j) ˇ jint (ˇv ) = f (v) − a(u, v) + jint (u), Φ , jint (ˇv ) − {τn }(i j) , jint (ˇv ) ∀v ∈ V. (34) ϵn ϵ˜n Note that by Nistche’s approach, see (32) or (33), the Lagrange multiplier has been condensed out in the last equation. Hence, the saddle point problem (9) can be summarized as a standard weak form with a single unknown. The resulting method derived from (34) is identical to the penalty method with an additional contact pressure term on the right-hand side. The corresponding contact integral terms read then by means of the gap linearization (30): ( ∫ ⟩ ( )T ) 1 ⟨ (i j) 1 (i j) (i j) (i j) ˇ jint (ˇv ) = ˇ Γ (i j) · ∇Φint ∇Φint j (u), JuK · Jˇv KΓ (i j) dγ , (35) int int ϵn int ϵn Γc(i j) ( ) ∫ ⟩ 1 ⟨ (i) (i j) ⟩ ⟨ 1 (i) (i j) (i j) Φ , jint (ˇv ) + {τn }(i j) , jint (ˇv ) = Φ + {τn }(i j) Jˇv KΓ (i j) · ∇Φint dγ . (36) (i j) int ϵ˜n ϵ ˜ n Γc (i j)

Finally, the contact stress {τn }(i j) in (36) is the result of a further projection onto the surface normal ∇Φint along (i j) Γc : ( )T (i j) (i j) {τn }(i j) = ∇Φint {τ }(i j) ∇Φint , (37) where τ is a constant approximation of the Cauchy stress σ . 3.5.2. Active set strategy Because of (34), our proposed Nitsche’s approach requires the approximated contact pressure τn . In contrast to the standard penalty method, the violation of the non-penetration constraint is not a sufficient criterion for the active (i j) set strategy from Section 3.3.2. In general, an appropriate active set is defined as the set of points A(i j) on Γc where the KKT conditions (3) are satisfied. Because only the negative part of the expressions in (32) is considered,

288

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

this set can be identified as: A(i j) := {x ∈ Γc(i j) | λ(in j) (x) < 0}.

(38)

Again, we assume contact between bodies denoted with Ω (i) and Ω ( j) . Now we take the expression from (32) for (i j) (i j) λn and substitute g(uˇ (i) , uˇ ( j) ) · ∇Φint for gn (uˇ (i) , uˇ ( j) ), where g(uˇ (i) , uˇ ( j) ) equals the right-hand side of (26) with tω = 1. This substitution yields: ) 1 ( (i) ˇ (i j) ) + Φω( j) (x + {u} ˇ (i j) ) . λ(in j) (x) := {τn (x)}(i j) + Φω (x + {u} ϵn (i j)

ˇ (i j) gives the updated locations of the intermediate contact surface. The choice For x ∈ Γc the expression x + {u} (i j) for Γc in (22) as contact surface is motivated by the arguments in Section 3.3.2 and has to be understood as a superset of the actual contact surface, which in fact is approximated by A(i j) . In total, we have to find for all contact pairs Ω (i) and Ω ( j) an active set A(i j) . 4. Discretization on a regular grid In this section the focus is on the numerical solution of the variational problem (34) in the context of the unfitted boundary approximation. We assume that the bounding box B consists of uniform cells and is aligned according to the Cartesian grid axes. A voxel is defined as Vlmn := [(l − 1)h 1 , lh 1 ] × [(m − 1)h 2 , mh 2 ] × [(n − 1)h 3 , nh 3 ] ,

(39)

where h k > 0 denotes the spatial step sizes for k = 1, . . . , d, and l, m, n ∈ N. Note that we can easily transfer to two dimensional problems by some minor adjustments. The bounding box B is then equal to the union of all voxels, ⋃ N1 ⋃ N2 ⋃ N3 i.e. B = l=1 m=1 n=1 Vlmn . In order to simplify notations, we denote the grid coordinate with θ := (l, m, n). Given a voxel Vθ from definition (39), we continue the decomposition of B with a subdecomposition of the voxels: Depending on the dimension and substructure, Vθ is then split into a certain number of simplexes, either triangles in 2D or tetrahedrons in 3D. We emphasize that we do not use a subdivision like in a quad- or octree-decomposition, but it might be considered as a preprocessing step to resolve more details in Vθ [23]. A simplex S is defined as the set of all convex combinations of given nodes s1 , . . . , sd+1 : S := {x ∈ B | ∃ςk : 0 ≤ ςk ≤ 1 with

d+1 ∑

ςk = 1, such that x =

d+1 ∑

ςk sk }.

k=1

k=1

Note that we only consider simplexes in shape of triangles and tetrahedrons, therefore there are d + 1 nodes in a simplex. A further decomposition of Vθ is then defined by NS simplexes: N

T (Vθ ) := {S1 , . . . , S NS }, such that Vθ =

S ⋃ ˙

Sk .

k=1

The number NS depends on the dimension d and decomposition pattern [39], see Fig. 7a. Here, the union of simplexes in the last equation is disjoint in the sense that two different simplexes share at most one face, edge or node. Furthermore we demand Θ(S) ⊊ Θ(Vθ ) for all S ∈ T (Vθ ) where Θ(·) denotes the set of nodes’ coordinates θ for either a simplex or a voxel. Concerning the bounding box, the voxels Vθ are identified as elements of B resulting in a structured discretization into quadrilateral or hexahedral elements. Clearly, it is possible to use each simplex in all voxels as elements, too, but we keep the previously mentioned choice of elements. Similarly to how voxels are indexed, grid positions are given by x θ = ((l − 1)h 1 , (m − 1)h 3 , (n − 1)h 3 ). Other quantities like level set functions Φ or displacements u at grid nodes ϑ are then denoted with Φϑ := Φ(x ϑ ) or uϑ := u(x ϑ ), respectively. We also introduce standard shape functions ϕϑ from the finite element method. Within each element Vθ of B, the ϕϑ fulfill the partition of unity property and supports the evaluation of uh at a x within Vθ , i.e. ∑ ∑ ϕϑ (x) = 1, and uh (x) = uϑ ϕϑ (x). (40) ϑ∈Θ(Vθ )

ϑ∈Θ(Vθ )

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

289

Fig. 7. Surface representation given by subtessellation and computing isosurfaces.

4.1. Surface reconstruction For the surface reconstruction, we use the level set approach from Section 3.1 in a discrete manner. In each simplex S ∈ T (Vθ ) we introduce linear shape functions ϕϑS like in (40) and approximate the surface or zero level set within S by solving the linear equation ∑ ΦhS (x) := Φϑ ϕϑS (x) = 0. ϑ∈Θ(S)

Certainly, we consider this approach only in simplexes whose edges intersect with the zero level set. In particular, one checks if there are grid coordinates ϑ − ∈ Θ(S) with Φϑ − ≤ 0 and ϑ + ∈ Θ(S) with Φϑ + > 0. Then, between each x ϑ + and x ϑ − , representing an edge in S, the cut points x c define a linear basis for the reconstructed surface points x in S satisfying the above equation. They can be easily computed by x ϑ + Φϑ − − x ϑ − Φϑ + xc = , (41) Φϑ − − Φϑ + such that ΦhS (x c ) = 0. Depending on the spatial dimension and on the positions of the cut points, the surface within S is a line, a triangle or a quadrilateral with x c as vertices, see Fig. 7b for an example in 3D. One can establish specific rules on how to construct these for the sake of a computer implementation [40]. Note that the zero contour surface L 0 (Φ) does not intersect simplexes with constant signs of level set values at their nodes. 4.2. Final linear system of equations Focusing on the expressions in (10), we proceed with the approaches from the previous section. Despite the unconventional, unfitted discretization ansatz on regular grids, the standard procedure for assembling the final system can be applied. Certainly, we ignore all elements which are voids, i.e. whose nodes have only positive level set values. Then, one has to distinguish which elements are uncut or cut by checking the signs of level set values: If Vθ ⊂ Ω , then this element has only negative level set signs at its nodes. Otherwise, heterogeneous level set signs at nodes indicate cut points between the zero level set and grid edges. Note that cut elements may contain portions of multiple domains Ω (i) . In uncut elements exact volume integration in Vθ can be applied to the operators in (10), see [2,19]. Otherwise, i.e. Vθ ∩ {Φh ≤ 0} ⊊ Vθ , we pursue a integration within the subdecomposition from the previous section. We ˜ = {(x q , ωq ) | x q ∈ S} ˜ for a simplex S. ˜ denote an arbitrary quadrature rule with a set of nodes and weights Q(S) Then the quadrature for any continuous function f in a cut element Vθ is accomplished by ∫ ∑ ∑ ∑ ˜ f (x) dx ≈ |S| ωq f (x q ). (42) Vθ ∩{Φh <0}

S∈T (Vθ ) S∈Ξ ˜ (S)

˜ (x q ,ωq )∈Q(S)

290

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

The set Ξ (S) := T (Υ− (S) ∪ Υc (S)) denotes a further subdecomposition of a simplex S, where Υ− (S) := {x ϑ | ϑ ∈ Θ(S) : Φϑ ≤ 0} contains all nodes with negative level set sign and Υc (S) all cut points (41) of S. For an illustration of these sets, Fig. 7b depicts an exemplary cut case. Surface integrals can be evaluated similarly to (42): One must only replace all computations concerning volumes with surface equivalents and further Ξ (S) with T (Υc (S)) because these cut points span the reconstructed plane Vθ ∩ {Φh = 0}. For more details about surface reconstruction in the context of finite elements, we refer to [41,42], and for integration approaches beyond this one [13,43]. The assembly of the total system of equation can now be realized by applying the integration rules for both uncut and cut elements Vθ to the operators in (34). The stiffness matrix A originates from the bilinear form a(uh , v h ) and the right-hand side vector b from f (v h ). The remaining terms are connected to Nitsche’s approach: (i j) (i j) The numerical integration of ⟨ jint (uˇ h ), jint (ˇv h )⟩ in (35) results in the structural matrix CCT . The remaining two (i j) (i j) (i) terms ⟨Φh , jint (ˇv h )⟩ and ⟨{τn }(i j) , jint (ˇv h )⟩ in (36) are then assembled into bΦ and bτn , respectively. Note that all expressions require projections (24) of quadrature points x q , in particular integer grid coordinates θ such that P (i) (x q ) ∈ Vθ ∩ {Φh(i) = 0}, analogously for P ( j) . Clearly, the assembly of these expressions has to be considered for all contact pairs (i, j) ∈ I × I. The constant Cauchy stress τ is reconstructed at grid nodes θ by stress recovery methods, i.e. by extrapolation from quadrature points. Although there are sophisticated approaches available in the context of non-conforming finite elements, like the least squares approximation [18,44], we use a simple volume average of stress results: ∫ 1 τθ = σ (uh ) dx, |Eθ | Eθ where Eθ := supp(ϕθ ) ∩ {Φh ≤ 0}. Note that the stress approximation of τ in voxels via (40) can be expressed by a linear operator S since we assumed linear elasticity at the beginning of this paper. A further involvement of Φint in (37) is then expressed by SΦ which corresponds to the term bτn . Defining NΦ as the integer number such that E(θk ∩ ){Φh(i) ≤ 0} ̸= ∅ for all k = 1, . . . , NΦ and i = 1, . . . , NΩ , we ⋃ NΦ define the discrete displacements vector as uh = uθk k∈Θ(E) , where E := k=1 Eθk . Finally, the complete discrete system reads: ) ( 1 1 (43) A + CCT uh = b − SΦ uh − bΦ , ϵn ϵ˜n where uh is a known solution for the displacements used for the construction of τ . Replacing the term bτn with SΦ uh motivates to regard (43) as a step in a fixed point method, see line 18 in Algorithm 2. Concerning convergence, the estimation of the spectral radius is presented for this scheme in [38]. The penalty parameter ϵn is supposed to be very small, i.e. ϵn ≪ 1, in order to effect a penalization in the system. We follow the approach from [4] and set it in relation to the Young’s modulus E: ϵn(i j) :=

hϵn0 , {E}(i j)

(44)

which can be further adjusted by ϵn0 for each contacting pair. For our experiments, we choose ϵn0 from the interval [10−1 , . . . , 10], unless otherwise stated. Furthermore, this factor is kept as a constant during the iterations. Finally, we present in Algorithms 1 and 2 pseudocodes for the initialization (focus on creation of level set functions and contact detection) and the routine for solving the contact problem, respectively. Expressions in capital letters denote methods which were mentioned in previous sections and sets in subscripts mean a spatially restrictive application or evaluation. Indeed, one can recognize that in Algorithm 2 the main effort (between lines 7 and 16) is expended for the active set search. 5. Numerical examples In this section we present several numerical test results realized by our proposed method. Initially, we consider simple contact problems in 2D and 3D whose analytical solutions are available to study the accuracy and effectiveness and show numerical convergence. Furthermore, we present more complex tests with multiple contacts in 3D to demonstrate the current range of application.

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

291

Algorithm 2 Routine for Contact Problems in Quasi-Static Deformation Processes 1:

Call Initialization and Contact Detection from Algorithm 1

2:

Assemble A and b

3:

k ← 0 and

4:

repeat

5:

u0h

▷ Body-wise assembly, see Section 4.2

← uh

▷ Initial displacements from Alg. 1

k ←k+1

6:

Ak ← ∅

7:

for all (i, j) ∈ IC do

8:

for l = i, j do ← EXTRAPOLATE(uk−1 , Φ (l) )|Ωc (i j) uk−1 h h|Ω (l)

9:

Φω(l)

10:

← LSM(Φ

(l)

, uk−1 h )|Ωc (i j)

▷ Application of level set method to (15)

11:

end for

12:

{τn } ← STRESSRECOVERY(uk−1 h )|Ω (i j)

13:

λn ← {τn }|Γ (i j) +

14:

(i j)

15:

A

ω,c

1 (Φω (i) ϵn

( j)

k

▷ Update of active set at quadrature points

A ←A ∪A end for

17:

Assemble CCT , bΦ and SΦ )−1 ( ( b − SΦ uk−1 − ukh ← A + ϵ1n CCT h

19:

▷ Explicit multiplier in current configuration

ω,c

(i j)

16:

18:

c

+ Φω )|Γ (i j)

← {x q | λn (x q ) ≤ 0}

k

▷ To ensure feasible motion field for LSM

until Ak = Ak−1 and ∥ukh −

uk−1 h ∥

▷ Based on current active set Ak 1 b ϵ˜n Φ

)

▷ Linear system from (43)

< TOL

For the simulations we consider bi- or trilinear shape functions. In all cases we applied a preconditioned conjugate gradient (CG) method in order to solve (43). As preconditioner we used the factor from an incomplete Cholesky decomposition of the matrix on the left-hand side of (43). As drop tolerance we chose 10−4 and 10−3 for two and three dimensional problems, respectively. For more details about solving this kind of discrete system of equations we refer to [45,46]. 5.1. Hertz problem 5.1.1. Two dimensional case We begin with the Hertz problem in 2D. The plain strain condition is assumed for this experiment, therefore a circular cylinder is pressed on a plate with homogeneous Dirichlet boundary conditions on three sides, see Fig. 8a. The cylinder is initially in contact at one point with the plate and is only stressed by a force F on its very top. There are analytical solutions available from literature, for instance from [4,47,48]. These formulas give explicit expressions for the contact stress σn and contact radius a. Note that each bodies’ level set function Φ (i) can be given analytically. To avoid a rigid body motion of the cylinder, two degree of freedoms (DOFs) are fixed in horizontal direction along the symmetry line. We consider a force F = 250 N for this experiment. For comparison we repeat the computation with F = 1000 N and explain differences. Furthermore, we assume for the external forces b(1) = b(2) = 0 N/mm3 in (2). Instead of a point load in Fig. 8a, we spread the force F symmetrically along the Neumann boundary Γn(1) with |Γn(1) | = 1 mm and obtain (1) an equivalent traction ˆt = F/|Γn(1) | where the cylinder has a radius of 10.018 mm. Concerning the analytical solution formulas [47], note that the plate can be considered as a cylinder with an infinite radius. For this experiment we choose the following set of material parameter: The cylinder Ω (1) has a Young’s modulus E (1) = 7000 MPa and a Poisson’s ratio ν (1) = 0.3, and for the plate E (2) = 106 MPa and ν (2) = 0.45. Through this choice of materials, a relatively soft cylinder is pressed on a almost rigid plate in this numerical test.

292

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

Fig. 8. Hertz problem in 2D.

Fig. 9. Convergence study.

Now we compare the theoretical solution for the contact pressure σn with the ones computed by our method, in particular τn(i) for i = 1, 2. For this purpose, Fig. 8b shows three graphs: One depicts the analytical result and the remaining two the pressures on the two contacting bodies. First, a good match of all three graphs means a good agreement with the theory. Second, the solutions for the two bodies almost coincide which implies that force balance along the contact area Γc(1,2) has been recovered. For most of our chosen grid sizes the resulting contact radii a were in acceptable agreement with the analytical solutions, so we omit to give details about that. We now examine the convergence w.r.t. contact stress σn . For this purpose, we consider the maximal stress on the contact surface Γc(1,2) and investigate its behavior when the grid size h decreases. In case of the smaller force F = 250 N the graph in Fig. 9b indicates first convergence and second that the limit (247.28 MPa) is close to the analytical solution (246.44 MPa). In case of the large force F = 1000 N, convergence is sill observable, but the limit (497.45 MPa) differs slightly from the theoretical result (492.88 MPa). One possible reason is that the

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

293

Table 1 Comparison between different penalty parameters for h = 0.1 mm with F = 250 N. ϵn0

Av. CG steps

Max. pressure [MPa]

Contact radius [mm]

Max. penetration [mm]

1 · 10−3

1328.7 35.83 37.95 41.84 45.62 48.26 51.65

320.35 246.69 245.27 245.19 245.26 241.45 202.77

1.26 0.62 0.61 0.62 0.66 0.76 0.99

8.59 · 10−4 2.91 · 10−4 2.89 · 10−4 2.95 · 10−4 3.55 · 10−4 1.12 · 10−3 6.28 · 10−3

1 · 10−2 0.1 1 10 100 1000

assumption of small deformations for the analytical results is not valid for F = 1000 N. Nevertheless, the deviation is about one percent and, therefore, negligible. Regarding convergence w.r.t. the displacements, we follow the procedure from [32,49] where in the latter source a similar study with the same unfitted finite element approach has been published. A reference solution uh,ref is computed on a very fine grid (h = 0.015 mm), since there is no analytical solution for u available. The outcome of our investigation can be viewed in Fig. 9a: It shows the plots of the relative error E rel = ∥uh − uh,ref ∥/∥uh,ref ∥ in the energy and H 1 (Ω )-norm as functions of the grid size. Obviously, the error decreases almost linearly in the energy norm. Comparing with the results from [3,50], this convergence behavior is optimal. The total convergence rate in H 1 (Ω )-norm appears to be similar. In contrast to the energy norm, the development in the H 1 (Ω )-norm does not show a uniform degradation. Because of results from other publications we assume that this convergence behavior is linked to the unfitted discretization approach, in particular the approximation of the contact surface. For stable convergence recovery we refer to [10,11,51,52]. Furthermore, we study the influence of the penalty factor ϵn in Eq. (43). Without the pressure term SΦ uh on the right-hand side, the method would be identical to the penalty method which has been extensively examined in [4]. Table 1 summarizes the results: We vary the parameter ϵn0 in (44) and observe that for the smallest value the solver’s performance for (43) suffers in form of many CG iterations. In addition, both contact pressure and radius are inaccurate. In general, a very small penalty parameter leads to ill conditioned systems, too. Concerning CG iterations, one can conclude that smaller values for the penalty factor cause fewer solver steps (unless ϵn0 is not too small) and smaller penetrations. Compared to the analytical solution (a = 0.645 mm), the contact areas are estimated as too small, which gets better with increasing the parameter. However, in many cases the stress results are in good agreement with the analytical solution (max |σn | = 246.44 MPa). Hence, our method benefits from the modified Nitsche’s approach (32). As observation, it is advisable to take a moderate choice for ϵn0 . For this study of the Hertz problem we used ϵn0 = 10. 5.1.2. Three dimensional case This experiment should mainly show that our new method is applicable to 3D problems and still provide accurate enough results. Figs. 10a and 10b illustrate the setting for the Hertz problem in 3D. The opaque surface in Fig. 10b is the intermediate surface from (21). The contacting bodies (a hemisphere as replacement for the full cylinder) are (1,2) colored according to the distance field Φint . As load case we set the external body force of the hemisphere b(1) equal to the gravitational force bg = (0, 0, g), where g = −9.81 N/mm3 . An equivalent force in negative z-direction is then given by F = gV (1) = 2gπr 3 /3, where r is the radius of the sphere (we set r = 3.36 mm) and V (1) is the volume of the hemisphere. These data yield then F = 779.37 N. As boundary conditions, we choose homogeneous Dirichlet conditions on all sides of the box, except the contact side, and two DOFs along the symmetry line of the hemisphere are fixed analogously to the 2D case. Both bodies have the same material parameters from the previous study, so E (1) = 7000 MPa, ν (1) = 0.3, E (2) = 106 MPa and ν (2) = 0.45. From [47] analytical solutions to the three dimensional case is known. Plugging in all geometrical and material data, we obtain the maximal contact stress max |σn | = 920.85 MPa and the contact radius a = 0.63 mm as reference values. As experiment, we compare our numerical results with the analytical ones. For h = 0.05 mm we obtained with our method maximal contact stresses max |τn(1) | = 913.75 MPa (hemisphere) and max |τn(2) | = 910.67 MPa (plate) and a contact radius ah = 0.66 mm. Hence, the contact area is well approximated by our level set approach. As

294

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

Fig. 10. Setting for Hertz problem in 3D.

Fig. 11. Clipped results for Hertz problem in 3D.

a small drawback, however, a larger difference between the contact pressures of the two bodies is recognizable, compared to the two dimensional case. However, the deviation for max |τn(1) | from max |σn | is about two percent. Besides, Figs. 11a and 11b show physically reasonable results for the stress in vertical direction. The stress results around the contact surfaces indeed look identical to the two dimensional case. Fig. 11b demonstrates also that penetration was prevented in a realistic way. The black lines are deformed grid edges (diagonal lines are only used for visualization) or reconstructed surfaces within elements. 5.2. Problems with multiple contacts Now, we turn to problems in 3D with more contacts. We present only numerical results for stresses and give some remarks about the solver’s performance. This section should primarily show how far the application range of our method is. For all subsequent numerical tests we assume zero body forces, so b(i) = 0 for all i ∈ I. 5.2.1. Linked rings This example treats rings or tori linked into a chain. But instead of the whole chain with all its links, we focus on a reduced model consisting of two halves and one complete ring. Fig. 12a shows a sketch of the setting, projected onto 2D. A box in dashed lines clips the mentioned objects into the resulting initial configuration. Regarding boundary conditions, Fig. 12b depicts another sketch where constant Dirichlet conditions are applied on all cut sides of the half rings. For a better illustration, bodies within the dashed box from Fig. 12a are rotated by 90 degrees along

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

295

Fig. 12. Contacts in a chain.

Fig. 13. Contact surface and active set.

the longitudinal axis. Here, we set uˆ = ±(0.06 cm, 0, 0) according to the sketch. As in the examples before, we fix two DOFs within the complete ring in order to prevent rotations around the longitudinal axis. These DOFs are indicated by crosses in Fig. 12a. For this experiment we choose a grid size h = 0.02 cm, such that the bounding box B contains 60 × 50 × 50 voxels, according to the directions along the x-, y- and z-axis. The rings have an outer radius of 0.3980 cm and an inner radius, defining the thickness, of 0.1 cm. Fig. 13a shows the zero level sets of all bodies in contact, colored by (i j) (i j) the corresponding Φint . As before, the opaque surface represents the intermediate level set Γint for a single collision (i j) (for better presentation we omit the other contact surface). Note that, on the complete torus, Φint is depicted w.r.t. only one collision. For all bodies we choose the same material parameters, in particular E (i) = 104 N/cm2 and ν (i) = 0.3 for i = 1, 2, 3. According to the boundary conditions, all bodies are stretched outwards. Obviously, two contact spots can be localized, each between a half ring and the complete torus. We choose a smaller penalty parameter (ϵn0 = 10−1 in (44)) than in the previous experiments. As maximal penetration we obtain then a value of 2.8 × 10−3 cm, on ( j) (i j) both intermediate surfaces. For all tests, we measured the penetration by simply evaluating Φω(i) + Φω on Γω,c , (i) ( j) where Ω and Ω are in contact. Regarding the amount of iterations and solver’s performance, we observe that our method needed 8 steps in loop 4 of Algorithm 2. Each step in this loop corresponds to an update of the active sets. The system has about 90 000 unknowns solved with approximately 109 CG steps, averaged over all outer iterations. In (38) we defined the active set strategy, we apply on the level of the quadrature points. From the graph in Fig. 13b one can first notice that both collisions need approximatively the same amount of quadrature points in each update of the active sets. Admittedly, we started with rather large active sets of quadrature points, but the amount decreases fast, which leads to the second observation: In fact, we could have stopped the search after the third iteration since then the active

296

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

Fig. 14. Stress results in deformed configuration.

sets have basically not changed. Apart from changes in the active sets A(i j) , our algorithm aborts the computations if the difference between solutions from subsequent iterations is below a chosen tolerance, see Algorithm 2. Finally, we present the results for stresses in the deformed configuration in Figs. 14a and 14b. There, the scene has been clipped longitudinally in the middle in order to view the stresses near the contacts better. In the figure on the left-hand side, the red spots indicate high values which can be found near the two contacts as well as on portions of the bodies which suffer most from being stretched or compressed. The rather elliptical form of the middle ring in the last figures demonstrates that our method provides good results for a realistic deformation. 5.2.2. Contact between fibers As last example, we arrange four beams or fibers in a network with rotational symmetry: Each fiber has contact with two other beams such that all together build a stable structure. The setting in voxel discretization can be viewed in Fig. 15a. All fibers have a length of 10.2 mm and a radius of 0.48 mm. For the depiction in Fig. 15a we chose a grid size h = 0.2 mm resulting in a fiber’s diameter of 5 voxels. So we use here a coarser grid than in the examples before, in order to demonstrate that our algorithm still provides realistic results. All fibers have a Young’s modulus E (i) = 104 MPa and a Poisson’s ratio ν (i) = 0.3 for i = 1, . . . , 4. As boundary conditions we choose again non-homogeneous Dirichlet conditions on the ends of the beams. A sketch can be found in Fig. 15b. The known displacements are depicted by arrows. As in the previous section, the sketch presents the scene as 2D projection onto two of the beams (back and front face). These conditions in the sketch are then applied on all four projections. In this test case we set uˆ = ±(0, 0, 0.5 mm). Regarding contact, Fig. 16a shows the reconstructions of all fibers by their zero level sets. In order to determine the signed distance functions, we use the analytical expression for cylinders. Approximating the surfaces of the fibers, the zero contours are colored in red, the colors turn to blue in the interiors. Again, the opaque surfaces are (i j) parts of all intermediate surfaces L 0 (Φint ). We adjust the penalty factor in (44) by setting ϵn0 = 10−1 . With all given parameters our method needed 14 iterations in the loop 4 of Algorithm 2, whereby after 5 steps the active set was finally determined. For solving the linear system (43) (∼ 24000 unknowns) the procedure used 106 CG steps per iteration in average. Moreover, we observed a maximal penetration of 9.7 × 10−3 mm. In Fig. 16b one can examine the results for the von Mises stress in the deformed configuration. The bending of each fiber around the other contacting beams is clearly recognizable. Moreover, the stress peaks appear where expected. Hence, our proposed method could realistically simulate a complex problem with multiple contacts between thin, structural components.

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

297

Fig. 15. Setting for the structure.

Fig. 16. Reconstruction by level sets and stress results.

6. Conclusion In summary, this paper deals mainly with two topics: The application of the finite element method on voxel-based meshes and a simultaneous treatment of contact between solid bodies. For the former we consider the established ideas from the immersed boundary [11] or non-conforming finite element method [8], both on Cartesian grids. For the latter issue we presented a novel approach by exploiting the level set framework [13,14]. Our numerical tests show that enforcing the contact constraints on the intermediate surface (21) gives accurate results. Although using a piece-wise linear approximation of this surface, we were able to verify, e.g. through the Hertz problem, that our approach provides results in good agreement with analytical solutions. Besides, the level set approach enabled us to expand our contact algorithm from 2D to 3D problems without any extra effort. In combination with Nitsche’s approach our method achieves, first, excellent contact stress results while preventing interpenetrations and, second, convergence when the grid size decreases. Although our method currently reaches at most first order convergence due to bi- and trilinear shape functions, similar work [38] indicates that higher order is possible, too. Regarding contact search and updating the active set, the level set method proved to be robust in all numerical tests and showed no complications, even for moderate large deformations. From our experience, a further great benefit is that the application of the LSM does not negatively affect the computing time of the actual contact routine, presented in Algorithm 2. As mentioned in the introduction, we aim to integrate our method for contact in the analysis of microstructures based on digital data [1]. Our last numerical test indicates that our contact algorithm works fine for simple, porous

298

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

structures built by thin fibers on a coarse grid discretization. Nevertheless, tangential interactions due to friction are currently not considered in our approach. For this purpose, one has to apply the friction law, for instance, by Coulomb in order to distinguish between stick and slide portions on the (intermediate) contact surface. More complicated problems like the compression or tensile loading of woven or non-woven fabrics are associated with large tangential motions and deformations of bodies. For those issues, we propose the same Eulerian surface dynamics from Section 3.2 in order to update both the bodies’ boundaries and contact locations. A similar approach has been published in the work [53]. Acknowledgments We would like to thank Kevin T. Chu (Serendipity Research) and Masa Prodanovic (University of Texas at Austin) for their open access to their Level Set Method Library (LSMLIB), which we implemented in our routines. Download link and further information can be found in [54]. Furthermore, we appreciate the support of the High Performance Center Simulation- and Software-Based Innovation in Kaiserslautern, Germany, as well as the German Research Foundation (DFG), in particular the projects AN 341/9-1 and VO 1487/16-1. References [1] GeoDict. https://www.math2market.com/. (Accessed: 05.09.18). [2] P. Wriggers, Computational Contact Mechanics, second ed., Springer Berlin Heidelberg, 2006. [3] B. Wohlmuth, Variationally consistent discretization schemes and numerical algorithms for contact problems, Acta Numer. 20 (2011) 569–734. [4] N. Kikuchi, J. Oden, Contact Problems in Elasticity - A Study of Variational Inequalities and Finite Element Methods, SIAM, 1988. [5] A. Popp, M. Gitterle, W. Gee, W.A. Wall, A dual mortar approach for 3D finite deformation contact with consistent linearization, Internat. J. Numer. Methods Engrg. 83 (2010) 1428–1465. [6] Z. Dostál, D. Horák, D. Stefanica, A scalable FETI-DP algorithm with non-penetration mortar conditions on contact interface, J. Comput. Appl. Math. 231 (2) (2009) 577–591. [7] R. Kornhuber, R. Krause, O. Sander, P. Deuflhard, S. Ertel, A monotone multigrid solver for two body contact problems in biomechanics, Comput. Vis. Sci. 11 (1) (2008) 3–15. [8] A. Hansbo, P. Hansbo, An unfitted finite element method, based on nitsche’s method, for elliptic interface problems, Comput. Methods Appl. Mech. Engrg. 191 (47–48) (2002) 5537–5552. [9] T. Belytschko, C. Parimi, N. Moës, N. Sukumar, S. Usui, Structured extended finite element methods for solids defined by implicit surfaces, Internat. J. Numer. Methods Engrg. 56 (4) (2003) 609–635. [10] J. Haslinger, Y. Renard, A new fictitious domain approach inspired by the extended finite element method, SIAM J. Numer. Anal. 47 (2009) 1474–1499. [11] T. Rüberg, F. Cirak, J. Manuel, G. Aznar, An unstructured immersed finite element method for nonlinear solid mechanics, Adv. Model. Simul. Eng. Sci. 3 (2016) 22. [12] M. Schneider, D. Merkert, M. Kabel, FFT-Based homogenization for microstructures discretized by linear hexahedral elements, Internat. J. Numer. Methods Engrg. 109 (10) (2016) 1461–1489. [13] S. Osher, R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, in: Applied Mathematical Sciences, vol. 153, Springer-Verlag, 2003. [14] J. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Sciences, Cambridge University Press, 1999. [15] P. Hansbo, Nitsche’s method for interface problems in computational mechanics, GAMM-Mitt. 28 (2) (2005) 183–206. [16] P. Wriggers, G. Zavarise, A formulation for frictionless contact problems using a weak form introduced by Nitsche, Comput. Mech. 41 (2008) 407–420. [17] M. Juntunen, R. Stenberg, Nitsche’s method for general boundary conditions, Math. Comp. 78 (267) (2009) 1353–1374. [18] M. Tur, J. Albelda, O. Marco, J.J. Ródenas, Stabilized method of imposing dirichlet boundary conditions using a recovered stress field, Comput. Methods Appl. Mech. Engrg. 296 (2015) 352–375. [19] D. Braess. Finite Elements: Theory, Fast. Solvers, Fast Solvers and Applications in Solid Mechanics, third ed., Cambridge University Press, 2007. [20] J. Nocedal, S. Wright, Numerical Optimization, second ed., Springer-Verlag New York, 2006. [21] F. Chouly, P. Hild, A nitsche-based method for unilateral contact problems: Numerical analysis, SIAM J. Num. Anal. 51 (2013) 1295–1307. [22] E. Brun, A. Guittet, F. Gibou, A local level-set method using a hash table data structure, J. Comput. Phys. 231 (6) (2012) 2528–2536. [23] C. Min, F. Gibou, A second order accurate level set method on non-graded adaptive cartesian grids, J. Comput. Phys. 225 (1) (2007) 300–321. [24] C. Min, F. Gibou, Robust second-order accurate discretizations of the multi-dimensional heaviside and Dirac delta functions, J. Comput. Phys. 227 (22) (2008) 9686–9695.

A. Leichner, H. Andr¨a and B. Simeon / Computer Methods in Applied Mechanics and Engineering 352 (2019) 276–299

299

[25] M. Mirzadeh, A. Guittet, C. Burstedde, F. Gibou, Parallel level-set methods on adaptive tree-based grids, J. Comput. Phys. 322 (2016) 345–364. [26] J.A. Sethian, A fast marching level set method for monotonically advancing fronts, Pnas 93 (4) (1996) 1591–1595. [27] J.A. Bærentzen, on the Implementation of Fast Marching Methods for 3D Lattices. Technical Report, Technical University of Denmark, 2001. [28] C. Min, On reinitializing level set functions, J. Comput. Phys. 229 (8) (2010) 2764–2772. [29] T.D. Aslam, A partial differential equation approach to multidimensional extrapolation, J. Comput. Phys. 193 (1) (2004) 349–355. [30] D. Adalsteinsson, J.A. Sethian, The fast construction of extension velocities in level set methods, J. Comput. Phys. 148 (1) (1999) 2–22. [31] D. Peng, B. Merriman, S. Osher, H. Zhao, M. Kang, A PDE-based fast local level set method, J. Comput. Phys. 155 (2) (1999) 410–438. [32] F. Chouly, R. Mlika, Y. Renard, An unbiased Nitsche’s approximation of the frictional contact between two elastic structures, Numer. Math. 139 (3) (2018) 593–631. [33] R.A. Sauer, L.D. Lorenzis, A computational contact formulation based on surface potentials, Comput. Methods Appl. Mech. Engrg. 253 (2013) 369–395. [34] R.A. Sauer, L.D. Lorenzis, An unbiased computational contact formulation for 3D friction, Internat. J. Numer. Methods Engrg. 101 (October 2014) (2015) 251–280. [35] S.W. Chi, C.H. Lee, J.S. Chen, P.C. Guan, A level-set enhanced natural kernel contact algorithm for impact and penetration modeling, Int. J. Numer. Methods Eng. 102 (2015) 839–866. [36] D. Durville, Contact-friction modeling within elastic beam assemblies: an application to knot tightening, Comput. Mech. 49 (2012) 687–707. [37] T.W. McDevitt, T.A. Laursen, A mortar-finite element formulation for frictional contact problems, Internat. J. Numer. Methods Engrg. 48 (10) (2000) 1525–1547. [38] M. Tur, J. Albelda, J.M. Navarro-Jiménez, J.J. Rodenas, A modified perturbed Lagrangian formulation for contact problems, Comput. Mech. 55 (2015) 737–754. [39] T. Apel, N. Düvelmeyer, Transformation of hexaedral finite element meshes into tetrahedral meshes according to quality criteria, Computing 71 (2003) 293–304. [40] C. Min, F. Gibou, Geometric integration over irregular domains with application to level-set methods, J. Comput. Phys. 226 (2) (2007) 1432–1443. [41] T.P. Fries, O. Samir, Higher-order accurate integration of implicit geometries, Internat. J. Numer. Methods Engrg. 106 (2016) 323–371. [42] M. Moumnassi, S. Belouettar, É. Béchet, S.P. Bordas, D. Quoirin, M. Potier-Ferry, Finite element analysis on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces, Comput. Methods Appl. Mech. Engrg. 200 (5–8) (2011) 774–796. [43] P. Smereka, The numerical approximation of a delta function with application to level set methods, J. Comput. Phys. 211 (1) (2006) 77–90. [44] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates, Part 1: The recovery technique, Internat. J. Numer. Methods Engrg. 33 (7) (1992) 1331–1364. [45] Z. Dostál, On preconditioning and penalized matrices, Numer. Linear Algebra Appl. 6 (2) (1999) 109–114. [46] M. Ferronato, C. Janna, G. Gambolati, Mixed constraint preconditioning in computational contact mechanics, Comput. Methods Appl. Mech. Engrg. 197 (45–48) (2008) 3922–3931. [47] A. Konyukhov, R. Izi, Introduction To Computational Contact: A Geometrical Approach, John Wiley & Sons, Inc., 2015. [48] V.L. Popov, Contact Mechanics and Friction, Springer-Verlag, Berlin Heidelberg, 2010. [49] M. Fabre, J. Pousin, Y. Renard, A fictitious domain method for frictionless contact problems in elasticity using Nitsche’ s method, SMAI J. Comput. Math. (2016). [50] S. Hüeber, B.I. Wohlmuth, An optimal a priori error estimate for nonlinear multibody contact problems, SIAM J. Numer. Anal. 43 (1) (2005) 156–173. [51] E. Burman, S. Claus, P. Hansbo, M.G. Larson, A. Massing, CutFEM: Discretizing geometry and partial differential equations, Internat. J. Numer. Methods Engrg. 104 (2015) 472–501. [52] S. Loehnert, A stabilization technique for the regularization of nearly singular extended finite elements, Comput. Mech. 54 (2) (2014) 523–533. [53] E. Vitali, D. Benson, Kinetic friction for multi-material arbitrary Lagrangian eulerian extended finite element formulations, Comput. Mech. 43 (6) (2009) 847–857. [54] K.T. Chu, Prodanovic. M., Level Set Method Library (LSMLIB). http://ktchu.serendipityresearch.org/software/lsmlib/index.html (Accessed: 02.08.18).