A displacement‐driven approach to frictional contact mechanics

The contribution at hand aims at a numerical method to represent contact between two or more different solid bodies. The proposed computational contact algorithm is based on a displacement‐driven approach to describe the geometrical constraint resulting from the impenetrability of the solid bodies in contact. A solution method of contact between one deformable and one rigid body is assumed for the initial computationally fast implementation of the underlying equations. Elastic or inelastic material behavior can be considered for the deformable body. Furthermore, the proposed approach is capable of finite deformations. Moreover, a node‐to‐segment (NTS) formulation is utilized for contact element discretization. Various numerical examples are presented to demonstrate the numerical method's applicability for different contact mechanics tasks. First, the presented examples verify the displacement‐driven approach by investigating its convergence behavior and consistency. Then, comparisons of solutions by the proposed formulation to the analytical solution and the penalty method are provided to validate its accuracy. Moreover, the stick and slip condition of large sliding motion in a frictional contact problem is successfully depicted.


INTRODUCTION
Contact mechanics is essential for describing the behavior of two or more solid bodies, which undergo particular motion that brings the considered body's initial admissible state in the reference configuration to occupy entirely or partially the same spatial domain.The considered solid bodies are sufficiently described as a set of particles or material points.Since the present contribution aims at a general continuum mechanical setting for the deformation behavior of the solid bodies, the material points can be parameterized.This parameterization enables an assignment of specific properties, such as stresses.A reference configuration is a domain corresponding to the initial position of material points of one specific solid body.The solid bodies might be deformed due to external loads, such as prescribed displacements, external forces, or contact, with or without friction.The Lagrangian description of motion, where a body and its material Furthermore, the NTS algorithm may fail to detect the contact for specific geometries.Non-physical penetration, which occurs if curved surfaces are considered, can be treated by extending the linear NTS algorithm to a higher order element. 18The inadmissible state can also be overcome by extending the NTS algorithm to the segment-to-segment (STS) strategy. 19,20Despite having a more cumbersome formulation, the STS strategy offers a more physical result.The development of the segment-to-segment method has become prevalent in recent years, leading to a more sophisticated method based on domain decomposition techniques, 21 called the mortar method.3][24] However, further discussion of this method is outside the scope of this article.
This research proposes a novel contact method based on the displacement-driven approach, incorporating the linear node-to-segment contact element principle.A relatively straightforward contact formulation that delivers a more stable computation to a general contact problem is expected.No penalty-like parameter is required in this model, even though an additional geometrical constraint is introduced.To the best of the authors' knowledge, the description of NTS combined with a geometrical displacement-driven approach for the frictional contact problem has not yet been presented.The solution method at hand can be applied and extended to solve various general computational contact problems.Moreover, finite deformations with elastic or inelastic material descriptions might be incorporated.In order to deliver a comprehensive understanding of the proposed contact formulation, the current work is structured in five sections as follows.The first section provides an overall description of the contribution.Then, in Section 2, the contact kinematics of the developed approach is introduced.The developed framework of the displacement-driven contact formulation is given in Section 3 with respect to a discretized FE setting.Section 4 demonstrates the applicability of the developments to different numerical examples.Additionally, the convergence behavior and the robustness of the computational solution method are verified.Finally, a conclusion and outlook are presented in Section 5.

KINEMATICS OF TWO BODIES IN THREE-DIMENSIONAL CONTACT
This contribution evaluates a three-dimensional frictional contact problem not limited to small strain but also covers large deformations and sliding.Even though the self-contact problem 25 is recently included in many finite element codes, contact between two bodies is considered in this displacement-driven contact formulation for a more general explanation.In the NTS algorithm, one body is selected as the "slave" while the other is identified as the "master".Then, contact between the slave node and the master segment is evaluated.In the further explanations, equations are written with the superscript that identifies the body in contact, where i = 1, 2 represent the slave and master bodies for the node-to-segment contact element, respectively.This numbering rule is consistent with the standard nomenclature used in contact mechanics.Figure 1 illustrates a continuum setting of two bodies Ω (i) undergoing finite deformations, where Ω (i) is composed of material point P (i) ∈ Ω (i) with the open sets Ω (1) ⊂ R 3 and Ω (2) ⊂ R 3 are the bodies in the reference configuration.The reference position is introduced by X (i) = (P (i) , t 0 ), where  (i) is a mapping of P (i) ∈ Ω (i) onto X ⊂ R. The current configuration at current time t is represented by Ω (1) t ⊂ R 3 and Ω (2) t ⊂ R 3 .Observing a material point P (i) , while Ω (i) moves, leads to the definition of the current position x (i) = (P (i) , t).Motion of the bodies from reference position X (i) to current position x (i) is defined by the displacement vector u (i) , which is expressed as u (i) = x (i) − X (i) .The deformation gradient F (i) is obtained by taking the derivative (i) . ( The potential contact surfaces are defined as boundary Ω (i) , which can be divided into three disjoint boundary sets that satisfy this following Equation (2) for each body i where Γ (i) u , Γ (i)  refer to the Dirichlet boundary, where displacements are prescribed and the Neumann boundary, where tractions are prescribed, respectively.In addition, since contact is expected, Γ (i) c is the contact boundary, where the contact constraints will be enforced.The counterparts of these boundaries in the current configuration are denoted as  (i) u ,  (i)  , and  (i) c .Furthermore, the traction and motion boundary conditions are defined for the reference configuration in a standard form as Schematic depiction of two bodies in contact from reference to current configuration.
where P (i) denotes the first Piola-Kirchhoff stress tensor and N (i) represents the component of the outward normal to Γ (i)  with respect to the reference configuration.Moreover, the prescribed traction on the Neumann boundary is denoted by t(i) 0 , while û(i) is the prescribed displacements applied on the Dirichlet boundary.Once the nodes of the slave surface Γ (1) c and the nodes forming the master surface Γ (2) c are defined, the event of contact is detected in the current configuration by measuring their relative position.Contact is evaluated via normal gap g N using Equation (4) based on the closest point projection as g N (X, t) = (x (1) (X (1) , t) − x (2) (X (2) , t)) ⋅ n (2) , ( where x (1) and x (2) are the current positions of a slave node and a master node, respectively.The corresponding master node x (2) must be the closest to the referred slave node, while n (2) denotes the outward unit normal of the master surface in the current configuration  (1) c .The unit normal of master surface n (2) is computed using the cross-product of two vectors representing the edge of the considered surface, which are connected to node x (2) (X (2) , t).This computation is adequate since a linear surface element is employed.To avoid gap function computation for all surfaces, the master node x (2) (X (2) , t), where the outward unit normal of the surface is computed, can be specified as the closest master node to the currently evaluated slave node.The calculation for the outward unit normal of the master surface is written as n (2) i,k = (x (2) left − x (2) i ) × (x (2) right − x (2) i ) ||(x (2) left − x (2) i ) × (x (2) right − x (2) i )|| Referring to Figure 2 and Equation (5), x (2) i is the position of node i, where the surface outward unit normal is evaluated, x (2) left and x (2) right in the node of surface edge that is connected to node i.The left and right nodes must be defined consistently so that the cross-product points outward.In more detail, Figure 3 illustrates contact checking between two bodies, where the cube acts as the master body and the tetrahedral as the slave body.The inadmissible state detection of the slave node at the position x (1) is carried out in this case.Therefore, the unit normal of three surfaces, namely S 1 , S 2 , and S 6 , is evaluated in node 1 of the master surface with position x (2) , with n (2) i,1 , n (2) i,2 , and n (2) i,6 , as the outward unit normal, respectively.Only master surfaces that fulfill the requirement of n (1) j ⋅ n (2) i,k < 0 ( 6 )   are used to evaluate the geometrical penetration condition.Equation (6) indicates that the unit normal of the slave node and the master surface is pointing in the opposite direction, which in Figure 3 is satisfied by only one surface n i,6 (x (2) ).
The outward unit normal at the slave node is evaluated using the outward unit normal of each surface connected to the corresponding slave node based on nodal averaging where n (1) j is the unit normal of slave node j and n j,k is the normal of surface element k evaluated at slave node j, which is illustrated in Figure 4.The non-penetration condition is violated once the value of g N is negative.Each slave node, that F I G U R E 4 Outward unit normal computation of the slave node based on nodal averaging.leads to the negative value condition is assigned as the active set node.Next, the initial contact displacement has to be computed and applied to this set of nodes.However, if multiple master surfaces fulfill this requirement, a master surface, which gives minimum distance, is used.Initial contact displacement d C is assumed to be the displacement required to bring slave node j from the inadmissible state to the geometrical non-penetration condition based on a geometric evaluation computed as Geometrical penetration checking between the slave node and the master surface is presented in Algorithm 1.Since the analysis is limited to the contact case between rigid and deformable bodies, only iteration for the deformable body, which is assigned as the slave body, is required.Afterwards, contact analysis with the displacement-driven contact approach is performed if an inadmissible state is detected.

DISPLACEMENT-DRIVEN CONTACT FORMULATION
The main idea of the proposed contact formulation is to obtain a geometrical non-penetration condition by enforcing a displacement load to nodal points in contact.The identification of nodal points in contact and the initial contact displacement is already explained in Section 2. These nodal points are called contact nodes, while the rest of the other nodes are classified as structure nodes.As an initial displacement is applied, a friction law is enforced on these contact nodes to acquire an equilibrium condition that satisfies the surface condition.The contact problem formulation is derived from three sets of residual equations, which are the residual of the structure nodes, the frictional contact nodes, and the additional geometrical constraints, which are also applied to the contact nodes.These residuals are represented by R s , R f , and R g , respectively.Then, these nonlinear sets of equations are coupled and consistently linearized with respect to the displacement increment.Linearization leads to an over-determined system of equations, resulting in a non-square matrix that cannot be directly inverted to solve the system of equations.Hence, a singular value decomposition technique is used to find the solution.Finally, the displacements are updated using the computed displacement increments once the decomposition is carried out, and the residuals are evaluated based on the current displacements.This iterative procedure is conducted until numerically acceptable residual values are obtained.Once the equilibrium conditions are achieved, the contact check is carried out to ensure the geometrical non-penetration of the solution.

Residual of structural node
As previously introduced, the residual of the structural nodes R s represents part of the continuum body that is not involved in the contact boundary.The kinetic equilibrium of this body in the spatial configuration is represented by the balance of linear momentum, which is derived from the global form, assuming the rate of the linear momentum is equal to the sum Algorithm 1. Geometric non-penetration checking for the current step t n+1 collect nodes and their connectivities for  (1) c and  (2) c state determine pair of master node closest to slave node for j = slave nodes (1) do for i = master surface nodes (2) do compute distance between nodes ‖x (1) j − x (2) i ‖ find potential contact pairs which fulfill min‖x (1) j − x (2) i ‖ end for i end for j for j = slave nodes (1) do compute outward unit normal on slave node n (1) j in Equation ( 7) for i = potential master surface nodes (2) do for k = master surface connected to node i do compute normal of master surface k on master node i, n (2) i,k in Equation ( 5) if Equation ( 6) then compute g N j in Equation ( 4)  (9)   where b(x, t) are the mass specific body forces and t are the surface traction forces.From the balance of mass, the Cauchy stress theorem t(x, t; n) ∶= (x, t)n and the Gauss theorem, the local form of the spatial linear momentum balance can be expressed as which represents the kinetic equilibrium condition of any deformation process in the spatial description.Galerkin's method is applied to Equation (11) to be used within the finite element method by applying a virtual test function u s , which is constant in time and satisfies u s ∶= {u s (x) ∈ Ω t |u s = 0∀x ∈  u }.Integrating this equation leads to Using Gauss theorem (see Equation 10), the symmetry condition  =  T , and the divergence theorem div(u) = div  ⋅ u +  ∶ grad(u), Equation ( 12) can be rewritten as which includes the external forces work in the body Ω t and the interia terms, b and  üs , respectively, as well as the external traction forces (u) ⋅ n = u ⋅ (n) = u ⋅ t, which work on the surface .Henceforward, Ω and Ω t are divided into the finite domains Ω E and Ω t E , where holds.Therefore, Equation ( 13) must be discretized using, for example, a linear ansatz function N(), where  is the local coordinate of an element.The displacement field u s = x s − X s is also discretized according to u s = N()d .In index notation, the displacements and the gradient of the displacement field are written as Then, the test function and its gradient in index notation are expressed by Thus, the discretized residual R s on node I for a standard displacement element under consideration of Equation (13), expressed in index notation, can be formulated as The internal and external force contributions are depicted within Equation (18).Note that the external contribution t is only evaluated if a traction load is prescribed at the surface  E of the element under consideration.Additionally, the inertia term can be neglected since the applicability of this solution is limited to the quasi-static case.

Residual of contact nodes
Contact is considered in the proposed formulation by applying a boundary constraint to the node in contact.The applied boundary constraint enforces a displacement load required to divert the condition from an inadmissible state to a geometrical non-penetration state.This displacement is defined in Section 2, specifically in Equation ( 8).The geometrical non-penetration condition is guaranteed by enforcing this condition to the node in contact.However, the contact node position does not describe the correct position since the overall body reaction under contact loading is not yet considered.A Newton iteration is implemented to attain equilibrium, which guarantees that the contact displacement meets both the geometric and the equilibrium state requirements.Consequently, a different condition on the contact node is required to allow displacement correction on the contact node.Hence, a friction law is included in the formulation to release the contact node from the enforced displacement load.Frictional contact behavior is considered within the friction law equation and permits contact nodes to adjust to the correct contact displacement in the equilibrium condition, as shown in Figure 5. F c is the reaction force on the nodal point of contact boundary due to the contact, called the contact force.These contact force components parallel and perpendicular to the surface are given by F I G U R E 5 Contact force calculation in the displacement-driven contact formulation.

F I G U R E 6
Notation for contact force decomposition on surface local coordinate in 3D space.
where vector n and t are directional unit vectors in normal and tangential directions of the master surface, respectively.The dyadic product of vector n with itself is the directional tensor N, and the dyadic product of vector t with itself is expressed as tensor T. Note that there are two axes parallel to the surface in the three-dimensional space.Here, T represents the dyadic product resulting from the sum of these two axes as depicted in Figure 6.During the iteration, F c is computed as the contribution of the internal force in elements projected to the connected contact node due to contact displacement.In the frictionless contact problem, F c will be equal to F N in the equilibrium state while it depicts the resultant of the normal and tangential contact forces in the presence of friction.Coulomb's friction law is enforced at the boundary constraint for the sake of presentation and without loss of generality, which is F I G U R E 7 2D illustration of normal force and tangential force on a node of a deformable body in frictional contact.
The coefficient of friction  is a positive scalar value that depends on the surface properties of the materials in contact.
The force perpendicular to surfaces in contact is known as the normal force, denoted as F N .F N defines the net force compressing two parallel surfaces together as depicted in Figure 7. f denotes the friction exerted by each surface  (i) c on the other, which value may vary from zero to maximum equal to ||F N ||.Since friction is a reaction force on the surface due to motion, friction always works in the opposing direction to the motion relative to each surface.The relative motion direction is represented by Next, the residual of contact nodes is formulated by enforcing Coulomb's friction law in Equation ( 20) to the nodes in contact.Therefore, to obtain equilibrium at the contact node, the mapping of the contact force perpendicular to the surface F T must be equal to the friction force in value of normal force F N multiplied by the friction coefficient  as illustrated in Figure 7. Therefore, the residual of the contact node is defined by where the scalar term |||| = √  ⋅ .Here,  is computed as A ij b j , where A represents an arbitrary tensor, and b denotes an arbitrary vector.The residual of the contact node can be expressed in index notation as In contrast to the residual of the structure node in Section 3.1, Equation ( 22) applies only to the slave node I of the surface in contact, and the formulation is derived directly using the equilibrium of forces at the nodal point.It is apparent that the initial contact displacement load induces the residual of the structure nodes and the contact nodes.Therefore, this residual should be evaluated prior to entering the contact equilibrium iteration, as shown in Algorithm 2, to check if carrying out the iteration scheme is necessary.In particular cases, the initial displacements might be extremely small, causing the residual norm to fall below zero tolerance.This indicates that the equilibrium is maintained, and the loading step can proceed.

Residual of geometrical constraint
In the attempt to solve contact issues by implementing boundary constraints on contact nodes, it is possible to encounter unstable conditions.In such cases, the contact node may oscillate between being situated inside and outside the master surface as it iterates for a displacement that fulfills both the equilibrium and geometric non-penetration conditions.An additional geometrical restraint is applied to the contact node to address this issue and improve numerical stability during iteration.This additional restraint allows for displacements only on the corresponding master surface and prevents any out-of-plane displacements.The geometrical constraint to the contact node is provided by x G(x (1) ) = Δu(x (1) ) ⋅ n(x (2) ), Algorithm 2. Initial checking for entering the Newton scheme for current step t n+1 Requirement: (1) a set of nodes in contact {C} (2) initial displacement to master segment where Δu and n are illustrated in Figure 8.Note that the displacement increment is used here rather than the displacement field to enforce this geometrical restraint in each iteration.The geometrical constraint G in this formulation is applied to the deformable body, which is the slave node since the derivation at hand is currently limited to contact cases between rigid and deformable bodies.Similar to Section 3.2, where the residual is directly derived on the nodal point, the residual of the geometrical constraint for each contact node I can be instantly written in index notation as In contrast to the other two residuals, which are vectors, the geometrical constraint residual is a scalar value.On the downside, this restriction precludes any changes to the contact node set that may happen during the equilibrium iteration.
The contact algorithm necessitates an additional check to consider the reopening condition.Further details regarding the investigation of the contact node-set modification can be found in Section 3.6.

Linearization of residuals
The numerical solution in the finite element analysis is performed by the Newton solution scheme.It is essential to ensure that numerical stability and quadratic convergence behavior are satisfied concurrently by achieving a physically correct result.Hence, a consistent linearization of the three sets of the residual equations is carried out.The residual on each node R I is formulated as The displacements as the nodal unknowns are separated into two parts, namely, the displacements of the structure nodes and contact nodes, which can be expressed by Using the formal structure of linearization for vector function f expressed by where Df(x) is evaluated by finally, the linearized set of residuals can be written as For the sake of simplicity, the partial derivatives of the three residuals for each displacement unknown are expressed as

Linearization of structural node residual
The derivative of the structure node residual equation (see Equation 18) with respect to the displacement of the structure node is derived as follows Accordingly, the structure node residual equation derivative with respect to the displacement of the contact node is derived likewise.Here, Equation ( 31) is valid for constitutive formulations that are linearizable with respect to an algorithmic tangent C.However, in certain cases of inelasticity at finite deformations, additional terms must be added due to a consistent linearization of internal variables, compare Fleischhauer et al. 26

Linearization of contact node residual
Now, the contact node residual equation, which is provided by Equation ( 21), has to be derived with respect to the unknown displacement of contact node S IJ 22 .Since the slave node and master surface pair are defined in the penetration checking, they are assumed constant within the iteration.Therefore, T and N are constant with respect to u. Considering the residual is formulated directly at a node, the derivative is presented as Similar to the formulation of the structure node residual in the previous section, the residual of the contact node is expressed in index notation as It can be seen that in Equation (33), the node index I is only given to the contact force F c because the directional tensors T and N are defined by the master surface pair.Moreover, a derivative of the force with respect to the displacement results in a nodal stiffness shown in Equation (31).Hence, one can view this equation as a modification of node stiffness triggered by the occurrence of contact.Furthermore, although the equation seems complicated, the computation is relatively straightforward since most variables are either constant or defined by geometrical conditions.In the same manner, , which represents the derivative of contact node residual with respect to structure node displacement, is formulated.

Linearization of residual from additional geometrical constraint
In contrast to the quite lengthy partial derivative of the contact node residual described in the previous section, the derivative of the geometrical constraint residual is expressed neatly as Note that n is a normal unit vector for a master surface of the corresponding slave node I, as explained in Equation (23).
The formulation of the geometrical constraint residual is derived directly considering the displacement increment of the contact node.Hence, the derivative of this residual equation with respect to the structure node is zero, as shown in Equation ( 35)

Singular value decomposition
As mentioned previously, an iterative Newton strategy is utilized to solve the problem in terms of the finite element method.Thus, the update of the unknown variables u s and u f is obtained by setting the residual R to zero, which leads to LinR equal zero as follows Rewriting the equation in matrix format reads which is of the form Ax = b that is usually solved by x = A −1 b.However, since the residual derivative S is not a square matrix, this matrix cannot be directly inverted.Hence, the authors refer to Schwarz 27 to obtain the unknown vector Δu where the singular value decomposition of the considered matrix is used.First, the residuals and the residual derivatives are calculated as previously explained.Then, the singular value decomposition of the m × n matrix S is computed as where U and V are unitary matrices with the size of m × m and n × n, respectively, which represent rotations.V * is the Hermitian conjugate of V, which is equal to V T , if matrix S is real.The same rule is applied to U.  is an m × n scaling matrix, which contains the eigenvalue of S in the diagonal terms.Once the decomposition is obtained, the term VV T is added to the left-hand side term, where it is equal to I since V is a unitary matrix.Next, both the left and right-hand sides of the equation are multiplied by the term U T .Thus, the current problem can be expressed as Here, x, b, and A are substituted by the variables of the current problem, where they represent the unknown displacement increment vector Δu, the residual vector R, and the residual derivative with respect to the increment of the solution S. Referring to Equation (38) and considering that in the contact case, matrix S is a real matrix,  can be expressed as  = U T S V.Then, substituting variables of the current problem to Equation (39) leads to Finally, the unknown displacement increment Δu can be calculated by rearranging Equation (40) as Then, the iterative nodal solution of can be computed.Next, the norm of the displacement increment is evaluated as Algorithm 3. Contact displacement evaluation for the current step t n+1 at iteration k in Equation ( 22) and R g (1)  i in Equation ( 24) for j = contact nodes do compute R f (1) i u (1) j in Equation (33) and R g(1) i u (1) j in Equation (34) end for j for j = structure nodes do compute R f (1) i u (1) j in Equation ( 33) end for j end for i for i = structure nodes do compute R s (1) i in Equation ( 18) for j = structure nodes do compute R s (1) i u (1) j in Equation (31)   end for j for j = contact nodes do compute R s (1) i u (1) j in Equation ( 31) end for j end for i assemble Lin R in Equation (36) solve for Δu (k) using SVD in Equation (41) go to Algorithm 4 Furthermore, the norm of the residual is computed for the updated displacement in the same manner, written as Convergence is achieved when both the displacement increment and residual norms approach zero or are lower than the determined zero tolerance.The Newton iteration in the proposed contact method is presented in Algorithm 3.

Contact node evaluation
The geometric non-penetration condition checking to obtain the contact node as well as the initial contact displacement enforced to it has been explained in Section 2. The pseudocode of this geometric checking algorithm can be found in Algorithm 1. Furthermore, the governing equations of the frictional displacement-driven contact, their derivatives, and the computation of iterative solutions have been previously described in this section.The Newton scheme for equilibrium iteration is provided in Algorithms 2 and 3. Next, the change in the contact node set during each analysis step must be considered to complete the proposed contact method.Since the out-of-plane movement is restricted by Equation ( 23), the reopening condition is evaluated after the displacement increment is obtained by evaluating the normal force as the result of the change of displacement.
Taking from the definition of normal force, which is a force exerted by the surface due to contact occurrences that work in the normal direction to prevent solid bodies from penetrating, then the reaction force working on the slave node in contact should point in the same direction as the normal of the surface, as depicted in Figure 7. Hence, the reopening condition on each contact node due to the contact displacement at the current iteration occurs if the requirement of Algorithm 4. Equilibrium check and contact node evaluation for the current step t n+1 at iteration k Requirement: (1) a set of nodes in contact {C} (2) displacement increment Δu (k) at {C} update u (k+1) n+1 = u (k) n+1 + Δu (k) compute contact force F c due to u (k+1) (2) ⋅ (n (2) ⊗ n (2) )F c ||(n (2) ⊗ n (2) is not met.Note that F c is the contact force on the contact node while n (2) represents the outward unit normal of the master surface corresponding to the current contact node.Equation (45) shows that during the iteration, if the reopening should occur, the unit normal force will be pointing in the opposite direction to the master surface outward unit normal.This condition means that instead of preventing penetration, the contact force working on the contact node pulls the node to stay in contact.Therefore, the node should be released from the contact node set.When the contact not set is changed, it affects the residual equation for both the structure and contact nodes.Consequently, it is necessary to perform a new calculation of residuals and their derivatives.If this is the case, then equilibrium iteration must be carried out from the beginning.This contact node evaluation and the equilibrium checking for each iteration at a certain load step are shown in Algorithm 4.
Next, the geometric non-penetration condition checking is repeated once the acceptable convergence level is achieved to guarantee the admissible state.In case penetration occurs when equilibrium is achieved, the step provided in Algorithm 2 located outside the iteration loop is performed again.Finally, the analysis can proceed to the subsequent loading step when the contact displacement and contact force fulfill both the convergence level and the admissible state.

NUMERICAL EXAMPLES
Several numerical examples of different contact cases are provided to demonstrate numerical stability, reliability, and accuracy of the proposed contact formulation.The numerical tolerance used in this formulation are 10 −13 and 10 −16 for the residual and geometric checking, respectively.

Elastic punch on a rigid plate
The analysis of the frictionless contact problem between an elastic punch and a rigid foundation is carried out.Zavarise and De Lorenzis 17 utilized this simulation to assess the performance of the NTS contact element in the penalty method.
In this analysis, uniform and non-uniform meshes are used.The regular NTS contact element failed the contact patch test with maximum and minimum errors of 50% and 12%, respectively.This problem arises from the constant penalty parameters assigned to the slave nodes, resulting in non-uniform gaps.This condition induces non-uniform contact pressure because it is a function of contact area, gap, and penalty parameter in the penalty method.Therefore, the same setup using a 50 × 50 mm 2 block is simulated to see the behavior of the NTS within the displacement-driven contact algorithm.Figure 9 illustrates the material properties, loading, and boundary conditions of this numerical example.Uniform and non-uniform discretizations are employed, as shown in Figure 10.The displacement distribution in x-and y-direction is shown in Figure 11.It can be seen that the displacement contours match reasonably well for all discretizations.Furthermore, an exact constant contact pressure in the opposite direction is expected since a constant loading pressure of 0.8 MPa is applied to the deformable block.The contact distribution on the contact surface is provided in Figure 12.A slight pressure inconsistency is observed at both the uniform and the non-uniform mesh.A maximum error of 0.89% with respect to the expected value is seen at the uniform mesh.On the other hand, the non-uniform mesh gives a higher maximum error of 1.41%.However, the minimum error of both discretizations is 0.03% with respect to the exact solution.This error value is considerably smaller than the previously mentioned result.This improvement in contact pressure arises due to an exact geometric non-penetration ensured by the contact displacement in the proposed algorithm.In this case, the contact force on each node is associated with the reaction force due to the uniform contact displacement.Therefore, employing a general NTS concept within the displacement-driven contact algorithm delivers a significantly more accurate contact pressure.

Rigid punch on an elastic foundation
A contact problem between a rigid punch and an elastic foundation, as proposed by Papadopoulos 28 , is recreated to show the presented contact model behavior in a singular problem that arises in the flat punch indentation case.Two different cases of this problem configuration are evaluated to validate the trustability of this contact algorithm.The first case involves displacement loading on a frictionless surface, where the stability of the contact algorithm are studied by

F I G U R E 12
Elastic punch on a rigid plate: Contact zone pressure distribution.
refining the mesh.The second scenario involves testing the friction feature of the proposed contact method through quasi-static frictional full sliding contact.For validation purpose, there is a straightforward final equation available for the analytical solution of this frictional contact case.

Frictionless contact
A block of a rigid flat indentor is driven into an elastic foundation by a vertical displacement v = 0.05mm.The elastic parameter for the foundation, the loading, and the boundary condition used in this contact problem are described in Figure 13A.Note that the analysis is carried out with an additional restraint in the out-of-plane direction to force a plain strain condition.In this analysis, a 8-node brick-type finite element is used to discretize the deformable foundation, where the behavior of four different mesh densities of the elastic foundation on the contact area will be compared.Figure 13B illustrates the discretization where five nodes of the deformable body are in contact.The other discretizations are not provided since the visibility of the finer mesh utilized for the contact area is compromised when applying mesh refinement levels of two, four, and eight times.As mentioned before, the node-to-segment contact element is incorporated into the contact formulation.Thus, in this example, the mesh sensitivity of the contact pressure in the contact interface is evaluated since more nodes are in contact if a finer mesh is employed.
Figure 14 shows the pressure at the deformable body surface, which is in contact with the rigid block due to the indentation and the zoomed-in part.The result shows a considerably reasonable pressure behavior, particularly on the middle area of the contact zone as the finer mesh is utilized.Furthermore, mesh refinement successfully fixes the numerical error around the singularity at the end of the contact zone.As an outcome, a more accurate pressure behavior due to flat indentation is obtained, indicated by a steeper gradient at both ends of the curve.In addition, the convergence study of the contact force with respect to mesh refinement shows that larger mesh density gives less discrepancy with respect to the force with a descending gradient.This behavior is captured in Figure 15.Therefore, it can be concluded that the algorithm delivers a reliable result satisfactorily.The final deformation presented by displacement in the y-direction is depicted in Figures 16.The singular problem at the edge of rigid flat punch is also depicted in Figure 17 as the Cauchy stress component  YY distribution of the elastic foundation.The stress contour shows the singularity problem at the edge of the rigid flat punch, where stress concentration is observed.

Frictional full sliding contact
A similar contact case with the same dimensions, material properties, and boundary conditions is analyzed in this example.The largest mesh density of the elastic foundation used in the previous example, with 40 nodes in contact zone,

F I G U E 14
Frictionless contact: Contact pressure distribution.

F I G U R E 15
Frictionless contact: Convergence study in terms of contact force.

F I G U R E 16
Frictionless contact: Final displacement of the elastic foundation in y-direction.

F I G U R E 17
Frictionless contact: Final Cauchy stress distribution of the elastic foundation  YY .
is employed.However, in this case, friction between surfaces is considered, which is set to  = 0.3, and a different loading scenario is applied.First, an external 3 N force establishes the contact.Then a tangential displacement is applied in addition to the force at the top of the rigid flat punch as described in Figure 18.Note that the applied force is considerably small to guarantee small deformations of the elastic foundation.Furthermore, the slow sliding is ensured to neglect the inertial effect.Thus, the quasi-static assumption is held, yielding that the proposed algorithm can solve the problem.This example is carried out to ensure the capability of the proposed contact algorithm in frictional contact, which is proven by comparing the result to an available analytical solution for frictional full-sliding contact. 29The contact pressure is computed using the analytical solution as where p(x) denotes the contact pressure at position x, P denotes the applied force in normal direction, and a denotes the half length of the contact zone, which is equal to the half dimension of the rigid block in contact with the elastic foundation in the flat punch case depicted in Figure 18.In this example, 2a is equal to 6mm, as illustrated in Figure 13. is a function of the friction coefficient  and Poisson's ratio , which is written as The comparison is depicted in Figure 19.The figure includes both an overall view and a zoomed-in part.The numerical result matches the analytical solution but not in the edge area, which can be caused by discretization and numerical rounding.Moreover, the dotted line denotes the analytical contact pressure for a contact scenario without friction to see the change of pressure due to friction.As the block slides in positive x-direction, the friction causes the pressure on the leading side to be lower than the pressure on the trailing side.This outcome is consistent with previous research. 30,31The comparison to the frictionless case demonstrates that friction does not significantly impact the contact traction when a slight indentation is reassured to maintain a small deformation on the elastic foundation.

Hertzian contact
Despite the moderately complex derivation, the renowned Hertzian contact theory leads to a relatively simple final solution. 32The Hertzian equation provides an analytical solution that is suitable for practical use or to validate a

F I G U R E 18
Frictional full sliding contact: Loading scenario.

F I G U R E 19
Frictional full sliding contact: Comparison between analytical and numerical contact pressure distribution.
numerical analysis for certain contact cases that fulfill the Hertzian contact limitation.Therefore, in this example, an analysis of the frictionless contact problem between a linear-elastic isotropic cylinder (E = 300 MPa, v = 0.25) and a rigid flat surface is performed.A uniform pressure p = 0.8 MPa is applied at the top of the cylinder, which is considerably small, to ensure the small strain condition.The loading and boundary conditions are designed to meet the Hertzian contact case so that a proper comparison can be accomplished.In this example, only a quarter of the cylinder is modeled for the numerical analysis due to symmetry, as shown in Figure 20.The numerical analysis results of the displacement-driven contact algorithm are validated using the analytical solution for the Hertzian cylinder contact case.3][24] The analytical contact zone per unit length of cylinder a and the normal traction p c equations 30 read where R * is the representative radius between two bodies in contact computed as where R 1 and R 2 denote the radius of cylinders one and two, respectively.On the other hand, E * represents the modulus of elasticity between two bodies in contact which is expressed by where E 1 and E 2 are the moduli of elasticty, while  1 and  2 are the Poisson's ratio of the two bodies.The numerical analysis of the Hertzian cylinder contact scenario is examined using two distinct elements.The first includes linear tetrahedral elements with one integration point.The second involves fully integrated linear brick elements.The cylinder body is discretized, with a higher density applied to the contact zone.The initial element size on the contact zone is 0.2mm, called mesh 1.Then, mesh refinement is applied on the contact zone for both types of discretization to study the influence of spatial discretization concerning accuracy.Then, smaller element sizes of 0.1 and 0.05mm are applied, named mesh 2, and mesh 3, respectively.Figure 21 shows the mesh created using tetrahedrons and hexahedrons for mesh 1.The finer discretizations are not provided due to visibility concerns.Both elements successfully establish a quadratic convergence behavior, which indicates that the formulation is derived correctly and consistently.Figures 22 and 23 show the convergence feature for an arbitrary load step for the tetrahedral and the hexahedral mesh, respectively.Note that the two lines illustrate the change of nodes in the contact set during the equilibrium iteration, which is allowed by the algorithm as described in Section 3.6.The blue dashed line indicates the equilibrium point of log(10 −13 ), the highest value of the logarithmic residual deemed acceptable.Next, the contact force due to the applied pressure is checked.Due to the absence of friction, the total contact force is equal to 8 N.This reaction force is equivalent to the total force created  by the applied pressure load but in the opposite direction, as expected.Figure 24 depicts the normal traction on the contact zone using tetrahedral elements.Some spatial oscillations are observed in the numerical solution when using tetrahedrons.Based on his behavior, it is apparent that the tetrahedral element with a single Gauss point is unable to match the analytical solution, even with a finer discretization in the contact zone.In contrast, the linear hexahedral element fits the analytical result sufficiently even for the coarse mesh, except for the outermost of the contact zone, see Figure 25.The two results indicate that the proposed algorithm is able to produce results that satisfy the analytical solution by utilizing the suitable element.Furthermore, a spatial convergence with an h-refinement measured in terms of the L 2 -norm of displacement error is studied.A minimum of three points is demanded for this study.The requirement is fulfilled using the three meshes of both elements.The error norm is calculated by where ū is the nodal displacement in the contact zone resulting from the numerical analysis of the proposed contact algorithm, and u is the nodal displacement in the contact zone obtained from the analytical solution using Equation (48) at the point where ū is evaluated.The error value is plotted with respect to the number of elements in one length unit [1∕mesh size] using a logarithmic scale as presented in Figure 26.It can be seen that finer elements reduce displacement error with a smaller gradient.Based on the three data points, this slope is expected to be lower as the mesh density increases.Overall, this suggests that at some point, with further mesh refinement, there will be a constant displacement error.This behavior implies that the displacement, as a body reaction due to the same loading and boundary condition, remains unchanged.In conclusion, the spatial convergence study confirms that the contact algorithm is able to produce a reliable result.

F I U R E 24
Hertzian contact: Results for several discretizations using tetrahedral elements.

F I G U R E 25
Hertzian contact: Results for several discretization using hexahedral elements.

F I G U R E 26
Hertzian contact: Spatial convergence study in terms of L 2 -norm of displacements.
Furthermore, the analytical contact zone obtained by the Hertzian equation is a = 0.798mm.The difference between the analytical and the numerical solution at the outermost contact zone can be reduced using the finer mesh.This problem is not discussed further since the other parts conform with the analytical solution.Nevertheless, the results are adequate to conclude that the displacement-driven contact model result successfully captures the correct behavior of the frictionless cylinder Hertzian contact problem.The final distribution of the Cauchy stress component  YY and displacement in the y-direction are shown in Figure 27A,B, respectively.The stress distribution analysis reveals that the stress decreases to zero toward the end of the contact zone.This expected behavior contrasts with the response in the previous example, where a flat punch is used.

4.4
Frictional sliding contact The following numerical example, inspired by Fischer, 22 aims to validate the stick/slip behavior and to observe the numerical stability of the proposed frictional contact model.Contact between a hyperelastic neo-Hookean cube (E = 100MPa, v = 0.3) and a rigid plate is analyzed.As Coulomb's friction law is employed in the formulation, the friction in the interface is represented by a constant coefficient of friction.Three different coefficients of friction are used, which are  = 0.0,  = 0.1, and  = 0.4, to see how friction influences the response of the deformable body during contact loading.First, a displacement of u y = −0.1mm is applied at the top of the block to push the body into the rigid surface.Then, an additional displacement load u x = 17.5mm is applied at the top of the compressed body to slide it along the rigid surface.The rigid body surface has a curvy geometry to give different contact settings to the deformable body along the sliding motion.The rigid surface is fixed, as well as the top surface of the block, where the displacement load is applied.
A more detailed description of loading and geometry is illustrated in Figure 28.The simulation is performed using linear brick elements for the discretization of the deformable body and the surface of the nonlinear rigid surface as pictured in Figure 29.
The numerical result is provided in terms of reaction forces for the three-dimensional Cartesian space instead of normal force and tangential force since the normal and tangential direction changes during the analysis.The reaction force of the deformable body is computed at the top of the block.Figure 30 depicts the reaction force in the x-direction for each coefficient of friction.It is plotted for the displacement load in the x-direction, which is measured from the point x equal zero as shown in Figure 28.Here, the influence of the friction coefficient on the reaction force can be seen clearly since the block is slid in the x-direction.As expected, the larger the friction coefficient is applied to the contact interface, the larger the reaction force is obtained.The reaction increases even more when the block is in the inclination area.In contrast, Figure 31, which depicts the reaction force in the y-direction, shows slight force reductions

F I G U R E 28
Frictional sliding contact: Geometry and loading.

F I G U R E 29
Frictional sliding contact: Discretization of the deformable body and the nonlinear rigid surface.

F I G U R E 30
Frictional sliding contact: Reaction force in x-direction for different coefficients of friction using the displacement-driven contact approach.
as the coefficient of friction increases, particularly from the middle of inclination to the middle of the declination area.This change can be explained as the change of the active contact zone due to friction.Due to friction, the stress at the front area of the contact interface is higher compared to the rear, resulting in an increased normal force at the front area.Since the relation between friction and normal force is linear, the force required to move the deformable body relative to the rigid surface also rises.At this point, the force due to the displacement load is inadequate, resulting in the stick condition at the front area.Then, the normal force at the rear contact interface reduces to balance the force or even lifted at some point, reducing the reaction force in y-direction.However, as the displacement at the top continues, a larger tangential force is produced.At a certain point, it might be larger than the resistance due to friction, which causes the slip condition.
To evaluate the result from the proposed algorithm, the same task is analyzed using the well-known penalty method.Figure 32 compares the reaction force in x-, y-, and z-direction for frictionless contact, analyzed using the penalty method (P) and the displacement-driven contact approach (DD).It can be concluded that the result from the contribution at hand agrees with the penalty method for the frictionless case.However, a slightly different reaction force in y-direction is observed.The penalty parameter is set to 10 5 MPa for this analysis to maintain numerical stability.However, a small penetration occurs at the surface, in the order of 10 −5 mm.While for the displacement-driven contact approach, the gap between contact surfaces falls under the tolerance value of 10 −13 mm.Furthermore, the comparison of contact responses for frictional contact is carried out.However, for the frictional case under a nonlinear contact surface, non-convergence is F I G U R E 31 Frictional sliding contact: Reaction force in y-direction for different coefficients of friction using the displacement-driven contact approach.

F I G U R E 32
Frictional sliding contact: Comparison between penalty and displacement-driven contact method for frictionless case.observed in the computation using the penalty method, and the analysis has to be terminated.In addition, the instability is observed to present earlier as a larger coefficient of friction is used, as can be seen in Figures 33 and 34.The simulation for friction coefficient  = 0.1 is stopped at sliding displacement u x = 3.33mm.Moreover, the simulation for friction coefficient  = 0.4 is terminated earlier at sliding displacement u x = 2.74mm, at the beginning of the inclination as depicted in Figure 35.
Furthermore, compared to the smooth reaction force by the penalty method, the proposed algorithm shows a non-smooth behavior, particularly at the displacements of 5 to 10 mm.The oscillation of the reaction force in x-direction is most enunciated for  = 0.1 and less noticeable for  = 0.4, as presented in Figure 30.The behavior of the algorithm used can justify the reason for this behavior.First, the closest point projection is utilized, which results in the displacement normal to the surface to achieve the geometric non-penetration condition.Then, the contact displacement on the surface or perpendicular to the normal of the surface is iterated within the Newton scheme.The singular value decomposition solves the contact displacement increment in each iteration, which gives the final tangential displacement in the equilibrium condition.This technique provides the closest result to the real solution, which might yield a slight difference depending on the problem characteristics.In this numerical example, the displacement normal to the surface mainly occurs in y-direction during block movement, making the reaction force in this direction (see Figure 31) smoother than in x-direction.
For a friction coefficient of  = 0.1, where the oscillation is most enunciated, the block experiences both sticking and sliding conditions.During the iteration, this condition is checked by the residual of the contact node.Furthermore, change in tangential displacements leads to alteration in forces, resulting in different states of stick and slip.Defining   tangential displacements is more challenging under these circumstances.When the friction coefficient is set to  = 0.4, the deformation pattern of the moving body indicates that the contact nodes frequently remain in a "stick" condition, where the restriction of tangential displacements is more apparent.As a result, a smoother curve is obtained.On the other hand, in the frictionless case, the deformation pattern during the movement implies that all the nodes in contact always slide.The deformation patterns, which show stick and slip for the three simulations at several load steps, are shown in Figures 36 and 37.The current position x 1 , x 2 , x 3 , and x 4 represent the initial position, where the block is compressed, the beginning of inclination, the top of inclination, and the declination zone, respectively.Here, the influence of friction on the behavior of the deformable body can be clearly observed.For instance, in the frictionless contact case, the node of the block surface can move relative to the rigid surface (see Figure 36A).As friction rises, this relative displacement decreases, as shown in Figure 36B.In Figure 36C, stick condition is observed, denoted by zero displacement in the contact zone.Additionally, in Figure 37I, the rear part of the deformable body is lifted due to the considerable friction, while this phenomenon does not occur when a smaller friction coefficient is used as depicted in Figure 37G,H.To summarize, despite the non-smooth behavior that requires improvement, the displacement-driven contact model achieves good numerical stability for this particular case.In addition, the proposed algorithm accounts for the effect of friction on the behavior of the deformable body, such as the stick and slip condition, as well as the reopening.Furthermore, the simulation effectively meets both the equilibrium and the geometrical non-penetration requirements.

CONCLUSIONS
The proposed contact method provides a displacement-driven approach for solving a frictional contact problem of solid bodies under quasi-static deformation.The solution is currently limited to interaction between a deformable and a rigid body.However, the applicability is not restricted to a specific constitutive model and a small deformation case.Moreover, an elastic or inelastic material description might be incorporated.The contact solution is implemented into a finite element method where the node-to-segment contact element formulation is incorporated for the contact boundary discretization.The displacement-driven approach means that the direction and the magnitude of the initial contact force at the contact boundary are computed according to the required displacement to ensure the geometrical non-penetration conditions as a prescribed value which is then solved in an iterative manner.The proposed solution method is derived straightforwardly from the three residual equations, which are the balance of linear momentum, friction law, and geometrical constraint to ensure a geometrical non-penetration condition of the interacting bodies.To solve the contact problem within a Newton solution scheme wihtin the FEM, linearization of these residual equations with respect to the unknown displacement field is carried out consistently to provide numerical stability and quadratic convergence behavior.
Several numerical examples are solved using the proposed algorithm, including the patch test, mesh study, validation using analytical solutions, and comparison to another contact method.The results indicate a considerable degree of numerical stability, accuracy, and reliability.Moreover, the stick and slip condition at large sliding in frictional contact is well depicted.Furthermore, the structure of this solution method provides the advantage of developing the contact model relatively easily.Therefore, extending the current contact formulation for the solution method of contact problems between deformable bodies is necessary to cover more general contact mechanics tasks.

F I G U R E 2 F I G U R E 3
Outward unit normal computation of the master surface.Graphical illustration of an inadmissible state detection between three-dimensional bodies.
0 then go to Algorithm 3 else go to the next load step end if F I G U R E 8 Illustration for the geometrical constraint in contact.

F I G U R E 9
Elastic punch on a rigid plate: Material properties, loading, and boundary conditions.(A) (B) F I G U R E 10 Elastic punch on a rigid plate: (A) Uniform and (B) non-uniform meshing.U R E 11 Elastic punch on a rigid plate: Displacement contour for (A) uniform and (B) non-uniform mesh.

21
Hertzian contact: Discretization of the quarter-cylinder using (A) structured hexahedral and (B) random tetrahedral elements.F I G U R E 22 Hertzian contact: Convergence behavior for tetrahedral discretization in terms of residual.F G U R E 23Hertzian contact: Convergence behavior for hexahedral discretization in terms of residual.

27
Hertzian contact: Final state of (A) Cauchy stress component  YY and (B) displacements in y-direction for mesh 3 discretization using 8-node brick elements.

F
I G U R E 33Frictional sliding contact: Reaction force in x-direction for different coefficients of friction using the penalty method.

F
I G U R E 34Frictional sliding contact: Reaction force in y-direction for different coefficients of friction using the penalty method.
Frictional sliding contact: Final deformation pattern for (A)  = 0.1 and (B)  = 0.4 using the penalty method.

F
I G U R E 36Frictional sliding contact: Deformation pattern for x-direction at certain load steps for each coefficient of friction.
Frictional sliding contact: Deformation pattern for y-direction at certain load steps for each coefficient of friction.