A finite volume penalty‐based implicit procedure for the treatment of the frictionless contact boundaries

This article presents a new implicit coupling procedure for mechanical contact simulations using an implicit cell‐centred finite volume method. Both contact boundaries are treated as Neumann conditions, where the prescribed contact force is calculated using a penalty law, which is linearised and updated within the iterative solution procedure. Compared to the currently available explicit treatment, the implicit treatment offers better efficiency for the same accuracy. This is achieved with the proposed implicit linearisation, which replaces the explicit under‐relaxation of the contact force. The proposed procedure, intended for frictionless contact of Hookean solids, can handle non‐conformal contact interface discretisations and faces in partial contact. The accuracy and efficiency of the implicit approach are compared with the explicit procedure on four benchmark problems, where it is shown that the proposed method can significantly improve efficiency and robustness.

The application of the finite volume method to mechanical contact problems began with the seminal works of Tropša et al., 11 and Jasak and Weller, 12 in which contact constraints were enforced via the Dirichlet-Neumann coupling procedure. Although such coupling was efficient, due to flaws in its extension to the three-dimensional problems, Cardiff et al. 13 redirected further development to the penalty method and Neumann-Neumann boundary treatment. The penalty-based contact procedure, initially developed for small-deformation frictionless contact problems, was later extended to large-deformation frictional contact problems by Cardiff et al. 7 and to lubrication by Škuric et al. 14 In these works, the contact force is calculated pointwise, and surface-to-surface interpolation 15 is used to transmit contact stresses. To eliminate the need for surface-to-surface interpolation, Batistić et al. 16 proposed a penalty-based finite volume version of the segment-to-segment algorithm, where contact force is calculated by integrating the normal gap distribution. These contact algorithms treat contact boundaries explicitly, that is, the explicitly calculated contact force is subsequently under-relaxed and updated within the iterative segregated solution procedure. To improve efficiency, Scolaro et al. 17 proposed an implicit treatment developed using an analogy with the Tuković et al. 18 procedure for multi-material interfaces. Implicit treatment is achieved by introducing additional constraints for contact traction and the normal component of the displacement field at the contact interface. The calculation of the normal gap is performed for face centres, where they introduced a sigmoid blending function for treatment of the contact area edge. Extension to the non-conformal discretisation at the contact interface is achieved by employing Arbitrary Mesh Interface (AMI) interpolation. 19 The aforementioned works are related to the cell-centred grid arrangement and implicit solution algorithm. It should be mentioned that there are other finite volume approaches in which the development of contact modelling is lately gaining more attention. An example are works from Bessenov et al., 20 Rucnie et al., 21 and Haider et al. 22 in which the authors propose contact algorithms intended for explicit finite difference and vertex-centred finite volume frameworks.
This work presents an extension of work previously published in References 13 and 16 by introducing a new implicit penalty approach for the Neumann contact boundaries. Augmentation of the contact force using the under-relaxation technique is replaced with linearisation of the penalty force; consequently, it is expected to achieve the same level of accuracy as the explicit approach but with an improved convergence rate for a wide range of penalty factors. This is beneficial since achieving an acceptable convergence rate in explicit contact procedures can be troublesome due to the adjustment of both the penalty factor and the under-relaxation factor. The proposed procedure can handle contact cases with partial contact and non-conformal discretisation at the contact interface. It is considered an ideal candidate for further extension towards cases with material and geometric nonlinearities (large deformations) since the penalty method allows for higher robustness due to controlled violation of contact conditions. In other words, there is no kinematic constraint on contact boundaries, which is the main difference from the implicit procedure proposed by Scolaro et al. 17 Also, no additional user-defined parameter can affect efficiency and accuracy, compared to Reference 17 where the blending coefficient is introduced.
The remainder of this article is organised as follows: Section 2 summarises the mathematical description of the considered problem, followed by Section 3, in which the numerical method is presented. Next, a discussion on discretisation and solution procedure is given, and the focus is kept on the proposed approach for contact boundary treatment. In Section 4, the primary features of the proposed procedure are demonstrated using four numerical examples. Finally, concluding remarks are given in Section 5.

MATHEMATICAL MODEL
In the case of contact between two isothermal, linear elastic bodies without inertial and body forces, the conservation of linear momentum reads as follows: where is the Cauchy stress tensor, is the del operator, and n is outward pointing unit normal; Ω denotes the continuum body volume whereas Γ denotes the body surface. Due to the small strain assumption, all quantities are connected to the reference geometry. Using Hooke's law, the Cauchy stress tensor is expressed in terms of displacement vector u: where I denotes the second-order identity tensor; and are the first and second Lamé coefficient, respectively.
At the portion of the boundary where contact occurs, non-adhesive frictionless contact is mathematically described using the Hertz-Signorini-Moreau conditions: where subscript f denotes face centre value and n f = f ∕| f | is the unit face normal. The first term on the right-hand side of Equation (5) represents the orthogonal contribution, which is treated implicitly. The second term, which is added to account for face non-orthogonality, is calculated explicitly with the face displacement gradient linearly interpolated between adjacent cell-centre values: where superscript * denotes the most recent solution and f x is the interpolation factor calculated as the ratio of distances f x = fN∕PN. The explicit value of the cell-centre gradient in Equation (6) is calculated using a weighted least square interpolation. 9 The second right-hand side term in Equation (4), which represents the divergence of stress, as well as the subtracted Laplace term, are calculated explicitly using face displacement gradient from Equation (6). Their value, added to the source term of the resulting system of equations, is updated using fixed-point/Picard iterations until convergence is met. Due to the different computational stencils used for explicit and implicit calculation of the linear Laplace terms, the converged solution will contain numerical diffusion, which tends to smear out possible solution oscillations. More details regarding the discretisation procedure can be found in the literature; for example, in Jasak and Weller 1 or Demirdžić and Muzaferija. 9

Contact boundary treatment
The portion of the discretised boundaries at which contact interaction is expected to occur is treated using the penalty-based Neumann boundary condition. The prescribed value of the boundary gradient depends on the solution from the adjacent boundary in contact; this coupling of the contact boundaries can be numerically resolved in two ways: • Explicit: the contact traction and displacement gradient depends solely on the spatial configuration of the contact boundaries from the previous iteration of the solution procedure; • Implicit: the contact traction and displacement gradient are a function of the unknown displacement solution of the current iteration, that is, contact traction depends on the displacement of the boundaries in contact.
Both approaches should lead to the same solution; however, a superior convergence rate is expected from the implicit coupling. The convergence rate of the currently-available explicit approach can be seen as the main drawback when solving static force-loading contact problems, as noted by Scolaro et al. 17 The implicit contact procedure proposed here is expected to overcome this limitation.
In addition to the fact that the implicit approach should allow the use of non-conformal discretisation at the contact interface and be extendable to scenarios with material and geometric nonlinearities; the following facts must also be kept in mind: • The computational points are located at the cell-centres, meaning the boundary displacements (on which contact tractions depend) are unknown. • The solution procedure uses fixed-point iterations; thus, only the final converged solution satisfies the momentum balance for control volumes. In other words, any intermediate solution of the displacement field can be far from the final solution, especially at the beginning of the iterative procedure. • The contact area is a priori unknown; therefore, boundary faces can change their contact status during the iterative procedure. Also, contact is mathematically described with contact conditions which are non-linear functions.
In the subsequent analysis, standard notation from the computational contact mechanics literature is used, and contact boundaries are denoted as slave and master. The presented approach is symmetric in the sense that the same set of equations are derived for both the master and the slave side, consequently, to generalise the presentation, subscript i is used for the considered boundary and subscript k for the adjacent. Below, the equations are derived concerning subscript i, which can represent the master or the slave surface.

3.3.1
Penalty method and normal gap calculation By applying the penalty method in a linear form, the normal contact pressure is linearly dependent on the normal gap: p n (g n ) = n ⟨−g n ⟩ = { n g n g n ⩽ 0, contact 0 g n > 0, no contact (7) where ⟨ • ⟩ are Macaulay brackets, n is the user-defined penalty factor (penalty stiffness) and g n is the normal gap distance (penetration). The penalty factor is initially estimated using the relationship proposed by Hallquist et al.: 26 where K b is the boundary bulk modulus, Γ b is cell contact area and Ω b the volume of the boundary cell in contact. f scale is the scaling factor, typically set to 1, but can be corrected depending on the resulting penetration. The current approach calculates an average penalty factor for the entire contact interface using Equation (8).
At the contact boundary, the normal contact traction t n , is calculated by integrating the contact pressure distribution: By defining contact pressure using the penalty method, that is, Equation (7), normal contact traction is written as a function of normal gap distance: The integration of the normal gap distance is accomplished using the finite volume version of the segment-to-segment method 16 in which the face integration area Γ bi is decomposed into j segments. Each segment has a linear variation of the normal gap, allowing the exact integration of face pressure distribution. 27 The segments are represented as line or surface elements, depending on problem dimensions, see Figure 2.
In accordance with the segment-to-segment method, the integral in Equation (10) is decomposed as follows: where j is the segment area vector and g nj is the segment gap distribution. A projection rule is established to construct segments, and the master surface normal is adopted as the projection direction. Figure 3 summarises the segment F I G U R E 2 2D and 3D representation of segment domain (coloured in red). construction steps. Each pair of contact faces is projected onto the auxiliary plane p defined by the master face m.
Projected polygons s and m are subjected to a polygon clipping algorithm to find their intersection s ∩ m, that is, segment j. The distribution of the normal gap across the segment is continuous and linear; by using the mean value theorem, Equation (11) can be written as follows: where g nj is the average value of the segment normal gap, which is evaluated by decomposing the segment into a finite number of triangles q, as illustrated in Figure 4. For all segment vertices, including the centre point, the normal gap is evaluated in the direction of the master normal vector. These values are used to compute the average value by summing the contribution of each triangle: In the case of partial overlap between faces m and s, special care must be taken, as integration is performed only over the area with a non-negative gap value.
Since the sum of the segment areas corresponds to the overall face area ∑ | j | = | bi |, Equation (12) can be also written as: where g ni is the average value of face normal gap. By multiplying Equation (14) with the face area, it is clear that contact force in the segment-to-segment method is, in fact, determined using the occluded penetration volume. This allows symmetric contact treatment since occluded penetration volume is uniquely defined. 29 In other words, surface-to-surface interpolation is avoided, as values for the slave surface are simultaneously defined when integrating the master contact surface. Compared to the pointwise gap calculation, the integration of contact pressure allows an accurate estimate of contact pressure for faces in partial contact. Moreover, by defining the projection rule (ray-tracing projection), other problems related to pointwise gap calculations are circumvented 30 (undefined orthogonal projection and undetected penetration). The integration of the contact pressure in the segment-to-segment approach can be characterised as explicit because integration is performed on the last spatial configuration of the contact surfaces. More details regarding implementing the segment-to-segment approach can be found in Batistić et al. 16

Contact traction linearisation
In the case of explicit contact coupling, contact traction is explicitly calculated using the last spatial configuration of the contact surfaces. Contact traction (t n ) * bi , calculated using Equations (10)- (14), is augmented using the under-relaxation technique to promote convergence: where n is the under-relaxation factor and n is the accumulated normal contact traction. Instead of the explicit accumulation of contact traction, this work proposes contact traction linearisation. Linearisation is performed according to the Taylor expansion, concerning the normal boundary displacement of both contact surfaces: Taking into account that partial derivatives of normal contact traction are equal to the penalty factor n , Equation (16) is rewritten as follows: where the term is added to account for faces in partial contact and denotes the relative contact area. In the case of conformal contact, is zero or unity, whereas for non-conformal contact, it can take a value in the range [0, 1]. The projection tensor n * bi n * bi is added to take into account only the normal component of displacement and is calculated with the spatial configuration of contact surfaces from the previous iteration. The first term represents the explicit value of contact traction, calculated with the normal gap value obtained using the segment-to-segment approach, whereas the second term represents the implicit correction, which tends to zero at convergence. Note that Equation (17) is derived assuming that the boundaries will remain in contact. If not, numerical adhesion is obtained, which is subsequently removed in the next iteration.

3.3.3
Conformal discretisation at the contact interface Firstly, the equations for a conformal discretisation are derived, while the extension to non-conformal discretisation is discussed in Section 3.3.5. Figure 5 shows two contact surfaces with different material properties and with orthogonal cells adjacent to the boundary. A similar approach proposed by Tuković et al. 18 and Scolaro et al. 17 is employed, with the main difference being that the kinematic constraint equation is replaced by the proposed equation for contact traction linearisation. This is required because imposing a kinematic impenetrability constraint at the contact interface is contradictory to the penalty method. From the over-relaxed momentum balance, 1 represented by Equation (4), one can write the equation for the traction vector at the boundary: F I G U R E 5 Conformal discretisation at the contact interface, orthogonal cells.
where superscript * is used to denote the value calculated using the solution from the previous iteration of the solution procedure. By rearranging Equation (18), the implicit boundary normal gradient is expressed as follows: By following the discretisation for the internal faces, the normal gradient at the boundary can be also expressed using the boundary and the cell-centre displacement value: Using equation (20), the boundary normal gradient in equation (19) is replaced by the unknown displacement of the boundary and cell-centre: At the contact interface, the continuity of normal contact traction, that is, action-reaction principle, is written as: Combining Equation (22) with Equation (21) it is possible to derive an equation for boundary displacement in terms of unknown cell-centre displacements and boundary displacements of the adjacent boundary. Here, for the purpose of subsequent analysis, it is necessary to derive the equation for boundary displacement of boundary u bk : where e is a substitution for the explicit part of the equation: Next, using Equation (23), the linearised traction vector defined by Equation (17) is rewritten in a way that dependency on adjacent boundary displacement u bk is removed: By substituting the linearised traction vector from Equation (25) into Equation (21) it is possible to obtain an equation for boundary displacement solely in terms of unknown cell-centre displacements: where G bi and N bi are: Finally, the equation for the normal boundary gradient is obtained by substituting Equation (26) back into Equation (20): Implicit coupling between contact boundaries is achieved by the term (ii), representing the coupling term. The coupling term is added to the off-diagonal positions of the solution matrix. In contrast, the term (i) is added to the corresponding central coefficient placed at the solution matrix's diagonal. Note that the same form of the equation is prescribed on both boundaries, which results in the symmetry of the coupling term. Accordingly, the symmetric structure of the resulting solution matrix, caused by the discretisation of the displacement Laplacian, will be preserved. The remaining terms are updated during the solution procedure and included in the right-hand side solution vector of the resulting system of equations. As contact evolves during the solution procedure, for cells in contact, diagonal and off-diagonal coefficients are updated.
In the case of the segregated solution procedure, the contribution from terms (i) and (ii) cannot be treated entirely implicitly. Therefore, deferred correction is used, and the resulting tensor coefficients are decomposed into the diagonal and off-diagonal parts. The diagonal part is directly included in the diagonal coefficients of the left-hand side matrix. In contrast, the off-diagonal part is multiplied with the available solution and added to the right-hand side solution vector. Such manipulation is inevitable for the segregated solution procedure; when the off-diagonal (explicit) part dominates, it is expected to reduce the robustness of the approach; accordingly, to ensure robustness in such cases, the prescribed boundary gradient is under-relaxed: where is the under-relaxation factor, whose value is set to = 0.4 in all cases presented here. Through testing, it has been concluded that a lower value is unnecessary and that a value of 0.4 is optimal to ensure robustness in all cases. The impact of applying under-relaxation on efficiency is investigated and discussed in the next section. It should be noted that in the case of implementation within a block-coupled solution procedure (e.g., Cardiff et al. 31 and Das et al. 32 ), it is possible to treat terms (i) and (ii) entirely implicitly, meaning that under-relaxation of the boundary gradient would not be needed. The derived equation for the boundary normal gradient can also be applied to contact between a deformable body and a rigid surface, its form is obtained by zeroing terms that refer to the adjacent rigid boundary. After each iteration of the solution procedure, Equation (26) is used to update/match the boundary displacement to the correspondingly applied normal gradient. In the case of boundary non-orthogonality, the resulting boundary displacement is corrected using the correction term explained next.

Boundary non-orthogonality correction
Since finite volume discretisation facilitates the application of automatic mesh generation, in many cases, boundary non-orthogonality is inevitable. Figure 6 shows a conformal contact interface with non-orthogonal cells at both contact surfaces.
To preserve second-order accuracy at the boundary, the non-orthogonal correction is performed by adding the explicit correction term, u corr bi , to the Equation (26). The correction term is calculated using the available gradient at the cell-centre: where k bi = d bi − ni n bi . The same procedure is followed for the non-conformal interface discretisation as well as for the correction of the normal boundary gradient, used to calculate stress at the boundary explicitly.

Handling non-conformal discretisation at the contact interface
To extend the procedure to handle non-conformal contact interface discretisation, required cell-centred and face-centred values from the adjacent boundary are reconstructed by employing Generalised Grid Interface interpolation 15 (GGI). For each cell on the considered boundary i, a virtual cell on the adjacent (shadow) boundary k is constructed, as shown in Figure 7.
Cell-centred and face-centred values of the virtual cell, are interpolated using area-weighted interpolation: virt where k is the interpolated variable of the adjacent boundary and k→i is the weighting factor calculated as the ratio between intersection area and overall area: By constructing a virtual cell on the adjacent boundary and using its values, conformal discretisation is mimicked, and the same set of equations is still valid. . Here should be pointed out that interpolation accuracy is not of prime importance since interpolation is required for the implicit correction in Equation (17), which tends to zero for the converged solution.

F I G U R E 6
Non-orthogonal cells sharing the conformal contact interface.

F I G U R E 7
Non-conformal discretisation at the contact interface, corresponding virtual cell of cell P i .

Solution procedure
After the discretisation procedure, the resulting linear algebraic equation for each CV has the following form: where a P is the diagonal tensorial coefficient, b P is the source vector and a N the neighbour tensorial coefficient. For CVs sharing the contact boundary, an additional tensorial coefficient a Pc and a Nc are added to the CV central and neighbouring coefficients. Assembling Equation (33) for each control volume results in the system of equations: where [A] is sparse and weakly diagonal dominant square coefficient matrix, {u} represents the unknown solution vector and {b} is the source vector. It is possible to use either the segregated/decoupled approach or the block-coupled one to solve this system of equations. Since the adopted underlying discretisation requires fixed-point iterations to resolve explicitly discretised terms and inter-equation coupling, the segregated approach is a logical choice. However, in that case, tensorial coefficients are decoupled into the diagonal and off-diagonal parts. Only the diagonal part is inserted in the coefficient matrix, as explained in Section 3.3.3 for the contact boundaries.
The following steps summarise the solution procedure for one loading/displacement increment: (1) Calculate the cell-centre gradient u * using the least square method. 9 Interpolate the gradients from cell-centres to face centres using Equation (6). Without contact modelling (step 2), the above-summarized solution procedure is an extension of the solution procedure used by other authors, see References 1 and 3. The implementation of the contact detection procedure (step 2.2) is divided into two phases. In the first phase, contact faces are encapsulated in object-oriented bounding boxes, and over them, detection is performed to find potential contact pairs, that is, faces potentially in contact. In the second stage, each pair of faces is projected onto a plane, and the separating axis theorem filters out false pairs. Since contact detection is performed on each iteration, computational costs are kept low by conducting the first phase only in the vicinity of the previous solution; for further information, see Reference 16.

NUMERICAL EXAMPLES
The proposed procedure is implemented using the open-source library foam-extend, a community-driven fork of the open-source toolbox, OpenFOAM. To test the proposed procedure, four different examples are considered. Examples are solved using one displacement or load increment with a targeted relative residual of 1 × 10 −6 . To examine the accuracy, obtained solutions are compared with results from the explicit-coupling approach and available analytical solutions. A primary focus is investigating possible efficiency gain compared to explicit contact treatment. To ensure a fair comparison between explicit and implicit contact treatment, explicit contact treatment is carried out with the optimal adjustment of contact traction accumulation, that is, the value of the under-relaxation factor is iteratively found to achieve the best possible efficiency.

Patch test
To confirm that a uniform pressure field can be transmitted across a non-conformal contact interface, a patch test proposed by Crisfield 33 is considered. The lower block is fixed at the bottom surface, whereas the upper block has a prescribed displacement Δ = 0.01 m at its top surface (Figure 8). Both blocks have the same dimensions, 1 m width and 0.5 m height, and the same material properties, Young's modulus E = 1 × 10 6 Pa and zero Poisson's ratio = 0. A different spatial discretisation is used for both blocks: the lower block has 25 CVs, and the upper 64 CVs. Under plane strain conditions, the analytical solution to the problem is: 33 In case > 0, Equation (35) can be used if a slip condition is prescribed on the bottom surface of the lower block and not fixed completely, as here. Three values of the penalty factor with different orders of magnitude are considered. Figure 9 shows the relative error field for the highest penalty factor used. One can see that the error distribution is uniform, meaning that the contact stress is correctly transmitted. The relative error is shown in percentages and is defined as: Table 1 depicts the required number of outer iterations to reach convergence in implicit and explicit contact coupling cases. Their efficiency is comparable when the explicit contact is optimally adjusted. Otherwise, if explicit contact is not adjusted correctly, it cannot converge or requires significantly more iterations. This can be seen in Figure 10, where the dependency of the outer iteration number on the under-relaxation factor is shown.

Cylindrical punch
In this example, a filleted punch is pressed into a flat-topped cylindrical foundation. Problem geometry is similar to the one proposed by NAFEMS; 34 however, a much larger fillet radius is taken to reduce peak contact pressure at the transition edge between the punch flat and filleted part. Furthermore, both bodies have the same material properties (E = 100 GPa, = 0.3) to allow the use of the semi-analytical solution. 35,36 From the analytical solution, the required punch top surface displacement and force can be calculated for the targeted ratio b∕a = 0.7 where a is the contact half-width and b is the radius of the punch flat zone (10 mm). Accordingly, this example is solved as a force-loading and displacement-loading case. For the force-loading case, calculated force F = 44657 N is prescribed as uniform pressure on the punch top surface. For the displacement-loading case, the prescribed displacement is Δ = 0.32115 mm. The employed axisymmetric computational mesh and problem geometry are shown in Figure 11. The problem is solved with three penalty factors of different orders of magnitude. Figure 12 shows normalised contact pressure distribution at the foundation contact surface. As expected, implicit and explicit contact coupling gives the same pressure distribution that agrees well with the analytical solution. The difference between results for different penalty factors is most visible at the pressure peak (the point where the fillet begins), thus the higher value of the penalty factor can capture pressure distribution better.

F I G U R E 11
Cylindrical punch: problem geometry (dimensions in mm) and axisymmetric computational mesh (4050 CVs).

F I G U R E 12
Cylindrical punch: normalised contact pressure distribution at the foundation contact surface. Table 2 depicts the required number of iterations for the implicit and explicit contact procedures. Results show that optimally adjusted explicit contact can be efficient for low values of penalty factor; however, such values cannot produce a satisfactory level of accuracy, see Figure 12. On the other hand, implicit contact maintains approximately the same number of iterations and is more efficient than explicit contact for higher values of penalty factor. Moreover, the same convergence behaviour is obtained in the case of the force-loading scenario, meaning that the implicit coupling can overcome the significant drawback of the explicit counterpart. Nonetheless, a force-loading scenario requires imposing initial geometrical interference, that is, an initial overlap, to avoid unrestrained rigid body motion.

Cylinders with parallel axes
To test how the level of implicitness affects efficiency, a commonly used example of contact between two identical infinite half-cylinders is chosen, 37 see Figure 13A. Three penalty factors, orders of magnitude different, are used, and simulations were performed for different spatial configurations of the cylinders. Each spatial configuration is obtained by rotation around origin 0 by angle Θ. The cylinders' deformation is the same since they are identical. As a result, the contact area is flattened, and contact is predominantly happening in the direction of axis z ′ . The same solution is obtained for each rotated configuration regardless of the rotation. However, the level of implicitness of the contact coupling differs for the segregated solution procedure, and a degradation of efficiency is expected when axis z ′ does not coincide with axis z.  Figure 14 shows the required number of iterations regarding the spatial configuration of the contact interface. For angles, Θ = 0 • and Θ = 90 • , when the implicitness of the contact treatment is the highest, efficiency is the same irrespective of the choice of the penalty factor. This is not the case for rotated configurations, and efficiency degradation is higher for higher values of the penalty factor. Dashed lines represent the solution if Equation (29) is under-relaxed with the under-relaxation factor of 0.8 which is twice as much compared to the proposed value of 0.4. One can see that the usage of under-relaxation reduces efficiency only when the implicitness is lower. This reduction in efficiency is proportional to the under-relaxation factor used and if proper under-relaxation is not applied convergence cannot be reached for such cases. An example is the simulation with n Δ∕E = 3.12 which cannot converge for rotation angles 15 • < Θ < 75 • with higher value of under-relaxation factor. Nonetheless, even when under-relaxation is applied, the resulting efficiency is hard to reach using explicit contact coupling. It should be noted once again that only the segregated solution procedure suffers from this, and in the case of the block-coupled solution procedure terms (i) and (ii) in Equation (19) will be completely included in the solution matrix, thus implicitness can remain the same irrespective of spatial configuration.
To check also the accuracy, the obtained solution is compared with the Hertz analytical solution, defined as: 38 where subscripts 1 and 2 are used to denote the contact body, p 0 is peak pressure, a is contact half-with, R c is effective curvature radius and E c is contact modulus. The example is solved using displacement loading, thus force F in Equation (37) is taken from the solution obtained with the highest penalty factor. The distribution of stresses along the z-axis is defined as: 38 As shown in Figure 15A,B the distribution of contact pressure and interior stresses fits the analytical solution. There are no oscillations and solutions from explicit and implicit treatment match; even the lower cylinder, which has 4 CVs per contact half-width, can capture the distribution of stresses accurately.

Contact between a hemisphere and a deformable/rigid block
This final example considers the common Hertz example of contact between a hemisphere and a block. The problem geometry and computational mesh are shown in Figure 16A,B. To confirm that the implementation can handle 3D problems and different material properties, the sphere and block material properties are set to E s = 210 GPa, s = 0.3, and E b = 70 GPa, b = 0.33. Two scenarios are considered, one in which the block is deformable and the second in which the block is rigid. Force and displacement loads are also considered, and the hemisphere top surface has either prescribed displacement or pressure. For the displacement loading case, the prescribed displacement is Δ = 0.01 m, and for the force loading case a pressure of 24.883 MPa is prescribed. The tangential displacement of the hemisphere's top surface for force loading is set to zero. The relationship between displacement and force is calculated using the Hertz analytical solution: 39 where effective curvature radius R c is equal to the hemisphere radius and contact modulus E c is calculated using Equation (37). For the case n Δ∕E s = 10, Figure 17A-C is showing zz , von Mises stress Eq and hydrostatic pressure p normalised with peak contact pressure p 0 . As expected, the obtained stress fields are smooth and continuous across the contact interface. Additionally, for the case n Δ∕E s = 1, the normalised cell-centre distributions of zz along the z-axis is given in Figure 17B. One can see no differences between the results from the implicit and explicit contact approaches. Also, the As in the previous test cases, three values of penalty factors with different orders of magnitude are considered, and results are summarised in Table 3. Again, the results confirm that the implicit procedure maintains efficiency when the contact direction is dominantly aligned with a Cartesian axis and as the penalty factor is increased. In other words, greater accuracy is achievable for the same efficiency. However, this is not the case for explicit contact, which requires significantly more iterations for higher accuracy. When the average relative gap is below 1%, the implicit procedure shows a speed-up of at least ten times. Although speed-up should be measured using wall clock time, it is not necessary here since the cost of each outer iteration is approximately the same; thus, the ratio between overall iteration numbers can faithfully depict speed-up.

CONCLUSIONS
This article presents a new procedure for penalty-based implicit coupling of contact boundaries using the cell-centred finite volume method. Contact boundaries are treated as Neumann boundaries and contact force, calculated by explicit integration of contact pressure distribution, is linearised and updated within a segregated solution procedure. As a result, compared to explicit coupling, the same accuracy can be obtained but at higher efficiency, which is the primary motivation behind this work. The proposed procedure is tested on four examples, and the following conclusions can be derived: • Overall solution efficiency directly depends on the contact coupling type. The results show that using the implicit contact coupling significantly improves efficiency, particularly when compared to a poorly adjusted explicit coupling. In the explicit coupling method, the user must adjust both the under-relaxation factor and the penalty factor value, affecting accuracy and efficiency; in contrast, the implicit coupling procedure only requires the specification of an under-relaxation factor, where a value of 0.4 is robust and efficient for all cases presented. • With implicitly coupled contact boundaries, it is possible to solve static force-loading contact problems. This is a troublesome task for explicitly coupled boundaries due to the imbalance of prescribed load and contact load during convergence and requires significant under-relaxation. • Since the implementation is done within a segregated solution procedure, the level of implicitness (efficiency) directly depends on the spatial orientation of the contact surfaces. In the case of high implicitness, the procedure maintains the same efficiency for different values of penalty factors. • In case of lower implicitness of contact treatment, the prescribed boundary gradient requires moderate under-relaxation to ensure robustness. Under-relaxation would be avoided in the case of implementation within a block-coupled solution procedure.
Future work will focus on extending the implicit procedure to contact problems with material and geometric nonlinearities and the implicit inclusion of the tangential (frictional) component of contact traction.