Improved area regularization technique for penalty-method-based node-to-segment contact analysis

Among numerous approaches that have been presented for computational contact analysis, the node-to-segment (NTS) approach is widely adopted in industrial applications, particularly with the penalty method. However, the pure penalty-method-based NTS approach fails the patch test even with the area regularization (AR) technique. There have been relatively few solutions for this problem compared to those for Lagrange-multiplier-method-based NTS approaches. This paper presents an improved AR technique and attempts to make the penalty-method-based NTS approach pass the patch test. The proposed approach improves the accuracy without the need to introduce additional geometry to the conventional penalty-method-based NTS approach. The entire competence area of a slave node is not directly utilized as a penalty weighting factor. Instead, a part of the competence area is extracted in accordance with the projection from a master node to a slave segment. In a complementary manner, a two-pass algorithm is used for the selected master nodes to achieve segment-wise force/moment equivalence. The results obtained using the proposed approach are introduced in a two-dimensional frictionless context. It is observed that the proposed approach passes the patch test. Furthermore, a significant improvement in accuracy compared to the conventional NTS and NTS-AR approaches is observed. We aim to extend the approach to three-dimensional applications in the future work.

The contact algorithms for a discretized interface are divided into several categories. The node-to-node (NTN) algorithm [3,4] defines the contact constraints between the nodes of slave and master bodies. Conventionally, its application has been restricted to a matching mesh with infinitesimal deformation. To increase the application scope, the NTN approach with a variable node-element technique has been developed for a nonmatching mesh with large deformation [5][6][7].
In addition to purely node-based approaches, Hughes et al. [8] and Hallquist [9] proposed the node-to-segment (NTS) approach for contact simulation. In this approach, a nonpenetration condition between a slave node and its projection on a master segment is defined. It enables the contact simulation of a nonmatching mesh with large deformation and sliding [10]. The classical one-pass NTS approach projects a slave node onto a master segment in a single direction. However, it fails the patch test, as reported by Taylor and Papadopoulos [11]. Several methods have been proposed to resolve this issue, including the two-pass algorithm, combinations of lower/higher-order interpolation functions [12], virtual slave node technique [13], and the contact domain method [14,15].
In contrast to node-wise approaches, the segment-tosegment (STS) approach involves the numerical integration over contacting segments to discretize a contact interface. Since the introduction of the contact segment concept by Simo et al. [16], the STS approach has been extended to large-deformation situations [17]. This approach passes the patch test and satisfies the stability requirements [18]. The continuous contact constraints were reported to provide better accuracy compared to the node-wise constraints [19]. Recent studies on contact algorithms have focused on the STS approach, which is based on the mortar method. The mortar method is based on the domain decomposition idea [20]. This method enforces traction continuity along an entire interface using a weak formulation. Numerous studies [21][22][23][24][25][26][27][28] have examined the mortar contact algorithm and its implementation in commercial software [29,30].
This study primarily focuses on the NTS approach. Although the mortar method is a well-established technique, the NTS approach is still widely used for practical and industrial applications. The computation cost of the node-wise projection that is integrated in the NTS approach is typically smaller than that of the integration point-wise projection used by the STS and mortar approaches. Although the NTS approach does not yield high accuracy, it is indispensable for nonlinear explicit analysis that requires a considerable computational cost. In addition, the lower-order interpolation function [12], which is the main drawback of the conventional NTS approach, is suitable for the linear-order solid/shell elements combined with single-point integration [29,31,32].
Another concern of the NTS approach is the formulation of contact constraint enforcement. Although the Lagrange multiplier method provides exact contact enforcement, the penalty method is widely used in practical situations [29]. The penalty method permits slight penetration of a slave body into a master body. Subsequently, the contact force is enforced by multiplying the penetration with a penalty parameter [33]. Despite the minor violation, this approach is simple and does not introduce additional variables such as the Lagrange multiplier. Regardless of whether the penalty or Lagrange multiplier method is employed, the conventional NTS approach still does not pass the patch test. Moreover, this problem cannot be resolved via mesh refinement. Studies that are aimed at improving the accuracy [12,14,15] are not purely based on the penalty method, except the work of Zavarise and De Lorenzis [13]. Most studies that consider the patch test criterion [18,34,35] use the Lagrange multiplier method.
The penalty parameter weighting may be used to improve the accuracy of the penalty-method-based NTS approach. The area regularization (AR) technique [36] employs the area of competence of a slave node as a weighting factor. Lee et al. [37] showed that the NTS approach with the AR technique (NTS-AR) passes the three-dimensional patch test when all slave nodes are projected onto an area with a single master element. They performed a three-dimensional explicit analysis to a collision scenario between a stiffened plate panel and simplified collider. The agreement between the predictions that are obtained using the NTS-AR and STS approaches was better than that between the conventional NTS and STS approaches. LS-DYNA, which is a commercial finite element software, enables penalty parameter weighting by area or mass, which is recommended for metal forming simulations [38,39]. However, Zavarise and De Lorenzis [13] showed that the NTS-AR approach fails the patch test when the projection of slave nodes is associated with multiple master segments. Furthermore, the accuracy of the two-pass NTS-AR approach was low. The authors of Ref. [13] added virtual slave nodes to each slave segment, in addition to the AR technique, for the NTS approach to pass the patch test. Then, the area of competence was projected onto the master segments for each virtual slave node. Finally, the concentrated contact force of the virtual slave nodes was transformed to (i) the uniform pressure on the projected area and (ii) concentrated force on the master nodes. Because the modified algorithm passed the patch test, it was referred to as the virtual slave node-to-segment technique passing the patch test (VTS-PPT).
In this paper, an improved AR technique is proposed for penalty-method-based NTS contact analysis. The entire competence area of a slave node is not directly utilized as a penalty weighting factor. In contrast, a part of the competence area is extracted by considering the geometry of a contiguous master segment. In a complementary manner, a two-pass algorithm is used for the selected master nodes, which utilizes the remaining competence area as the weighting factor. The performance of the proposed NTS analysis is demonstrated for a two-dimensional frictionless contact situation. The achievement of local force/moment equivalence is illustrated by performing various patch tests. Finally, an application of the proposed approach to various numerical examples are presented. The proposed approach improves accuracy without the need to introduce additional geometry to the conventional penalty-method-based NTS approach.

Theoretical background
This section introduces the theoretical background of frictionless NTS contact analysis based on the penalty method. The inaccuracy of conventional NTS contact analysis is demonstrated through a patch test example. The details of frictionless NTS contact analysis are presented in a twodimensional context with a linear finite element.

Contact formulation with penalty method
Let B α denote two solid bodies in contact, as illustrated in Fig. 1. Each solid body is enclosed by surfaces α , where superscript α ∈ {sl, ms} indicates the slave and master bodies. The initial position of both solid bodies, x α 0 , is displaced to the current position, x α = x α 0 +u α , by deformation, where u α denotes displacement. c = sl c = ms c denotes the contact surface between the solid bodies.
Prior to determining the contact surface, c , the contravariant convective coordinate, ξ 1 , is employed to express the master surface coordinates, ρ = ρ ξ 1 . Then, a distance function, d ξ 1 = x sl − ρ ξ 1 2 , is defined between the fixed location, x sl , at the slave surface, sl c , and an arbitrary location, ρ, at the master surface, ms c . The contact pair,x ms = ρ ξ 1 , which is the closest point projection of x sl on ms c , is obtained by solving the distance minimization of Eq. 1.ξ 1 is the convective coordinate ofx ms , which is the solution of Eq. 1.
ρ 1 = dρ/dξ 1 is the covariant tangent vector at the master surface. A normal gap function, g N , is defined as follows: Here, n denotes the surface normal unit vector atx ms , and it points outward from the master solid body.
In the penalty method, the normal contact traction, t N , is typically expressed, as shown below: N is the penalty parameter. H (−g N ) is the Heaviside function, which is zero for g N > 0 and unity for g N ≤ 0. In contrast to exact contact enforcement, the penalty method permits penetration, i.e., a negative gap function, g N . An appropriate penalty parameter must be selected because small or large penalty parameters lead to nonphysically large penetration or a numerically ill-conditioned problem, respectively.
The contact constraints contribute to the variation formulation, δW = δW int +δW ext +δW c , where δW , δW int , δW ext , and δW c indicate the total, internal, external, and the contact virtual work, respectively. For frictionless situations, δW c is formulated as follows: A more detailed explanation of contact formulation is provided in Refs. [33,40].

Contact element description
The NTS contact element is defined as a single pair of a slave node and master segment. A description of this process is shown in Fig. 2, where S at position x sl is the slave node, and M 1 and M 2 at positions x ms 1 and x ms 2 are two master nodes comprising the linear master segment, respectively.
In the discretized field, the virtual work in Eq. 4 is approximated, as shown below: The elemental contact force vector, F ce , is expressed using Eq. 7.
A detailed derivation and review of the conventional NTS approach are presented in Ref. [41]. The linearization of the contact virtual work is explained in "Appendix A".
The elemental contribution of Eq. 7 is assembled into the global equation. For instance, an equation for nonlinear static analysis can be written as follows:  [42]. The Newton-Raphson iteration scheme is employed with the force convergence criterion of R T R/F T ext F ext < 1 × 10 −12 for nonzero F ext and the energy criterion of U T R/ U T 0 R 0 < 1 × 10 −20 , otherwise. Here, U is a global displacement increment vector. U 0 and R 0 indicate the values of U and R that are obtained using the initial Newton-Raphson iteration, respectively.

Limitation of conventional approach
The limitation of the conventional NTS approach is demonstrated using a simple patch test. Figure 3 illustrates patch test Case A, wherein the contact surface is uniformly discretized and all the master nodes do not intersect the projection of the slave segment. The other patch test cases are consecutively shown in a row for the necessary explanation in the remaining sections. In Case A, the uniformly distributed pressure p 0 = 1 Pa acts on the upper rectangular block, which presses the lower block. At the interface between the two blocks, the upper and lower surfaces represent the slave and master surfaces, respectively. The material properties of the solid body are an elastic modulus of E = 10 6 Pa and a Poisson's ratio of ν = 0. The magnitude of the analytical contact pressure is p = p 0 along the entire contact surface. In addition, the analytical penetration between the slave and master segments is uniform. In the conventional NTS approach, equal penalty parameters are assigned to each slave node, as shown in Fig. 4. This approach cannot predict uniform contact pressure. As explained by Zavarise and De Lorenzis [13], a constant penalty parameter cannot produce uniform contact pressure for uniform penetration. In contrast, uniform penetration cannot be obtained for uniform contact pressure. Therefore, local force and moment equivalence cannot be achieved for the individual contact segment. Figure 5 shows the numerical results obtained, which accompanies the nonuniform contact penetration illustrated in Fig. 6. The two-pass algorithm also experiences such nonuniformity. It should be noted that the contact pressure shown in Fig. 5 strongly depends on the penalty parameters even though the problem is not ill conditioned. This is because for the nonequivalent local force and moment, the

Improvement over constant penalty parameter
An AR technique that varies the penalty stiffness was proposed to overcome the disadvantages of the conventional NTS approach [36]. Rather than considering only a single node as the slave component, as shown in Fig. 2, the AR technique extends the standpoint to the geometry of neighboring slave segments. Figure 7 shows the physical interpretation of the AR technique using a slave surface that comprises two segments, S 1 S and SS 2 . In Fig.7a, the distributed penalty parameter of the constant stiffness,ˆ N , is assigned along the area of segments. Then, as shown in Fig. 7b, the equivalent Patch test Case A with the conventional NTS approach: deformed configuration ( N = 2 × 10 6 N/m). Scale factor for the displacement is 2 × 10 5 Fig. 7 AR technique: a equally distributed penalty parameter, b equivalent nodal discrete penalty parameter discrete penalty stiffness, N , of the slave node is obtained by multiplying the competence area andˆ N . For each slave node, the competence area, i.e., penalty weighting factor w, is the sum of half of the areas of the neighboring slave segments.
The AR technique can be employed by substituting N e = N w e into the contact force vector, F ce , in Eq.7. An additional linearization contribution for the consistent tangent stiffness is derived in "Appendix A.2", which is neglected in this study.
In Case A, the penalty parameter depends on each slave node when the AR technique is employed. Therefore, the NTS-AR approach can predict uniform contact pressure for uniform penetration, as described in Fig. 8.
The numerical results shown in Fig. 9 demonstrate the advantage of the NTS-AR approach over the conventional NTS approach. The accuracy of the NTS-AR approach is maintained for patch test Case B in Fig. 10, wherein the contact surface is discretized in a biased manner. The relevant numerical results are shown in Fig. 11. Except for the biased discretization, Case B is similar to Case A in that all master nodes are matched to the slave node.

Limitation of AR technique
The NTS-AR approach fails the patch test when the projection of the slave segment spans multiple master segments. This further results in two typical adverse effects, i.e., nonequivalent moment and the mismatch of the competence area, which were also reported by Zavarise and De Lorenzis [13]. Figure 12 shows patch test Case C, wherein one of the master nodes is located in the middle of the slave segment projection. Figure 13 shows the contact nodal force and pressure when the NTS-AR approach is used with uniform penetration. Although the slave surface is under uniform pressure, the moment equivalence is not sufficed for the individual master segment, as shown in Fig. 14. This nonequivalence stems from the misaligned loading point of Patch test Case C with NTS-AR approach: a penalty parameter distribution, b contact nodal force, and c pressure distribution for the uniform penetration the slave node to the master segments. Thus, the addition of a penalty stiffness acting point or a change in its magnitude is required to achieve local moment equivalence.
The patch test Case D in Fig. 15 shows a more fatal situation compared to Case C, which is the mismatch of the competence area. As illustrated in Fig. 16, the area of the left master segment is too narrow for the projection of the competence area. Consequently, the contact force is skewed toward the left master segment, and a combination of force/moment nonequivalence occurs for the individual master segment. Therefore, the geometry of the contiguous master segment should be considered to overcome this limitation.

Improved AR technique
In this section, we describe a modified AR technique to resolve the aforementioned limitations. The numerical procedure of the modified technique is presented in detail.

Procedure
The proposed approach comprises three components: slaveto-master projection, master-to-slave projection, and penalty parameter estimation.
The slave and master surfaces are shown in Fig. 17. The slave and master surfaces are discretized by nodes S i and M i , where i ∈ {1, . . . , 5}. Here, S * i and M * i are the projections of nodes S i and M i , respectively.

Slave-to-master projection
The procedure begins with the slave-to-master projection in Algorithm 1. Here, E α n and E α s denote the list of contact nodes and segments, respectively, where α ∈ {sl, ms}.Ẽ α n ⊂ E α n indicates the list of contact nodes with a contact pair.ξ 1,α I is the list of convective coordinates for contact segment I ∈ E α s . Furthermore,ξ 1 i I indicates the convective coordinate of the projection from the i-th node to segment I .
For each slave node, the contact-pair master segment is sought by the closest point projection, which is similar to the conventional NTS and NTS-AR approaches. If a contact pair is found, the convective coordinate of the projection from the slave node to the master segment is determined. The main difference is that the information of the projection on the master segment is utilized for the penalty parameter estima-  3. Thus, the proposed approach arranges the obtained convective coordinate in a segmentwise manner, as shown in Line 7 of Algorithm 1. This is an additional step in the node-wise storage of the conventional NTS approach, wherein a single convective coordinate is usually allocated for each slave node. Therefore, Algorithm 1 can be easily implemented in existing NTS-based contact analysis with minor modifications to the computational storage strategy.

Master-to-slave projection
The master-to-slave projection presented in Algorithm 2 is performed to consider the geometry of the contiguous master segment. i 1I and i 2I denote the list of two-end nodes for segment I corresponding to the convective coordinates of ξ 1 = −1 and 1, respectively. Pair (i) is the contact-pair segment for node i. InBetween (I 1 , I 2 ) denotes the list of the intermediate nodes between segments I 1 and I 2 .
It can be inferred that this part differs from the conventional approaches. For each slave segment, whether the projections of the two slave nodes are on different master segments is examined. If this condition is satisfied, the slave segment becomes the projection target, and the closest point projection is applied from the intervening master node to the target segment. Such a predetermination of the projection target helps eliminate the steps required to find the closest segment. If the projection is enclosed inside the slave segment, the master node and slave segment constitute an additional contact element. In Fig. 17, nodes M 2 , M 3 , and M 4 are selected as the master nodes to generate projections M * 2 , M * 3 , and M * 4 , respectively. In addition, the proposed algorithm considers a special situation wherein one of the slave nodes contains a contact pair and the other does not. Although this makes the proposed algorithm more complicated, it primarily consists of conditional statements to account for the general contact situation.
In contrast to the slave-to-master projection, wherein segment search is required, the proposed master-to-slave projection is applied to a predetermined slave segment. Thus, the computational cost is considerably lower than that of slaveto-master projection.

Penalty parameter estimation
The final procedure is Algorithm 3 for the penalty parameter estimation. Here, g i N e is the normal gap function between the i-th node and its contact-pair segment. Neighbor (i) denotes the list of neighboring segments for the i-th node, where the number of entries is either one or two. Further, l I is the area of segment I , and F i I ce is the elemental contact force vector defined by the contact element of slave node i and master segment I . As the proposed method involves a two-pass algorithm, i and I are obtained from the slave and master surfaces.

Algorithm 3 Penalty parameter estimation
Input : finite element representation, contact pair information,Ẽ α n , ξ 1,α For each active contact pair, the penalty parameter is computed according to the information of the segment-wise projection,ξ 1,α I , obtained using Algorithms 1 and 2. The penalty weighting factor for each contact element,w e , is determined, along with the master-to-slave element. The weighting factor is a competence area for the nodes wherein neighboring slave segments include no projection, e.g., nodes S 4 and S 5 in Fig. 17, similar to the conventional NTS-AR approach. Otherwise, the weighting factor is the sum of the modified areas, or the sum of half of the area between the nearest projection and the target node for the projected segment and half of the area of the nonprojected segment. The classification of the area on the basis of the projection and node enables the approximation of moment equivalence within a pair of node-projection sections, e.g., S 1 M * 2 and S * 1 M 2 in Fig. 17. After the penalty parameter, N e =ˆ Nwe , is obtained, the elemental contact contribution is computed using Eq. 7 for each contact element, without modifying the original formulation of the NTS approach.
The additional linearization contribution described in "Appendix A.3" should be considered if consistent tangent stiffness is required. A brief numerical investigation on the impact of additional linearization is provided in Sect. 5.2. Except that, similar to the NTS-AR approach, such additional linearization contribution is neglected in this study. Despite this inconsistency, the number of Newton-Raphson iterations is the same for the conventional NTS, NTS-AR, and proposed approaches for the given patch test.
The cost of penalty parameter estimation is almost negligible compared to that of the contact force vector. In the latter, the additional cost compared to the conventional onepass NTS approach increases proportionally with the number Fig. 18 Patch test Case C with the proposed approach: a penalty parameter distribution, b contact nodal force, and c pressure distribution for the uniform penetration of active contact elements between the master node and the slave segment.

Improvement over conventional approaches
Prior to the numerical simulation, the theoretical results of the proposed approach are examined based on the assumption that the penetration is uniform. Figure 18 shows the contact pressure distribution for Case C predicted by the proposed approach. The master-to-slave penalty parameter (red) rearranges the load distribution and resolves the moment nonequivalence that is induced in the NTS-AR approach in Fig. 13. In addition, the competence area mismatch in Case D is resolved, as illustrated in Fig. 19. This is because the distributed penalty parameter of the upper slave node in Fig. 16 is truncated by the proposed method before it spans multiple master segments. The master node handles the truncated penalty parameter in a complementary manner. Figures 20 and 21 show that the proposed approach numerically passes the patch test. The approach predicts an exactly uniform pressure distribution mainly because local force and moment equivalence is satisfied for the slave and master segments. As the proposed approach passes the patch test, the accuracy is expected to improve via mesh refinement for general applications. In addition, unexpected excessive pressure, which is a problem in conventional approaches, is reduced in the proposed approach. The NTS-AR approach, which is more accurate than the conventional NTS approach, as shown in Figs. 9 and 11, does not ensure better accuracy when it fails the patch test.

Further comparison with other refinements
The proposed approach is further compared with the virtual slave node-to-segment approach, i.e., VTS-PPT, developed by Zavarise and De Lorenzis [13]. Both approaches are derived from the penalty-method-based NTS approach; they employ the AR technique and pass the flat patch test. The virtual slave node technique does not project the existing slave node; instead, it utilizes two virtual slave nodes per slave segment. The virtual slave node is located in the middle of its competence area to resolve local moment nonequivalence. In addition, it projects the competence area on multiple master segments with a series of force transformations, as described in Sect. 1. This procedure aims to reduce the competence area mismatch. However, the virtual slave node technique is indispensable for modifying the original NTS formulation of contact force estimation because of the employment of additional nodes.
In contrast, the proposed approach employs an existing master node to resolve local nonequivalence. It adopts the Fig. 22 Insufficiency of the proposed method original NTS formulation, which also indicates enhanced compatibility with existing NTS-based contact analysis. In addition, the advantages of the two-pass approach, including the detection of master-to-slave penetration, which cannot be achieved by the one-pass approach, are obtained. Furthermore, the proposed approach is expected to prevent over-constraint or locking, which are the main disadvantages of the two-pass approach, by projecting the selected master nodes. These properties are investigated using the numerical examples shown in Sect. 5.4.
However, the limitation of the proposed approach may be attributed to the master-to-slave projection, as shown in Fig. 22. The two-end nodes of the slave segment are projected on different master segments. Node M 2 becomes an intermediate master node between slave node projections S * 1 and S * 2 . In this case, projection M * 2 from M 2 is not enclosed by target slave segment S 1 S 2 . Therefore, the slave segment area is not divided by the projection. For such a segment, the proposed approach does not utilize M 2 . Instead, it divides S 1 S 2 without the projection, similar to the conventional NTS-AR approach.
As the aforementioned situation frequently occurs for an almost node-matching configuration, such selective replacement of the conventional NTS-AR approach is reasonable. Nonetheless, the one-pass virtual slave-to-segment approach helps improve the conventional NTS and NTS-AR approaches without master-to-slave projection. Otherwise, the STS and mortar approaches can be used; however, their computational cost is high.

Additional tasks toward three-dimensional analysis
In this paper, the proposed approach is investigated in a twodimensional context. The discussions on the requirements for the three-dimensional extension are significant. The overall procedure is similar to that in Sect. 4.1. However, the strategy to determine the penalty weighting factor which satisfies the local force and moment equivalence in three dimensions must be developed. In two-dimensions, the contact surface consists of the line segments, wherein the equivalence is attained by dividing the line between projected point and node in half. In contrast, the three-dimensional contact deals with the additional geometry, surface element. Therefore, specific methodology that divides the surface element based on the projection information while satisfying the local force/moment equivalence should be developed. In addition, contact between the various types of surface elements, e.g., triangular and quadrilateral elements, should be considered. The three-dimensional extension of the proposed approach is currently being studied, and it will be presented in the future.

Additional patch test
The patch test Case E, which originates from the study of Taylor and Papadopoulos [11], is examined. As shown in Fig. 23, the smaller upper rectangular block presses the larger lower block. Except the contact interface, uniform pressure p 0 = 1 Pa acts on the upper surface of both solid bodies. Both blocks are discretized uniformly, and comprised of the same material properties: E = 10 6 Pa and ν = 0. The discretized configurations are symmetric to the vertical (X = 0) axis. At the interface between the two blocks, the upper and lower surfaces represent the slave and master surfaces, respectively. Figure 24 shows the corresponding numerical results. Because the contact surfaces are loaded by external pressure as well as the contact pressure, their resultant pressure is displayed instead of the contact pressure only. As shown, the proposed approach passes the patch test as well as VTS-PPT. In contrast, the conventional NTS, NTS-AR, and their two-pass approaches fail the patch test as reported by Refs. [11,13].
The modified Case E, referred to as Case F in Fig. 25, is finally investigated in terms of patch test. It was introduced by Zavarise and De Lonrenzis [13] for the validation of VTS-PPT approach. Case F is the same as Case E except that the lower block is discretized in a biased manner, which invokes the competence area mismatch. As in Case E, the discretized configurations are symmetric to the vertical (X = 0) axis. As shown in Fig. 26, the proposed and VTS-PPT approaches still pass the patch test despite the magnified discrepancy of the other approaches.

Internally pressurized cylinder
An example of an internally pressurized cylinder is employed to verify whether the numerical prediction converges to the analytical solution with mesh refinement. As illustrated in Fig. 27, the quarter portion of two hollow cylinders is mod-  eled with a symmetric boundary condition. The inner and outer radii of the inner cylinder are R 1 = 1 m and R 2 = 1.2 m, respectively. The inner radius of the outer cylinder is the same as R 2 , and the outer radius is R 3 = 1.4 m. The same material properties are applied to both cylinders, i.e., E = 1 Pa and ν = 0.3. The inner cylinder is internally pressurized by pressure p 0 = 0.005 N/m, and it presses the outer cylinder outward. At the interface between the two cylinders, the inner and outer surfaces are designated as the slave and master surfaces, respectively.
The mesh refinement level, N r , varies in the range 0-5. For each level, the inner cylinder is discretized by 10 · 2 N r + 1 × 2 · 2 N r quadrilateral elements. The outer cylinder consists of 10 · 2 N r − 1 × 2 · 2 N r elements to prevent the node-matching configuration, except at the nodes of the symmetric boundary. Figure 28 shows the mesh refinement from N r = 0 to 1. Each contact surface is equally discretized. The penalty parameter for the conventional NTS approach is N = N T S N = 1 N/m. To impose a similar penalty parameter, the distributed penalty parameter, N = N T S N /l sl , is utilized for the NTS-AR, VTS-PPT, and proposed approaches. Here, l sl is the slave segment area. First, the convergence property of the entire discretized field is investigated in terms of the error energy norm, ε E , defined in Eq. 9 [43].
denotes the volume of the entire solid body. σ ex and ε ex denote the stress and engineering strain vectors of the analytical solution, respectively, which are in the Voigt notation. σ h and ε h denote the vectors that are obtained through numerical prediction. The analytical solution is described in Ref. [44]. Figure 29 indicates a decrease in ε E according to the average mesh size, h. For a certain range of mesh refinement levels, linear convergence is observed for all approaches. From N r = 0 (h = 0.21 m) to 2 (h = 0.05 m), the VTS-PPT approach exhibits the lowest error energy norm. The proposed approach show better accuracy for the remaining regime, although the difference in ε E is small.
The improved accuracy of the proposed approach is demonstrated by localizing the aspect to the contact surface. Figures 30 and 31 show the contact pressure on the contact surface according to the mesh refinement level. Here, θ is the azimuthal angle in polar coordinates. When θ is 15 • -75 • , the conventional NTS and NTS-AR approaches predict contact pressures that are almost identical to each other, regardless of the discretization. In the vicinity of θ = 45 • , a nonuniform low pressure is observed for both approaches, which is not resolved even by the finest discretization. In the Dirichlet boundary condition regime (θ = 0 • and 90 • ), the NTS approach produces a distinctly excessive pressure, which increases with the mesh refinement level. Although the NTS-AR approach alleviates this overestimation, the nonuniformity still persists. The VTS-PPT approach resolves the inaccuracy for those regimes. However, another nonunifor-mity occurs in the adjacent regimes, which is retained for the finest refinement.
In contrast, the proposed approach predicts the contact pressure with an improved convergence compared to the analytical solution. It almost eliminates pressure nonuniformity for N r = 4. Subsequently, as the mesh is refined, the numerical prediction gradually approaches the analytical solution and yields excellent accuracy. The current example explains the meaning of the patch test results. Mesh refinement allows for the convergence of the contact pressure to the exact solution only when the contact element passes the patch test.
It should be noted that the error energy norm of the conventional approaches is low (Fig. 29) despite their inaccuracy in contact pressure prediction. This is because the relative proportion of the contact contribution to the entire field decreases with the mesh refinement level; Fig. 32 shows the highlighted scenario. Altogether, the stress distributions obtained using the conventional and proposed approaches are almost identical. However, in the localized aspect, the conventional NTS approach shows stress discontinuity across the contact surface. Although the discontinuity is alleviated by the NTS-AR approach, nonuniform stress still persists along the azimuthal direction. In contrast, the proposed approach produces a continuous and azimuthally uniform stress distribution. Also, the VTS-PPT predicts more precise results compared to the conventional NTS and NTS-AR approaches despite the slight nonuniformity.
The impact of the consistent tangent stiffness is investigated through the current example. The tangent stiffness of the conventional NTS approach is consistent and symmetric for the frictionless contact. However, the NTS-AR, VTS-PPT, and the proposed approaches require the additional stiffness terms to achieve consistency. Although such terms improve the convergence property, their asymmetry deteriorates the computational efficiency per each Newton-Raphson iteration. It is due to that the linear algebraic solver for a general asymmetric matrix should be introduced, which requires far more computational cost than that for a symmetric positive definite matrix. The additional stiffness terms are derived in "Appendices A.2 and A.3" for the NTS-AR and proposed approaches, respectively. Further, Ref. [13] explained the effects of the consistency for the VTS-PPT approach, of which inconsistent version exhibited the slow convergence even for the simple patch test. Figure 33 depicts the force convergence behavior in terms of the Newton-Raphson iteration numbers. It is shown that the consistent tangent stiffness results in a slight improvement in the convergence of the NTS-AR and proposed approaches, and that the inconsistent version shows sufficiently fast convergence behavior. Therefore, it would be sufficient to employ only the inconsistent stiffness terms for the proposed approach, which are symmetric and simple to construct.

Hertzian contact
The contact between two convex surfaces is examined using the Hertzian contact situation illustrated in Fig. 34. Two half-cylinders with the same radii (R = 1 m) and material properties (E = 1 Pa and ν = 0.2) are approximated using the symmetric boundary condition, and they collide with each other. The bottom of the lower cylinder is vertically constrained, and the top of the upper cylinder experiences distributed pressure p 0 = 0.0025 N/m. At the interface between the two cylinders, the upper and lower surfaces are designated as the slave and master surfaces, respectively. The original situation consists of a fully circular cylinder with a concen- Fig. 33 Internally pressurized cylinder: convergence behavior for the discretization with a N r = 2 and b N r = 5 trated force at the top [45]. Nevertheless, the approximation to the half configuration with equivalent pressure has been employed in numerous studies, and the contact approaches have been effectively demonstrated [15,22,28].
The analytical solution yields quantities that are related to the contact region. The contact area, l ex , maximum contact pressure, p ex max , and contact pressure distribution, p ex (X ), are expressed as follows [45]:  Figures 36 and 37 show the contact pressure distribution along the contact surface. Near the center, the conventional NTS and NTS-AR approaches predicted a highly fluctuating pressure distribution, and this becomes more apparent for the finer discretization. Although the fluctuation decreases as the distance from the center increases, it still results in a node-wise varying discrepancy compared to the analytical solution. As shown in Fig. 38, such local nonsmoothness eventually causes stress discontinuity across the contact surface whose magnitude alternates according to the position. In contrast, the VTS-PPT and the proposed approaches predict a smooth contact distribution on the slave and master surfaces, as shown in Figs. 36 and 37. Further, compared to the analytical solution, the discrepancy decreases as the mesh refinement level increases. This improvement is illustrated in Fig. 38, wherein a smooth stress distribution is preserved with continuity across the interface.

Largely bent beam
Because the proposed approach employs the two-pass concepts, the locking-free behavior is investigated. Figure 39 illustrates the example, which was presented by Puso and Laursen [22] and employed in numerous studies [15,28,46]. As shown, two slender solid bodies with the same size and material properties (E = 1 Pa and ν = 0) press each other through the ambient pressure of p 0 = 0.1 N/m. At the tip plane, the applied moment induces large bending. At the interface between the two bodies, the upper and lower surfaces are designated as the slave and master surfaces, respectively. Similar to the approach in the previous Hertzian contact example, two sets of mesh (N r = 0 and 1) are discretized for the investigation. As shown in Fig. 40, the lower body is equally discretized by 10 · 2 N r × 1 · 2 N r quadrilateral elements. The upper body is also uniformly discretized by 10 · 2 N r + 1 × 1 · 2 N r elements, except the biased discretization at each ends. For N r = 0, distributed penalty parameters ofˆ N = 4, 8, 16, 32, 64, and 128 N/m 2 are employed for the simulation. For the halved contact segments, N r = 1, twice larger parameters are used for the regularization. If the over-constraint problem exists, the ambient pressure would induce the excessive impenetration condition at the interface, thus resulting in locking for the applied moments. Figure 41 shows the vertical displacement, v, obtained through the numerical simulation. Figures 42 and 43 display the deformed configurations for N r = 0 and 1, respectively. The results are compared with that of the conformal discretization, which consists of 10 · 2 N r × 2 · 2 N r uniform elements. Except forˆ N = 64 and 128 N/m 2 , the proposed approach with N r = 0 simulates the large bending, i.e., locking-free behavior. However, the vertical displace-  locking. Therefore, it can be inferred that the proposed contact approach is locking free, given the moderate magnitude of distributed penalty parameter. Further, the mesh refinement for the proposed approach suppresses the impact of distributed penalty parameter.
It is worth noting that the analogous situation can be observed in the numerical study for the contact domain method [15]. Here, overestimation and locking occurred for the coarser discretization, and the behavior was highly dependent on the stabilization parameter. Then, such parameter dependency and locking were relieved by the mesh refinements. For the contact domain method, the stabilization parameter determined the amount of elemental elimination of discrete Largrange multipliers.

Half hollow cylinder with foundation
Another large deformation contact is considered, and a corresponding example was presented by Fischer and Wriggers [47]. As shown in Fig. 44, a half hollow cylinder is pressed by prescribed vertical displacement, v 0 = −40 mm, and it presses an almost rigid foundation. The cylinder consists of two cylinders of the equal thickness. The material properties are E = 100 and E = 1 GPa for the inner and outer cylinders with same Poisson ratio of ν = 0.3. The foundation consists of the material property of E = 100 TPa and ν = 0. The outer surface of cylinder and the upper surface of foundation are designated as the slave and master surfaces, respectively.
The discretization also follows Ref. [47], wherein the cylinder and foundation are uniformly discretized by 64 × 2 and 52 × 10 quadrilateral elements, respectively. The penalty parameter for the conventional NTS approach is N = N T S N = 100 kN/mm. Further, the distributed penalty parameter for AR-related approaches is employed asˆ N =   Figures 45 and 46 show the deformed configuration that is predicted by the numerical analysis. The mid-point of the cylinder initially contacts with the foundation. Then, the contact regime moves to both sides, and the mid-point is displaced vertically upward. All of the approaches depict such behavior without significant difference from each other. It can be concluded that the proposed approach is capable of predicting the large deformation similar to other approaches.
In addition, the final displacement of the cylinder midpoint is compared, as shown in Table 1. Fischers and Wriggers [47] showed that the penalty method predicted larger vertical displacement than the Lagrange multiplier method. They discussed that this was because the excessive force gradient was created by the penalty parameter, between In those aspects, despite the small underprediction from the other penaltymethod-based approaches, the proposed approach partially suppresses such gradient by distributing the contact contribution, from slave-to-master elements to the master-to-slave elements.

Conclusions
This paper presents an improved AR technique to overcome the limitations of the conventional penalty-method-based NTS and NTS-AR approaches. The inaccuracy of conventional approaches is primarily attributed to the estimation of the penalty parameter without considering the master surface geometry. The proposed method attempts to utilize the mas- ter surface geometry by adding master-to-slave projection.
The penalty weighting factor of the slave node is modified such that it is half of the area confined within the contiguous master-to-slave projection. The selected master nodes are employed in a complementary manner to achieve local force and moment equivalence. The improvement of the proposed NTS-AR approach is validated by the flat patch test, which includes a slave segment that spans multiple master segments.
The limitation of the proposed approach is that the penalty parameter estimation is an approximation for a curved contact surface. In addition, its applicability is restricted to linear-order finite elements. Nonetheless, the accuracy of the proposed approach is significantly higher than that of the conventional approaches in numerical examples with curved surfaces. In addition, the algorithm in the proposed approach maximally utilizes the environment of the original NTS formulation. Based on these outcomes, it is expected that the proposed approach will be compatible with the existing NTS framework.
The proposed approach is based on the two-dimensional static context with implicit analysis. However, the recent industrial applications of the penalty-method-based NTS approach focus on three-dimensional nonlinear explicit analysis. Thus, the research on the proposed approach is currently being extended, and an efficient master-to-slave projection strategy for three-dimensional situations is being investigated.

Declarations
Competing interests The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adap-tation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A: Linearization of contact virtual work
For the contact surface discretized by NTS elements, the linearization, δW c , of the contact virtual work, δW c , in Eq. 4 is approximated as follows:

A.1 Conventional NTS approach
For the conventional NTS approach, δW N ce is zero because the penalty parameter, N e , is constant. The remaining linearization contributions are expressed as follows [33,40,41]:

A.2 Area regularization (AR) technique
When the AR technique is used, δW N ce in Eq. 2 is no longer be eliminated because the penalty parameter varies according to the weighting factor. Therefore, its contribution should be added if the consistent tangent stiffness matrix is required. To derive the penalty stiffness increment, N e =ˆ N e w e , the penalty weighting factor, w e , and its increment, w e , are formulated as follows:  Figure 47 shows the slave node, S, with the neighboring segment, S i S. N sl se is the number of slave segments in the neighborhood of slave node S of the e-th contact element.
x sl ie is the position of the other end node, S i , for each neighboring slave segment, S i S. τ sl ie is the unit tangent vector from x sl ie to x sl e . Note thatK N ce,i is defined between δx T e and û ie , which results in an asymmetric tangent stiffness matrix.

A.3 Proposed modified AR technique
In the proposed improved NTS-AR approach, the penalty weighting factor,w e , and its increment, w e , differ from w e and w e in Eqs. A10 and A11. Figure 48 illustrates this scenario. For convenience, it shows the case in which the slave node, S, at position x sl e corresponds to ξ 1,sl = 1 in terms of the convective coordinate for the neighboring slave segment, S i S. Here, S i S is defined by positions x sl ie and x sl e . In addition, it is assumed that all neighboring slave segments are projected using master-to-slave projection. Then,w e and w e are formulated as follows: x sl e − x sl ie 2 ρ sl 1,ie · x ms ie − x sl ie (A16) Among the projections to slave segment S i S,ξ sl ie indicates the convective coordinate that is the closest to slave node S. In addition,x sl ie is its physical position.x ms ie is the position of the projecting master node,M i . ρ sl 1,ie is the covariant tangent vector at slave segment S i S. m 11,sl ie = 1/ ρ sl 1,ie · ρ sl 1,ie is the corresponding contravariant metric component. Hereafter, the double sign becomes positive and negative when slave node S corresponds to ξ 1,sl = −1 and ξ 1,sl = 1, respectively, in terms of the convective coordinate for S i S. The inverted double sign becomes negative and positive when slave node S corresponds to ξ 1,sl = 1 and ξ 1,sl = −1, respectively. Substituting N e =ˆ N e w e into Eq. A2, δW N ce is expressed as follows: ce,i defines the connection between δx T e and ũ ie . For the neighboring slave segment with no master-toslave projection,K N ,A ce,i andK N ,B ce,i are replaced withK N ce,i , similar to the conventional NTS-AR approach. Additionally, the role of the slave and master bodies is inverted for the master-to-slave contact elements.