1 Introduction
During the last couple of years, ”digital twins” appear to be one of the key concepts towards digitalization. They provide computeraided assistance for realworld models allowing to monitor their state at specified positions or make prediction for the the future states [10]
. Model order reduction realizes digital twins by enabling realtime execution. It reduces the degrees of freedom while keeping the essential physical properties of the model. Moreover it enables the reusability of such models for other scenarios. An extensive literature is devoted to model order reduction applied for different kinds of partial differential equations, see e.g.,
[35, 42, 5, 18, 16, 9]. Model reduction methods are typically categorized between datadriven methods [40, 38, 44] or physicsbased methods such as balanced truncation [8, 12, 33], Krylov subspace methods [3, 29, 13, 7, 45, 11] and modal reduction [5].In this paper we focus on reduction approaches for dynamic contact problems in linear elasticity. Such problems are naturally nonlinear due to the unknown moving contact interface, see, e.g., [20, 31, 46, 21] for introduction to contact mechanics and [32] for first results on reduction approaches in this matter. In [4]
, a dynamic problem with a linear contact condition is discussed and both the displacements and the Lagrange multipliers are reduced using methods such as the singular value decomposition or the nonnegative matrix compression. In contrast, a reduction scheme for linear nodetonode contact problems was introduced in
[36], which reduces only the primal displacements. Such contact conditions allow, however, only sticking at the contact and loosing the contact in the normal direction .Mechanical contact problems represent variational problems under unilateral constraints [46, 15], which are commonly solved by the penalty method [28] or the augmented Lagrange method [22] after the discretization. Despite their advantages, these methods turn out to be inefficient for designing reduction schemes. Moreover, the contact shape is usually not known a priori. This is a major challenge for reduction schemes, since the reduced model needs to account for the contact impact.
In this work we propose a novel reduction method for the class of nodetosegment contact problems as described in, e.g., [46]. This discretization technique of the contact condition is widespread in the field of contact mechanics since it can handle large deformations and allows more flexibility at the contact zone such as sliding, see [48, 47, 49]. Solvers for such contact problems have been integrated into nonlinear finite element commercial software [2, 23]. For early implementations of nodetosegment techniques we refer to [26, 27].
The novel reduction approach is purely physicsinformed, i.e., no trajectories of the full model are required. The projection matrix can be computed based on the system matrices and the external load vector. Once the reduced model is established within an offline procedure, we turn to the adjoint space of the reduced displacement, i.e., the space of the dual Lagrange multipliers. In contrast to
[36], the timedependent dual problem represents a Nonlinear Complementarity Problem (NCP). The latter is approximated by a sequence of Linear Complementarity Problems (LCPs) which can be solved within fixedpoint iterations for each time step. In terms of convergence, a quadratic convergence rate for the fixedpoint iterations is guaranteed due the general theoretical framework [1, 37, 30]. Combining with the fact that the sequential LCPs stem from the reduced space and the number of the constraints is small, the convergence rate provides major savings with respect to execution time.More specifically, the model order reduction is essentially based on a partitioning of the the overall displacement degrees of freedom in the flavor of the CraigBampton method [41], which separates the interior nodes, i.e., the slave nodes, from the contact nodes, i.e., the master nodes. Due to the nonlinear constraints, the system stiffness matrix turns out to depend on the Lagrange multipliers. However, the partitioning of the nodes provides linear system blockmatrices for the slave nodes allowing an accurate reduction. After the reduced slave nodes are computed, the nonlinear blockmatrices of the master nodes are utilized for the static condensation restoring the coupling between the master and slave nodes. Thus the partitioning provides essential improvement in the accuracy for recovering the Lagrange multipliers and the contact shape as well. To put it in a nutshell, our proposed reduction scheme combines Krylov subspace methods, adjoint methods for the linearized problem and the CraigBampton partitioning technique. The overall methodology is general in nature, but so far it has been worked out in detail and tested in the planar case only.
Our current contribution presents the new reduction algorithm in the following steps: In Section 2 the variational formulation of dynamic contact problems with quadratic constraints is derived. The common nodetosegment technique is described in Section 3. Furthermore, Section 4 and Section 5 discuss the solving procedure of the contact problem. Instead of solving the underlying NCP, the corresponding linearized system (LCP) is sequentially solved within fixedpoint iterations. The resulting reduction method is outlined in Section 6. The performance of the approach is highlighted by means of two numerical examples in Section 7. Finally, a conclusion is drawn in Section 8.
2 The Generic Formulation
In this section we first outline the governing equations for dynamic contact in the semidiscretized form. Afterwards, the derivation of the variational formulation is briefly sketched.
2.1 The Generic Semidiscretized Model
In this paper we restrict ourselves to finite element methods for a comprehensive contact treatment. The contact is discretized by the nodetosegment technique that allows the sliding of a node over segments defined on the contact area. As explained below, such an approach leads to a quadratic inequality in the displacement vector that is combined with the equations of structural dynamics. In this way, a generic model for the system of semidiscretized equations can be derived that reads
(1a)  
(1b)  
(1c)  
(1d) 
where denotes the vector of nodal displacement variables, the system mass and stiffness matrix, respectively, and the external loads. Moreover,
is a thirdorder tensor,
stands for the constraint matrix and for a given initial clearance vector. In (1b) and (1d) the quadratic term is defined as with for . Analogously, the inequality constraints are to be interpreted componentwise and stand for the discretized nonpenetration conditions of all variables in the contact interface where the vector of Lagrange multipliers enforces the constraint in the dynamic equation (1a). The multiplier , standing for the discretized contact pressure, is always positive and its inner product with the constraints satisfies a complementarity condition. Note also that(2) 
which means that this term can be viewed as a nonlinear contribution to the stiffness matrix that depends on the Lagrange multipliers.
In this paper, we will introduce a physicsbased model order reduction method for the semidiscretized equations (1) that projects the nodal variables onto a much smaller space while explicitly preserving the constraints and the Lagrange multipliers. For this purpose, we assume a frictionless, adhesivefree normal contact in combination with small deformation theory and a linearelastic material. Those assumptions will lead to (1) if

the Lagrange multiplier method is used to enforce the nonpenetration condition

the finite element method is applied for the spatial discretization

the nodetosegment contact technique is used for the discretization of the contact zone.
We do not consider approaches like the Augmented Lagrangian method or the Nitsche method. The same holds for additional friction effects. In contrast to [36], the nodetosegment provides quadratic inequalities as constraints. This technique will be sketched in greater detail in Section 3.
For a contact problem with bodies that fits into the framework described so far, the mass matrix and the stiffness matrix consist of diagonal block matrices that stem from the discretization of each individual body. The constraint tensors and will then typically exhibit a sparse structure where each row stands for a nodetosegment contact condition. For further details of various contact models and discretization schemes we refer to [31, 46]. Before we proceed further with the numerical approaches, we briefly summarize the setting of the classical variational formulation of contact problems in the dynamic case.
2.2 Strong Formulation of the Dynamic Contact Problem
We introduce some notational preliminaries. The initial undeformed solid or is a disjoint union of bounded connected subdomains i.e., with for On we consider a frictionless, adhesivefree normal multibody contact problem and a linearelastic material describing the displacement field and the contact pressure where is the spatial variable and stands for the temporal variable. The restriction of the displacement field on the corresponding domain we denote as The boundary of each elastic body is decomposed into where the latter represents the contact interface. The contact problem in strong form is then given by a dynamic boundary value problem (BVP): For all and for the corresponding initial data find such that
(3)  
where the constant denotes the mass density of the body , the volume force and the surface traction. The stress tensor and the linearized strain tensor are given by
(4)  
(5) 
with Young’s modulus and Poisson’s ratio Furthermore, the boundary value problem (3) is subject to the contact conditions
(6) 
with where is the gap function described as
(7) 
and is the contact pressure on , respectively. Note, that is the outer normal vector on and is the projection point of to the body in the outer normal direction, fulfilling
(8) 
On the one hand the gap function (7) simply measures the distance between and But on the other hand depending on the underlying geometry, (7) can lead to quite complex expressions. Moreover, in case of the nonuniqueness of (8), the gap function (7) can turn into a nonsmooth function. For simplicity, we assume that there exist functionals and such that for all it holds
(9) 
i.e., the outer normal depends on the displacement vector affinelinearly.
In [36] we considered a nonpenetration condition allowing only normal contact such that no movement in the tangential direction was possible, see also [46]. In particular the gap function was linear. In this work, however, using the assumption (9) we want to generalize the constraint to a quadratic condition allowing also sliding movement such that a point can move along the contact area.
2.3 Weak Formulation of the Dynamic Contact Problem
In order to state the weak formulation of the contact problem, we introduce the function spaces
(10)  
(11) 
where is the set of admissible displacements and is decomposed bodywise such that for all it holds and on . The set is the convex cone of the Lagrange multipliers defined on which is the union of all Furthermore, the following abstract notation is introduced:
(12) 
The nonpenetration condition in weak form with test function is in general given by
(13) 
In particular for a twobody problem, i.e. the righthand side of (13) reads:
(14) 
Note, that we provide the weak formulation solely for the twobody problem. The extension to the general multibody case is straightforward. Using (14) we define the trilinear form on
(15) 
the bilinear form on
(16) 
and the linear form on
(17) 
After these preparations, the weak form of the dynamic contact problem is stated as follows: For each find the displacement field and the contact pressure such that
(18)  
The system (18) is discretized with respect to the spatial variable by applying the standard Galerkin projection with basis functions and nodal variables to the displacement field.
Moreover, the integral over in the weak dynamic equation can be approximated by a discrete sum
(19)  
with the constants where is the area around node Putting finally as discrete pressure variable, the discretized constraint with the thirdorder tensor the matrix and contact offset lead to a quadratic inequality of the form (1). Note that only a few displacements are involved for the contact condition and hence most of the entries of and are zero.
3 The NodeToSegment Contact Condition
In this section we study the contact condition using the nodetosegment technique in more detail. As already indicated above it turns out that this condition is quadratic in terms of the displacements. In the following we assume a 2D setting where linear finite elements are used for the discretization. Let
(20) 
be two neighbouring contact nodes. They define a contact segment that an be written as
(21) 
Furthermore, we assume that after the discretization a nodesegment correspondance is established at the contact zone. Let be the corresponding node of the segment . The distance of the segment to the node is given by (see Figure 1(a))
(22) 
Similarly as for the gap function in (7), the distance function in (22) can lead to a quite complex expression. Therefore we define an angular function representing an easier alternative to (22). For a node let be the nearest segment to and let denote the normal vector orthogonal to the line containing . Then the angular function is given by (see Figure 1(b))
(23) 
Note, that for two nonzero vectors the equality always holds, where denotes the angle enclosed by and
For (23) three different scenarios can be distinguished:

, which means that lies in the half of the line with positive normal.

, which means that lies on the line containing the segment .

, which means that lies in the half of the line with negative normal.
The normal to depends on the end nodes of and is given by
(24) 
In order to derive a formulation for (23) in terms of the displacements and also to keep the notation concise, we write with the vectors standing for the initial position of the node and for the corresponding displacement. The label runs through Based on this notation, the function (23) can be reformulated as
(25) 
The function defined in (25) is at most of quadratic order in terms of the displacements. Furthermore, by assembling the corresponding terms of the summation within the equation (25) in a proper way we can determine the global constraint matrices , the global vectors and and the scalars . For the nodal displacement vector and contact nodes the contact condition reads
(26) 
In a more compact form (26) can be stated as
(27) 
with
(28) 
The latter represents a quadratic contact condition derived from the nodetosegment discretization technique.
4 Solving Static Contact Problems
A common approach to solve mechanical contact problems is the augmented Lagrangian method [39, 22, 46]. The method is very flexible to different descriptions of the contact region. However, here we apply a dual approach, allowing us to turn to the adjoint problem resulting into a NCP. In the context of model order reduction for constrained problems dual approaches have a great benefit, which will be highlighted in the following sections. First we address the dual formulation and the solving procedure of a NCP in case of static contact problems.
4.1 NCP for the static contact problem
In the following we consider the quadratic contact condition
(29) 
Then, the variational formulation of the mechanical contact problem starts by considering the energy expression
(30) 
The corresponding KKTconditions read
(31)  
Solving for the displacements in the KKTconditions leads to the representation
(32) 
4.2 Solving a sequence of LCPs for the static contact problem
A wellknown and successful technique of solving NCPs is to consider their linear approximation, representing an LCP system, for which an LCP solver can be applied repeatedly. This method has been introduced in [30] where the problem statement is given in a generalized framework. More about iterative techniques for solving variational inequalities can be found in [37, 1]. A general form of an LCP is defined as follows:
(36) 
LCP problems may be solved by applying either FischerBurmeister type algorithms [19] or
more efficiently Lemke’s algorithm (see [17], 4.4.5). A detailed discussion of Lemke’s algorithm as well as a much more comprehensive study of Linear Complementarity Problems and can be found in [17]. Specifically, we use a Python implementation of the algorithm [34].
The core idea, which below is also used for our reduction technique, is to exchange the original NCP problem
(37) 
by the linear approximation
(38) 
where
(39) 
with
(40) 
represents the quadratic constraint function. The system (38) is solved by an LCP solver (e.g., the Lemke method) until a carefully chosen stopping criterion depending on the initial starting point is satisfied. Then is updated by and the LCP problem is solved again. Finally, this procedure yields an approximation for the solution vector of the nonlinear problem (37).
5 Solving Dynamic Contact Problems
After the static contact problems have been addressed, we turn our attention to dynamic contact problems. Similar to the static case, the key idea for solving the timedependent contact problems is to replace the nonlinear dual system by its linear approximation.
5.1 Notational Convention
For subsequent use, we introduce a new notation for our solution variables. By we denote where stands for the th step of the discretized time grid. The Lagrange multiplier carries two indices, where the superscript denotes the corresponding fixedpoint iteration of the LCP. If it is clear from the context that the time instance is kept fixed, we make use of the variable
5.2 NCP for the dynamic contact problem
The dynamical model combined with the KuhnTacker conditions (31) leads to the system
(41a)  
(41b)  
(41c)  
(41d) 
All four conditions in (41) must be fulfilled in each time step. The balance of momentum (41a) is discretized in time by the implicit Euler scheme, a straightforward and unconditionally stable method, also known as the first order backward differentiation formula. To this end, we replace the second order derivative by the finite difference
(42) 
where is the time stepsize. Note that for the velocity it holds a fact that is hidden behind (42). Furthermore, due to the presence of the the second order time derivative, the implicit Euler leads to a twostep method. Since (42) possesses first order of accuracy only, other integration schemes can be applied instead, e.g. the generalized method, which is very common in structural mechanics, see [43]. In this work, however, we mainly focus on designing a reduction scheme and therefore, we continue working with the implicit Euler method. Inserting (42) into equation (41) and keeping in mind the abbreviation (33), we obtain
(43) 
Assuming that the previous two time steps are known, (43) is solved for , which leads to
(44) 
We remark that the twostep time discretization requires initial values and to start. If and as initial displacement and velocity are given, one can compute by an explicit Euler step and then continue with the twostep formula (43). The initialization of the Lagrange multiplier is shortly addressed in Section 6. If the inequality constraints were replaced by the equality constraints , the index of the resulting differentialalgebraic equation would equal 3, which means that special care must be taken for the time integration [14, 25, 43]. Inserting (44) into the constraints (41b)(41d) we obtain an NCP problem similar to (35). After computing the Lagrange multiplier the displacements are given by (44). Note that the NCP has to be solved in each time step.
5.3 Solving a sequence of LCPs for the dynamic contact problem
After having introduced the transient model, we address the corresponding solution procedure in more detail. The core idea of the algorithm is the generalization of the iterative solver for the static problems to an iterative solver for the dynamic problemsby substituting the displacement vector in (34) by
(45) 
The quadratic constraint function is defined as
(46) 
The function (46) defines for each time step the adjoint problem of (41), which represents an NCP problem, i.e.,
(47) 
Instead of directly solving (47), we prefer to solve the linear approximation of (47) within fixedpoint iterations for each time step. For this purpose, we approximate the nonlinear constraint function (46) by its linearization. As soon as the linearized constraint system comes into play, we make use of the abbreviation see also Section 5.1. This leads to the linearized system
(48)  
Note that at the time step it holds for which is assumed to be known. For each the linearization is performed at Furthermore, we write and such that the system (48) can be transformed into the LCP problem
(49) 
The system (49) is repeatedly solved within a fixedpoint iteration, until the error between and is small enough and the system converges, yielding the final for a certain
According to the general theoretical framework [1, 37, 30], quadratic convergence is guaranteed when iteratively solving (48). This means that only a few iterations suffice to obtain a solution for (47) at each time step. Within this work we do not pursue the theoretical aspects any further, but rather show a very good agreement with the theory based on our computational examples, see Section 7.
5.4 Computation of the Jacobian
In order to proceed with (48) we need to compute the Jacobian matrix of the function
(50) 
with
(51) 
The derivative matrix is given by
(52) 
Due to the chain rule for the partial derivative of
with respect to it holds:(53) 
In turn, the partial derivative of with respect to results in
(54) 
Finally, inserting (54) into (53), for we obtain
(55) 
For the sake of a more compact representation we introduce the term
(56) 
leading to the matrix
(57) 
Using the new abbreviation the partial derivatives can be computed as follows
(58) 
Finally, the Jacobian matrix can be stated as
(59) 
6 Reduction of the Contact Problem
In many real world applications the numerical integration of constrained systems such as (41) is not possible in real time due to the large number of degrees of freedom. Model order reduction strategies [6] introduce a reduced state with , where is defined by
(60) 
6.1 Computation of ROM
One way to obtain the reduction matrix for the system (41a
) is to use modal reduction. Setting up the eigenvalue problem of (
41a)(61) 
and taking the first eigenvectors, the matrix may be defined by
(62) 
A preferable technique may be the Krylov subspace methods [3, 42]. The subspace is defined by
(63) 
The Krylov base may be computed by the Arnoldi algorithm, which delivers an orthonormal base of the subset. Inserting the reduction (60), where is defined by (62) or (63) into the differential equations system (41a) and multiplying by , one obtains the reduced system
(64) 
where
(65)  
(66) 
are the reduced mass and stiffness matrices. Similar to the full dimensional case, we apply implicit Euler to the reduced system (64), which leads to
(67) 
Note that (67) involves only operations with matrices and vectors in smaller dimensions or . Although , the acting force has only a few entries and for the Lagrange multiplier it holds with .
6.2 Contact Treatment
After the general reduction method has been outlined, we introduce a more accurate way of choosing the basis vectors for the reduced space of the displacements. The underlying idea for increasing the accuracy of a reduction approach for contact problems is to partition the nodal variables into, so called, master and slave nodes. One of the most prominent partitioning methods is the Guyan reduction, also called static condensation, see