FINITE ELEMENT IMPLEMENTATION OF GEOMETRICALLY NONLINEAR CONTACT

The OOFEM finite element software has been recently updated to include contact algorithms for small strain applications. In this work, we attempt to extend the contact algorithms to large strain problems. Reviewing the current code and comparing it with approaches encountered in literature, we arrive at a specific algorithmic solution and integrate it into the current code base. The current code is explained, the necessary extensions are derived and documented, and the algorithmic changes are described. Tests confirm the functionality and quadratic rate of convergence of the proposed implementation.


Introduction
Having originated with Heinrich Hertz in 1881 [1], contact mechanics has been a fringe discipline of structural mechanics for most of the 20 th century, owing mainly to its limited options for analytical analysis. Only with the advent of numerical methods, such as the finite element method (FEM), has development in this field sprung forward [2]. This surge in new advances is only exacerbated by the ever more rapid evolution of computational hardware capabilities experienced today; ever more complicated contact problems necessitating ever more complicated methods previously unimaginable due to hardware constraints are coming into focus [2,3].
The tool used for the purposes of this article is the finite element computational software OOFEM [4]. OOFEM is a free, open-source C++ code developed mainly at the Czech Technical University in Prague. It is based on object-oriented design with emphasis on variability and extendability. Currently, OOFEM's finite element code supports a number of modules for solution of various types of numerical problems. Of interest for our purposes is the sm module intended for structural mechanics.
Basic introduction of contact algorithms into the OOFEM software has been achieved in [5]. This implementation, however, has certain limitations, among them chiefly the lack of support for large strain, i.e., geometrically non-linear problems. In this work, we strive to remove this particular limitation and demonstrate a simple method of extending existing basic penalty-based formulation to handle large strain problems.

Basic equation
Without loss of generality, let us consider two deformable bodies B γ , γ = 1, 2. Moreover, let Ω γ 0 and Ω γ represent the regions occupied by the body in the reference and deformed configuration, respectively. Position of a point in the reference configuration is denoted by X γ , and a position x γ of a point in the deformed configuration is obtained be means of the one-to one deformation mapping ϕ γ as The basic equations for each body under consideration are the following: • Dirichlet boundary conditions • Neumann boundary conditions • Constitutive equation Here, we use the so-called Saint-Venant Kirchhoff material model. Note that the contact implementation is not affected by the choice of constitutive law and extension to inelastic law is straightforward.

FEM Implementation of Geometrically Nonlinear Contact
In the equations above, F is the deformation gradient, E is the Green-Lagrangian strain tensor, Div is divergence operator in the reference configuration, P is the first Piola-Kirchhoff stress,B are prescribed body forces in the reference configuration, u are prescribed displacements,N is a unit normal vector in the reference configuration, Γ u and Γ σ are the Dirichlet and Neumann boundaries which are mutually disjoint, i.e., Γ γ u ∩ Γ γ σ = ∅, and Γ γ u ∪ Γ γ σ = Γ γ , λ γ and µ γ are Lame's parameters.
Moreover, contact conditions introduce essentially a new set of boundary conditions into the task. Note that only normal contact is considered in what follows. The core issue of contact problems can be described by means of the contact constraints in the form of Karush-Kuhn-Tucker optimality conditions [6]: The first condition represents the fact that contact forces t c are only compressive, the second condition is characterized by a gap function g(u) and express a non-penetration condition. The last equation is simply complementarity between two distinct and mutually exclusive states for the system; either contact occurs, and therefore the gap function g c is zero, while a non-zero compressive contact force is present, or the contact force is zero and the gap function is larger than zero.

Penalty Method
Due to the constraint conditions in form of ineqaulities, see (7), weak form of the above introduced formulation is represented by a variational inequality, see [2]. Several methods how to handle the contact constraints are possible, e.g., the Lagrangian multiplier method, the augmented Lagrangian method, or the Nitsche method [3,6]. Currently in OOFEM, the Lagrangian multiplier method and the penalty method are implemented [5]. It is the penalty method which we have selected for the purpose of large strain extension, as its formulation is more straightforward and the necessary code complexity is notably lower. For the penalty method, the weak form can be writ- where δ represents variation and W c is the penaltybased contact contribution along contact interface Γ c defined as where H is the Heaviside function and p is the penalty parameter. The main disadvantages of the penalty method is the that the achieved solution is imprecise (it would only be precise for p → ∞), but moreover, any improvement in precision acts directly against the stability of the computation, as larger penalty parameters introduce larger members into computational stiffness matrices, which disturbs most numerical solvers [2]. Those disadvantages are balanced by the relatively simple formulation of the method and, in comparison with the Lagrangian multiplier method, by simple resulting equation systems, as no new variables are introduced into the computation.

Node-to-Segment Discretization
An implementation of contact problems into the framework of the finite element method necessitates selecting a distinctive approach at contact discretization. As the physical world is described by FEM in the form of finite elements and connecting nodes, contact can only be considered between those. Typical contact discretizations include nodeto-node, node-to-segment, segment-to-segment and other approaches [3]. At current time, node-to-node and node-to-segment discretizations are implemented in OOFEM [5]. We are only extending the node-tosegment discretization for large strains, as the nodeto-node discretization is by definition inadequate for this purpose.
In node-to-segment discretization, a single node comes into contact with a segment. Various objects may be understood under the term "segment", most prominently a set of element boundaries (element edges in 2D domains or element surfaces in 3D domains) or an arbitrarily defined rigid analytical surface.
Contact detection in node-to-segment discretization necessitates determining three things [3,6]: (1.) The gap function g c , essentially a projection of the node onto the contact segment, which determines whether there is contact or not; analogical to the penetration function in the analytical approach to contact problems. (2.) The normal vector in the deformed configuration n, which determines the orientation of the contact surface at the contact point. Obtained for most segments in OOFEM as vector perpendicular to the contact segment's surface at the contact point [5]. (3.) The extended N-matrix N * , which ensures allocation of the external forces and stiffness to the proper degrees of freedom in the FEM formulation. When the segment in question is an element boundary, the matrix takes the form [5] where I (d×d) is a square unit matrix of the appropriate dimension, which corresponds to degrees of freedom of the contacting node, and N e (x c ) is the

FEM Formulation
To integrate contact constraints into the FEM framework, considerations have to be made for the contact forces and energy discussed above within the principles of FEM discretization. For the penalty method, contact energy is expressed as Note that the energy contact contribution above is calculated only in the contact point, which can be interpreted as Lobatto integration of (9), see, e.g., [6] for more details.
Variation of this part of energy yields It is worth noting that the term pg(u) represents the contact force. Variation of the gap function can be further expressed as where δd c are the variational nodal displacements of relevant FEM nodes and n is the normal vector in deformed configuration. The external forces of contact are therefore To obtain the associated tangential stiffness matrix. the internal forces have to be subjected to differentiation with respect to nodal displacements d c [5]. For small strains, it is generally assumed that the direction of contact does not change throughout the analysis, and therefore n remains independent of d c [5,6].

Extension to Large Strain Problems
Now we stand before a question: which of the previously described considerations change if we consider large strains? The mechanics of introducing contact constraints to FEM remain the same, as well as the definition of contact energy. Two critical differences appear. Firstly, it can no longer be assumed that contact follows the same direction throughout the whole computation. The normal vector has to be reassessed after each computational step and also the stiffness matrix has to reflect this. Secondly, large deformations result in change of contact point on the segment, which is reflected in the evaluation of matrix N * .
The normal vector is deformation dependent and can be expressed geometrically [2] n = t 1 × t 2 (16) where × is the vector cross product and t i are tangential vectors of the element boundary in the deformed configuration (in the case of an element edge, i.e., contact formulation in 2D, one is the tangent and the other is the e 3 , that is, a unit vector perpendicular to the plane).
For the new expression of stiffness, we elect to utilize a formula presented in [2]. Specifically for the case of node-to-2D-segment penalty method, it reads where B * is the extended B-matrix. Note that the first member of the stiffness matrix corresponds with the previously derived term for the small strain case. For a linear edge in 2D, the extended matrices are defined as where ξ is the parametric coordinate on the element edge.

Class Structure
The OOFEM code has a strictly object-oriented structure [4]. The integration of contact into this structure has been mainly accomplished in [5] and only slight additions were necessary for the large strain extension. Figure 1 displays the relevant parts of OOFEM class structure together with the classes used for contact.
The general idea of class structure centers around two main classes: ElementEdgeContactSegment and Node2SegmentPenaltyContact.
ElementEdgeContactSegment is derived from the general ContactSegment class and represents element edges in 2D. When initialized from an OOFEM input file, the class receives a set of element boundaries. Internal processes ensure that contact is considered always with the edge which is closest to the contacting node [5].
Node2SegmentPenaltyContact is a subclass of ActiveBoundaryCondition, an OOFEM feature vol. 30/2021 FEM Implementation of Geometrically Nonlinear Contact which allows boundary conditions to contribute by themselves to the internal forces, external forces, and stiffness matrix. The utilization of this feature avoids creating any form of a "contact element", otherwise prevalent among FEM contact implementations [3]. The contact boundary condition is likewise initialized from the OOFEM input file with information on the size of the penalty parameter, a set of nodes and a set of one or more contact segments. Any node is always checked for contact with any segment [5].

Operations
The assembly of external forces and stiffness matrix takes place in the Node2SegmentPenaltyContact class, specifically in its overridden methods assembleVector() and assemble().
For the determination of gap g c , the contact segment method computePenetration() is invoked, which computes it from a purely geometrical consideration, considering the node to have penetrated the segment when its deformed coordinates are on a different side of the segment than the undeformed ones.
The contact segment is also responsible for the computation of the tangent and normal vectors. The tangent vector is determined as where B * seg is the segment part of the extended Bmatrix and x i , y i coordinates of the segment nodes. Meanwhile, the normal vector is determined as for large strains and from the Element::edgeEvalNormal() function for small strains (which returns the undeformed normal vector).
Neither of those vectors is normalized. The normalization of the normal vector happens in the contact condition, when the N v , B v and T v matrices are determined according to (18).
In the assembly of stiffness, the contact conditions asks the contact segment whether the large strain part shall be assembled. This is accomplished by the ContactSegment::hasNonLinearGeometry() function and based on the setting of the element in contact.

Demonstration
To demonstrate the functionality of our approach, we select a simple computational example in a 2D plane strain idealization. The mesh used can be seen in Figure 2a.
A double-clamped beam is positioned above a cantilever. The length of the beam is 5 m and length of the cantilever is 3 m. Both are 0.25 m high with a gap of 0.25 m between them. The top of the beam is loaded in every step by an uniform distributed load of 8000 kN/m, causing vertical deformation and thus    subsequently contact between the two bodies. Loading is performed in 200 uniform loading steps. Both beam and cantilever are modelled by rectangular plane strain elements (OOFEM element Quad1PlaneStrain). In total, the task mesh consists of 805 elements. Refinement of the mesh has been accomplished with the help of a custom script in the MATLAB software [8]. All elements are made of the same Saint-Venant Kirchhoff material with the parameters λ = 100 MPa and µ = 150 MPa.
Contact conditions are introduced between the lower edge of the beam and the upper edge of the cantilever below. On the latter interface, two boundary conditions are defined, one pairing the beam's nodes with the cantilever's element edges and the other vice versa. This ensures that no clipping caused by crude or otherwise imperfect meshing can take place. Figure 2b shows the deformed state of the mesh and the vertical component of the stress tensor in step 4 of 200, when the beam and the cantilever first come into contact. Likewise, Figure 2c shows the same in step 200 at the end of the analysis. To render those vol. 30/2021 FEM Implementation of Geometrically Nonlinear Contact figures, the ParaView open-source software was used [9]. Values of stress in the visualizations are given in Pa.
Note that between steps 4 and 200, the contact surface moves significantly to the left both on the beam and the cantilever. This was accomplished by alternating activation of the two opposite boundary conditions on the interface, which managed to subsequently introduce contact between neighboring nodesegment pairs, while the contact gap between the original pairs opened up again.
The main objective of this test was to prove quadratic convergence of the new approach to large strain contact stiffness matrix.
The vast majority of the 200 steps computed managed to converge in 6 iterations or less. The only exception is step 4, in which the contact occurs for the first time. Table 1 shows error evolution in a typical computational step (step 14 was chosen). Moreover, the dependence of error (norm of residual forces) on the number of iterations during the equilibrium Newton-Raphson process for steps 4, 14, and 164 is visualized in Figure 3. Thanks to the semilogarithmic scale used, it is apparent that approximate quadratic convergence has been achieved.

Conclusions
The paper describes an implementation of a two dimensional large strain node-to-segment penaltybased approach to contact mechanics into an opensource finite element code OOFEM. The emphasis is given on proper design of the hierarchy of classes with a focus on modularity and extensibility towards different methods for handling contact constraints, contact with friction, or three dimensional formulation. The functionality of the presented approach is successfully verified on a numerical example and a quadratic rate of convergence of the Newton-Raphson iteration procedure is demonstrated.