Continuous gap contact formulation based on the screened Poisson equation

We introduce a PDE-based node-to-element contact formulation as an alternative to classical, purely geometrical formulations. It is challenging to devise solutions to nonsmooth contact problem with continuous gap using finite element discretizations. We herein achieve this objective by constructing an approximate distance function (ADF) to the boundaries of solid objects, and in doing so, also obtain universal uniqueness of contact detection. Unilateral constraints are implemented using a mixed model combining the screened Poisson equation and a force element, which has the topology of a continuum element containing an additional incident node. An ADF is obtained by solving the screened Poisson equation with constant essential boundary conditions and a variable transformation. The ADF does not explicitly depend on the number of objects and a single solution of the partial differential equation for this field uniquely defines the contact conditions for all incident points in the mesh. Having an ADF field to any obstacle circumvents the multiple target surfaces and avoids the specific data structures present in traditional contact-impact algorithms. We also relax the interpretation of the Lagrange multipliers as contact forces, and the Courant–Beltrami function is used with a mixed formulation producing the required differentiable result. We demonstrate the advantages of the new approach in two- and three-dimensional problems that are solved using Newton iterations. Simultaneous constraints for each incident point are considered.


Introduction
Considerable interest has recently emerged with the enforcement of essential boundary conditions [1] and contact detection algorithms [2] using approximate distance functions (ADFs). Realistic contact frameworks for finite-strain problems make use of three classes of algorithms that are often coupled (see [3]): 1. Techniques of unilateral constraint enforcement. Examples of these are direct elimination (also reducedgradient algorithms), penalty and barrier methods, and Lagrange multiplier methods with the complementarity condition [4]. Stresses from the underlying discretization are often used to assist the normal condition with Nitsche's method [5].
3. Discretization methods. Of concern are distance calculations, estimation of tangent velocities and general discretization arrangements. In contemporaneous use are (see [3]) node-to-edge/face, edge/face-to-/edge/face and mortar discretizations.
In terms of contact enforcement and friction algorithms, finite displacement contact problems are typically addressed with well-established contact algorithms, often derived from solutions developed by groups working on variational inequalities and nonsmooth analysis. However, in the area of contact discretization and the related area of contact kinematics [11], there are still challenges to be addressed in terms of ensuring the complete robustness of implicit codes. One of the earliest papers on contact discretization was by Chan and Tuba [12], who considered contact with a plane and used symmetry to analyze cylinder-to-cylinder contact. Francavilla and Zienkiewicz [13] later proposed an extension to node-to-node contact in small strains, allowing for direct comparison with closed-form solutions.
The extension to finite strains requires further development, and the logical development was the node-to-segment approach, as described in the work of Hallquist [14]. Although node-to-segment algorithms are known to entail many defficiencies, most of the drawbacks have been addressed. Jumps and discontinuities resulting from nodes sliding between edges can be removed by smoothing and regularization [15]. Satisfaction of patch-test, which is fundamental for convergence, can be enforced by careful penalty weighing [16,17]. It is well known that single-pass versions of the node-to-segment method result in interference and double-pass can produce mesh interlocking, see [18,19]. This shortcoming has eluded any attempts of a solution and has motivated the development of surface-tosurface algorithms. One of the first surface-to-surface algorithms was introduced by Simo, Wriggers, and Taylor [11].
Zavarise and Wriggers [20] presented a complete treatment of the 2D case and further developments were supported by parallel work in nonconforming meshes, see [21]. A review article on mortar methods for contact problems [22] where stabilization is discussed and an exhaustive detailing of most geometric possibilities of contact was presented by Farah, Wall and Popp [23]. This paper reveals the significant effort that is required to obtain a robust contact algorithm. An interesting alternative approach to contact discretization has been proposed by Wriggers, Schröder, and Schwarz [24], who use an intermediate mesh with a specialized anisotropic hyperelastic law to represent contact interactions. In the context of large, explicit codes, Kane et al. [25] introduced an interference function, or gap, based on volumes of overlapping, allowing non-smooth contact to be dealt in traditional smooth-based algorithms.
In addition to these general developments, there have been specialized algorithms for rods, beams, and other structures. Litewka et al. [26] as well as Neto et al. [27,28], have produced efficient algorithms for beam contact.
For large-scale analysis of beams, cables and ropes, Meier et al. [29] have shown how beam specialization can be very efficient when a large number of entities is involved.
Recently, interest has emerged on using approximate distance functions [30,31,32,2] as alternatives to algorithms that heavily rely on computational geometry. These algorithms resolve the non-uniqueness of the projection, but still require geometric computations. In [30], Wolff and Bucher proposed a local construction to obtain distances inside any object, but still require geometric calculations, especially for the integration scheme along surfaces. Liu et al. [31] have combined finite elements with distance potential discrete element method (DEM) in 2D within an explicit integration framework. A geometric-based distance function is constructed and contact forces stem from this construction. In [32], the analysis of thin rods is performed using classical contact but closed-form contact detection is achieved by a signed-distance function defined on a voxel-type grid. In [2], a pre-established set of shapes is considered, and a function is defined for each particle in a DEM (discrete element method) setting with a projection that follows. In the context of computer graphics and computational geometry, Macklin et al. [33] introduced an element-wise local optimization algorithm to determine the closest-distance between the signed-distance-function isosurface and face elements. Although close to what is proposed here, no solution to a partial differential equation (PDE) is proposed and significant geometric calculations are still required.
In this paper, we adopt a different method, which aims to be more general and less geometric-based. This is possible due to the equivalence between the solution of the Eikonal equation and the distance function [34]. It is worth noting that very efficient linear algorithms are available to solve regularized Eikonal equations, as discussed by Crane et al. [35]. The work in [35] provides a viable solution for contact detection in computational mechanics. Applications extend to beyond contact mechanics and they provide a solution for the long-standing issue of imposing essential boundary conditions in meshfree methods [1].
We solve a partial differential equation (PDE) that produces an ADF (approximate distance function) that converges to the distance function as a length parameter tends to zero. The relation between the screened Poisson equation (also identified as Helmholtz-like equation), which is adopted in computational damage and fracture [36,37] and the Eikonal equation is well understood [38]. Regularizations of the Eikonal equation such as the geodesics-inheat [35] and the screened Poisson equation are discussed by Belyaev and Fayolle [39]. We take advantage of the latter for computational contact mechanics. Specifically, the proposed algorithm solves well-known shortcomings of geometric-based contact enforcement: 1. Geometric calculations are reduced to the detection of a target element for each incident node. 3. The contact force direction is unique and obtained as the gradient of g(x).
4. Since the Courant-Beltrami penalty function is adopted, contact force is continuous in the normal direction.
Concerning the importance of solving the uniqueness problem, Konyukhov and Schweizerhof [40] have shown that considerable computational geometry must be in place to ensure uniqueness and existence of node-to-line projection.
Another important computational effort was presented by Farah et al. [23] to geometrically address all cases in 3D with mortar interpolation. Extensions to higher dimensions require even more intricate housekeeping. When compared with the geometrical approach to the distance function [32,2], the algorithm is much simpler at the cost of a solution of an additional PDE. Distance functions can be generated by level set solutions of the transport equation [41].
The remainder of this paper is organized as follows. In Section 2, the approximate distance function is introduced as the solution of a regularization of the Eikonal equation. In Section 3, details concerning the discretization are introduced. The overall algorithm, including the important step-control, is presented in Section 4. Verification and validation examples are shown in Section 5 in both 2D and 3D. Finally, in Section 6, we present the advantages and shortcomings of the present algorithm, as well as suggestions for further improvements.

Approximate distance function (ADF)
Let Ω ⊂ R d be a deformed configuration of a given body in d-dimensional space and Ω 0 the corresponding undeformed configuration. The boundaries of these configurations are Γ and Γ 0 , respectively. Let us consider deformed coordinates of an arbitrary point x ∈ Ω and a specific point, called incident, with coordinates x I . For reasons that will become clear, we also consider a potential function φ (x I ), which is the solution of a scalar PDE. We are concerned with an approximation of the signed-distance function. The so-called gap function is now introduced as a differentiable function g [φ (x I )] such that [3]: If a unique normal n (x I ) exists for x I ∈ Γ, we can decompose the gradient of g [φ (x I )] into parallel ( ) and n (x I ) on those points. In the frictionless case, the normal contact force component is identified as f c and contact conditions correspond to the following complementarity conditions [4]: The vector form of the contact force is given by Assuming that a function g [φ (x I )], and its first and second derivatives are available from the solution of a PDE, we approximately satisfy conditions (2) by defining the contact force as follows: where κ is a penalty parameter for the Courant-Beltrami function [42,43]. The Jacobian of this force is given by: Varadhan [44] established that the solution of the screened Poisson equation: produces an ADF given by −c L log [φ (x)]. This property has been recently studied by Guler et al. [38] and Belyaev et al. [34,39]. The exact distance is obtained as a limit: which is the solution of the Eikonal equation. Proof of this limit is provided in Varadhan [44]. We transform the ADF to a signed-ADF, by introducing the sign of outer (+) or inner (−) consisting of Note that if the plus sign is adopted in (7), then inner points will result in a negative gap g(x). The gradient of and the Hessian of g(x) is obtained in a straightforward manner: Note that c L is a square of a characteristic length, i.e. c L = l 2 c , which is here taken as a solution parameter. Using a test field δφ(x), the weak form of (5) is written as where dA and dV are differential area and volume elements, respectively. Since an essential boundary condition is imposed on Γ such that φ (x) = 1 for x ∈ Γ, it follows that δφ(x) = 0 on Γ and the right-hand-side of (10) vanishes.

Discretization
In an interference condition, each interfering node, with coordinates x N , will fall within a given continuum element.
The parent-domain coordinates ξ for the incident node x I also depends on element nodal coordinates. Parentdomain coordinates are given by: and it is straightforward to show that for a triangle, with similar expressions for a tetrahedron [45]. The continuum element interpolation is as follows: where N K (ξ) with K = 1, . . . , d + 1 are the shape functions of a triangular or tetrahedral element. Therefore, (11a) can be written as: In (13), we group the continuum node coordinates and the incident node coordinates in a single array x C = x N x I T with cardinality #x C = d (d + 2). We adopt the notation x N for the coordinates of the continuum element. For triangular and tetrahedral discretizations, ξ I is a function of x N and x I : The first and second derivatives of a N with respect to x C make use of the following notation: Source code for these functions is available in Github [45]. A mixed formulation is adopted with equal-order interpolation for the displacement u and the function φ. For a set of nodal displacements u N and nodal potential values φ N : where in three dimensions where ξ I = a N (x C ). First and second derivatives of N K [a N (x C )] are determined from the chain rule: For the test function of the incident point, we have For linear continuum elements, the second variation of φ (ξ) is given by the following rule: Since the gradient of φ makes use of the continuum part of the formulation, we obtain: where j is the Jacobian matrix in the deformed configuration. The element residual and stiffness are obtained from these quantities and available in Github [45]. Use is made of the Acegen [46] add-on to Mathematica [47] to obtain the source code for the final expressions.

Algorithm and step control
All nodes are considered candidates and all elements are targets. A simple bucket-sort strategy is adopted to reduce the computational cost. In addition, Step control is required to avoid the change of target during Newton-Raphson iteration. The screened Poisson equation is solved for all bodies in the analyses. Figure 1 shows the simple mesh overlapping under consideration.
The resulting algorithm is straightforward. Main characteristics are: • All nodes are candidate incident nodes.
• All elements are generalized targets.
The modifications required for a classical nonlinear Finite Element software (in our case SimPlas [48]) to include this contact algorithm are modest. Algorithm 1 presents the main tasks. In this Algorithm, blue boxes denote the contact detection task, which here is limited to: 1. Detect nearest neighbor (in terms of distance) elements for each node. A bucket ordering is adopted. 2. Find if the node is interior to any of the neighbor elements by use of the shape functions for triangles and tetrahedra. This is performed in closed form.
3. If changes occurred in the targets, update the connectivity table and the sparse matrix addressing. Gustavson's algorithm [49] is adopted to perform the updating in the assembling process.
In terms of detection, the following algorithm is adopted: 1. Find all exterior faces, as faces sharing only one tetrahedron.
2. Find all exterior nodes, as nodes belonging to exterior faces. ii. Calculates the projection on the element (ξ I ) and the corresponding shape functions N (ξ I ).
iii. If 0 ≤ N K (ξ I ) ≤ 1 then e is the target element. If the target element has changed, then flag the solver for connectivity update.
Since the algorithm assumes a fixed connectivity

Numerical tests
Numerical examples are solved with our in-house software, SimPlas [48], using the new node-to-element contact element. Only triangles and tetrahedra are assessed at this time, which provide an exact solution for ξ I .

Mathematica
[47] with the add-on Acegen [46] is employed to obtain the specific source code. All runs are quasistatic and make use of a Neo-Hookean model. If C is the right Cauchy-Green tensor, then

Patch test satisfaction
We employ a corrected penalty so that the contact patch test is satisfied in most cases. This is an important subject for convergence of computational contact solutions and has been addressed here with a similar solution to the one discussed by Zavarise and co-workers [17,16].
We remark that this is not a general solution and in some cases, our formulation may fail to pass the patch test. Figure 2 shows the effect of using a penalty weighted by the edge projection, see [16]. However, this is not an universal solution.

Two-dimensional compression
We begin with a quasi-static two-dimensional test, as shown in Figure 3, using three-node triangles. This test consists of a compression of four polygons (identified as part 3) in Figure 3 by a deformable rectangular punch (part 1 in the same Figure). The 'U' (part 2) is considered rigid but still participates in the approximate distance function (ADF) calculation. To avoid an underdetermined system, a small damping term is used, specifically 40 units with L −2 MT −1 ISQ dimensions. Algorithm 1 is adopted with a pseudo-time increment of ∆t = 0.003 for t ∈ [0, 1].
For h = 0.020, h = 0.015 and h = 0.010, Figure 4 shows a sequence of deformed meshes and the contour plot of φ(x).
The robustness of the algorithm is excellent, at the cost of some interference for coarse meshes. To further inspect the interference, the contour lines for g(x) are shown in Figure 5. We note that coarser meshes produce smoothedout vertex representation, which causes the interference displayed in Figure 4. Note that g(x) is determined from Using the gradient of φ (x), the contact direction is obtained for h = 0.02 as shown in Figure 6. We can observe that for the star-shaped figure, vertices are poorly represented, since small gradients are present due to uniformity of φ (x) in these regions. The effect of mesh refinement is clear, with finer meshes producing a sharper growth of reaction when all four objects are in contact with each other. In contrast, the effect of the characteristic length l c is not noticeable.
In terms of the effect of l c on the fields φ(x) and g(x), Figure 7 shows that information, although φ(x) is strongly affected by the length parameter, g(x) shows very similar spatial distributions although different peaks. Effects of h and l c in the displacement/reaction behavior is shown in Figure 8. The mesh size h has a marked effect on the results up to h = 0.0125, the effect of l c is much weaker.  g(x) g(x) g(x)

Three-dimensional compression
In three dimensions, the algorithm is in essence the same. Compared with geometric-based algorithms, it significantly reduces the coding for the treatment of particular cases (node-to-vertex, node-to-edge in both convex and non-convex arrangements). The determination of coordinates for each incident node is now performed on each tetrahedra, but the remaining tasks remain unaltered. We test the geometry shown in Figure 9 with the following objectives: • Assess the extension to the 3D case. Geometrical nonsmoothness is introduced with a cone and a wedge.
• Quantify interference as a function of l c and κ as well as the strain energy evolution.
Deformed configurations and contour plots of φ (x) for this problem are presented in Figure 9, and the corresponding CAD files are available on Github [45]. A cylinder, a cone and a wedge are compressed between two blocks.
Dimensions of the upper and lower blocks are 10 × 12 × 2 consistent units (the upper block is deformable whereas the lower block is rigid) and initial distance between blocks is 8 consistent units. Length and diameter of the cylinder are 7.15 and 2.86 (consistent units), respectively. The cone has a height of 3.27 consistent units and a radius of 1.87.
Finally, the wedge has a width of 3.2, a radius of 3.2 and a swept angle of 30 degrees. A compressible Neo-Hookean law is adopted with the following properties: • Blocks: E = 5 × 10 4 and ν = 0.3.  Figure 11. It is noticeable that l c has a marked effect near the end of the compression, since it affects the contact force.

Two-dimensional ironing benchmark
This problem was proposed by Yang et al. [50] in the context of surface-to-surface mortar discretizations. Figure 12 shows the relevant geometric and constitutive data, according to [50] and [51]. We compare the present approach with the results of these two studies in Figure 13. Differences exist in the magnitude of forces, and we infer that this is due to the continuum finite element technology. We use finite strain triangular elements with a compressible neo-Hookean law [52]. The effect of l c is observed in Figure 13. As before, only a slight effect is noted in the reaction forces. We use the one-pass version of our algorithm, where the square indenter has the master nodes and targets are all elements in the rectangle. Note that, since the cited work includes friction, we use here a simple model based on regularized tangential law with a friction coefficient, µ f = 0.3. 0 2e+08 4e+08 6e+08 8e+08 1e+09

Three-dimensional ironing benchmark
We now perform a test of a cubic indenter on a soft block. This problem was proposed by Puso and Laursen [18,19] to assess a mortar formulation based on averaged normals. The frictionless version is adopted [19], but we choose the most demanding case: ν = 0.499 and the cubic indenter. Relevant data is presented in Figure 14. The rigid 1 × 1 × 1 block is located at 1 unit from the edge and is first moved down 1.4 units. After this, it moves longitudinally 4 units inside the soft block. The soft block is analyzed with two distinct meshes: 4 × 6 × 20 divisions and 8 × 12 × 40 divisions. Use is made of one plane of symmetry. A comparison with the vertical force in [19] is performed (see also [18] for a clarification concerning the force components). We allow some interference to avoid locking with tetrahedra. In [19], Puso and Laursen employed mixed hexahedra, which are more flexible than the crossed-tetrahedra we adopt here. Figure 15 shows the comparison between the proposed approach and the mortar method by Puso and Laursen [19]. Oscillations are caused by the normal jumps in gradient of φ(x) due to the classical C 0 finite-element discretization (between elements). Although the oscillations can be observed, the present approach is simpler than the one in Puso and Laursen.  [19].
We introduced a discretization and gap definition for a contact algorithm based on the solution of the screened Poisson equation. After a log-transformation, this is equivalent to the solution of a regularized Eikonal equation and therefore provides a distance to any obstacle or set of obstacles. This approximate distance function is smooth and is differentiated to obtain the contact force. This is combined with a Courant-Beltrami penalty to ensure a differentiable force along the normal direction. These two features are combined with a step-control algorithm that ensures a stable target-element identification. The algorithm avoids most of the geometrical calculations and housekeeping, and is able to solve problems with nonsmooth geometry. Very robust behavior is observed and two difficult ironing benchmarks (2D and 3D) are solved with success. Concerning the selection of the length-scale parameter l c , which produces an exact solution for l c = 0, we found that it should be the smallest value that is compatible with the solution of the screened Poisson equation. Too small of a l c will produce poor results for φ (x). Newton-Raphson convergence was found to be stable, as well as nearly independent of l c . In terms of further developments, a C 2 meshless discretization is important to reduce the oscillations caused by normal jumps and we plan to adopt the cone projection method developed in [7] for frictional problems.