# Coupling surrogate neural networks and classical numerical solvers for non-smooth mechanical dynamics with contact

Mohammad Issa, Xavier Merlhiot, Anders Thorin

► To cite this version:

Mohammad Issa, Xavier Merlhiot, Anders Thorin. Coupling surrogate neural networks and classical numerical solvers for non-smooth mechanical dynamics with contact. 2025. hal-05385069

HAL Id: hal-05385069

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

Preprint submitted on 27 Nov 2025

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 Attribution 4.0 International License
# Coupling surrogate neural networks and classical numerical solvers for non-smooth mechanical dynamics with contact

Mohammad ISSA*†

Xavier Merlhiot*

Anders Thorin*

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Hybrid formulation and time-stepping framework</td><td>4</td></tr><tr><td>2.1</td><td>Nonsmooth dynamics and Moreau-Jean scheme</td><td>4</td></tr><tr><td>2.2</td><td>Neural network surrogate model</td><td>6</td></tr><tr><td>2.3</td><td>Coupled system and surrogate integration</td><td>6</td></tr><tr><td>2.3.1</td><td>Moreau-Jean scheme for the coupled system</td><td>7</td></tr><tr><td>2.3.2</td><td>Hybrid MJ formulation as an NCP</td><td>8</td></tr><tr><td>2.3.3</td><td>SOR – Prox hybrid iteration</td><td>9</td></tr><tr><td>2.4</td><td>Position-level integrator (XDE)</td><td>11</td></tr><tr><td>2.5</td><td>Algorithmic implementation and convergence</td><td>12</td></tr><tr><td>3</td><td>Numerical results</td><td>13</td></tr><tr><td>3.1</td><td>Harmonic oscillator-mass contact</td><td>14</td></tr><tr><td>3.2</td><td>Beam-mass impact simulation</td><td>15</td></tr><tr><td>4</td><td>Conclusion</td><td>18</td></tr></table>

<footnote>*Université Paris-Saclay, CEA, List, F-91120, Palaiseau, France</footnote>

<footnote>†Corresponding author: mohammad.issa2@cea.fr</footnote>
Abstract

Machine learning surrogates enable data-driven modeling of mechanical systems with nonlinear and non-smooth dynamics, yet their use typically remains restricted to smooth regimes without intermittent contact. We introduce a hybrid framework that couples a recurrent neural surrogate with a classical non-smooth contact solver for multibody systems with intermittent contact. The surrogate replaces the flexible subdomain while remaining fully coupled through the standard contact resolution loop. Within each implicit time-stepping iteration, the network is queried to provide the flexible body's kinematic response, its contribution enters the contact conditions without any matrix assembly or local matrix inversion for the flexible block, so that the solver structure and complementarity formulation are preserved.

The framework is exemplified using surrogates as gated recurrent unit (GRU) that maps applied contact forces and a hidden state to the next configuration and velocity, trained offline on physically consistent trajectories. Two variants of the coupled scheme are considered: the velocity-level Moreau–Jean (MJ) integrator and an in-house position-level integrator (XDE).

The framework is validated on a harmonic oscillator–mass contact and a beam–mass impact, where the beam is replaced by the surrogate. In both cases, the hybrid solver reproduces the nonlinear contact responses with accuracy comparable to full finite element simulations, demonstrating that the surrogate can be inserted non-intrusively inside the contact solver without degrading physical consistency. This extends the use of surrogates to configurations involving intermittent contact.

**Keywords:** Non-smooth dynamics, unilateral contact, nonlinear complementarity problem (NCP), surrogate modeling, recurrent neural network

# 1 Introduction

The simulation of mechanical systems subjected to unilateral constraints has long been a central challenge in computational mechanics. Non-smooth behavior such as contact, impact, and friction gives rise to complex dynamics that require efficient and robust numerical methods.

A foundational contribution was made by Moreau [10, 11], who introduced differential inclusions to describe systems with unilateral constraints and dry friction along a first nonsmooth time-stepping scheme. A generalization of this scheme by Jean [8] led to the so-called Moreau–Jean time-stepping scheme which is widely recognized for handling impacts and frictional contacts among systems of rigid and deformable bodies. These two schemes both handle the unilateral constraints associated with the normal direction of contact at the velocity level, with the advantage that impact laws (such as Newton's impact law) can be easily supported, but with the drawback that the unilateral constraints in position may not be strictly respected (the approximation being of first order with respect to the time step [1]). To tackle this problem, several modified variants of the Moreau-Jean scheme have been proposed ([1, 13, 15]), one of the central ideas being to rely on an adaptation of the GLL (Gear-Gupta-Leimkuhler) stabilization method [6] to unilateral constraints in the context of a first-order nonsmooth time stepping scheme. Alternatively, if one accepts to support only perfectly plastic impact laws, it is possible to
resort to the simpler approach outlined in [15] and referred to by Studer as *time-stepping by preconditioning*.

Despite their robustness, time-stepping schemes coupled with finite element (FE) discretizations [16] can be computationally demanding for flexible subsystems with fine meshes, due to large numbers of degrees of freedom and possibly small admissible time steps. This increases the computational cost per time step and can make large-scale simulations impractical.

Machine learning provides an attractive alternative for such subdomains. Neural Networks (NNs) have demonstrated strong capability to approximate nonlinear systems efficiently while offering constant-time inference [7]. NNs have been successfully employed as surrogates for complex mechanical responses, for instance in deep reduced-order models for nonlinear microstructure dynamics [5], graph-based neural models for parametrized PDEs [4], and recurrent surrogates for time-dependent systems [3]. These advances motivate the integration of learned surrogates into classical time-stepping frameworks for non-smooth dynamics.

<figure><img src="image_2.png" /><figcaption>Figure 1: Controlled robot interacting with a soft body. The present contribution offers a method to couple efficient time-integration methods for the robotic arm with a neural-network surrogate for the soft body.</figcaption></figure>

In many engineering assemblies, only a subset of components exhibits complex constitutive behavior, while others can be considered perfectly rigid. In these cases, the flexible subdomain dominates the computational cost, motivating the replacement of its FE model by a surrogate while retaining full compatibility with classical non-smooth contact solvers.

We therefore approximate the flexible subdomain response using a neural network trained offline on physically consistent trajectories. During simulation, the surrogate is queried within the contact-resolution loop of an event-capturing nonsmooth time-stepping scheme. Crucially, the surrogate contribution enters the complementarity conditions associated to contact without assembling stiffness or mass matrices for the flexible block, and without local matrix inversion. The overall solver (Gauss-Seidel/SOR-Prox) structure and NCP formulation remain unchanged, enabling a non-intrusive replacement of the flexible FE subdomain while preserving coupling with the rest of the system.

In the present implementation, the surrogate is instantiated as a gated recurrent unit (GRU) network that maps contact forces and a hidden state to the next configuration and velocity. The hybrid methodology is demonstrated using two integrator variants: the velocity-level Moreau-Jean (MJ) and the position-level XDE schemes.

In this work, our objective is to introduce a hybrid framework that couples a neural
network surrogate for the computationally expensive part with a classical multibody system through intermittent contact interactions governed by nonsmooth contact laws, as illustrated in Figure 1. The approach preserves the essential features of non-smooth dynamics—event capturing, unilateral constraints, and implicit time stepping—while eliminating the need to assemble or invert the matrices associated with the flexible block. The methodology is validated on two benchmarks (harmonic oscillator–mass and beam–mass impact), demonstrating accurate reproduction of nonlinear contact responses compared to full FE simulations. This demonstrates that machine-learning-based reduction techniques can be reliably extended to dynamic scenarios involving strong nonlinearities and intermittent impacts.

The paper is organized as follows. Section 2.1 introduces the theoretical background of non-smooth dynamics and recalls the Moreau–Jean scheme, which forms the basis of the proposed coupling strategy. Section 2.3 presents the hybrid formulation combining the neural-network surrogate with the classical time-stepping framework, including the SOR–Prox contact solver. Section 2.4 introduces the variant of the Moreau–Jean scheme that is implemented in the XDE Physics software [9] and which relies on a treatment of unilateral contact constraints at position level similar to the approach proposed by Studer in [15]. Algorithmic details of the coupled solver and convergence refinements, including the adaptive relaxation strategy, are summarized in Section 2.5. Section 3 reports numerical experiments.

## 2 Hybrid formulation and time-stepping framework

We recall here the classical formulation of nonsmooth mechanical systems with unilateral contact and its discretization via the Moreau–Jean scheme, which will form the basis of the coupled surrogate–solver framework introduced later.

### 2.1 Nonsmooth dynamics and Moreau–Jean scheme

Consider a mechanical system whose configuration is represented by the generalized coordinates $x \in \mathbb{R}^d$ and velocity $v = \dot{x}$. The dynamics of a system subject to perfect unilateral constraints can be expressed in strong form [2, 11] as

$$m(x)\dot{v} = f(x, v, t) + J_c^\top(x)\lambda, \quad (1a)$$

$$v = \dot{x}, \quad (1b)$$

$$0 \le g(x) \perp \lambda \ge 0, \quad (1c)$$

where $g = [g_1, \dots, g_p]^\top$ groups the $p$ constraint functions, $J_c(x) = \frac{\partial g}{\partial x}(x)$ is the constraint Jacobian, and the multiplier $\lambda(t) \in \mathbb{R}^p$ represents the contact reaction forces associated with the active constraints. The complementarity condition (1c) means that for all $l \in \{1, \dots, p\}$, we have

$$\left\{ \begin{array}{l} g_l(x(t)) \ge 0, \\ \lambda_l(t) \ge 0, \\ g_l(x(t)) \cdot \lambda_l(t) = 0. \end{array} \right.$$
109 **Weak measure-based formulation.** Impacts introduce discontinuities in velocity, preventing a classical derivative. Following Moreau [11, 14], one adopts a weak measure-based formulation, leading to the second-order Moreau process, expressed here with Newton's impact law:

$$m(x) dv = f(x, v^+, t) dt + J_c^\top(x) dr, \qquad (2a)$$

$$v^+ = \dot{x}^+, \qquad (2b)$$

$$U = J_c(x) v, \qquad (2c)$$

$$-dr_l \in N_{T_{\mathbb{R}_+}(g_l(x))}(U_l^+ + e_l U_l^-), \quad l = 1, \dots, p, \qquad (2d)$$

where $U(t)$ is the gap velocity, $dt$ is the Lebesgue measure on $\mathbb{R}$, $dv$ is the differential measure associated with jumps in velocity, $dr \in \mathbb{R}^p$ is a differential measure capturing constraint reaction impulses, and $e_l$ is the coefficient of restitution. The superscripts “-” and “+” denote pre- and post-impact values, respectively (left and right time limits at the impact instant).

Because impacts induce instantaneous jumps in velocity, the reaction $dr$ decomposes into a regular part and a sum of instantaneous impulses at the impact instants:

$$dr = \lambda(t) dt + \sum_{i \in I} p_i \delta_{t_i},$$

where $\lambda(t)$ represents the regular reaction and $p_i \in \mathbb{R}^p$ are the impulsive reactions, which are called percussions, at the impact times $t_i$.

According to Moreau's viability lemma [12], the contact conditions (2d) can also be expressed in impulsive form, linking velocities and impulses at impact. For an impact time $t_i$ with impulse $p_i \in \mathbb{R}^p$. We have, for each contact $l$,

$$g_l(x(t_i)) = 0 \Rightarrow 0 \le U_l^+(t_i) + e_l U_l^-(t_i) \perp (p_i)_l \ge 0, \qquad (3)$$

and $g_l(x(t_i)) > 0 \Rightarrow (p_i)_l = 0$. In particular, a nonzero impulse occurs only at active contact and enforces $U_l^+ + e_l U_l^- = 0$ in the normal direction. This complementarity form is the foundation of the time-stepping schemes of Moreau-Jean introduced next.

128 **Moreau-Jean time-stepping scheme.** Let $(x_n, v_n)$ be the approximation of $(x(t_n), v^+(t_n))$.

129 The Moreau-Jean (MJ) integrator discretizes (2) over $]t_n, t_{n+1}]$ using a modified $\theta$-method [8, 11]

$$\begin{aligned}
m(x_{n+\theta})\Delta t_n k_{v,n} &= \Delta t_n f(x_{n+\theta}, v_{n+\theta}, t_{n+\theta}) + J_c^\top(x_{n+1}) P_{n+1}, \\
k_{x,n} &= v_{n+\theta}, \\
v_{n+1} &= v_n + \Delta t_n k_{v,n}, \\
x_{n+1} &= x_n + \Delta t_n k_{x,n}, \\
U_n &= J_c(x_n)v_n, \\
U_{n+1} &= J_c(x_{n+1})v_{n+1},
\end{aligned} \qquad (4)$$

for $l \in \{1, \dots, p\}$, if $\bar{g}_{l,n+1} \le 0$, then $0 \le U_{l,n+1} + e_l U_{l,n} \perp P_{l,n+1} \ge 0$, otherwise $P_{l,n+1} = 0$,
with

$$ \Delta t_n = t_{n+1} - t_n, \quad t_{n+\theta} = t_n + \theta \Delta t_n, \quad v_{n+\theta} = v_n + \theta \Delta t_n k_{v,n}. $$

A contact l is treated as active if a predicted gap based on the free evolution of the system satisfies

$$ \bar{g}_{l,n+1} = g(x_{n+1,free}) \le 0, $$

where $x_{n+1,free}$ is defined as the value that would be taken by $x_{n+1}$ if $P_{n+1}$ vanished. Consequently one sees that, after selection of the active contacts thanks to a free prediction step, the contact laws are enforced at the velocity level as a nonlinear complementarity problem (NCP).

## 2.2 Neural network surrogate model

We model the flexible subdomain (0) by a neural surrogate trained offline to reproduce its dynamics subjected to external contact forces. The reference dynamics solved during data generation are governed by:

$$ \begin{aligned} m^{(0)}(x^{(0)}, t) \dot{v}^{(0)} &= f^{(0)}(x^{(0)}, v^{(0)}, t) + f_{\text{contact}}^{(0)}(t), \\ \dot{x}^{(0)} &= v^{(0)}, \end{aligned} \qquad (5) $$

where $m^{(0)}(x^{(0)}, t)$ is the mass matrix, $f^{(0)}(x^{(0)}, v^{(0)}, t)$ collects internal elastic and damping forces, and $f_{\text{contact}}^{(0)}(t)$ denotes the external contact forces. The neural network parameters are obtained through supervised training on trajectories generated by this system.

The temporal dependencies is captured in the form of a recurrent neural network (RNN), with an internal state encoding the memory of past dynamics $h_n$. At each time step, the surrogate must realize a discrete-time map of the form

$$ (f_{\text{contact},n}^{(0)}, h_n) \xrightarrow{\text{RNN}} (x_{n+1}^{(0)}, v_{n+1}^{(0)}, h_{n+1}), \qquad (6) $$

Without exclusion of other architectures that follows this input-output relationships, this map was chosen in this work as a Gated Recurrent Unit (GRU) layers. This structure ensures that the surrogate can be queried inside the time-stepping scheme using only the current contact impulse and its own internal state, without the need to pass the full flexible configuration to the solver. Its hidden state encodes the memory of past dynamics (positions and velocities), enabling the prediction of nonlinear responses with reduced input dimensionality.

Figure 2 illustrates the surrogate architecture: a stack of GRU blocks propagates the hidden state across time, while the final linear layer produces the configuration-level outputs. At each time step, the GRU updates its hidden state.

## 2.3 Coupled system and surrogate integration

We consider a system composed of two subdomains coupled through unilateral contact: a flexible part (0) and a rigid or FEM part (1). The generalized coordinates and velocities are partitioned as:

$$ x = \begin{bmatrix} x^{(0)} \\ x^{(1)} \end{bmatrix}, \quad v = \begin{bmatrix} v^{(0)} \\ v^{(1)} \end{bmatrix}. $$
<figure><img src="image_3.png" /><figcaption>Figure 2: GRU-based surrogate model.</figcaption></figure>

The Moreau-Jean update over $]t_n, t_{n+1}]$, which hence has block-diagonal form for the smooth contributions, reads

$$ \begin{pmatrix} m^{(0)}(x_{n+\theta}^{(0)}) & 0 \\ 0 & m^{(1)}(x_{n+\theta}^{(1)}) \end{pmatrix} \Delta t_n \begin{pmatrix} k_{v,n}^{(0)} \\ k_{v,n}^{(1)} \end{pmatrix} = \Delta t_n \begin{pmatrix} f^{(0)}(x_{n+\theta}^{(0)}, v_{n+\theta}^{(0)}, t_{n+\theta}) \\ f^{(1)}(x_{n+\theta}^{(1)}, v_{n+\theta}^{(1)}, t_{n+\theta}) \end{pmatrix} + \begin{pmatrix} J_{c,0}^{\top} \\ J_{c,1}^{\top} \end{pmatrix} P_{n+1}, \quad (7) $$

where:

$$ J_{c,i} = \frac{\partial g}{\partial x^{(i)}}(x). $$

The contact velocity is $U = J_{c,0}v^{(0)} + J_{c,1}v^{(1)}$, and the two subdomains interact only through the contact impulse $P_{n+1}$. Indeed, if $P_{n+1} = 0$, the two subsystems evolve independently under their own dynamics. When contact is active, the same contact impulse enters each subdomain in (7) through the respective generalized impulses $J_{c,0}^{\top} P_{n+1}$ and $J_{c,1}^{\top} P_{n+1}$.

While the integrator uses the impulse $P_{n+1}$, a potential neural surrogate of block (0) would be driven by the force applied during the step. We therefore should supply to the RNN the step-averaged contact force

$$ f_{\text{contact},n}^{(0)} = \frac{1}{\Delta t_n} J_{c,0}^{\top} P_{n+1}. \quad (8) $$

This idea is actually the starting point of the coupling approach presented in this work. In other terms, the coupling strategy that will be detailed hereafter is purely algebraic and is done through the contact impulse: block (1) enters through the classical Moreau-Jean terms in (7), while block (0) contributes through the surrogate response to the force $f_{\text{contact},n}^{(0)}$. This structure will be exploited to replace the expensive flexible block (0) by a neural surrogate without modifying the contact solver or complementarity structure.

### 2.3.1 Moreau-Jean scheme for the coupled system

We now write the Moreau-Jean update when the flexible subdomain (0) is replaced by the neural surrogate, while the rigid/FEM subdomain (1) is integrated with the classical time-stepping scheme. The surrogate is queried at each time step using the contact
impulse, through the interface relation (8), whereas the contact law and the global complementarity structure remain unchanged. The update over $[t_n, t_{n+1}]$ then reads:

$$ 
\begin{aligned}
(x_{n+1}^{(0)}, v_{n+1}^{(0)}, h_{n+1}) &= \text{RNN}(f_{\text{contact},n}^{(0)}, h_n), \\
m^{(1)}(x_{n+\theta}^{(1)}) \Delta t_n k_{v,n}^{(1)} &= \Delta t_n f^{(1)}(x_{n+\theta}^{(1)}, v_{n+\theta}^{(1)}, t_{n+\theta}) + J_{c,1}^{\top}(x_{n+1}) P_{n+1}, \\
k_{x,n}^{(1)} &= v_{n+\theta}^{(1)}, \\
v_{n+1}^{(1)} &= v_n^{(1)} + \Delta t_n k_{v,n}^{(1)}, \\
x_{n+1}^{(1)} &= x_n^{(1)} + \Delta t_n k_{x,n}^{(1)}, \\
v_{n+\theta}^{(1)} &= v_n^{(1)} + \theta \Delta t_n k_{v,n}^{(1)}, \\
x_{n+\theta}^{(1)} &= x_n^{(1)} + \theta \Delta t_n k_{x,n}^{(1)}.
\end{aligned} 
\quad (9) $$

The free prediction used for contact activation is

$$ 
\begin{aligned}
x_{n+1, \text{free}}^{(0)} &= \text{RNN}(0, h_n), \\
k_{x,n, \text{free}}^{(1)} &= k_{x,n|P_{n+1}=0}^{(1)}, \\
x_{n+1, \text{free}}^{(1)} &= x_n^{(1)} + \Delta t_n k_{x,n, \text{free}}^{(1)}, \\
\bar{g}_{n+1} &= g(x_{n+1, \text{free}})
\end{aligned} 
\quad (10) $$

The contact velocity is

$$ 
\begin{aligned}
U_n &= J_{c,0}(x_n) v_n^{(0)} + J_{c,1}(x_n) v_n^{(1)}, \\
U_{n+1} &= J_{c,0}(x_{n+1}) v_{n+1}^{(0)} + J_{c,1}(x_{n+1}) v_{n+1}^{(1)}.
\end{aligned} 
\quad (11) $$

A contact $l$ is treated as active if a predicted gap based on the free velocity satisfies $\bar{g}_{l,n+1} \le 0$. For each $l \in \{1, \dots, p\}$:

$$ \bar{g}_{l,n+1} \le 0 \Rightarrow 0 \le U_{l,n+1} + e_l U_{l,n} \perp P_{l,n+1} \ge 0, \quad \text{else } P_{l,n+1} = 0. $$

This restates the classical MJ contact law within the coupled system; it is precisely this complementarity condition that will incorporate the surrogate response in the next section.

### 2.3.2 Hybrid MJ formulation as an NCP

We adopt the following strategy for each Moreau-Jean step:

(i) **linearization of subdomain (1).** The dynamics of subdomain (1) are linearized once per time step about the impulse-free predictor computed with $P_{n+1} = 0$. In particular, on $[t_n, t_{n+1}]$ we suppose that the mass matrix is constant and independent of position $m^{(1)}(x^{(1)}) = m_1$, and that the regular forces are affine functions $f^{(1)}(x^{(1)}, v^{(1)}, t) = f_1 - K_1x^{(1)} - B_1v^{(1)}$.

(ii) **subdomain (0) surrogate.** Subdomain (0) is kept nonlinear, and its dynamics is evaluated by a trained RNN surrogate queried with the current contact load,

$$ (x_{n+1}^{(0)}, v_{n+1}^{(0)}, h_{n+1}) = \text{RNN}(\frac{1}{\Delta t_n} J_{c,0}^{\top} P_{n+1}, h_n), $$

so that $v_{n+1}^{(0)}$ depends nonlinearly on $P_{n+1}$ via the contact load.
203 (iii) gap computation. The gap is approximated as affine in generalized position x:

$$g(x) = g_{01} + J_c x = g_{01} + (J_{c,0} \quad J_{c,1}) \begin{pmatrix} x^{(0)} \\ x^{(1)} \end{pmatrix}.$$

204 **Proposition 1** (MJ NCP in hybrid surrogate formulation). Consider a Moreau-Jean time step over $[t_n, t_{n+1}]$ with parameter $\theta \in [0, 1]$. Under the strategy above, block (1) satisfies
205
206

$$\hat{m}_1 \Delta t_n k_{v,n}^{(1)} = \Delta t_n \hat{f}_{n,1} + J_{c,1}^{\top} P_{n+1}, \quad (12)$$

207 with

$$\hat{m}_1 = m_1 + \theta \Delta t_n B_1 + \theta^2 \Delta t_n^2 K_1, \quad \hat{f}_{n,1} = f_1 - K_1(x_n^{(1)} + \theta \Delta t_n v_n^{(1)}) - B_1 v_n^{(1)}.$$

208 Moreover, the contact update is governed by the nonlinear complementarity problem

$$0 \le c_n^{(1)} + G_0^{(1)} P_{n+1} + \Psi_{\text{MJ}}(P_{n+1}) \perp P_{n+1} \ge 0, \quad (13)$$

209 with surrogate contribution

$$\Psi_{\text{MJ}}(P) = J_{c,0}(e v_n^{(0)} + v_{n+1}^{(0)}(P)), \quad v_{n+1}^{(0)}(P) = \text{RNN}\left(\frac{1}{\Delta t_n} J_{c,0}^{\top} P, h_n\right),$$

210 and

$$c_n^{(1)} = J_{c,1}(e v_n^{(1)} + v_n^{(1)} + \Delta t_n \hat{m}_1^{-1} \hat{f}_{n,1}), \quad G_0^{(1)} = J_{c,1} \hat{m}_1^{-1} J_{c,1}^{\top}.$$

211 Sketch of proof. Starting from the block form of the MJ scheme, substituting the affine approximation for $f^{(1)}$, and the constancy of $m^{(1)}$, expanding $x_{n+\theta}^{(1)}$, $v_{n+\theta}^{(1)}$, one isolates $k_{v,n}^{(1)}$. After rearrangement, we obtain (12). The MJ contact law provides $0 \le U_{n+1} + e U_n \perp P_{n+1} \ge 0$ with $U = J_{c,0}v^{(0)} + J_{c,1}v^{(1)}$. Split the contact velocity, the MJ identity $J_{c,0}(e v_n^{(0)} + v_{n+1}^{(0)})$ still holds, but $v_{n+1}^{(0)}$ is now provided directly by the surrogate as a function of $P_{n+1}$. For the linearized block (1), the standard MJ algebra yields $J_{c,1}(e v_n^{(1)} + v_{n+1}^{(1)}) = c_n^{(1)} + G_0^{(1)} P_{n+1}$ with the stated $c_n^{(1)}$ and $G_0^{(1)}$. $\square$

215 In the proposed hybrid formulation, no stiffness or mass matrix is assembled for block (0), the surrogate directly provides $v_{n+1}^{(0)}$ as a function of $P_{n+1}$, and the problem is an NCP that can be handled by the Gauss-Seidel metho. During each SOR-Prox iteration, only the FEM block (subdomain (1)) requires classical linear solving, while surrogate contributes nonlinearly through the same NCP structure. This yields a full NCP formulation that blends numerical and data-driven mechanics without altering the contact resolution loop.

216 **2.3.3 SOR – Prox hybrid iteration**

217 At each time step, the contact update reduces to a nonlinear complementarity problem in the impulse $P_{n+1}$ of the generic form

$$0 \le G_0^{(1)} P_{n+1} + c_n^{(1)} + \Psi(P_{n+1}) \perp P_{n+1} \ge 0,$$
where $G_0^{(1)}P_{n+1} + c_n^{(1)}$ is the classical contribution of the rigid block (1), and $\Psi(P_{n+1})$ is the nonlinear surrogate contribution evaluated at the current impulse estimate. The NCP is solved using a Successive Over-Relaxation with proximal projections (SOR – Prox), widely used for nonsmooth contact dynamics.

The proximal projection onto a closed convex set $C \subset \mathbb{R}^p$ with respect to a positive definite metric $R$ is defined by

$$ \mathbf{prox}_C^R(x) = \arg \min_{y \in C} \|x - y\|_R, \quad \|x\|_R^2 = x^\top R x, $$

which yields the fixed-point form

$$ P_{n+1} = \mathbf{prox}_C^R \left( P_{n+1} - R^{-1} \left( G_0^{(1)} P_{n+1} + c_n^{(1)} + \Psi(P_{n+1}) \right) \right). \quad (14) $$

For computational efficiency we use block-diagonal metrics,

$$ R = \frac{1}{r_1} \mathbb{I}_1 \oplus \cdots \oplus \frac{1}{r_p} \mathbb{I}_p, $$

so that the update decouples component-wise:

$$ P_{l,n+1} = \mathbf{prox}_{C_l} \left( P_{l,n+1} - r_l \left( \sum_{j=1}^{p} (G_0^{(1)})_{lj} P_{j,n+1} + c_{l,n}^{(1)} + \Psi_l(P_{n+1}) \right) \right), \quad (15) $$

where $\Psi_l$ denotes the $l$-th component of the surrogate term.

<caption>**Algorithm 1:** SOR – Prox with relaxed Gauss-Seidel for contact NCP</caption>

**Input:** $P^0$, relaxation $\omega \in (0, 1]$

**Output:** $P_{n+1}$

1. Initialize $P_{l,n+1}^0 \leftarrow P_{l,n}$ if $\bar{g}_{l,n+1} \le 0$, else 0;

2. **repeat**

3. Evaluate surrogate once to compute $\Psi(P^\nu)$;

4. **for** $l = 1, \dots, p$ **do**

5. $\xi_l \leftarrow P_{l,n+1}^\nu - r_l \left( \sum_{j<l} (G_0^{(1)})_{lj} P_{j,n+1}^{\nu+1} + \sum_{j \ge l} (G_0^{(1)})_{lj} P_{j,n+1}^\nu + c_{l,n}^{(1)} + \Psi_l(P^\nu) \right)$;

6. $P_{l,n+1}^{\nu+1} \leftarrow (1-\omega)P_{l,n+1}^\nu + \omega \mathbf{prox}_{C_l}(\xi_l)$;

7. **end**

8. **until** converged;

9. **return** $P_{n+1}$;

In classical settings, one sets $r_l = 1/G_{ll}$ [15]. In our hybrid formulation, only $G_0^{(1)}$ is available, as the surrogate bypasses the FEM assembly of subdomain (0). Choosing $r_l = 1/G_{0,ll}^{(1)}$ neglects the missing block and may reduce stability. In practice, convergence is recovered by tuning the relaxation $\omega$, with small values (e.g. $\omega = 0.1$) for stiff hybrid dynamics. Thus, the SOR – Prox sweep incorporates the surrogate response, retaining the structure of classical schemes while avoiding assembly of the flexible block.
## 2.4 Position-level integrator (XDE)

In contrast to the velocity-level Moreau-Jean scheme, which enforces contact constraints by correcting velocities through impulses, the XDE variant enforces the Signorini law directly at the position level. Throughout this section we assume fully implicit Euler ($\theta = 1$) and zero restitution ($e = 0$), since a rebound velocity cannot be imposed consistently in a purely positional formulation. Active constraints are detected geometrically at $t_n$ via

$$g_n = \phi(x_n^{(0)}, x_n^{(1)}),$$
and the position-level Signorini condition is imposed at $t_{n+1}$:

$$0 \le g_{n+1}(x_{n+1}^{(0)}, x_{n+1}^{(1)}) \perp P_{n+1} \ge 0.$$

Hence, the impulse $P_{n+1}$ is computed to enforce non-penetration in configuration space, without relying on a free velocity prediction.

**Proposition 2 (XDE NCP in hybrid surrogate formulation).** Under the assumptions $\theta = 1, e = 0$, and $g$ smooth near $x_n$, the XDE integrator leads to the NCP

$$0 \le P_{n+1} \perp G_0^{(1)} P_{n+1} + c_n^{(1)} + \Psi_{\text{XDE}}(P_{n+1}) \ge 0, \quad (16)$$

with

$$\Psi_{\text{XDE}}(P) = J_{c,0}(x_{n+1}^{(0)}(P) - x_n^{(0)}), \quad x_{n+1}^{(0)}(P) = \text{RNN}\left(\frac{1}{\Delta t_n} J_{c,0}^{\top} P, h_n\right),$$

and

$$G_0^{(1)} = \Delta t_n J_{c,1} \hat{m}_1^{-1} J_{c,1}^{\top}, \quad c_n^{(1)} = g_n + \Delta t_n J_{c,1} v_n^{(1)} + \Delta t_n^2 J_{c,1} \hat{m}_1^{-1} f_n^{(1)}.$$

*Proof.* To obtain a tractable update over $[t_n, t_{n+1}]$, the gap is linearized around the known configuration $x_n$:

$$g_{n+1} \approx g_n + J_{c,0}(x_{n+1}^{(0)} - x_n^{(0)}) + J_{c,1}(x_{n+1}^{(1)} - x_n^{(1)}), \quad (17)$$

where $J_{c,i} = \partial g/\partial x^{(i)}|_{x_n}$ for $i = 0, 1$. The subdomain (1) is integrated implicitly and linearized once per time step around a predictor obtained without contact impulse. This yields an affine displacement response to the impulse $P_{n+1}$:

$$x_{n+1}^{(1)} - x_n^{(1)} \approx \Delta t_n v_n^{(1)} + \Delta t_n^2 \hat{m}_1^{-1} f_n^{(1)} + \Delta t_n \hat{m}_1^{-1} J_{c,1}^{\top} P_{n+1},$$

so that, after left-multiplication by $J_{c,1}$, we can precompute once per step the classical terms $G_0^{(1)}$ and $c_n^{(1)}$. The flexible subdomain (0) is replaced by the RNN surrogate, that is fed by the contact force $f_{\text{contact},n}^{(0)} = \frac{1}{\Delta t_n} J_{c,0}^{\top} P_{n+1}$, and returns the next configuration $x_{n+1}^{(0)} = \text{RNN}(f_{\text{contact},n}^{(0)}, h_n)$. Substituting into (17), we obtain the NCP (16). $\square$

During the SOR-Prox iterations, the surrogate is queried with $\frac{1}{\Delta t_n} J_{c,0}^{\top} P^{\nu}$ to compute $\Psi_{\text{XDE}}(P^{\nu})$. After convergence on $P_{n+1}$, a final RNN call gives $(x_{n+1}^{(0)}, v_{n+1}^{(0)})$, and block (1) is updated by the implicit step. The XDE scheme thus embeds the surrogate at the displacement level, preserving the Gauss-Seidel contact structure while enforcing non-penetration geometrically and avoiding velocity prediction.
## 2.5 Algorithmic implementation and convergence

* **Phase 1: Free evolution.** If no active contact is detected, the surrogate RNN is queried using the current state and zero contact load to evolve block (0), while a standard implicit $\theta$-method computes the free response of block (1).

* **Phase 2: Contact resolution.** When contact constraints are active, the nonlinear complementarity problem (NCP) is solved by SOR-Prox method. At each iteration, the impulse estimate is updated, and the RNN is re-evaluated accordingly until convergence. Both integrators share the same SOR-Prox contact solver; they differ only in (i) how the active set and coefficients $(c_n^{(1)}, G_0^{(1)})$ are formed, (ii) whether a free-velocity prediction is used (MJ only), and (iii) the surrogate contribution $\Psi(\cdot)$.

* **Phase 3: Final update.** Once converged, the final impulse $P_{n+1}$ updates positions and velocities of both subsystems.

Algorithm 2 summarizes the procedure.

**Algorithm 2:** Hybrid time step with SOR – Prox

**Input:** integrator ∈ {MJ, XDE}; $x_n, v_n$; trained RNN; $J_{c,0}, J_{c,1}$; $\Delta t_n$; $\theta$; $f_{n,1}$; tolerances $\varepsilon$; relaxation $\omega$

**Output:** $P_{n+1}, x_{n+1}, v_{n+1}, h_{n+1}$

1. Precompute $\hat{m}_1$, and $\hat{f}_{n,1}$

2. **if integrator = MJ then** // velocity-level

3. Build $\bar{g}_{n+1}$ from free-velocity forecast; mark active if $\bar{g}_{l,n+1} \le 0$

4. $c_n^{(1)} \leftarrow J_{c,1}(e v_n^{(1)} + v_n^{(1)} + \Delta t_n \hat{m}_1^{-1} \hat{f}_{n,1})$

5. $G_0^{(1)} \leftarrow J_{c,1} \hat{m}_1^{-1} J_{c,1}^{\top}$

6. Define $\Psi(P) \leftarrow J_{c,0}(e v_n^{(0)} + v_{n+1}^{(0)}(P))$, where $v_{n+1}^{(0)}(P) = \text{RNN}(\frac{1}{\Delta t_n} J_{c,0}^{\top} P, h_n)$

7. **else** // position-level

8. Detect active constraints geometrically from $g_n = \phi(x_n^{(0)}, x_n^{(1)})$

9. $c_n^{(1)} \leftarrow g_n + \Delta t_n J_{c,1} v_n^{(1)} + \Delta t_n^2 J_{c,1} \hat{m}_1^{-1} f_{n,1}$

10. $G_0^{(1)} \leftarrow \Delta t_n J_{c,1} \hat{m}_1^{-1} J_{c,1}^{\top}$

11. Define $\Psi(P) \leftarrow J_{c,0}(x_{n+1}^{(0)}(P) - x_n^{(0)})$, where $x_{n+1}^{(0)}(P) = \text{RNN}(\frac{1}{\Delta t_n} J_{c,0}^{\top} P, h_n)$

12. **end**

13. Initialize $P^0$ (warm-start on active set, else 0)

14. $P_{n+1} \leftarrow \text{SOR} - \text{Prox}(P^0, G_0^{(1)}, c_n^{(1)}, \Psi, \omega, \varepsilon)$

15. Update $(x_{n+1}, v_{n+1}, h_{n+1})$ with converged $P_{n+1}$; perform one final RNN call

16. **return** $P_{n+1}, (x_{n+1}, v_{n+1}, h_{n+1})$
**Error monitoring for convergence** For unilateral contact constraints, convergence is checked by monitoring violations of the Signorini conditions. The residual at iteration $\nu$ is defined as

$$err = \max_l \begin{cases} -\gamma_l, & \gamma_l < 0, \\ \min(\gamma_l, R_{ll}P_{l,n+1}^{(\nu)}), & \gamma_l > 0, \end{cases}$$

where $\gamma = GP_{n+1}^{(\nu)} + c_n$ is the constraint residual, $P_{n+1}^{(\nu)}$ is the current impulse iterate, and $R_{ll}$ is the local effective mass. This criterion ensures that both complementarity conditions and non-penetration are satisfied at convergence.

**Adaptive relaxation strategy** In contact simulations, high-frequency oscillations after impact can slow or even destabilize SOR-based contact solver when the relaxation parameter $\omega$ is not well chosen. A constant small $\omega$ is safe but unnecessarily slows convergence in well-behaved regimes. Instead, we employ an adaptive strategy: $\omega$ is reduced only when oscillatory residuals or stagnation are detected, allowing aggressive relaxation elsewhere while damping problematic oscillations in the presence of stiff contacts or rapid beam vibration.

At each sweep, the current Gauss-Seidel error err is monitored. If the error history exhibits oscillatory behavior (e.g., A-B-A motif), or stagnates over multiple iterations, $\omega$ is reduced (typically multiplied by a factor $\eta < 1$) but never below a threshold $\omega_{min}$.

<caption>**Algorithm 3:** Adaptive update of the relaxation parameter $\omega$</caption>

Input: $\omega_{init}$, $\omega_{min}$, $\eta < 1$, $\varepsilon_{osc}$, $\varepsilon_{tol}$, $k$
Output: Updated $\omega$ per SOR sweep
1 $\omega \leftarrow \omega_{init}$;
2 for each SOR sweep do
3 Record current error err and update rolling history [$e_1, e_2, e_3$];
4 if $|e_1 - e_3| < \varepsilon_{osc}$ and $|e_2 - e_3| > \varepsilon_{osc}$ // oscillation then
5 $\omega \leftarrow \max(\eta \omega, \omega_{min})$;
6 end
7 else if $|err - err_{last}| < \varepsilon_{tol}$ for k sweeps // stagnation then
8 $\omega \leftarrow \max(\eta \omega, \omega_{min})$;
9 end
10 end

Typical values are $\omega_{init} = 0.5$, $\omega_{min} = 10^{-4}$, $\eta = 0.7$, $\varepsilon_{osc} = 10^{-10}$, with stagnation detected after $O(10)$ sweeps with negligible error change. This adaptive scheme preserves stability in challenging contact regimes while avoiding unnecessary damping in smooth phases, improving overall robustness.

## 3 Numerical results

We evaluate the proposed coupling strategy on two test cases, see Fig. 3. The first one is an academic one-dimensional problem: a damped harmonic oscillator in intermittent contact with a free mass. This setup provides a controlled scenario to verify three aspects: the ability of GRU-based surrogate to reproduce nonlinear dynamic behaviors of the
<figure><img src="image_4.png" /><figcaption>Figure 3: Two systems used to exemplify the coupling method. In blue, the subsystem simulated by means of a NN surrogate model.</figcaption></figure>

flexible subsystem, the convergence of Gauss-Seidel contact iterations in the presence of the surrogate, and the integration of an RNN surrogate within a classical nonsmooth time-stepping solver. The second test involves the impact between a free mass and a flexible beam, where the nonlinear response of the beam is replaced by a neural surrogate.

All simulations consider frictionless unilateral contact. Jacobians and gaps are linearized at each time step. For the position-level integrator (XDE) we restrict to ($e = 0, \theta = 1$).

## 3.1 Harmonic oscillator–mass contact

We first consider contact between a damped oscillator and a free mass. The oscillator (subdomain (0) with mass $m^{(0)} = 1$, damping $c = 0.1$, stiffness $k = 3$) is replaced by a GRU surrogate, while the free body (subdomain (1)) has mass $m^{(1)}$. This test assesses the surrogate’s ability to reproduce the oscillator dynamics and to remain stable inside a nonsmooth solver.

Both the velocity-level MJ integrator with restitution coefficient $e = 0$ and the position-level XDE integrator are tested and compared against a reference solution computed with the same integrator but without a neural surrogate.

The surrogate consists of three stacked GRU layers with 128 hidden units each, trained to predict the oscillator’s position and velocity ($x^{(0)}$, $v^{(0)}$) from the applied contact-force sequence. The training dataset is generated by simulating oscillator trajectories under random windowed impulsive loadings. Once trained, the surrogate reproduces the free evolution with high accuracy, with only a negligible transient at the beginning due to hidden-state initialization.

A representative contact case is reported with $x_0^{(0)} = 0, v_0^{(0)} = 1, x_0^{(1)} = 3, v_0^{(1)} = -1.5, \Delta t_n = 10^{-2}$, and $e = 0$.

Figures 4–5 compare the MJ and XDE coupling formulations to the reference simulation. Both formulations accurately reproduce the reference positions, velocities, and contact impulses. Quantitatively, both MJ and XDE integrators closely track the reference trajectory. Over the two bodies, the mean displacement error is $5 \times 10^{-2}$ m (MJ) and $6 \times 10^{-2}$ m (XDE), with the maximum deviation remains below $1.7 \times 10^{-1}$ m. Consistently with theory, MJ allows a small bounded penetration of order $10^{-2}$ m, while XDE maintains non-penetration within numerical tolerance ($|g_{min}| \approx 10^{-7}$ m).
<figure><img src="image_5.png" /><figcaption>Figure 4: Harmonic oscillator-mass contact: comparison of the coupled solver (MJ and XDE) with the reference simulation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Gap behavior highlights the integrator difference: velocity-level MJ vs. position-level XDE.</figcaption></figure>

These results confirm that the surrogate can replace the flexible subdomain without altering the global nonsmooth behaviour: the coupled scheme retains the expected contact-law properties and the reconstructed trajectories remain close to the FEM reference. The GRU surrogate accurately reproduces the nonlinear dynamics across contact activation and release; small initial transients are attributed to hidden-state cold-start and vanish after a few iterations.

## 3.2 Beam-mass impact simulation

We now consider a problem in which a flexible cantilever circular beam (Reissner type) interacts with a lumped point mass through a unilateral contact constraints. The beam lies in the $x-z$ plane, is clamped at $(0, 0)$ and uniformly discretized into $N-1$ finite elements with $N$ nodes. Each node carries three degrees of freedom: the horizontal displacement $x$, the vertical displacement $z$, and the rotation $\theta_y$ about the $y$-axis. The free mass has two translational degrees of freedom $(x, z)$ and is initially located above the beam. Neither gravity nor damping is applied.

The material and geometric properties of the beam are: Young's modulus: $E = 6.9 \times 10^{10}$ Pa, Poisson's ratio: $\nu = 0.33$, density: $\rho = 2700$ kg/m$^3$, section radius: $r = 5 \times 10^{-3}$ m, length: $L = 0.5$ m, number of nodes: $N = 10$.
In the coupled formulation, the beam dynamics is replaced by a neural surrogate, while the lumped mass is integrated with the nonsmooth scheme as a classical rigid body.

**Surrogate training setup.** A dataset of trajectories is generated by simulating the full beam-mass system under multiple initial conditions of the lumped mass. At each time step we store:

* beam states: nodal positions $(x, z)$ and rotations $(\theta_y)$, and their time derivatives $(\dot{x}, \dot{z}, \dot{\theta}_y)$.

* contact force split across the two adjacent nodes $(i, i+1)$ of the active contact segment with weights $(1-\lambda)$ and $\lambda$ with $\lambda \in [0, 1]$ according to the contact location, projected onto each node's local frame.

The surrogate is trained to map the contact force history to the time evolution of the beam DOFs (excluding the clamped node). A three-layer GRU with 256 hidden units per layer is used. Although no timing comparison is reported here (low DOF case), inference is constant-cost per call and thus does not scale with mesh size, unlike FEM solves.

**Contact handling.** At each iteration of the coupled solver, a contact manager locates the closest beam segment and computes the quantities required for contact resolution: the signed gap $g(x)$, the outward normal $n(x)$, and the Jacobians $J_{c,0}$ (w.r.t. beam DOFs) and $J_{c,1}$ (w.r.t. mass DOFs).

The tangent vector of the segment $[p_i, p_{i+1}]$ is $t = p_{i+1} - p_i$ and the outward unit normal is $n = \frac{[-t_z, t_x]}{\|[t_z, t_x]\|}$. Let $q$ be the mass center and $c = p_i + \alpha(p_{i+1} - p_i)$ the contact point with $\alpha \in [0, 1]$, then the signed gap to the beam surface is $g(x) = n^\top(q - (c + rn))$.

Only the segment nodes $(i, i+1)$ contribute to $J_{c,0}$, with partitioning $(1-\alpha, \alpha)$, rotated into the local frame by $R(\theta_k)n = \begin{bmatrix} \cos\theta_k & -\sin\theta_k \\ \sin\theta_k & \cos\theta_k \end{bmatrix} \begin{bmatrix} n_x \\ n_z \end{bmatrix}$, while $J_{c,1} = n$ for the mass. The resulting Jacobians are sparse and assembled as in Algorithm 1.

<caption>**Algorithm 1** Contact Jacobian Assembly for Beam-Mass System (2D)</caption>

1: **Input:** Beam node positions $\{p_k\}$, nodal rotations $\{\theta_k\}$, mass position $q$
2: Find closest segment $[p_i, p_{i+1}]$ and compute $\alpha \in [0, 1]$
3: Compute tangent $t = p_{i+1} - p_i$ and normal $n = \frac{[-t_z, t_x]}{\|[t_z, t_x]\|}$
4: **Node i:** $(J_{c,0})_i \leftarrow -(1-\alpha) R(\theta_i)n$
5: **Node i+1:** $(J_{c,0})_{i+1} \leftarrow -\alpha R(\theta_{i+1})n$
6: **Mass:** $J_{c,1} \leftarrow n$

**Numerical results.** We consider $\Delta t_n = 10^{-4}$ s and an initial mass position $x^{(1)} = (0.49, 0.007)$, with downward velocity $v^{(1)} = -0.2$. both the MJ and the XDE integrators are run under identical conditions and compared against a FEM reference solver. Figure 6 reports the contact dynamics. The evolution of the contact impulse $P_{n+1}$ for both integrators is consistent and matches the reference, with discrepancies vanishing when averaged over a few time steps. The residual errors of the SOR-Prox iterations remain
<figure><img src="image_7.png" /><figcaption>Figure 6: Beam–mass: contact dynamics (MJ/XDE vs FEM reference).</figcaption></figure>

well below the tolerance threshold. As expected, MJ admits a small bounded penetration, whereas XDE keeps the gap strictly non-negative within solver precision.

The vertical motion of the two nodes nearest the contact (Nodes 8 – 9) together with the mass is shown in Figure 7. The surrogate predictions under both integrators closely match the reference solution, including the post-impact rebound. Extending the comparison to several beam nodes (Figure 8) shows similarly good agreement. A short transient near the clamped is attributable to GRU hidden-state and vanishes after a few time steps, but the deviation remains within acceptable error bounds. Rotational degrees of freedom are likewise reproduced with high accuracy, including during contact activation (Figure 9).

<figure><img src="image_8.png" /><figcaption>Figure 7: Vertical displacements of Nodes 8–9 and the mass (MJ/XDE vs FEM reference).</figcaption></figure>

Quantitatively, the comparison is restricted to the contact-relevant displacements (beam nodes 8–9 and the rigid mass). The maximum deviation from the FEM reference remains below $4 \times 10^{-4}$ m for MJ and below $3.3 \times 10^{-4}$ m for XDE. Again, MJ exhibits a bounded penetration ($\sim 10^{-4}$ m), while XDE preserves $g \ge 0$ within tolerance.

The fact that both translational and rotational DOFs are captured from force histories alone demonstrates the suitability of RNN surrogates as drop-in replacement within the coupled nonsmooth solvers. The surrogate captures the beam deformation over the entire contact event, including activation and release phases.

Overall, both MJ and XDE integrators reproduce the FEM reference with comparable accuracy. MJ provides efficient velocity-level enforcement with a small bounded
<figure><img src="image_9.png" /><figcaption>Figure 8: Vertical displacements of selected nodes (MJ/XDE vs FEM reference).</figcaption></figure>

penetrations, while XDE guarantees exact satisfaction of the Signorini condition at the configuration level. The surrogate can therefore replace the flexible subsystem without altering the global non-smooth behaviour or the structure of the contact solver.

# 4 Conclusion

We proposed a hybrid framework that couples neural-network surrogates with classical time-stepping solvers for nonsmooth mechanical systems under unilateral constraints. A trained recurrent surrogate replaces the flexible subdomain inside the standard contact-resolution loop. The surrogate contribution enters the complementarity conditions without assembling or inverting any matrix block for the flexible part, thereby preserving both the SOR-Prox solver structure and the global NCP formulation.

The approach was instantiated to both velocity-level (Moreau-Jean) and position-level (XDE) integrators, and validated on harmonic-oscillator and beam-mass benchmarks. The results show accurate reproduction of nonlinear dynamics, stable Gauss-Seidel convergence, and consistent enforcement of unilateral constraints.

Overall, this work demonstrates that learned surrogates can be integrated into established nonsmooth dynamics solvers in a non-intrusive way, preserving modularity while avoiding FEM assembly/inversion for the surrogate subdomain. Future work will address
<figure><img src="image_10.png" /><figcaption>Figure 9: Rotations $\theta^{(0)}$ of selected nodes (MJ/XDE vs FEM reference).</figcaption></figure>

frictional contact, scalability to large systems, and physics-consistent surrogate training. The current surrogate does not enforce energy consistency or reciprocity; incorporating such constraints could further improve fidelity.

References

[1] V. Acary. “Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb’s friction”. In: *Computer Methods in Applied Mechanics and Engineering* 256 (2013), pp. 224–250. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2012.12.012.

[2] M. Anitescu and F. A. Potra. “A time-stepping method for stiff multibody dynamics with contact and friction”. In: *International Journal for Numerical Methods in Engineering* 55.7 (2002), pp. 753–784. DOI: https://doi.org/10.1002/nme.512.

[3] S. Conti, Z. Guo, and A. Manzoni. “Multi-fidelity surrogate modeling using long short-term memory networks”. In: *arXiv preprint arXiv:2208.03115* (2022).

[4] N. Franco, S. Fresca, F. Tombari, and A. Manzoni. “Deep learning-based surrogate models for parametrized PDEs: handling geometric variability through graph neural
networks”. In: *Journal of Computational Physics* 493 (2023), p. 112882. DOI: 10.1016/j.jcp.2023.112882.

445

446

[5] S. Fresca, G. Gobat, A. Fedeli, A. Frangi, and A. Manzoni. "Deep-learning-based reduced order models for the real-time simulation of the nonlinear dynamics of microstructures". In: *International Journal for Numerical Methods in Engineering* 123.20 (2022), pp. 4749-4777. DOI: 10.1002/nme.7006.

447

448

449

450

[6] C. Gear, B. Leimkuhler, and G. Gupta. "Automatic integration of Euler-Lagrange equations with constraints". In: *Journal of Computational and Applied Mathematics* 12-13 (1985), pp. 77-90. ISSN: 0377-0427. DOI: https://doi.org/10.1016/0377-0427(85)90008-1.

451

452

453

454

[7] I. Goodfellow, Y. Bengio, and A. Courville. *Deep Learning*. MIT Press, 2016.

455

[8] M. Jean. "The non-smooth contact dynamics method". In: *Computer Methods in Applied Mechanics and Engineering* 177.3 (1999), pp. 235-257. ISSN: 0045-7825. DOI: https://doi.org/10.1016/S0045-7825(98)00383-1.

456

457

458

[9] X. Merlhiot. "On some industrial applications of time-stepping methods for nonsmooth mechanical systems: issues, successes and challenges". In: *Euromech Colloquium [516] — Nonsmooth contact and impact laws in mechanics*. Grenoble, France, July 2011.

459

460

461

462

[10] J. J. Moreau. "Application of convex analysis to the treatment of elastoplastic systems". In: *Applications of Methods of Functional Analysis to Problems in Mechanics*. Ed. by P. Germain and B. Nayroles. Berlin, Heidelberg: Springer Berlin Heidelberg, 1976, pp. 56-89. ISBN: 978-3-540-38165-5.

463

464

465

466

[11] J. J. Moreau. "Unilateral Contact and Dry Friction in Finite Freedom Dynamics". In: *Nonsmooth Mechanics and Applications*. Ed. by M. J.J. and P. P.D. Vol. 302. International Centre for Mechanical Sciences (Courses and Lectures). Springer, 1988, pp. 1-82. DOI: 10.1007/978-3-7091-2624-0_1.

467

468

469

470

[12] J. Moreau. "Numerical aspects of the sweeping process". In: *Computer Methods in Applied Mechanics and Engineering* 177.3 (1999), pp. 329-349. ISSN: 0045-7825. DOI: https://doi.org/10.1016/S0045-7825(98)00387-9.

471

472

473

[13] S. Schoeder, H. Ulbrich, and T. Schindler. "Discussion of the Gear-Gupta-Leimkuhler method for impacting mechanical systems". In: *Multibody System Dynamics* 31.4 (May 2013), 477-495. ISSN: 1573-272X. DOI: 10.1007/s11044-013-9370-y.

474

475

476

[14] D. E. STEWART and J. C. TRINKLE. "AN IMPLICIT TIME-STEPPING SCHEME FOR RIGID BODY DYNAMICS WITH INELASTIC COLLISIONS AND COULOMB FRICTION". In: *International Journal for Numerical Methods in Engineering* 39.15 (1996), pp. 2673-2691. DOI: https://doi.org/10.1002/(SICI)1097-0207(19960815)39:15<2673::AID-NME972>3.0.CO;2-I.

477

478

479

480

481

[15] C. Studer. *Augmented Time-Stepping Integration of Non-Smooth Dynamical Systems*. 2009.

482

483

[16] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu. *The Finite Element Method: Its Basis and Fundamentals*. 7th. Butterworth-Heinemann, 2013.

484

485