# FEM approximation of dynamic contact problem for fracture under fluid volume control using generalized HHT-α and semi-smooth Newton methods

Victor A. Kovtunenko, Yves Renard

► To cite this version:

Victor A. Kovtunenko, Yves Renard. FEM approximation of dynamic contact problem for fracture under fluid volume control using generalized HHT-α and semi-smooth Newton methods. Applied Numerical Mathematics: an IMACS journal, 2025, 218, pp.148-158. 10.1016/j.apnum.2025.07.009. hal-05235234

HAL Id: hal-05235234

https://hal.science/hal-05235234v1

Submitted on 15 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Research Paper

# FEM approximation of dynamic contact problem for fracture under fluid volume control using generalized HHT-$\alpha$ and semi-smooth Newton methods

<figure><img src="image_2.png" /></figure>

Check for
updates

Victor A. Kovtunenko<sup>a,b, c, d, e</sup>, Yves Renard<sup>c, d</sup>

## ARTICLE INFO

MSC:

35L85

49M15

74R20

74S05

Keywords:

Fluid-driven fracture

Volume control

Non-penetrating crack

Elastodynamic contact

Semi-smooth Newton method

## ABSTRACT

A class of elastodynamic contact problems for fluid-driven cracks stemming from hydro-fracking application is considered in the framework of finite element approximation. The dynamic contact problem aims at finding a non-negative fracture opening and a mean fluid pressure which are controlled by the volume of pumped fracturing fluid. Well-posedness of the fully discrete variational problem is proved rigorously by using the Lagrange multiplier and penalty methods for the minimization problem subjected to both: unilateral and non-local constraints. Numerical solution of the dynamic nonlinear equation is computed in 2D experiments using the semi-smooth Newton and the generalized Hilber–Hughes–Taylor $\alpha$-method.

## 1. Introduction

The hydro-fracturing is widely used in petroleum engineering applications pumping out oil, natural gas, and geothermal water from reservoirs in the earth. Hydro-fractures are created by injecting a fracturing fluid at high pressure through a well bore and usually controlled by the injection rate into inlet holes providing the pumped fluid volume in time. Mathematical modeling of reservoirs with hydraulic fractures saturated by a Newtonian fluid can be found in the works by [1–4] and others. For analysis of crack evolution based on a method of balanced viscosity solution we cite [5]. In order to describe mass balance of fluid in the open crack in an averaged way, a volume constraint was suggested in [6]. In [7] we introduced complementarity conditions providing non-penetration between opposite faces of a fluid-driven crack within variational formulation of the poroelastic model. For the variational theory of elastostatic models of solids with non-penetrating cracks we refer the readers to [8–11], and to [12] for its numerical treatment.

In the framework of computational contact and impact mechanics (see [13]), well-posedness of elastodynamics contact problems remains in general an open question. It can be recovered or established by employing a viscous damping together with the regularization of contact conditions. Using penalization and regularization arguments, mathematical results on the existence of solution

<footnote><sup>a</sup> Department of Mathematics and Scientific Computing, Karl-Franzens University of Graz, NAWI Graz, Heinrichstr. 36, 8010 Graz, Austria</footnote>

<footnote><sup>b</sup> Lavrentyev Institute of Hydrodynamics, Siberian Division of the Russian Academy of Sciences, 630090 Novosibirsk, Russia</footnote>

<footnote><sup>c</sup> INSA Lyon, CNRS, Ecole Centrale de Lyon, Université Claude Bernard Lyon 1, Université Jean Monnet, ICJ UMR5208, LaMCoS UMR5259, Villeurbanne, 69621, France</footnote>

<footnote>* Corresponding author at: Department of Mathematics and Scientific Computing, Karl-Franzens University of Graz, NAWI Graz, Heinrichstr. 36, 8010 Graz, Austria. E-mail addresses: victor.kovtunenko@uni-graz.at (V.A. Kovtunenko), Yves.Renard@insa-lyon.fr (Y. Renard).</footnote>
of dynamic viscoelastic models were obtained in [14–16]. For solvability of the elastodynamic problem for non-penetrating cracks with modified contact law we cite [17]. The general theory of coupled systems of variational and hemivariational inequalities can be found in the works by [18,19] and co-authors. We refer the readers to [20] for relevant theories and uncertainty quantification in variational inequalities, and to [21] for a penalization treatment of optimization problems under volume constraints of non-local type.

For solution of complementarity problems, a semi-smooth Newton method was proposed, see [22], based on a generalized gradient of non-smooth merit functions. The Newton iteration converges locally super-linear, and globally monotone when implemented in the form of primal-dual active set (PDAS) algorithm. Generalized Newton's methods are used for solution of contact and frictional problems in mechanics, see [23,24]. The PDAS strategy was applied to the Signorini's contact problem with Coulomb friction in [25], to dynamic frictional contact problem in [26], and to non-penetrating crack problems in [27,28]. In [29] we study rigorously convergence of the semi-smooth Newton method and equivalent PDAS algorithm for a FEM approximation of dynamic two-body contact and crack problems.

The theory of finite element method (FEM) approximation of elastodynamic contact problems can be found in [30,31]. The space semi-discretization is ill-posed and standard time approximation schemes are not stable with respect to the mechanical energy when decreasing time step or increasing simulation time. To stabilize the FEM discretization, methods of restitution coefficient [32], mass redistribution [33,34] were suggested. Between them, Nitsche's method [30] is consistent and has better convergence properties than a penalty method. For temporal discretization, the Hilber–Hughes–Taylor method [35] called HHT-$\alpha$ is adopted in the literature, see its further developments in [36,37]. In [38] the schemes of Newmark family described by two weight parameters $\gamma$ and $\beta$ were considered as the particular case when $\alpha = 1$. The generalized HHT-$\alpha$ scheme with the vector parameter $\alpha$ was introduced in [39] because it is more high frequency dissipative.

In the present work, in Section 2 we formulate the elastodynamic problem for crack subjected to non-penetration and volume constraints. In Section 3 the corresponding FEM discretization in space and the generalized HHT-$\alpha$ scheme in time are introduced, and its variational solution is proved rigorously. We present the semi-smooth Newton iteration and equivalent PDAS algorithm for the penalized problem in Section 4. Its properties of locally super-lineal and globally monotone convergence will be validated numerically. Motivated by hydro-fracking, in Section 5 we simulate the extension-compression loop motion of a symmetric 2D isotropic body with a fluid-driven crack controlled by the pumped fluid volume.

## 2. Setting of the problem

We model a reservoir consisted of two layers separated by a single fracture by two deformable elastic bodies coming into contact along the common interface, see Fig. 1. In the undeformed configuration, let $i$-th body occupy the domain $\Omega^i$ in $\mathbb{R}^d$, $d = 2, 3$, where $i = 1$ stands for the first and $i = 2$ for the second body. For simplicity we suppose that the Lipschitz continuous boundary $\partial\Omega^i$ comprises two non-overlapping parts $\Gamma_D^i$ and $\Gamma_C^i$. Let the bodies be clamped on the Dirichlet boundary $\Gamma_D^i$ and touch each other at the common portion of the boundary $\Gamma_C^i = \Gamma_C^2 := \Gamma_C$ having a unit normal $\mathbf{n}$. We choose the direction such that $\mathbf{n}$ is inward $\Omega^i$ and outward $\Omega^2$ (see Fig. 1). Let us define the domain $\Omega = \Omega^1 \cup \Omega^2$ with the crack $\Gamma_C$ and the outer boundary $\Gamma_D := \Gamma_D^i \cup \Gamma_D^2$.

For a time $t \in [0, T)$, where the final time $T > 0$ is prescribed, we denote the cylinder $\Omega_T = (0, T) \times \Omega$ with two side surfaces $\Gamma_{DT} = (0, T) \times \Gamma_D$ and $\Gamma_{CT} = (0, T) \times \Gamma_C$. Under small deformation assumption, the displacement field $\mathbf{u}$ in $\Omega_T$ builds the linearized strain tensor as symmetric part of the gradient:

$$ \varepsilon(\mathbf{u}) = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T), $$

and the stress tensor field is determined from Hooke's law:

$$ \sigma(\mathbf{u}) = A\varepsilon(\mathbf{u}) \quad \text{in } \Omega_T. \qquad (2.1) $$

In (2.1) the fourth order symmetric tensor of elastic coefficients $A$ has the usual properties of uniform ellipticity and boundedness. For the volume force $f$ prescribed in $\Omega_T$ and constant density of the elastic material $\rho$, the elastodynamic problem consists in finding the displacement field $\mathbf{u}: [0, T) \times \Omega \to \mathbb{R}^d$ and acceleration $\ddot{\mathbf{u}}$ satisfying the constitutive equation (2.1) and the equation of motion:

$$ \rho\ddot{\mathbf{u}} - \nabla \cdot \sigma(\mathbf{u}) = \mathbf{f} \quad \text{in } \Omega_T, \qquad (2.2) $$

using the tensor contraction “$\cdot$”, which is supported by the initial conditions:

$$ \mathbf{u}(0, \cdot) = \mathbf{u}_0, \quad \dot{\mathbf{u}}(0, \cdot) = \dot{\mathbf{u}}_0 \quad \text{in } \Omega. \qquad (2.3) $$

In (2.3) the displacement $\mathbf{u}$ and velocity $\dot{\mathbf{u}}$ at $t = 0$ are prescribed by the initial fields $\mathbf{u}_0$ and $\dot{\mathbf{u}}_0$. Equations (2.1)–(2.3) are endowed with the homogeneous Dirichlet condition on the outer boundary surface:

$$ \mathbf{u} = \mathbf{0} \quad \text{on } \Gamma_{DT} \qquad (2.4) $$

and conditions on the contact surface $\Gamma_{CT}$ stated below.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Example crack problem in the undeformed 2D configuration.</figcaption></figure>

For arbitrary displacement field $\mathbf{v}$ we adopt its orthogonal decomposition into the normal $v_n$ are tangential $\mathbf{v}_\Gamma$ components:

$$ \mathbf{v} = v_n \mathbf{n} + \mathbf{v}_\Gamma, \quad v_n := \mathbf{v} \cdot \mathbf{n} \quad \text{on } \Gamma_C. $$

Defining the relative jump of fields across the contact surface:

$$ \llbracket \mathbf{u} \rrbracket = \mathbf{u} \Big|_{\Gamma_C^1} - \mathbf{u} \Big|_{\Gamma_C^2}, $$

the crack opening should be non-negative:

$$ \llbracket u_n \rrbracket = \llbracket \mathbf{u} \rrbracket \cdot \mathbf{n} \ge 0 \quad \text{on } \Gamma_{CT}, \qquad (2.5) $$

and the normal stress $(\boldsymbol{\sigma}(\mathbf{u})\mathbf{n}) \cdot \mathbf{n} =: \sigma_n(\mathbf{u})$ should be continuous:

$$ \llbracket \sigma_n(\mathbf{u}) \rrbracket = 0 \quad \text{on } \Gamma_{CT}. \qquad (2.6) $$

Complementary to (2.5) and (2.6) conditions for the contact force $\lambda$ are (see [7]):

$$ \sigma_n(\mathbf{u}) + P := \lambda \le 0, \quad \lambda \llbracket u_n \rrbracket = 0 \quad \text{on } \Gamma_{CT}, \qquad (2.7) $$

where the mean fluid pressure $P: [0, T) \mapsto \mathbb{R}$ is to find, and negative sign of $\lambda$ corresponds to compression. The friction-free contact implies zero tangential stress:

$$ \boldsymbol{\sigma}(\mathbf{u})\mathbf{n} - \sigma_n(\mathbf{u})\mathbf{n} := (\boldsymbol{\sigma}(\mathbf{u})\mathbf{n})_\Gamma = \mathbf{0} \quad \text{on } \Gamma_{CT}. \qquad (2.8) $$

Let a continuous function $A(t) > 0$ starting at $A(0) = 0$ be given. It describes a volume of fluid pumped into the fracture keeping the fluid volume balance:

$$ \int_{\Gamma_C} \llbracket u_n \rrbracket d\Gamma = A(t) \quad \text{for } t \in [0, T). \qquad (2.9) $$

Gathering together governing relations we validate (2.1)–(2.9) in function spaces. For this task we introduce the Hilbert space owing to the Dirichlet condition (2.4):

$$ \mathbf{V} = \{ \mathbf{v} \in H^1(\Omega)^d : \mathbf{v} = \mathbf{0} \text{ a.e. } \Gamma_D \}, $$

and the corresponding Bochner space

$$ \mathbf{W} = \{ \mathbf{v} \in L^2(0, T; \mathbf{V}), \dot{\mathbf{v}} \in L^2(0, T; L^2(\Omega)^d) \}. $$

Admissible displacements satisfying the non-penetration (2.5) build the convex cone:

$$ \mathbf{K} = \{ \mathbf{v} \in \mathbf{V} : \llbracket v_n \rrbracket \ge 0 \text{ a.e. } \Gamma_C \}. $$

Let the force $\mathbf{f} \in C([0, T]; L^2(\Omega)^d)$ and initial fields in (2.3) satisfy $\mathbf{u}_0 \in \mathbf{K}$ with $\llbracket u_{0n} \rrbracket = 0$ and $\dot{\mathbf{u}}_0 \in L^2(\Omega)^d$. For smooth $\mathbf{u}$ and $\mathbf{v}$ the Green's formula holds (see [10]):

$$ -\int_{\Omega} \nabla \cdot \boldsymbol{\sigma}(\mathbf{u}) \cdot \mathbf{v} \, d\mathbf{x} = \int_{\Omega} \boldsymbol{\sigma}(\mathbf{u}) : \varepsilon(\mathbf{v}) \, d\mathbf{x} - \int_{\Gamma_D} \boldsymbol{\sigma}(\mathbf{u})\mathbf{n} \cdot \mathbf{v} \, d\Gamma + \int_{\Gamma_C} \llbracket \sigma_n(\mathbf{u})v_n + (\boldsymbol{\sigma}(\mathbf{u})\mathbf{n})_\Gamma \cdot \mathbf{v}_\Gamma \rrbracket \, d\Gamma, \qquad (2.10) $$

where “ : ” implies the tensor double contraction. Testing (2.10) with $\mathbf{v} - \mathbf{u}$ after substitution of the equation of motion (2.2) and boundary conditions (2.6), (2.8) yields
$$ \int_{\Omega} ((\rho \ddot{\mathbf{u}} - \mathbf{f}) \cdot (\mathbf{v} - \mathbf{u}) + \sigma(\mathbf{u}) : \varepsilon(\mathbf{v} - \mathbf{u})) d\mathbf{x} = - \int_{\Gamma_C} \sigma_n(\mathbf{u}) [\![v_n - u_n]\!] d\Gamma \quad \text{for } t \in (0, T). \quad (2.11) $$

Integration of (2.11) by parts over time using initial conditions (2.3) and inequalities (2.5), (2.7) leads to the variational formulation of the problem (2.1)–(2.9):

$$ \left\{
\begin{array}{ll}
\text{Find } \mathbf{u} \in \mathbf{W}, P(t) \text{ and } \mathbf{u}(t, \cdot) \in \mathbf{K} \text{ for } t \in (0, T) & \\
\text{such that: } \mathbf{u}(0, \cdot) = \mathbf{u}_0, \quad \displaystyle\int_{\Gamma_C} [\![u_n]\!] d\Gamma = A(t) \text{ for } t \in [0, T), & \\
\displaystyle\int_{\Omega_T} (-\rho \ddot{\mathbf{u}} \cdot (\dot{\mathbf{v}} - \ddot{\mathbf{u}}) + \sigma(\mathbf{u}) : \varepsilon(\mathbf{v} - \mathbf{u})) d\mathbf{x} dt - \displaystyle\int_{\Gamma_{CT}} P[\![v_n - u_n]\!] d\Gamma dt & (2.12) \\
\geq \displaystyle\int_{\Omega} \rho \ddot{\mathbf{u}}_0 \cdot (\mathbf{v}(0, \cdot) - \mathbf{u}_0) d\mathbf{x} + \displaystyle\int_{\Omega_T} \mathbf{f} \cdot (\mathbf{v} - \mathbf{u}) d\mathbf{x} dt & \\
\text{for all } \mathbf{v} \in \mathbf{W}, \mathbf{v}(t, \cdot) \in \mathbf{K} \text{ and } \mathbf{v} = \mathbf{u} \text{ for } t \geq T - \zeta \text{ with } \zeta > 0. &
\end{array}
\right. $$

We note that, in the strong formulation P can be reduced from the system when inserting the fluid volume balance equation (2.9) and $\mathbf{v} = \mathbf{0}$ into (2.11) such that

$$ P = \frac{1}{A} \int_{\Omega} ((\rho \ddot{\mathbf{u}} + \mathbf{f}) \cdot \mathbf{u} - \sigma(\mathbf{u}) : \varepsilon(\mathbf{u})) d\mathbf{x} \quad \text{for } t \in (0, T). \quad (2.13) $$

Existence of a solution to the elastodynamic contact problem (2.12) in the function setting is open.

### 3. Discretization of the variational problem

Let $\mathcal{T}_h$ be a regular triangulation of the domain $\Omega$ of the mesh size $h = \max_{K \in \mathcal{T}_h} h_K$, where $h_K$ is the diameter of $K$. We introduce the finite-dimensional vector space $\mathbf{V}_h$ which is built by piecewise polynomials of degree $p \in \mathbb{N}$ on $\mathcal{T}_h$:

$$ \mathbf{V}_h = \{ \mathbf{v}_h \in C^0(\bar{\Omega})^d : \mathbf{v}_h|_K \in (\mathbb{P}_p(K))^d \text{ for all } K \in \mathcal{T}_h, \mathbf{v}_h = \mathbf{0} \text{ on } \Gamma_D \}. $$

The triangulation is supposed conformal to the subdivision of crack surfaces $\Gamma_C^{1h} \subset \Gamma_C^1$ and $\Gamma_C^{2h} \subset \Gamma_C^2$, which coincide. On the discrete set $\Gamma_C^h \subset \Gamma_C$ of $N_C^h$ nodal points $\mathbf{x}_h$ at the crack, we look for the nodal variables $\lambda_h(t)$ verifying complementarity conditions (2.5)–(2.7) on $\Gamma_{CT}^h = (0, T) \times \Gamma_C^h$. The space semi-discretized problem (2.12) in the mixed variational form reads:

$$ \left\{
\begin{array}{ll}
\text{Find } \mathbf{u}_h : (0, T) \mapsto \mathbf{V}_h, \lambda_h : (0, T) \mapsto \mathbb{R}^{N_C^h} \text{ and } P : (0, T) \mapsto \mathbb{R} \text{ such that:} & \\
\mathbf{u}_h(0, \cdot) = \mathbf{u}_{0h}, \quad \displaystyle\int_{\Gamma_C} [\![u_{hn}]\!] d\Gamma = A(t) \text{ for } t \in [0, T), & \\
\displaystyle\int_{\Omega_T} (-\rho \ddot{\mathbf{u}}_h \cdot \dot{\mathbf{v}}_h + \sigma(\mathbf{u}_h) : \varepsilon(\mathbf{v}_h)) d\mathbf{x} dt + \displaystyle\int_0^T \sum_{\Gamma_C^h} (\lambda_h - P)[\![v_{hn}]\!] (\mathbf{x}_h) dt & (3.1) \\
= \displaystyle\int_{\Omega} \rho \ddot{\mathbf{u}}_0 \cdot \mathbf{v}(0, \cdot) d\mathbf{x} + \displaystyle\int_{\Omega_T} \mathbf{f} \cdot \mathbf{v}_h d\mathbf{x} dt \quad \text{for all } \mathbf{v}_h \in C^0([0, T]; \mathbf{V}_h), & \\
[\![u_{hn}]\!] (\mathbf{x}_h) \ge 0, \quad \lambda_h \le 0, \quad \lambda_h [\![u_{hn}]\!] (\mathbf{x}_h) = 0 \quad \text{for } \mathbf{x}_h \in \Gamma_{CT}^h, &
\end{array}
\right. $$

where approximation of the initial displacement $\mathbf{u}_{0h} \in \mathbf{V}_h$ and velocity $\dot{\mathbf{u}}_{0h} \in \mathbf{V}_h$ has been used.

For given $N \in \mathbb{N}$ and time step size $\tau = T/N$, let us consider the uniform time discretization by points $t^m = m\tau$, $m = 0, \dots, N$, where $t^0 = 0$ and $t^N = T$. At times $t^m$ we discretize the displacement (respectively velocity and acceleration) by $\mathbf{u}_h^m$ (respectively $\dot{\mathbf{u}}_h^m$ and $\ddot{\mathbf{u}}_h^m$), the force $\mathbf{f}^m = \mathbf{f}(t^m, \cdot)$, the Lagrange multiplier $\lambda_h^m = \lambda_h(t^m)$, and the mean pressure $P^m = P(t^m)$. Moreover, we introduce intermediate time steps $m + \alpha_2$ weighted by parameter $\alpha_2 > 0$ in the following way:

$$ \mathbf{v}_h^{m+\alpha_2} = \alpha_2 \mathbf{v}_h^{m+1} + (1-\alpha_2) \mathbf{v}_h^m. \quad (3.2) $$

Let weights $\gamma \in [0, 1]$, $\beta \in [0, 0.5]$, and $\boldsymbol{\alpha} = (\alpha_1, \alpha_2) \in [0, 1]^2$. Using (3.2) and following [39], a generalized HHT-$\boldsymbol{\alpha}$ method applied to (3.1) has the following form for $m \ge 0$:
$$
\left\{
\begin{array}{l}
\text{Find } \mathbf{u}_h^{m+1}, \dot{\mathbf{u}}_h^{m+1}, \ddot{\mathbf{u}}_h^{m+1} \in \mathbf{V}_h, \lambda_h^{m+1} \in \mathbb{R}^{N_C^h}, \text{ and } P^{m+1} \in \mathbb{R} \text{ such that:} \\
\mathbf{u}_h^{m+1} = \mathbf{u}_h^m + \tau \dot{\mathbf{u}}_h^m + \frac{\tau^2}{2} \ddot{\mathbf{u}}_h^{m+2\beta}, \\
\dot{\mathbf{u}}_h^{m+1} = \dot{\mathbf{u}}_h^m + \tau \dot{\mathbf{u}}_h^{m+\gamma}, \\
\int_{\Gamma_C} \left[ \|u_{hn}^{m+\alpha_2}\| d\Gamma = A^{m+\alpha_2}, \quad \int_{\Omega} (\rho \dot{\mathbf{u}}_h^{m+\alpha_1} \cdot \mathbf{v}_h + \sigma(\mathbf{u}_h^{m+\alpha_2}) : \boldsymbol{\varepsilon}(\mathbf{v}_h)) \right] d\mathbf{x} \\
+ \sum_{\Gamma_C^h} (\lambda_h^{m+\alpha_2} - P^{m+\alpha_2}) \|v_{hn}\|(\mathbf{x}_h) = \int_{\Omega} \mathbf{f}^{m+\alpha_2} \cdot \mathbf{v}_h \, d\mathbf{x} \quad \text{for all } \mathbf{v}_h \in \mathbf{V}_h, \\
\left[ \|u_{hn}^{m+\alpha_2}\| (\mathbf{x}_h) \ge 0, \quad \lambda_h^{m+\alpha_2} \le 0, \quad \lambda_h^{m+\alpha_2} \left[ \|u_{hn}^{m+\alpha_2}\| (\mathbf{x}_h) = 0 \right] = 0 \quad \text{for } \mathbf{x}_h \in \Gamma_C^h,
\end{array}
\right.
\tag{3.3}
$$

where functions $\mathbf{v}_h$ associate a nodal basis in $\mathbf{V}_h$, and initial fields $\mathbf{u}_h^0 = \mathbf{u}_{0h}$, $\dot{\mathbf{u}}_h^0 = \dot{\mathbf{u}}_{0h}$. Note that the scheme is unconditionally stable for $\alpha_1 \ge \alpha_2 \ge 1/2$, a second order time integration scheme is obtained when $\gamma = 1/2 + \alpha_1 - \alpha_2$, and high-frequencies dissipation is ensured for $\beta = 1/4(1 + \alpha_1 - \alpha_2)^2$ (see [39]). Therefore, in our experiments we test the standard HHT-$\alpha$ scheme as $\alpha = (1, 0.9)$ and the generalized one as $\alpha = (1.1, 1)$, which both satisfy the above conditions with $\gamma = 0.6$ and $\beta = 0.3025$.

Reducing the unknown $\dot{\mathbf{u}}_h^{m+1}$, $\ddot{\mathbf{u}}_h^{m+1}$ by the identities:

$$
\left\{
\begin{aligned}
\dot{\mathbf{u}}_h^{m+1} &= \frac{\gamma}{\beta\tau}(\mathbf{u}_h^{m+1} - \mathbf{u}_h^m) + (1-\frac{\gamma}{\beta})\dot{\mathbf{u}}_h^m + \tau(1-\frac{\gamma}{2\beta})\ddot{\mathbf{u}}_h^m, \\
\ddot{\mathbf{u}}_h^{m+1} &= \frac{1}{\beta\tau^2}(\mathbf{u}_h^{m+1} - \mathbf{u}_h^m) - \frac{1}{\beta\tau}\dot{\mathbf{u}}_h^m + (1-\frac{1}{2\beta})\ddot{\mathbf{u}}_h^m
\end{aligned}
\right.
$$

from the implicit HHT scheme (3.3), it can be rewritten explicitly as follows:

$$
\left\{
\begin{array}{l}
\text{Find } (\mathbf{u}_h^{m+\alpha_2}, \lambda_h^{m+\alpha_2}, P^{m+\alpha_2}) \in \mathbf{V}_h \times \mathbb{R}^{N_C^h} \times \mathbb{R} \text{ such that:} \\
\mathbf{A}_\tau(\mathbf{u}_h^{m+\alpha_2}, \mathbf{v}_h) + \sum_{\Gamma_C^h} (\lambda_h^{m+\alpha_2} - P^{m+\alpha_2}) \|v_{hn}\|(\mathbf{x}_h) = \mathbf{F}_\tau^m(\mathbf{v}_h) \\
\text{for all } \mathbf{v}_h \in \mathbf{V}_h, \\
\qquad \int_{\Gamma_C} \|u_{hn}^{m+\alpha_2}\| d\Gamma = A^{m+\alpha_2}, \\
\qquad \|u_{hn}^{m+\alpha_2}\|(\mathbf{x}_h) \ge 0, \quad \lambda_h^{m+\alpha_2} \le 0, \quad \lambda_h^{m+\alpha_2} \|u_{hn}^{m+\alpha_2}\|(\mathbf{x}_h) = 0 \quad \text{for } \mathbf{x}_h \in \Gamma_C^h.
\end{array}
\right.
\tag{3.4}
$$

In (3.4) the bilinear operator $\mathbf{A}_\tau : \mathbf{V}_h \times \mathbf{V}_h \mapsto \mathbb{R}$ is given by

$$
\mathbf{A}_{\tau}(\mathbf{u}_{h}^{m+\alpha_2}, \mathbf{v}_{h}) = \int_{\Omega} \left( \frac{\rho \alpha_1}{\alpha_2 \beta \tau^2} \mathbf{u}_{h}^{m+\alpha_2} \cdot \mathbf{v}_{h} + \sigma(\mathbf{u}_{h}^{m+\alpha_2}) : \boldsymbol{\varepsilon}(\mathbf{v}_{h}) \right) d\mathbf{x}, \quad (3.5)
$$

and the linear operator $\mathbf{F}_\tau^m : \mathbf{V}_h \mapsto \mathbb{R}$ is defined as

$$
\mathbf{F}_{\tau}^{m}(\mathbf{v}_{h}) = \int_{\Omega} \left[ \mathbf{f}^{m+\alpha_2} + \frac{\rho \alpha_1}{\beta \tau^2} \left( \frac{1}{\alpha_2} \mathbf{u}_h^m + \tau \dot{\mathbf{u}}_h^m + \tau^2 \left( \frac{1}{2} - \frac{\beta}{\alpha_1} \right) \ddot{\mathbf{u}}_h^m \right) \right] \cdot \mathbf{v}_h \, d\mathbf{x}. \quad (3.6)
$$

The ellipticity and boundedness properties of the elasticity tensor in (2.1) yield the following uniform estimates. There exist constant $C_E, C_K, C_I > 0$ independent of the mesh size $h$ such that

$$
\|\boldsymbol{\varepsilon}(\mathbf{v}_h)\|_{0,\Omega} \le \|\nabla \mathbf{v}_h\|_{0,\Omega}, \quad \|\boldsymbol{\sigma}(\mathbf{v}_h)\|_{0,\Omega} \le C_E \|\nabla \mathbf{v}_h\|_{0,\Omega}, \tag{3.7}
$$

the Korn's and Poincaré's inequalities yield:

$$
\int_{\Omega} \sigma(\mathbf{v}_h) : \boldsymbol{\varepsilon}(\mathbf{v}_h) d\mathbf{x} \ge C_K \|\mathbf{v}_h\|_{1,\Omega}^2, \quad \|\mathbf{v}_h\|_{1,\Omega}^2 = \|\mathbf{v}_h\|_{0,\Omega}^2 + \|\nabla \mathbf{v}_h\|_{0,\Omega}^2, \tag{3.8}
$$

and for the quasi-uniform mesh $\mathcal{T}_h^i$ the inverse inequality holds:

$$
\|\mathbf{v}_h\|_{0,\Omega} \ge C_1 h \|\mathbf{v}_h\|_{1,\Omega}. \tag{3.9}
$$

With the help of (3.7)–(3.9) we prove the existence theorem.

**Theorem 3.1.** At each time step $m$, there exist one unique solution to the mixed variational problem (3.4) subjected to non-penetration and volume constraints.

**Proof.** For fixed $m$, let us define the feasible set:

$$
\mathbf{K}_h^m(\mathcal{A}) = \{ \mathbf{v} \in \mathbf{V}_h : \quad \llbracket v_{hn} \rrbracket \ge 0 \text{ on } \Gamma_C^h, \quad \int_{\Gamma_C} \llbracket v_{hn} \rrbracket d\Gamma = A^m \},
$$
which is convex and closed, and consider the variational inequality:

$$ egin{cases} 
\text{Find } \mathbf{u}_h^{m+a_2} \in \mathbf{K}_h^m(A) \text{ such that:} \\
\mathbf{A}_\tau(\mathbf{u}_h^{m+a_2}, \mathbf{v}_h - \mathbf{u}_h^{m+a_2}) \ge \mathbf{F}_\tau^m(\mathbf{v}_h - \mathbf{u}_h^{m+a_2}) & \text{for all } \mathbf{v}_h \in \mathbf{K}_h^m(A). 
\end{cases} 
\quad (3.10) $$

The operator $\mathbf{A}_\tau$ is bounded due to (3.7) and coercive owing to (3.8) and (3.9):

$$ \mathbf{A}_\tau(\mathbf{v}_h, \mathbf{v}_h) \ge \left( \frac{\rho \alpha_1 h^2}{\alpha_2 \beta \tau^2} C_1^2 + C_K \right) \|\mathbf{v}_h\|_{1,\Omega}^2. $$

The Lions–Stampachia theorem for variational inequalities justifies one unique solution to (3.10). It follows the mean pressure $P^{m+a_2} \in \mathbb{R}$ by formula (2.13):

$$ P^{m+a_2} := \frac{1}{\Lambda^m} (\mathbf{A}_\tau(\mathbf{u}_h^{m+a_2}, \mathbf{u}_h^{m+a_2}) - \mathbf{F}_\tau^m(\mathbf{u}_h^{m+a_2})), \quad (3.11) $$

and the Lagrange multiplier $\lambda_h^{m+a_2} \in \mathbb{R}^{N_C^h}$ using Green's formula (2.10) yields:

$$ \sum_{\Gamma_C^h} \lambda_h^{m+a_2} \llbracket \nu_{hn} \rrbracket(\mathbf{x}_h) := \mathbf{A}_\tau(\mathbf{u}_h^{m+a_2}, \mathbf{v}_h) - \mathbf{F}_\tau^m(\mathbf{v}_h) + P^{m+a_2} \sum_{\Gamma_C^h} \llbracket \nu_{hn} \rrbracket(\mathbf{x}_h), \quad (3.12) $$

for $\mathbf{v}_h \in \mathbf{V}_h$. Then (3.10)–(3.12) compose the solution triple to the problem (3.4). ☐

For solution of (3.4) we apply the penalty and semi-smooth Newton methods.

## 4. Semi-smooth Newton method of solution

For fixed parameters $r > 0$ and small $\delta > 0$ we penalize the problem and express the complementarity conditions using the minimum-based merit function:

$$ egin{cases} 
\text{Find } (\mathbf{u}_{h\delta}^{m+a_2}, \lambda_{h\delta}^{m+a_2}) \in \mathbf{V}_h \times \mathbb{R}^{N_C^h} \text{ such that:} \\
\mathbf{A}_\tau(\mathbf{u}_{h\delta}^{m+a_2}, \mathbf{v}_h) + \sum_{\Gamma_C^h} \left( \lambda_{h\delta}^{m+a_2} - \frac{1}{\delta} \left[ A^{m+a_2} - \int_{\Gamma_C} \llbracket \mathbf{u}_{h\delta n}^{m+a_2} \rrbracket d\Gamma \right] \right) \llbracket \nu_{hn} \rrbracket(\mathbf{x}_h) \\
= \mathbf{F}_\tau^m(\mathbf{v}_h) & \text{for all } \mathbf{v}_h \in \mathbf{V}_h, \\
\min(\llbracket \mathbf{u}_{h\delta n}^{m+a_2} \rrbracket, -r\lambda_{h\delta}^{m+a_2})(\mathbf{x}_h) = 0 & \text{for } \mathbf{x}_h \in \Gamma_C^h. 
\end{cases} 
\quad (4.1) $$

For the minimum function, a generalized gradient can be introduced by

$$ \nabla \min(\llbracket \mathbf{u}_{hn} \rrbracket, -r\lambda_h) = (\mathbf{1}_{A(\llbracket \mathbf{u}_{hn} \rrbracket, \lambda_h)}, -r\mathbf{1}_{I(\llbracket \mathbf{u}_{hn} \rrbracket, \lambda_h)}), \quad (4.2) $$

where the indicator function $1$ is used for the strictly active set:

$$ \mathcal{A}(\llbracket \mathbf{u}_{hn} \rrbracket, \lambda_h) = \{\mathbf{x}_h \in \Gamma_C^h : \quad (\llbracket \mathbf{u}_{hn} \rrbracket + r\lambda_h)(\mathbf{x}_h) < 0\}, \quad (4.3) $$

and for its complementary inactive set:

$$ I(\llbracket \mathbf{u}_{hn} \rrbracket, \lambda_h) = \{\mathbf{x}_h \in \Gamma_C^h : \quad (\llbracket \mathbf{u}_{hn} \rrbracket + r\lambda_h)(\mathbf{x}_h) \ge 0\}. \quad (4.4) $$

The function is semi-smooth (see [22]) on the finite set:

$$ egin{cases} 
\|\min(\llbracket \mathbf{u}_{hn} \rrbracket, -r\lambda_h) - \min(\llbracket \mathbf{v}_{hn} \rrbracket, -r\mu_h) \\ -\nabla \min(\llbracket \mathbf{u}_{hn} \rrbracket, -r\lambda_h) \begin{bmatrix} \llbracket \mathbf{u}_{hn} - \mathbf{v}_{hn} \rrbracket \\ \lambda_h - \mu_h \end{bmatrix} \|_{\infty, \Gamma_C^h} = 0 & (\llbracket \mathbf{u}_{hn} - \mathbf{v}_{hn} \rrbracket \ne 0) \\ 
\|\min(\llbracket \mathbf{u}_{hn} \rrbracket, -r\lambda_h) \begin{bmatrix} \llbracket \mathbf{u}_{hn} - \mathbf{v}_{hn} \rrbracket \\ \lambda_h - \mu_h \end{bmatrix} \|_{\infty, \Gamma_C^h} = 0 & (\llbracket \mathbf{u}_{hn} - \mathbf{v}_{hn} \rrbracket = 0) 
\end{cases} 
\quad (4.5) $$

Based on (4.2)–(4.5) we formulate the semi-smooth Newton method for iterative solution of the nonlinear system (4.1). Initializing with a guess $(\mathbf{u}_{h\delta}^{m+a_2,0}, \lambda_{h\delta}^{m+a_2,0}) \in \mathbf{V}_h \times \mathbb{R}^{N_C^h}$ at $k = 0$, at every $k \ge 0$ solve the linear system:
$$
\left\{
\begin{array}{l}
\text{Find } (\mathbf{u}_{h\delta}^{m+a_2,k+1}, \lambda_{h\delta}^{m+a_2,k+1}) \in \mathbf{V}_h \times \mathbb{R}^{N_C^h} \text{ such that:} \\
\mathbf{A}_\tau(\mathbf{u}_{h\delta}^{m+a_2,k+1}, \mathbf{v}_h) + \sum_{\Gamma_C^h} (\lambda_{h\delta}^{m+a_2,k+1} + \frac{1}{\delta} \int_{\Gamma_C} \|\mathbf{u}_{h\delta n}^{m+a_2,k+1}\| d\Gamma) \|\mathbf{v}_{hn}\|(\mathbf{x}_h) \\
= \mathbf{F}_\tau^m(\mathbf{v}_h) + \frac{1}{\delta} A^{m+a_2} \sum_{\Gamma_C^h} \|\mathbf{v}_{hn}\|(\mathbf{x}_h) \quad \text{for all } \mathbf{v}_h \in \mathbf{V}_h, \\
\|\mathbf{u}_{h\delta n}^{m+a_2,k+1}\|_A \mathcal{A}(\|\mathbf{u}_{h\delta n}^{m+a_2,k}\|, \lambda_{h\delta}^{m+a_2,k}) - r \lambda_{h\delta}^{m+a_2,k+1} \mathbf{I}_A(\|\mathbf{u}_{h\delta n}^{m+a_2,k}\|, \lambda_{h\delta}^{m+a_2,k}) = 0.
\end{array}
\right.
\tag{4.6}
$$

**Theorem 4.1.** At each time step *k*, the semi-smooth Newton iterate (4.6) has one unique solution. If the initial guess is chosen sufficiently close to the solution of problem (4.1), then the iteration sequence converges super-linearly by increasing *k*.

**Proof.** The unique solution to the linear system (4.6) is a consequence of Theorem 3.1.

Let us gather the terms in the following increment:

$$
\begin{align*}
I := \|\mathbf{u}_{h\delta n}^{m+a_2,k} - \mathbf{u}_{h\delta n}^{m+a_2}\|_A & \mathcal{A}(\|\mathbf{u}_{h\delta n}^{m+a_2,k}\|, \lambda_{h\delta}^{m+a_2,k}) - r(\lambda_{h\delta}^{m+a_2,k} - \lambda_{h\delta}^{m+a_2}) \mathbf{I}_A(\|\mathbf{u}_{h\delta n}^{m+a_2,k}\|, \lambda_{h\delta}^{m+a_2,k}) \\
& = -\min(\|\mathbf{u}_{h\delta n}^{m+a_2,k}\| + r\lambda_{h\delta}^{m+a_2,k}, 0) + \min(\|\mathbf{u}_{h\delta n}^{m+a_2}\| + r\lambda_{h\delta}^{m+a_2}, 0) \\
& \qquad + (\|\mathbf{u}_{h\delta n}^{m+a_2,k} - \mathbf{u}_{h\delta n}^{m+a_2}\| + r(\lambda_{h\delta}^{m+a_2,k} - \lambda_{h\delta}^{m+a_2})) \mathbf{I}_A(\|\mathbf{u}_{h\delta n}^{m+a_2,k}\|, \lambda_{h\delta}^{m+a_2,k}).
\end{align*}
$$

In virtue of (4.5) it follows the asymptotic estimate:

$$
\|I\|_{\infty} = o\left( \left\| \|u_{h\delta n}^{m+a_2,k} - u_{h\delta n}^{m+a_2}\| + r(\lambda_{h\delta}^{m+a_2,k} - \lambda_{h\delta}^{m+a_2}) \right\|_{\infty} \right). \quad (4.7)
$$

Subtracting the variational equations (4.1) and (4.6) we get:

$$
\mathbf{A}_\tau(\mathbf{u}_{h\delta}^{m+a_2,k+1} - \mathbf{u}_{h\delta}^{m+a_2}, \mathbf{v}_h) + \frac{1}{\delta} \int_{\Gamma_C} \|\mathbf{u}_{h\delta n}^{m+a_2,k+1} - \mathbf{u}_{h\delta n}^{m+a_2}\| d\Gamma \sum_{\Gamma_C^h} \|\mathbf{v}_{hn}\|(\mathbf{x}_h) = - \sum_{\Gamma_C^h} (\lambda_{h\delta}^{m+a_2,k+1} - \lambda_{h\delta}^{m+a_2}) \|\mathbf{v}_{hn}\|(\mathbf{x}_h) \quad \text{for all } \mathbf{v}_h \in \mathbf{V}_h. \tag{4.8}
$$

For $\mathbf{v}_h$ such that $\|\mathbf{v}_h\|_{1,\Omega} \le C_1 \|\|\mathbf{v}_{hn}\|\|_{\infty}$ with $C_1 > 0$, using the bounds (3.7) and $C_2 \|\|\mathbf{v}_{hn}\|\|_{0,\Gamma_C} \le \|\|\mathbf{v}_{hn}\|\|_{\infty}$ with $C_2 > 0$ we estimate
from (4.8) the dual norm:

$$
\| \lambda_{h\delta}^{m+a_2,k+1} - \lambda_{h\delta}^{m+a_2} \|_{\infty} \le \left( \frac{\rho\alpha_1}{\alpha_2\beta\tau^2} + C_E \right) C_1 \| \mathbf{u}_{h\delta}^{m+a_2,k+1} - \mathbf{u}_{h\delta}^{m+a} \|_{1,\Omega} + \frac{1}{C_2\delta} \| \mathbf{u}_{h\delta n}^{m+a_2,k+1} - \mathbf{u}_{h\delta n}^{m+a_2} \|_{\infty}. \quad (4.9)
$$

Testing (4.8) with $\mathbf{v}_h = \mathbf{u}_{h\delta}^{m+a_2,k+1} - \mathbf{u}_{h\delta}^{m+a_2}$ and using the bounds (3.8), (3.9) gives

$$
\begin{align}
& \left( \frac{\rho \alpha_1 h^2}{\alpha_2 \beta \tau^2} C_1^2 + C_K \right) \| \mathbf{u}_{h\delta}^{m+a_2,k+1} - \mathbf{u}_{h\delta}^{m+a_2} \|_{1,\Omega}^2 + \frac{C_2}{\delta} \| \| \mathbf{u}_{h\delta n}^{m+a_2,k+1} - \mathbf{u}_{h\delta n}^{m+a_2} \| \|_{0,\Gamma_C}^2 \nonumber \\
& \le \| I \|_{\infty} \left( \left\| \frac{1}{\Gamma} \| \mathbf{u}_{h\delta n}^{m+a_2,k} - \mathbf{u}_{h\delta n}^{m+a_2} \| \right\|_{\infty} + \| \lambda_{h\delta}^{m+a_2,k} - \lambda_{h\delta}^{m+a_2} \|_{\infty} \right). \tag{4.10}
\end{align}
$$

The estimates (4.8)–(4.10) yield the super-linear convergence for the iterates, when the initial guess $(\mathbf{u}_{h\delta}^{m+a_2,0}, \lambda_{h\delta}^{m+a_2,0})$ is sufficiently close to $(\mathbf{u}_{h\delta}^{m+a_2}, \lambda_{h\delta}^{m+a_2})$. $\square$

We realize the semi-smooth Newton iteration (4.6) as a primal-dual active set (PDAS) iterates of the active (4.3) and complementary inactive (4.4) sets.

**Algorithm 4.1.** PDAS

1. Choose some initial guess $\mathcal{A}_{h\delta}^{m+1,-1} \subset \Gamma_C^h$, for instance $\mathcal{A}_{h\delta}^{m+1,-1} = \emptyset$, and its complement $I_{h\delta}^{m+1,-1}$.

2. At every iterate $k \ge -1$ solve successively the linear system:

$$
\left\{
\begin{array}{l}
\text{Find } (\mathbf{u}_{h\delta}^{m+a_2,k+1}, \lambda_{h\delta}^{m+a_2,k+1}) \in \mathbf{V}_h \times \mathbb{R}^{N_C^h} \text{ such that:} \\
\mathbf{A}_\tau(\mathbf{u}_{h\delta}^{m+a_2,k+1}, \mathbf{v}_h) + \sum_{\Gamma_C^h} (\lambda_{h\delta}^{m+a_2,k+1} + \frac{1}{\delta} \int_{\Gamma_C} \|\mathbf{u}_{h\delta n}^{m+a_2,k+1}\| d\Gamma) \|\mathbf{v}_{hn}\|(\mathbf{x}_h) \\
= \mathbf{F}_\tau^m(\mathbf{v}_h) + \frac{1}{\delta} \mathcal{A}^{m+a_2} \sum_{\Gamma_C^h} \|\mathbf{v}_{hn}\|(\mathbf{x}_h) \quad \text{for all } \mathbf{v}_h \in \mathbf{V}_h, \\
\|\mathbf{u}_{h\delta n}^{m+a_2,k+1}\| = 0 \quad \text{on } \mathcal{A}_{h\delta}^{m+1,k}, \quad \lambda_{h\delta}^{m+a_2,k+1} = 0 \quad \text{on } I_{h\delta}^{m+1,k}.
\end{array}
\right.
\tag{4.11}
$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Displacement field $\mathbf{x}_h + \mathbf{u}_h^m$ for $m = 0, \dots, 8$ within extension-compression loop.</figcaption></figure>

3: Compute the complementary active and inactive sets:

$$ egin{cases} 
\mathcal{A}_{h\delta}^{m+1,k+1} = \{\mathbf{x}_h \in \Gamma_C^h : (\|\mathbf{u}_{h\delta n}^{m+a_2,k+1}\| + r\lambda_{h\delta}^{m+a_2,k+1})(\mathbf{x}_h) < 0\}, \\ 
\mathcal{I}_{h\delta}^{m+1,k+1} = \{\mathbf{x}_h \in \Gamma_C^h : (\|\mathbf{u}_{h\delta n}^{m+a_2,k+1}\| + r\lambda_{h\delta}^{m+a_2,k+1})(\mathbf{x}_h) \ge 0\}. 
\end{cases} 
\quad (4.12) $$

4: If $\mathcal{A}_{h\delta}^{m+1,k+1} = \mathcal{A}_{h\delta}^{m+1,k}$ and $\|\mathbf{u}_{h\delta}^{m+a_2,k+1} - \mathbf{u}_{h\delta}^{m+a_2,k}\|_{\infty} < \text{tol}$ then stop.

The first stopping condition in the step 4 of Algorithm 4.1 guarantees fulfillment of the complementarity conditions at $\Gamma_C^h$, whereas the second condition with a prescribed tolerance $\text{tol} > 0$ arises due to penalization of the volume constraint. In computer simulation $\text{tol} = 10^{-10}$ and the penalty parameter $\delta = 10^{-10}$ are set. In the next section we validate numerically properties of the convergence stated in Theorem 4.1.

## 5. Numerical experiment

For computer simulation, we choose the 2D geometry (given in [m]) drawn in Fig. 1:

$$ egin{cases} 
\Omega^1 = \{\mathbf{x} \in (0, 2.5) \times (1, 2)\}, & \Omega^2 = \{\mathbf{x} \in (0, 2.5) \times (0, 1)\}, \\ 
\Gamma_D = \{x_1 \in [0, 2.5], x_2 = 0, 2\} \cup \{x_1 = 0, 2.5, x_2 \in (0, 2)\}, \\ 
\Gamma_C = \{x_1 \in (0, 2.5), x_2 = 1\}. 
\end{cases} $$

For the isotropic solid, the material density $\rho = 2700$ [kg/m³], Young's modulus $E = 73$ [GPa] and Poisson ratio $\nu = 0.34$ yield the shear modulus $\mu = E/(2(1+\nu)) \approx 27$ [GPa]. We test the elastic body which is uniformly compressed with the force $\mathbf{f} = (0, -50)$ [kN] in $\Omega^1$, and $\mathbf{f} = (0, 50)$ [kN] in $\Omega^2$, such that the crack $\Gamma_C$ is closed.

The fluid volume $A$ [m²] pumped in the crack is prescribed to be linear:

$$ A(t) = \frac{0.5}{T} t \quad \text{for } t \in [0, T/2), \quad A(t) = \frac{0.5}{T} (T - t) \quad \text{for } t \in [T/2, T], $$

within extension and compression, respectively. The initial velocity $\dot{\mathbf{u}}_0 = \mathbf{0}$ and acceleration $\ddot{\mathbf{u}}_0 = \mathbf{0}$. We utilize the standard piecewise $\mathbb{P}_1$-polynomial FEM for $\mathbf{V}_h$. For illustration, in Fig. 2 the extension-compression loop during $T = 2.5$ [s] is shown at times $t^m = m\tau$, $m = 0, \dots, 8$ for grid points $\mathbf{x}_h + \mathbf{u}_h^m$ in the current configuration $\mathbf{x}_h \in \Omega$ with the mesh-size $h = 0.25$ and step-site $\tau = 0.3125$.

At each time step, the system matrix of the problem has the number of rows and columns which corresponds to the number of degrees of freedom (DOF). The DOF is proportional to the vector dimension, that is two in 2D, and inverse proportional to the mesh-size squared. For instance, for $h = 1/4$ in Fig. 2 the DOF number is 88. For isotropic elastic equations the system matrix has a sparse structure and forms both a P-matrix and a M-matrix. We find numerical solution to the linear system using a direct solver in Matlab that is the backslash command.

To examine stability of the dynamic contact problem, we consider the energy:

$$ E_{h\delta}^m = \frac{1}{2} \int_{\Omega} (\rho |\dot{\mathbf{u}}_{h\delta}^m|^2 + \sigma(\mathbf{u}_{h\delta}^m) : \varepsilon(\mathbf{u}_{h\delta}^m)) d\mathbf{x}, \quad m = 0, \dots, N. \quad (5.1) $$

The discrete energy $E_{h\delta}^m$ [kJ] in (5.1) computed by the Cranck–Nicolson (CN) scheme [$\gamma = 0.5, \beta = 0.25, \alpha = (1, 1)$] is compared with the standard HHT scheme [$\gamma = 0.6, \beta = 0.3025, \alpha = (1, 0.9)$] and the generalized HHT scheme [$\gamma = 0.6, \beta = 0.3025, \alpha = (1.1, 1)$]. In the corresponding plots (a)–(c) of Fig. 3, three different curves are related to different time steps $\tau \in \{1/25, 1/20, 1/15\}$ at fixed $h = 1/20$ (that is 1818 DOF). In Fig. 3 it can be seen that CN becomes unstable in the plot (a), spurious oscillations in the energy are reduced by the $\alpha_2$-damping in the plot (b), whereas the oscillations are much better suppressed by the $\alpha_1$-damping in the plot (c).

For solution we realize the PDAS Algorithm 4.1, where $r = 1$. Starting with $\mathcal{A}^{m+1,-1} = \emptyset$, the $k$-iterates at the time step $m = 89$ with $\tau = 1/72$ are demonstrated in Fig. 4. The spatial system has 10796 DOF for $h = 1/50$. The crack opening $\|\mathbf{u}_{h\delta n}^{m+1,k}\|$, Lagrange
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 3. The discrete energy $E_{h\delta}^m$ at times $t^m$ for selected $[\gamma, \beta, \alpha]$ when decreasing $\tau$.</figcaption></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 4. PDAS iterates k: crack opening $\llbracket u_{h\delta n}^{m+1,k} \rrbracket$, contact force $\lambda_{h\delta}^{m+1,k}$, active set $\mathcal{A}_{h\delta}^{m+1,k}$.</figcaption></figure>

multiplier $\lambda_{h\delta}^{m+1,k}$, and active set $\mathcal{A}_{h\delta}^{m+1,k}$ are depicted along $N_C^h = 90$ points at the discrete crack $\Gamma_C^h$. The algorithm converges in only 6 iterations. From the plots (a) and (c) of Fig. 4 we report monotone properties of the convergence for $\llbracket u_{h\delta n}^{m+1,k} \rrbracket$ and $\mathcal{A}_{h\delta}^{m+1,k}$, but not for $\lambda_{h\delta}^{m+1,k}$ in the plot (b).

For the problem stated under the uniform force $\mathbf{f} = (f_1, f_2)$, one can observe in Fig. 2 and Fig. 4 that the fracture opening occurs first time near the interface ends $x_1 = 0$ and $x_1 = 2.5$. In another test problem we impose the inhomogeneous vertical force $-f_2(\mathbf{x})$ in $\Omega^1$ and $f_2(\mathbf{x})$ in $\Omega^2$, where

$$f_2(\mathbf{x}) = (1 - x_1/1.25)^2. \qquad (5.2)$$

The corresponding extension-compression loop $\mathbf{x}_h + \mathbf{u}_h^m$ over the grid points $\mathbf{x}_h \in \Omega$ as $h = 0.25$ and $\tau = 0.3125$ is drawn in Fig. 5.

In Fig. 6 there are shown the crack opening $\llbracket u_{h\delta n}^{m+1,k} \rrbracket$, Lagrange multiplier $\lambda_{h\delta}^{m+1,k}$, and active set $\mathcal{A}_{h\delta}^{m+1,k}$ at the same parameters $m = 89$, $\tau = 1/72$ and $h = 1/50$, for $k = 0, \dots, 6$ iterates under the inhomogeneous vertical force from (5.2). Here the last three iterates are required to attain the tolerance with respect to the penalty approximation. In comparison with the previous test, in Fig. 6 the fracture opening occurs first time around the interface center $x_1 = 1.25$. For the relevant physical foundation by cohesion forces imposed at the interface we refer to [40,41].

### CRediT authorship contribution statement

Victor A. Kovtunenko: Investigation. Yves Renard: Conceptualization.

### Acknowledgement

V.A.K. is supported by the OeAD Scientific & Technological Cooperation (MULT 06/2023) financed by the Austrian Federal Ministry of Science, Research and Economy (BMFWF); the author acknowledges the financial support by the University of Graz. Y.R. is supported by the Ministry for Europe and Foreign Affairs (France) within Programme Danube 2023-2025 (Project PHC Danube n°49916UM).
<figure><img src="image_11.png" /><figcaption>Fig. 5. Displacement field $x_h + u_h^m$ for $m = 0, \dots, 8$ under inhomogeneous force.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 6. PDAS iterates $k: [u_{h\delta}^{m+1,k}], \lambda_{h\delta}^{m+1,k}$, and $A_{h\delta}^{m+1,k}$ under inhomogeneous force.</figcaption></figure>

References

[1] A.P. Bunger, E. Detournay, D.I. Garagash, Toughness-dominated hydraulic fracture with leak-off, Int. J. Fract. 134 (2005) 175–190, https://doi.org/10.1007/s10704-005-0154-0.

[2] C. Giovannini, L. Fumagalli, F.S. Patacchini, Predicting nonlinear-flow regions in highly heterogeneous porous media using adaptive constitutive laws and neural networks, J. Comput. Phys. (2025) 114093, https://doi.org/10.1016/j.jcp.2025.114093.

[3] A. Mikelić, M.F. Wheeler, T. Wickl, Phase-field modeling of a fluid-driven fracture in a poroelastic medium, Comput. Geosci. 19 (2015) 1171–1195, https://doi.org/10.1007/s10596-015-9532-5.

[4] V.V. Shelukhin, V.A. Baikov, S.V. Golovin, A.Y. Davletbaev, V.N. Starovoitov, Fractured water injection wells: pressure transient analysis, Int. J. Solids Struct. 51 (2014) 2116–2122, https://doi.org/10.1016/j.ijsolstr.2014.02.019.

[5] M.A. Efendiev, A. Mielke, On the rate-independent limit of systems with dry friction and small viscosity, J. Convex Anal. 13 (2006) 151–167.

[6] S. Almi, G.D. Maso, R. Toader, Quasi-static crack growth in hydraulic fracture, Nonlinear Anal. Theory Meth. Appl. 109 (2014) 301–318, https://doi.org/10.1016/j.na.2014.07.009.

[7] V.A. Kovtunenko, Poroelastic medium with non-penetrating crack driven by hydraulic fracture: variational inequality and its semidiscretization, J. Comput. Appl. Math. 405 (2022) 113953, https://doi.org/10.1016/j.cam.2021.113953.

[8] M. Bach, V.A. Kovtunenko, A.M. Khludnev, Derivatives of the energy functional for 2D-problems with a crack under Signorini and friction conditions, Math. Methods Appl. Sci. 23 (2000) 515–534, https://doi.org/10.1002/(SICI)1099-1476(200004)23:6<515::AID-MMA122>3.0.CO;2-S.

[9] A.M. Khludnev, J. Sokolowski, Modeling and Control in Solid Mechanics, Birkhäuser, Basel, 1997.

[10] A.M. Khludnev, V.A. Kovtunenko, Analysis of Cracks in Solids, WIT-Press, Southampton, Boston, 2000.

[11] A.M. Khludnev, V.A. Kovtunenko, A. Tani, Evolution of a crack with kink and non-penetration, J. Math. Soc. Jpn. 60 (2008) 1219–1253, https://doi.org/10.2969/jmsj/06041219.

[12] V.A. Kovtunenko, Numerical simulation of the non-linear crack problem with non-penetration, Math. Methods Appl. Sci. 27 (2004) 163–179, https://doi.org/10.1002/mma.449.

[13] T.A. Laursen, Computational Contact and Impact Mechanics, Springer, Berlin, 2003.

[14] O. Chau, A. Heibig, A. Petrov, A class of dynamic unilateral contact problems with sub-differential friction law, in: N. Daras, M. Rassias, N. Zographopoulos (Eds.), Exploring Mathematical Analysis, Approximation Theory, and Optimization, Springer, Cham, 2023, pp. 17–31.

[15] C. Eck, J. Jarušek, M. Krbec, Unilateral Contact Problems: Variational Methods and Existence Theorems, CRC Press, Boca Raton, 2005.

[16] A. Petrov, M. Schatzman, Mathematical results on existence for viscoelastodynamic problems with unilateral constraints, SIAM J. Math. Anal. 40 (2009) 1882–1904, https://doi.org/10.1137/070695101.

[17] T. Kashiwabara, H. Itou, Unique solvability of a crack problem with Signorini-type and Tresca friction conditions in a linearized elastodynamic body, Phil. Trans. R. Soc. A 380 (2022) 20220225, https://doi.org/10.1098/rsta.2022.0225.

[18] S. Migórski, A. Ochal, M. Sofonea, Nonlinear Inclusions and Hemivariational Inequalities: Models and Analysis of Contact Problems, Springer, New York, 2013.

[19] S. Zeng, A.A. Khan, S. Migórski, A new class of generalized quasi-variational inequalities with applications to Oseen problems under nonsmooth boundary conditions, Sci. China, Math. 67 (2024) 315–333, https://doi.org/10.1007/s11425-022-2069-0.

[20] J. Gwinner, B. Jadamba, A.A. Khan, F. Raciti, Uncertainty Quantification in Variational Inequalities: Theory, Numerics, and Applications, Chapman and Hall/CRC, New York, 2021.
[21] N. Aguilera, H.W. Alt, L.A. Caffarelli, An optimization problem with volume constraint, SIAM J. Control Optim. 24 (1986) 191–198, https://doi.org/10.1137/0324011.

[22] M. Bergounioux, K. Ito, K. Kunisch, Primal-dual strategy for constrained optimal control problems, SIAM J. Control Optim. 37 (1999) 1176–1194, https://doi.org/10.1137/S0363012997328609.

[23] R. Kučera, K. Motyčková, A. Markopoulos, J. Haslinger, On the inexact symmetrized globally convergent semi-smooth Newton method for 3D contact problems with Tresca friction: the R-linear convergence rate, Optim. Methods Softw. 35 (2019) 65–86, https://doi.org/10.1080/10556788.2018.1556659.

[24] Y. Renard, Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity, Comput. Methods Appl. Mech. Eng. 256 (2013) 38–55, https://doi.org/10.1016/j.cma.2012.12.008.

[25] S. Hüeber, G. Stadler, B.I. Wohlmuth, A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction, SIAM J. Sci. Comput. 30 (2008) 572–596, https://doi.org/10.1137/060671061.

[26] S. Abide, M. Barboteu, S. Cherkaoui, S. Dumont, Unified primal-dual active set method for dynamic frictional contact problems, Fixed Point Theory Algorithms Sci. Eng. 2022 (2022) 19, https://doi.org/10.1186/s13663-022-00729-4.

[27] M. Hintermüller, V.A. Kovtunenko, K. Kunisch, Generalized Newton methods for crack problems with non-penetration condition, Numer. Methods Partial Differ. Equ. 21 (2005) 586–610, https://doi.org/10.1002/num.20053.

[28] M. Hintermüller, V.A. Kovtunenko, K. Kunisch, A Papkovich-Neuber-based numerical approach to cracks with contact in 3D, IMA J. Appl. Math. 74 (2009) 325–343, https://doi.org/10.1093/imamat/hxp017.

[29] V.A. Kovtunenko, Y. Renard, Convergence analysis of semi-smooth Newton method for mixed FEM approximations of dynamic two-body contact and crack problems, J. Comput. Appl. Math. 471 (2025) 116722, https://doi.org/10.1016/j.cam.2025.116722.

[30] F. Chouly, P. Hild, Y. Renard, Finite Element Approximation of Contact and Friction in Elasticity, Birkhäuser, Cham, 2023.

[31] V. Linck, G. Bayada, L. Baillet, T. Sassi, J. Sabil, Finite element analysis of a contact with friction between an elastic body and a thin soft layer, J. Tribol. 127 (2005) 461–468, https://doi.org/10.1115/1.1866170.

[32] L. Paoli, M. Schatzman, A numerical scheme for impact problems II: the multidimensional case, SIAM J. Numer. Anal. 40 (2002) 734–768, https://doi.org/10.1137/S003614290037873X.

[33] H.B. Khenous, P. Laborde, Y. Renard, Mass redistribution method for finite element contact problems in elastodynamics, Eur. J. Mech. A Solids 27 (2008) 918–932, https://doi.org/10.1016/j.euromechsol.2008.01.001.

[34] F. Dabaghi, P. Krejčí, A. Petrov, J. Pousin, Y. Renard, A weighted finite element mass redistribution method for dynamic contact problems, J. Comput. Appl. Math. 345 (2019) 338–356, https://doi.org/10.1016/j.cam.2018.06.030.

[35] H.M. Hilber, T.J.R. Hughes, R.L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthq. Eng. Struct. Dyn. 5 (1977) 283–292, https://doi.org/10.1002/eqe.4290050306.

[36] P. Hauret, P.L. Tallec, Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact, Comput. Methods Appl. Mech. Eng. 195 (2006) 4890–4916, https://doi.org/10.1016/j.cma.2005.11.005.

[37] H. Huang, N. Pignet, G. Drouet, F. Chouly, HHT-$\alpha$ and TR-BDF2 schemes for dynamic contact problems, Comput. Mech. 73 (2024) 1165–1186, https://doi.org/10.1007/s00466-023-02405-9.

[38] T. Laursen, V. Chawla, Semismooth and semiconvex functions in constrained optimization, Int. J. Numer. Methods Eng. 40 (1997) 863–886, https://doi.org/10.1002/(SICI)1097-0207(19970315)40:53.0.CO;2-V.

[39] J. Chung, G.M. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method, J. Appl. Mech. 60 (1993) 371–375, https://doi.org/10.1115/1.2900803.

[40] H. Itou, V.A. Kovtunenko, N.P. Lazarev, Poroelastic problem of a non-penetrating crack with cohesive contact for fluid-driven fracture, Appl. Eng. Sci. 15 (2023) 100136, https://doi.org/10.1016/j.apples.2022.100136.

[41] H. Itou, V.A. Kovtunenko, K.R. Rajagopal, Poroelastic problem of a non-penetrating crack with cohesive contact for fluid-driven fracture, Math. Models Methods Appl. Sci. 35 (2025), https://doi.org/10.1142/S021820252550040X.