1 Introduction↩︎

The finite element method is a central computational tool of virtual prototyping and simulation-driven engineering design. It enables simulation-driven design before physical prototypes are built and motivates robust numerical methods for strongly nonlinear boundary value problems. Contact mechanics is a prominent example. Its difficulty arises from evolving contact zones, changing boundary conditions, large local deformations and localized force transfer, cf. [1]–[3]. These effects occur in metal forming, crash simulations and sealing, but also in soft robotics and contact-aided mechanisms, where contact becomes part of the desired response, cf. [4], [5]. Classical finite element contact formulations discretize potential contact interfaces explicitly and enforce the constraints by penalty, Lagrange multiplier, augmented Lagrangian or barrier-type methods, cf. [1], [2]. Surface-to-surface and mortar formulations extend this framework to non-matching contact meshes, cf. [6], [7]. Despite their maturity, these methods rely on contact surfaces, gap functions, active sets and contact search algorithms. For large sliding, self-contact, severe mesh distortion or topology optimization, where contact boundaries may evolve or may not be known in advance, this becomes a major algorithmic burden, cf. [5], [8], [9]. This motivates formulations in which contact emerges from a continuum model instead of being imposed on an explicit lower-dimensional interface. The third medium contact method follows this idea. Instead of enforcing contact constraints directly between two surfaces, the space between potentially contacting bodies is filled with a highly compliant fictitious medium. This third medium has only a negligible influence before contact, but stiffens strongly when compressed to nearly vanishing volume. Contact forces are then transmitted through the compressed medium, while explicit contact search and inequality constraints are avoided. The method is introduced in [10] for frictionless finite deformation contact using an isotropic-anisotropic material model for the intermediate medium. A related fictitious contact material combined with high-order finite elements is investigated in [11]. The extension to isogeometric analysis, including a study of material parameters and higher-order spatial convergence, is presented in [12]. An isogeometric-meshfree coupling strategy is proposed in [13], showing that the concept can also be transferred beyond standard finite element discretizations. A renewed interest in third medium contact emerges in density-based topology optimization. In this setting, the void phase naturally provides a region in which a fictitious contact medium can be placed, while contact boundaries are not known in advance. Internal contact modeling for finite strain topology optimization is introduced in [8]. This work also proposes the so-called HuHu regularization, which penalizes the Hessian of the displacement field in the third medium to control severe element distortions. The use of third medium contact for topology optimization of self-contacting structures is further developed in [5], where additional design requirements are introduced to improve the robustness of optimized contact-aided mechanisms. Internal contact is also used as a design mechanism for nonlinear elastic springs with tailored force-displacement responses in [14]. The concept is extended to contact-aided thermo-mechanical regulators in [15], where self-contact is used to tune heat transfer through switches, diodes and triodes. Frictional third medium contact is introduced in [16] by adding an anisotropic shear contribution and a crystal-plasticity-inspired slip mechanism that mimics Coulomb friction.

Despite these advances, the stabilization of the third medium remains a central issue. The medium must be soft enough to avoid forces before contact, but stable enough to withstand extreme compression, shear and mesh distortion during contact. The HuHu regularization is highly effective, but it penalizes all second-order deformation modes, including bending. To reduce this effect, the HuHu-LuLu regularization is introduced in [17]. There, a Laplacian contribution is subtracted from the Hessian contraction, which reduces the penalization of bending and quadratic compression while maintaining control over undesirable skew deformation modes. A deformation-gradient averaging regularization is proposed in [18]. This approach penalizes deviations of the deformation gradient at the integration points from an element-wise representative value and thereby controls spatial variations of the deformation gradient without additional degrees of freedom. A rotation-based regularization is introduced in [19]. It targets local changes of rotation, avoids penalizing stretch deformation modes and remains compatible with first-order elements without adding auxiliary fields. In [20] a thin layer is introduced between the solid and the third medium. The idea is to make the third medium even softer with a stiffening in the thin layer. By this the influence of the third medium on the deformation of the solids is reduced. Another active direction is the development of low-order-compatible formulations. Many gradient-based regularizations require second derivatives of the displacement field and therefore higher-order finite elements. This increases implementation effort and computational cost. First-order finite element formulations are introduced in [9], where auxiliary scalar fields approximate gradients of rotation-like or skew-symmetric deformation measures. This enables triangular, quadrilateral, tetrahedral and hexahedral low-order elements, but introduces additional unknown fields. A related thermo-mechanical formulation based on low-order ansatz spaces is proposed in [21], where the third medium carries both mechanical contact forces and heat flux without explicit interface conditions. A stabilization-free virtual element formulation for third medium contact is presented in [22], allowing polygonal discretizations while avoiding the classical virtual element stabilization difficulty in the presence of higher-order regularization terms. A three-dimensional third medium contact framework for hyperelastic contact and pneumatically actuated systems is developed in [23]. This formulation includes a three-dimensional regularization strategy and a pneumatic loading contribution. An efficient solution strategy for additional regularization fields is proposed in [24]. There, the fields are solved in a staggered manner at quadrature points by means of a neighbored-element method, which reduces the number of global unknowns while retaining low-order displacement approximations. The current state of third medium contact therefore provides several robust stabilization strategies. However, each strategy comes with specific numerical consequences. Some approaches require higher-order derivatives or higher-order finite elements. Others introduce additional global fields, rely on element-wise averaging, or use specialized discretizations. Low-order finite elements remain attractive since they are simple, robust and computationally cheap. At the same time, they require stabilization mechanisms that do not rely on the direct evaluation of second derivatives of the displacement field.

In the present work, an auxiliary-field stabilization for third medium contact formulations is introduced. A deformation-gradient-like field is added in the third medium and treated independently from the displacement field. A penalty contribution couples this auxiliary field to the physical deformation gradient, while the stabilizing term acts on the gradient of the auxiliary field rather than on the gradient of the deformation gradient itself. As a result, the formulation avoids the direct evaluation of second displacement derivatives and remains compatible with first- and second-order finite elements. Continuous and element-wise discontinuous auxiliary-field interpolations are considered in order to assess the influence of inter-element continuity on the regularization effect. Several benchmark problems involving large deformation, severe third-medium compression and progressive self-contact are used to evaluate the robustness and accuracy of the proposed formulation.

The remainder of this paper is organized as follows. First, existing third-medium contact stabilization strategies are reviewed and the present contribution is placed in context. The finite-deformation setting, the third-medium energy and the auxiliary-field stabilization are then introduced, followed by the finite element discretization for continuous and discontinuous auxiliary fields. A one-dimensional low-order example isolates the central mechanism of the method. Even though the deformation gradient is constant within each first-order element, a continuous auxiliary field transfers element-wise jumps into non-vanishing auxiliary-field gradients. This explains the effectiveness of gradient-type stabilization for low-order displacement approximations. The formulation is subsequently tested in two-block compression, self-contact within a box, a C-shaped self-contact problem and a three-dimensional extension. The paper closes with the main findings, limitations and possible extensions.

2 Overview of existing third-medium stabilization strategies↩︎

Several third-medium contact regularization strategies have been proposed to control excessive mesh distortion during severe compression and sliding. Tab. 1 summarizes the main developments relevant to the present work. Existing approaches either penalize higher-order displacement information directly, reduce the regularization to selected rotation-, skew- or volume-related measures, exploit special discretization technologies, or avoid additional unknowns by element-wise averaging. The formulation proposed here follows a different route. A deformation-gradient-like auxiliary field is introduced as an independently interpolated mixed variable. Its gradient provides the stabilizing contribution, while the penalty coupling to the physical deformation gradient controls the consistency of the auxiliary field. For discontinuous auxiliary interpolations, the additional degrees of freedom remain element-local and can be eliminated by static condensation.

4pt

3 Continuum mechanical foundations↩︎

Within this section, the continuum formulation of the proposed third medium contact approach is introduced. The physical solid domain is denoted by \(\B_{\mathrm{s}}\), while the space between potentially contacting surfaces is represented by the fictitious third medium domain \(\B_{\mathrm{tm}}\). Both domains are described within finite deformation kinematics. For the solid, a standard hyperelastic material model is employed. In the third medium, the same hyperelastic base energy is used in a strongly scaled form and supplemented by an additional stabilization contribution.

3.1 Boundary value problem, kinematics and weak form↩︎

All balance equations are formulated with respect to the reference configuration \(\B \subset \mathbb{R}^3\). Under static loading, the balance of linear momentum reads \[\Div{\bP} + \rho_0 \bb = \bzero \qquad \mathrm{in}\quad \B , \label{eq:BaMo}\tag{1}\] where \(\bP\) is the first Piola-Kirchhoff stress tensor, \(\rho_0\) is the mass density in the reference configuration and \(\bb\) denotes the body force per unit mass. For the boundary, the standard decomposition \[\partial\B = \partial\B_{u} \cup \partial\B_{t}, \qquad \partial\B_{u} \cap \partial\B_{t} = \emptyset\] is used, together with \[\bu=\bar{\bu} \quad\mathrm{on}\quad \partial\B_{u} \qquad\textrm{and}\qquad \bP\cdot\bN=\bar{\bt} \quad\mathrm{on}\quad \partial\B_{t}. \label{eq:BC}\tag{2}\] Here, \(\bar{\bu}\) is the prescribed displacement, \(\bar{\bt}\) is the prescribed traction and \(\bN\) denotes the outward unit normal in the reference configuration. A material point \(\bX\in\B\) is mapped to the current configuration by \[\bx = \Bvarphi(\bX,t) = \bX+\bu(\bX,t).\] Consequently, the deformation gradient, its Jacobian and the right Cauchy-Green tensor are given by \[\bF = \frac{\partial \bx}{\partial \bX} = \bI+\nabla\bu , \qquad J = \det\bF \qquad\textrm{and}\qquad \bC = \bF^{T}\cdot\bF . \label{eq:F}\tag{3}\] \(J\) measures the local volume change. The volume-preserving part of the deformation is described by the isochoric deformation gradient \[\widehat{\bF} = J^{-1/3}\bF .\] In the solid domain, the strain-energy density is denoted by \(\psi\). Here, a compressible Neo-Hookean model is used as \[\psi(\bC) = \frac{K}{2} \left(\ln J\right)^2 + \frac{\mu}{2} \left( J^{-2/3}\tr\bC - 3 \right), \label{eq:NeoHook}\tag{4}\] where \(K\) and \(\mu\) are the bulk and shear moduli, respectively. Stress measures follow from \[\bS = 2 \frac{\partial\psi(\bC)}{\partial\bC}, \qquad \bP = \bF\cdot\bS . \label{eq:PK}\tag{5}\] Multiplication of Eq. 1 with an admissible virtual displacement \(\delta\bu\) and integration over the reference domain yields the weak form \[\begin{align} \delta\Pi = \int_{\B} \frac{\partial\psi(\bC)}{\partial\bF} : \delta\bF \,\mathrm{d}V - \int_{\B} \rho_0 \bb \cdot \delta\bu \,\mathrm{d}V - \int_{\partial\B_t} \bar{\bt} \cdot \delta\bu \,\mathrm{d}A = 0 . \end{align} \label{eq:dynamic95weak95form}\tag{6}\] Here, \(\delta\bF=\nabla\delta\bu\) denotes the variation of the deformation gradient. Body forces \(\bb\) are neglected within the following.

3.2 Third medium energy↩︎

In the third medium, the corresponding strain-energy density \(\psi\) is used as a fictitious base material that is scaled by the small parameter \(\gamma\), while \(W^{\mathrm{tm}}\) denotes the total third-medium energy obtained by integration over \(\B_{\mathrm{tm}}\). This scaling keeps the influence of the third medium negligible before contact, whereas the underlying energy still increases strongly under severe compression and thereby provides the desired contact-barrier effect. For three-dimensional simulations, the third-medium energy is defined as \[W^{\mathrm{tm}} = \gamma \int_{\B_{\mathrm{tm}}} \left[ \frac{K}{2} \left(\ln J\right)^2 + \frac{\mu}{2} \left( J^{-2/3} \tr\bC - 3 \right) \right] \,\mathrm{d}V . \label{eq:NeoHookTM}\tag{7}\] Volumetric changes are controlled by the first term, while the second term describes the isochoric deformation of the third medium. For the three-dimensional computations, the full third-medium energy in Eq. 7 is retained. In the two-dimensional plane strain examples, however, a reduced contact-barrier energy is used. This choice follows the common interpretation of the third medium as a fictitious contact medium rather than as a physical material phase: its stiffness should remain negligible before contact, while the energy has to increase strongly when the intermediate region is compressed towards a nearly collapsed state, cf. [9], [10]. Within the present plane strain setting, the deformation is embedded in three dimensions by setting \(F_{33}=1\). The isochoric contribution then already provides the required barrier effect for strongly collapsing in-plane deformations, since it becomes singular as \(J \to 0\). Consequently, the volumetric contribution is omitted in the reduced two-dimensional model, leading to \[W^{\mathrm{tm}}_{\mathrm{red}} = \gamma \int_{\mathcal{B}_{\mathrm{tm}}} \frac{\mu}{2} \left( J^{-2/3}\operatorname{tr}\mathbf{C}-3 \right) \,\mathrm{d}V . \label{eq:Wtm95red}\tag{8}\] Similar reduced two-dimensional third-medium energies based on distortional or isochoric contributions have been used in [24]. It should be emphasized that Eq. 8 is not intended as a complete three-dimensional material law. It is used only as a reduced two-dimensional contact-barrier energy in the plane strain examples, whereas all three-dimensional simulations employ the full energy in Eq. 7 .

3.3 Deformation-gradient-based stabilization↩︎

A useful way to motivate the proposed stabilization is to consider a gradient-type control of the deformation gradient in the third medium. Such a term would penalize spatial variations of \(\bF\) and can formally be written as

\[W^{\mathrm{tm}}_{\nabla\bF} = \int_{\B_{\mathrm{tm}}} \frac{\alpha_r}{2} \left| \nabla \bF \right|^2 \,\mathrm{d}V , \label{eq:TMC95direct95gradF}\tag{9}\] where \(\alpha_r\) is a regularization parameter. Here, Eq. 9 is used only as a motivating reference form and is not meant to represent the specific structure of all existing third-medium regularizations. Instead, it highlights the quantity that the present formulation aims to control without evaluating second derivatives of the displacement field directly. For low-order finite elements, however, such a direct term is not suitable. With a first-order displacement interpolation, the deformation gradient is element-wise constant and its gradient cannot provide an effective intra-element stabilization. To avoid the direct evaluation of \(\nabla \mathbf{F}\), an additional deformation-gradient-like field \(\BTheta\) is introduced in the third medium. This field is interpolated independently from the displacement field and is weakly coupled to the physical deformation gradient by the penalty contribution as \[W^{\mathrm{tm}}_{p} = \int_{\B_{\mathrm{tm}}} \frac{p_\Theta}{2} \left| \BTheta-\bF \right|^2 \,\mathrm{d}V , \label{eq:TMC95penalty95theta95F}\tag{10}\] where \(p_\Theta\) is the penalty parameter. For sufficiently large \(p_\Theta\), the auxiliary field is driven towards the deformation gradient. Regularization is then applied to the spatial gradient of the auxiliary field, \[W^{\mathrm{tm}}_{r} = \int_{\B_{\mathrm{tm}}} \frac{\alpha_r}{2} \left| \nabla\BTheta \right|^2 \,\mathrm{d}V . \label{eq:TMC95grad95theta}\tag{11}\] In this way, a gradient regularization of \(\bF\) is approximated without computing \(\nabla\bF\) directly. Only first derivatives of the independently interpolated field \(\BTheta\) are required. For three-dimensional simulations, the stabilized third-medium energy becomes \[\begin{align} W^{\mathrm{tm}}_{\mathrm{stab}} &= W^{\mathrm{tm}} + W^{\mathrm{tm}}_{p} + W^{\mathrm{tm}}_{r} \\[2mm] &= \int_{\B_{\mathrm{tm}}} \bigg[ \gamma \left[ \frac{K}{2} \left(\ln J\right)^2 + \frac{\mu}{2} \left( J^{-2/3} \tr\bC - 3 \right) \right] + \frac{p_\Theta}{2} \norm{\BTheta-\bF}^2 + \frac{\alpha_r}{2} \norm{\nabla\BTheta}^2 \bigg] \,\mathrm{d}V . \end{align} \label{eq:TMC95total95energy95theta}\tag{12}\] For the reduced two-dimensional setting, \(W^{\mathrm{tm}}\) in Eq. 12 is replaced by \(W^{\mathrm{tm}}_{\mathrm{red}}\) from Eq. 8 . Hence, the proposed stabilization provides a low-order-compatible approximation of a deformation-gradient regularization while avoiding the direct computation of second derivatives of the displacement field.

3.4 Finite element discretization and static condensation↩︎

The proposed stabilization leads to a mixed finite element formulation in the third medium. The displacement field \(\bu\) and the deformation-gradient-like auxiliary field \(\BTheta\) are treated as independent unknowns. The physical deformation gradient \(\bF\) remains kinematically tied to the displacement field, while \(\BTheta\) serves as an additional stabilization field. Both fields are coupled weakly through the penalty contribution in Eq. 10 . The regularization acts on \(\nabla\BTheta\) and therefore avoids the direct use of \(\nabla\bF\). Consequently, the formulation mimics a deformation-gradient regularization without requiring second derivatives of the displacement field. Within an isoparametric finite element setting, the displacement field and the auxiliary field are interpolated independently. On element level, the approximations read \[\underline{\bu} = \underline{\IN}^{e} \underline{\bd}_{u}^{e}, \qquad \underline{\BTheta} = \underline{\IN}_{\Theta}^{e} \underline{\bd}_{\Theta}^{e} \label{eq:FE95interpolation}\tag{13}\] while associated gradient quantities are computed as \[\underline{\bF} = \underline{\bI} + \underline{\IB}_{u}^{e} \underline{\bd}_{u}^{e}, \qquad \nabla\underline{\BTheta} = \underline{\IB}_{\Theta}^{e} \underline{\bd}_{\Theta}^{e}. \label{eq:FE95gradients}\tag{14}\] A central feature of the formulation is that the interpolation order of \(\BTheta\) is not tied to the interpolation order of \(\bu\). The auxiliary field may therefore be approximated with the same polynomial degree as the displacement field, but lower-order approximations are equally possible. Since \(\BTheta\) serves as a numerical stabilization field rather than as an independent physical field, no physical continuity requirement is imposed a priori. This allows both continuous and element-wise discontinuous auxiliary interpolations to be considered within the same framework. In the present work, the interpolation pairs \(\mathrm{T}_{1}^{u}\mathrm{T}_{1}^{\Theta}\), \(\mathrm{T}_{2}^{u}\mathrm{T}_{2}^{\Theta}\) and \(\mathrm{T}_{2}^{u}\mathrm{T}_{1}^{\Theta}\) are compared with their discontinuous counterparts \(\mathrm{T}_{1}^{u}\mathrm{T}_{1}^{\Theta,d}\), \(\mathrm{T}_{1}^{u}\mathrm{T}_{0}^{\Theta,d}\), and \(\mathrm{T}_{2}^{u}\mathrm{T}_{1}^{\Theta,d}\). Here, the subscript denotes the polynomial degree, while the superscript identifies the approximated field. No additional superscript indicates a continuous auxiliary interpolation, whereas the superscript \(d\) denotes an element-wise discontinuous approximation of \(\BTheta\). For discontinuous auxiliary interpolations, all \(\BTheta\)-degrees of freedom remain local to the element. Static condensation can therefore be performed on element level, such that the global Newton system contains only displacement degrees of freedom. For continuous auxiliary interpolations, the corresponding degrees of freedom are globally coupled and are therefore assembled in the standard way. After linearization, the element contribution takes the block form \[\begin{bmatrix} \underline{\IK}^{e}_{uu} & \underline{\IK}^{e}_{u\Theta} \\[2mm] \underline{\IK}^{e}_{\Theta u} & \underline{\IK}^{e}_{\Theta\Theta} \end{bmatrix} \begin{bmatrix} \Delta\underline{\bd}^{e}_{u} \\[2mm] \Delta\underline{\bd}^{e}_{\Theta} \end{bmatrix} = - \begin{bmatrix} \underline{\IR}^{e}_{u} \\[2mm] \underline{\IR}^{e}_{\Theta} \end{bmatrix}. \label{eq:block95system}\tag{15}\]

For continuous auxiliary-field interpolations, both fields contribute to the global system, and Eq. 15 represents the element-level form of the coupled mixed problem. For discontinuous auxiliary-field interpolations, all auxiliary degrees of freedom are element-local. In this case, the mixed character is retained at element level, but the auxiliary unknowns can be eliminated before global assembly. For element-wise discontinuous auxiliary fields, elimination of the auxiliary increment yields the condensed displacement equation \[\underline{\IK}^{e,\mathrm{c}}_{uu} \, \Delta\underline{\bd}^{e}_{u} = - \underline{\IR}^{e,\mathrm{c}}_{u}, \label{eq:condensed95system}\tag{16}\] with \[\underline{\IK}^{e,\mathrm{c}}_{uu} = \underline{\IK}^{e}_{uu} - \underline{\IK}^{e}_{u\Theta} \left( \underline{\IK}^{e}_{\Theta\Theta} \right)^{-1} \underline{\IK}^{e}_{\Theta u} \quad\textrm{and}\quad \underline{\IR}^{e,\mathrm{c}}_{u} = \underline{\IR}^{e}_{u} - \underline{\IK}^{e}_{u\Theta} \left( \underline{\IK}^{e}_{\Theta\Theta} \right)^{-1} \underline{\IR}^{e}_{\Theta}. \label{eq:schur95complement}\tag{17}\] Once the global displacement increment is known, the auxiliary increment is recovered locally by \[\Delta\underline{\bd}^{e}_{\Theta} = - \left( \underline{\IK}^{e}_{\Theta\Theta} \right)^{-1} \left( \underline{\IR}^{e}_{\Theta} + \underline{\IK}^{e}_{\Theta u} \Delta\underline{\bd}^{e}_{u} \right) \label{eq:theta95increment95recovery}\tag{18}\] and the auxiliary degrees of freedom are updated according to \[\underline{\bd}^{e,it+1}_{\Theta} = \underline{\bd}^{e,it}_{\Theta} + \Delta\underline{\bd}^{e,it}_{\Theta}. \label{eq:DOF95update}\tag{19}\] Thus, the discontinuous variants retain a local auxiliary-field coupling without increasing the size of the global system, while the continuous variants provide the stronger inter-element regularization mechanism investigated below.

4 Numerical examples↩︎

Within this section the performance of the proposed auxiliary-field stabilization is assessed. The numerical study is organized around benchmark problems that test different aspects of the formulation, including the influence of the third-medium scaling, the penalty coupling, the gradient regularization and the interpolation order. First, two elastic blocks are pressed into contact to analyze the role of the third-medium parameters and the residual gap. Further verification is provided by the self-contact-within-a-box benchmark and the classical C-shaped boundary value problem, for which reference results are available in the literature. Most examples are formulated in two dimensions in order to isolate the relevant stabilization effects and allow a clear comparison of the different interpolation pairs. A three-dimensional version of the self-contact-within-a-box problem is included to demonstrate that the formulation is not restricted to plane settings.

4.1 \(\nabla\BTheta\)-regularization for continuous low-order interpolation↩︎

This example illustrates why a low-order discretization can still generate a non-vanishing gradient of the auxiliary field \(\BTheta\), and why this mechanism requires a continuous interpolation of \(\BTheta\). For clarity, the setting is reduced to a one-dimensional boundary value problem discretized by two linear finite bar elements. The total length is \(L=2\), with element lengths \(L_1=0.25\) and \(L_2=1.75\). Three geometrical nodes are obtained, and the displacement field is fully prescribed by \(u^1=0\), \(u^2=1/8\) and \(u^3=7/4\). Hence, the deformation field is fixed and only the auxiliary degrees of freedom associated with \(\BTheta\) remain unknown. Depending on the interpolation chosen for \(\BTheta\), different field distributions arise inside the elements and across the shared element boundary. Fig. 2 shows the prescribed displacement field and the resulting auxiliary-field distributions. Due to the linear displacement interpolation, the deformation gradient is constant within each element, but differs between the short left element and the longer right element. With a continuous interpolation of \(\BTheta\), the auxiliary field shares one degree of freedom at the common node. This shared value has to represent the deformation-gradient information of both adjacent elements and therefore acts as an inter-element transition value. As a result, \(\BTheta\) becomes piecewise linear and develops a non-zero gradient across the two elements, as shown in Fig. 2b). For an element-wise discontinuous interpolation, each element owns its local auxiliary degrees of freedom. The auxiliary field can then represent the constant deformation gradient of each element independently, without enforcing compatibility at the shared geometrical node. Consequently, the field becomes discontinuous at the element interface, as shown in Fig. 2c), and no continuity-driven gradient-based coupling between neighboring elements is generated. Thus, the continuous interpolation introduces a non-local stabilization mechanism: changes of the deformation gradient between adjacent elements are transferred into gradients of \(\BTheta\) and can be penalized by the regularization term in Eq. 11 . This observation is relevant not only for the present formulation, but also for related third-medium contact stabilizations that rely on auxiliary fields or gradient-type regularization measures.

4.2 Self-contact between two blocks↩︎

In this example the influence of the third-medium scaling parameter \(\gamma\) on the contact response is investigated. Two elastic blocks of size \(100\times50\) are separated by a third medium layer of the same in-plane size, as shown in Fig. [fig:BVPcontactBlocksRef]. Dirichlet boundary conditions are prescribed on the top and bottom boundaries. The lower boundary is fixed by \(\overline{\bu}=\bzero\), while the upper boundary is displaced according to two loading cases. A T\(_1^u\)T\(_1^\Theta\) discretization is applied.

In case 1, a purely vertical compression is applied with \(\overline{\bu}=[0,-60]^T\). In case 2, a combined horizontal and vertical compression is applied with \(\overline{\bu}=[10,-60]^T\). The solid blocks are modeled by the compressible Neo-Hookean material in Eq. 4 with \(K_{\mathrm{s}}=5/3\) and \(\mu_{\mathrm{s}}=5/14\). Inside the third medium, the same base energy is used with \(\gamma=\{10^{-3},10^{-4},10^{-5},10^{-8}\}\) and the stabilization parameters \(p_{\Theta}=1\) and \(\alpha_r=1\). This setup is used to assess how the stiffness scaling of the third medium affects the residual gap between both blocks. For the evaluation, two reference nodes are tracked at the right edge of the third medium. The quantity \(u_2^o\) denotes the vertical displacement of the upper reference node \(P(100,100)\), while \(u_2^u\) denotes the vertical displacement of the lower reference node \(P(100,50)\). The difference \(u_2^o-u_2^u\) is used as a scalar measure for the remaining gap between the two blocks.

Fig. [fig:contactPlotVertical] reports the displacement-gap response for case 1 and different values of \(\gamma\). At the beginning of the loading path, the curves are almost identical, which confirms the weak influence of the third medium before contact. After the gap is closed, the curves separate and the effect of the scaling parameter becomes visible.

The trend observed in Fig. [fig:contactPlotVertical] can be interpreted directly from the role of the scaling parameter \(\gamma\). This parameter controls the stiffness level of the fictitious third medium relative to the surrounding solid. Larger values of \(\gamma\) make the intermediate layer more resistant against compression and therefore lead to a more pronounced residual gap after contact has been established, see also Fig. [fig:contactBlocksVertDeformed]. Reducing \(\gamma\) weakens the third medium before contact and allows the two solid bodies to approach each other more closely, which results in a smaller residual gap. This improved contact closure, however, comes at the price of a more demanding nonlinear solution process, since the third medium has to sustain severe compression while carrying only a very small stiffness contribution away from contact. The increasing number of Newton iterations in Tab. [tab:convstudy] reflects this trade-off between contact accuracy and nonlinear robustness. Within the investigated range, all converged curves remain smooth in the contact regime, indicating a stable transition from free motion to contact force transmission.

This visual impression is quantified in Tab. [tab:convstudy], where the final residual gap and the total number of Newton iterations are listed for different penalty parameter combinations of \(\gamma\) and \(p_\Theta\). For \(p_\Theta=1\), reducing \(\gamma\) from \(10^{-3}\) to \(10^{-8}\) decreases the residual gap from \(1.4377\) to \(0.0015\). The same trend is observed for \(p_\Theta=10^{-1}\), where the gap values remain close to those obtained with \(p_\Theta=1\) for the converged cases. The number of Newton iterations increases as \(\gamma\) decreases, which reflects the stronger enforcement of the contact-like response through a softer third medium. For the largest value \(\gamma=10^{-3}\) and the smaller penalty \(p_\Theta=10^{-1}\), no converged solution is obtained.

For case 2, the upper block is subjected to a combined horizontal and vertical displacement with \(\overline{\bu}=[10,-60]^T\). Compared with the purely vertical loading in case 1, this load case additionally introduces tangential motion between the two blocks and therefore activates a more demanding shear-dominated deformation state in the third medium. Fig. [fig:contactPlotShear] shows the corresponding displacement-gap response for different values of \(\gamma\). The overall trend remains consistent with case 1: before contact, the curves are close to each other, while the influence of \(\gamma\) becomes visible once the gap is closed. In the contact regime, larger values of \(\gamma\) again lead to a larger residual gap, whereas smaller values allow a tighter closure of the interface. The zoom in Fig. [fig:contactPlotShear]b highlights that the response remains smooth even under combined compression and tangential loading. This indicates that the proposed stabilization remains robust when the third medium is exposed not only to compression, but also to pronounced shear deformation.

The final deformed configurations for case 2 are shown in Fig. [fig:contactBlocksVertDeformedShear]. The superposed reference configuration illustrates the horizontal shift of the upper block and the resulting distortion of the third medium layer. Compared with the vertical loading case, the deformation state is visibly more asymmetric since the third medium has to accommodate both normal compression and tangential motion. For larger values of \(\gamma\), the intermediate layer remains more pronounced and the two solid blocks stay further apart. For smaller values of \(\gamma\), the interface closes more tightly while the deformation remains regular. The comparison confirms that the proposed stabilization controls the third-medium distortion also in the presence of shear-dominated contact kinematics.

4.3 Self-contact within a box↩︎

Comparison of different interpolations. As a second benchmark, the proposed TMC formulation is tested with the classical self-contact-within-a-box problem shown in Fig. [fig:BoxT1T1]a) and b). A deformable solid frame encloses a third medium region inside the rectangular domain \([0,2]\times[0,0.5]\). Thus, the frame represents the physical solid, while the interior of the box is filled by the fictitious contact medium. For the solid material, the bulk and shear moduli are chosen as \(K=20\) and \(\mu=10\), respectively, following the benchmark setup in [9]. Boundary conditions are prescribed at three points: the lower left point is fixed by \(\overline{\bu}(0,0)=\bzero\), the lower right point is constrained in vertical direction by \(\overline{u}_2(2,0)=0\), and the upper loading point is displaced by \(\overline{\bu}(1,0.5)=[0,-1]^T\). Inside the third medium, the auxiliary-field stabilization introduced above is applied. An unstructured triangular mesh serves as discretization. To assess the role of the interpolation order, the continuous pairs \(\mathrm{T}_{1}^{u}\mathrm{T}_{1}^{\Theta}\), \(\mathrm{T}_{2}^{u}\mathrm{T}_{2}^{\Theta}\) and \(\mathrm{T}_{2}^{u}\mathrm{T}_{1}^{\Theta}\) are compared with the discontinuous variants \(\mathrm{T}_{1}^{u}\mathrm{T}_{1}^{\Theta,d}\), \(\mathrm{T}_{1}^{u}\mathrm{T}_{0}^{\Theta,d}\), and \(\mathrm{T}_{2}^{u}\mathrm{T}_{1}^{\Theta,d}\). Special attention is given to the sensitivity of the contact response with respect to the penalty parameter \(p_\Theta\) and the regularization parameter \(\alpha_r\). As accuracy measure, the residual gap between the upper and lower parts of the deforming frame is evaluated, where smaller values indicate a tighter enforcement of the third-medium contact response. The reference and deformed configurations of the self-contact-within-a-box benchmark are compared in Fig. [fig:BoxT1T1] for two penalty parameters. A coarse unstructured mesh with 354 elements is used in Figs. [fig:BoxT1T1]a) and c), whereas Fig. [fig:BoxT1T1]b) and d) shows a refined discretization with 3208 elements. Both parameter sets are chosen such that a stable contact response with a small residual gap is obtained. As expected for penalty-type couplings, the suitable magnitude of \(p_\Theta\) is problem-dependent and has to be balanced against the third-medium stiffness, the mesh size and the interpolation order. The configurations therefore serve as a first qualitative assessment of the proposed stabilization. A systematic parameter study is carried out in the following to identify robust and accurate choices of \(p_\Theta\) and \(\alpha_r\) for the different interpolation pairs, with the aim of keeping the residual gap between the upper and lower flanges small. The sensitivity study starts with the penalty parameter \(p_\Theta\), as summarized in Tab. [tab:sensitivity95study95ptheta95cont]. For the continuous interpolation pairs, all simulations converge for \(p_\Theta=10^{-1}\) over the considered range of \(\gamma\). This value is therefore selected for the subsequent study of the gradient penalty \(\alpha_r\), reported in Tab. [tab:sensitivity95study95alpha95r95cont]. For this parameter choice, all tested combinations converge, which indicates a robust coupling between the auxiliary field and the deformation gradient.

The parameter studies in Tabs. [tab:sensitivity95study95ptheta95cont], [tab:sensitivity95study95alpha95r95cont] and [tab:sensitivity95study95ptheta95disc] show several consistent trends. For the continuous auxiliary-field interpolations, \(p_\Theta=10^{-1}\) provides a robust parameter range for the present benchmark, since all tested values of \(\gamma\) reach the final load step and produce small residual gaps. Larger penalty parameters do not necessarily improve the contact response. Instead, they can deteriorate the nonlinear convergence behavior and may lead to premature termination of the load path. The regularization parameter \(\alpha_r\) mainly affects the residual gap up to a moderate value. Beyond this range, the additional improvement becomes small, indicating a saturation of the gradient-regularization effect for this benchmark. Comparing the interpolation pairs, the \(\mathrm{T}_{2}^{u}\mathrm{T}_{1}^{\Theta}\) approximation gives results that are very close to those obtained with \(\mathrm{T}_{2}^{u}\mathrm{T}_{2}^{\Theta}\), while using fewer auxiliary degrees of freedom. Thus, enriching the displacement field while keeping a lower-order auxiliary field appears to be an efficient choice for the considered self-contact problem. The \(\mathrm{T}_{1}^{u}\mathrm{T}_{1}^{\Theta}\) approximation remains attractive from a low-order perspective, but its residual gaps are larger than those of the quadratic displacement formulations. The discontinuous auxiliary-field interpolations in Tab. [tab:sensitivity95study95ptheta95disc] show a different behavior. For the quadratic displacement approximation, the \(\mathrm{T}_{2}^{u}\mathrm{T}_{1}^{\Theta,d}\) and \(\mathrm{T}_{2}^{u}\mathrm{T}_{0}^{\Theta,d}\) variants give almost identical results over the investigated parameter range. This indicates that the additional local auxiliary degrees of freedom of the discontinuous linear field do not translate into a stronger regularization effect for this benchmark. Both variants provide small residual gaps for moderate values of \(\gamma\), but the convergence deteriorates when the third medium becomes very soft. For example, at \(p_\Theta=10^{-1}\) both quadratic discontinuous variants reach the final load step down to \(\gamma=10^{-5}\), whereas smaller values of \(\gamma\) lead to premature termination of the load path. The \(\mathrm{T}_{1}^{u}\mathrm{T}_{1}^{\Theta,d}\) interpolation is less favorable in this example, since several parameter sets stop well before the final load step and the results are nearly insensitive to \(p_\Theta\). The element-wise constant \(\mathrm{T}_{1}^{u}\mathrm{T}_{0}^{\Theta,d}\) variant is more robust with respect to load completion for moderate penalty values, but it has to be interpreted as a limiting case because \(\nabla\BTheta\) vanishes within each element. Consequently, it does not provide a true gradient-driven stabilization and mainly reflects the effect of the local penalty coupling. Overall, the behavior of the discontinuous variants is consistent with the one-dimensional discussion in Sec. 4.1: element-local auxiliary fields are attractive from an implementation point of view, but the missing continuity-driven inter-element coupling weakens the gradient-type regularization mechanism compared with the continuous auxiliary-field formulations.

Comparison with deformation-gradient averaging - simple box. A comparison with established regularization strategies is essential to assess both the robustness and the computational efficiency of the proposed formulation. Recently, Faltus et al. [18] introduced a deformation-gradient averaging approach that avoids additional degrees of freedom and therefore does not enlarge the global system of equations. In their formulation, a representative deformation gradient \(\bar{\bF}\) is computed at the element center and enforced at the remaining integration points by a penalty contribution. In this way, variations of \(\bF\) inside an element are suppressed without explicitly regularizing \(\nabla\bF\). The method can therefore be interpreted as an element-wise deformation-gradient regularization. For a fair comparison, the same benchmark and the same \(Q_1\) displacement interpolation are used in both formulations. The formulation by Faltus et al. is evaluated with a \(Q_1^u\) discretization, while the present mixed formulation uses a \(Q_1^uQ_1^\Theta\) interpolation. Although the proposed approach introduces additional auxiliary degrees of freedom, it also provides a direct mixed representation of the deformation-gradient-like stabilization field. The parameters are kept fixed for all mesh refinements. For the proposed formulation, \(p_\Theta=10\), \(\alpha_r=1000\) and \(\gamma=10^{-6}\) are used. For the deformation-gradient averaging approach, the penalty parameter is chosen as \(\kappa_{\bar F}=0.5\) with \(\gamma_{\bar F}=10^{-6}\).

4pt

The results in Tab. 2 show the expected increase in system size caused by the additional auxiliary field of the proposed mixed formulation. For all mesh refinements, however, this larger system size is accompanied by a substantially more stable nonlinear load path. The proposed formulation consistently requires fewer load steps than the deformation-gradient averaging approach and also reduces the total number of Newton iterations for all considered meshes up to \(409600\) elements. At \(102400\) elements, the computation reaches the final configuration in \(25\) load steps and \(153\) iterations, compared with \(113\) load steps and \(868\) iterations for the formulation by Faltus et al. [18]. At this mesh level, the improved nonlinear robustness compensates for the additional algebraic cost and reduces the total time from \(183.108\,\mathrm{s}\) to \(144.586\,\mathrm{s}\). The largest discretization with \(409600\) elements represents an extremely fine mesh and leads to more than \(1.7\) million unknowns for the mixed formulation. At this scale, the cost of the additional auxiliary degrees of freedom becomes dominant and the total solution time increases to \(1895.320\,\mathrm{s}\), compared with \(803.366\,\mathrm{s}\) for the deformation-gradient averaging approach. Even in this demanding case, however, the proposed formulation reaches the final load in only \(30\) load steps, whereas the reference approach requires \(108\) load steps. Thus, the mixed auxiliary-field formulation trades a larger algebraic system for a considerably more robust load stepping behavior, with the overall efficiency depending on the balance between nonlinear robustness and the cost of the additional auxiliary unknowns.

4.4 Self-contact of a C-shaped box↩︎

The proposed TMC stabilization is further assessed with the classical C-shaped box benchmark shown in Fig. [fig:BVP95CShape]. A two-dimensional C-shaped hyperelastic solid with outer dimensions \(1000\times500\) and a wall thickness of \(100\) is considered. Inside the gap, a third medium fills the open space, while an additional third-medium layer of thickness \(12.5\) is placed at the right boundary to allow contact between the upper and lower arms. Along the left boundary, the solid is fully clamped. At the upper right corner \(P(1000,500)\), a prescribed displacement \(\bar{\bu}=(0,\bar{u}_y)^T\) drives the closing motion of the C-shape. Large rotations, severe compression of the third medium and progressive self-contact are therefore induced by a simple displacement-controlled loading. For the solid, the compressible Neo-Hookean material in Eq. 4 is used with \(K_{\mathrm{s}}=5/3\) and \(\mu_{\mathrm{s}}=5/14\). As a first check, the influence of the discretization is assessed with the continuous \(\mathrm{T}_{1}^{u}\mathrm{T}_{1}^{\Theta}\) interpolation. Fig. [fig:CShape95meshes] shows the deformation evolution for four increasingly refined meshes.

Since the penalty coupling is mesh-dependent, the parameter \(p_\Theta\) is adjusted for each discretization rather than kept fixed. This procedure should therefore not be interpreted as a classical mesh-convergence study with one fixed parameter set, but as a calibrated mesh sequence used to assess whether the same deformation mechanism can be obtained over a range of discretizations. For the four meshes shown in Fig. [fig:CShape95meshes], the values \(p_\Theta=\{0.001,0.005,0.01,0.05\}\) are used. Across all meshes, the same global deformation mode is obtained, i.e., the upper arm bends into the cavity, the third medium is squeezed into a narrow band, and self-contact develops along the inner side of the C-shape. With increasing mesh resolution, the deformation path becomes smoother and the compressed third-medium layer is represented more sharply. Qualitative agreement across all calibrated meshes indicates that the proposed stabilization is not tied to a particular mesh density.

A more detailed view of the deformation process is given in Fig. [fig:CShapeT1T1]. At early loading, the upper arm starts to rotate downward while the third medium remains broadly distributed inside the cavity. With increasing load, the gap closes and the third medium is progressively compressed between the approaching solid surfaces. At the final load step, deformation localizes along the contact region, while the third-medium mesh remains regular and no element collapse is observed. Even with the continuous \(\mathrm{T}_{1}^{u}\mathrm{T}_{1}^{\Theta}\) discretization, the large rotation and subsequent self-contact are captured robustly and the overall deformation pattern remains smooth and stable throughout the loading path.

Fig. [fig:CShape95LoadDisp] shows the vertical reaction force as a function of the prescribed vertical displacement \(u_y\). The reaction force is obtained from the vertical support reaction associated with the imposed displacement at the loading point. Over most of the loading range, the TMC curve follows the NTS reference closely. During the initial bending-dominated phase, both curves are almost indistinguishable and remain in good agreement after contact is established. Only in the magnified initial contact regime in Fig. [fig:CShape95LoadDisp]b, a small deviation becomes visible. Such a deviation is expected for a continuum contact regularization with finite third-medium stiffness. For engineering-scale simulations, the difference remains small compared with the overall force level, while the TMC formulation avoids explicit contact search and provides a smooth transition into self-contact.

Comparison with deformation-gradient averaging - C-box. A second comparison with the deformation-gradient averaging approach of Faltus et al. [18] is performed for the C-shaped benchmark. Here, the focus is not on computational cost, but on the robustness of the third-medium deformation under different interpolation choices. Fig. [fig:CShape95Faltus] shows the deformed configurations at the load parameter \(\lambda_{\mathrm{load}}=0.2\). \(\bar \bF\) is calculated at the center points of each individual element. A similar approach calculating the volume averaged \(\bar \bF=\frac{1}{V}\int_{\B}\bF\)dv reaches the same results.

The continuous auxiliary-field formulation provides the most robust response in this comparison. For both triangular and quadrilateral discretizations, the third medium remains confined to the closing gap and follows the deformation of the C-shaped structure without visible element collapse. Although the result of the element \(\mathrm{T}_{1}^{u}\mathrm{T}_{1}^{\Theta}\) is not displayed, the same behavior is observed in corresponding simulations.

The element-wise discontinuous variants behave fundamentally differently. Since \(\nabla\BTheta\) vanishes inside each element for discontinuous low-order auxiliary interpolations, the gradient regularization does not provide an inter-element stabilization mechanism. For linear triangles the difference between \(\boldsymbol{F}\) and \({\boldsymbol{\Theta}}\) is zero since the deformation gradient \(\boldsymbol{F}\) is constant in the element. Thus the regularization cannot work. Furthermore, it is observed that for quadratic triangles the regularization also does not work. This is due to the fact that the displacement \(\boldsymbol{u}_h\) is nearly a self-affine function for T\(_2\) triangle. We note, that a self-affine (linear) displacement field occurs in the quadratic triangle when the displacements at the mid points of the element are given by \[\bu_4=\frac{\bu_1+\bu_2}{2}\,,\quad \bu_5=\frac{\bu_2+\bu_3}{2}\,\quad and\,\,\, \bu_6=\frac{\bu_3+\bu_2}{2}. \label{eq:mid95nodes}\tag{20}\] The resulting linear displacement field yields a constant deformation gradient \(\boldsymbol{F}\) and thus, the direct stabilization term Eq. 9 is zero or the difference in Eq. 10 is zero and the regularization does not work. Furthermore, it can be shown by a Taylor series expansion that, e.g., for the displacements of a mid node \(4\) of the element at \(\bx_m\) which lies between the vertex nodes \(\bx_i\) and \(\bx_j\) that the following estimate holds with \(\bh=\frac{1}{2}(\bx_j -\bx_i)\) for the quadratic displacement field \[\begin{align} \bu(\bx_m+\bh) &= \bu(\bx_m) + \nabla\bu(\bx_m)\,\bh + \frac{1}{2} \nabla^2\bu(\bx_m) [\bh,\bh] + \mathcal{O}(\|\bh\|^3)\\ \bu(\bx_m-\bh) &= \bu(\bx_m) - \nabla\bu(\bx_m)\,\bh + \frac{1}{2} \nabla^2\bu(\bx_m) [\bh,\bh] + \mathcal{O}(\|\bh\|^3)\,. \end{align}\] Adding the two equations and solving for the mid point displacement yields \[\bu(\bx_m) = \frac{\bu_i+\bu_j}{2} - \frac{1}{2} \nabla^2\bu(\bx_m) [\bh,\bh] + \mathcal{O}(\|\bh\|^4)\] and therefore, to second-order accuracy the displacement at the mid nodes is provided by one half of the sum of the vertex nodes. This means that Eq. 20 is fulfilled up to second order accuracy. We note that with \(\mathbf{h} = \frac{h}{2}\mathbf{t}\) we obtain \[\frac{1}{2}\nabla^2\bu(\mathbf{x}_m)[\bh,\bh] = \frac{h^2}{8} \frac{\partial^2\bu}{\partial s^2}(\mathbf{x}_m).\] where \(\mathbf{t}\) is the unit tangent vector at the edge and \(s\) the coordinate along the edge.

For the quadratic triangular element, \(\bu_m\) is still an independent nodal displacement. It is generally not exactly the average of the end-node displacements. But the difference between \(\bu_m\) and the average is of order \(h^2\), so it vanishes under mesh refinement. It is actually zero for an affine (linear) displacement field where \(\nabla^2 \bu_h=\mathbf{0}\). Thus the effect of the regularisation functional \(W^{\textrm{tm}}_p\) is considerably small and goes to zero for finer meshes and does not work as regularization. Even changing the penalty parameter \(p_\Theta\) cannot prevent the observed element collapse and local self-penetration, and the parameter \(\alpha_r\) has no effective influence in this setting. The deformation-gradient averaging yields the same result.

Hence, for the triangular discretization, see Fig. [fig:CShape95Faltus] c), the third medium is pushed laterally out of the gap and forms an artificial extrusion at the right boundary, whereas the quadrilateral case remains stable. Overall, the comparison shows that the stabilizing effect of the proposed method relies on a continuous auxiliary-field interpolation when low-order gradient information is needed. Discontinuous variants remove this mechanism, while deformation-gradient averaging remains more dependent on the underlying element type.

4.5 Three-dimensional self-contact within a box↩︎

Finally, the self-contact-within-a-box benchmark is extended to a three-dimensional setting in order to demonstrate that the proposed formulation is not restricted to plane problems. Following the two-dimensional setup, a deformable solid frame encloses a third medium region, and the imposed displacement drives the structure into self-contact. A low-order \(\mathrm{T}_{1}^{u}\mathrm{T}_{0}^{\Theta,d}\) tetrahedral discretization with an element-wise auxiliary field is used in the third medium and selected material parameters \(\gamma=10^{-5}\), \(\alpha_r=1\) and \(p_\Theta=1\). Fig. [fig:BoxT1T0953D] shows the deformed configuration from two perspectives. During loading, the upper part of the frame bends into the box, the third medium is compressed between the approaching solid surfaces, and the deformation remains regular. Overall, the example confirms that the auxiliary-field stabilized third medium formulation carries over directly to three-dimensional finite deformation contact problems.

5 Conclusion and outlook↩︎

In this work, an auxiliary-field stabilized mixed finite element formulation for third medium contact at finite deformations is presented. The main idea is to introduce a deformation-gradient-like field \(\BTheta\) in the third medium and to couple it weakly to the physical deformation gradient \(\bF\) by a penalty contribution. Regularization acts on \(\nabla\BTheta\) rather than on \(\nabla\bF\) directly. In this way, gradient-type regularization is incorporated without evaluating second derivatives of the displacement field, which makes the approach applicable to first- and second-order finite elements.

A central observation is that the regularization effect depends strongly on the interpolation of the auxiliary field. For continuous auxiliary-field interpolations, neighboring elements share the nodal values of \(\BTheta\). Element-wise changes of the deformation gradient are therefore transferred into spatial variations of the auxiliary field, which are penalized by the regularization parameter \(\alpha_r\). Penalizing these spatial variations stabilizes the third medium and reduces the residual gap in the considered contact benchmarks. Element-wise discontinuous auxiliary fields are attractive from an implementation point of view, since their additional unknowns remain local to the element. However, the numerical results show that the missing continuity-driven inter-element coupling weakens the gradient-type stabilization mechanism compared with the continuous auxiliary-field formulations.

Across the numerical examples, large deformation contact, severe third-medium compression and progressive self-contact are captured in a stable manner. In the two-block benchmark, the scaling parameter \(\gamma\) controls the trade-off between residual gap and nonlinear solution effort, i.e., smaller values improve contact closure, while the Newton process becomes more demanding. The self-contact-within-a-box example confirms that moderate penalty parameters provide a robust coupling between \(\BTheta\) and \(\bF\), whereas overly large values of \(p_\Theta\) may deteriorate convergence. For quadratic displacement approximations, the \(\mathrm{T}_{2}^{u}\mathrm{T}_{1}^{\Theta}\) interpolation gives results close to the \(\mathrm{T}_{2}^{u}\mathrm{T}_{2}^{\Theta}\) formulation while using fewer auxiliary degrees of freedom. The C-shaped benchmark further shows that the proposed third medium formulation gives a force-displacement response close to a classical augmented-Lagrange node-to-segment contact formulation, while retaining the main advantage of third medium contact, since no explicit contact search, active-set strategy or predefined contact interface is required.

Overall, the formulation combines the smooth continuum character of third medium contact with a low-order-compatible stabilization mechanism. Contact forces are transmitted through a fictitious continuum region, and self-contact develops naturally as part of the deformation process. Such a continuum-based treatment is particularly attractive for problems in which potential contact zones are not known in advance, for example in strongly deforming structures, contact-aided mechanisms and metamaterials with internal self-contact.

Several aspects remain open for future work. A systematic parameter selection strategy for \(\gamma\), \(p_\Theta\) and \(\alpha_r\) is an important next step, especially for complex three-dimensional applications. Further work should address adaptive choices of the third-medium stiffness and more efficient solution strategies for large-scale problems. From an application point of view, the method appears especially promising for the analysis and design of mechanical metamaterials that undergo large shape changes and repeated internal contact.

Acknowledgement↩︎

The authors gratefully acknowledge the computing time granted by the Center for Computational Sciences and Simulations (CCSS) of the University of Duisburg-Essen and provided on the supercomputer magnitUDE at the Zentrum für Informations- und Mediendienste (ZIM).

References↩︎