An Efficient Model Order Reduction Scheme for Dynamic Contact in Linear Elasticity
AAn Efficient Model Order Reduction Scheme for DynamicContact in Linear Elasticity
D. Manvelyan ∗ , B. Simeon † , U. Wever ‡ February 9, 2021
Abstract
The paper proposes an approach for the efficient model order reduction of dynamic con-tact problems in linear elasticity. Instead of the augmented Lagrangian method that iswidely used for mechanical contact problems, we prefer here the Linear ComplementarityProgramming (LCP) method as basic methodology. It has the advantage of resulting in themuch smaller dual problem that is associated with the governing variational principle andthat turns out to be beneficial for the model order reduction. Since the shape of the contactzone depends strongly on the acting outer forces, the LCP for the Lagrange multipliers hasto be solved in each time step. The model order reduction scheme, on the other hand, isapplied to the large linear system for the displacements and computed in advance by meansof an Arnoldi process. In terms of computational effort the reduction scheme is very ap-pealing because the contact constraints are fully satisfied while the reduction acts only onthe displacements. As an extension of our approach, we furthermore take up the idea of theCraig-Bampton method in order to distinguish between interior nodes and the nodes in thecontact zone. A careful performance analysis closes the paper.
Keywords:
Dynamic contact, linear elasticity, model order reduction, linear complementarityprogram, digital twin technology
The model order reduction for problems in structural mechanics gains rapidly in importance dueto its relevance for simulative operation support [10]. Such support is possible if the simulationmodel runs in parallel to the operation and is synchronized by sensor values at discrete time ∗ Siemens AG, Corporate Technology, Munich , Germany, [email protected] † University of Kaiserslautern, Mathematics, Germany, [email protected] ‡ Siemens AG, Corporate Technology, Munich , Germany, [email protected] a r X i v : . [ m a t h . NA ] F e b oints. The term ”digital twin” has become popular for such a model of the real system. Inparticular, a digital twin allows to monitor the state of a system at any specified position andat any point in time (and thus, also predict the future states). Obvious benefits of digital twinsare, e.g., simplified inspection and service planning, lifetime prediction, advanced fault detectionas well as control and optimization during operation. The challenge, on the other hand, lies inproviding small-scale models that are real time capable and that still preserve the physical keyproperties of the specific application.There is a vast literature on model order reduction for different kinds of partial differentialequations, see [23, 4, 11, 7] for some examples. We concentrate here on reduction methods fordynamic contact problems in linear elasticity. Despite the linearity of the continuum mechanicsmodel for the displacements, such problems are nonlinear in nature due to the unknown movingcontact interface, see, e.g., [13, 17, 25] for background on contact mechanics. The reduced modelshould preserve here both the system dynamics in terms and the shape of the contact area. Notethat for linear elastic models without contact efficient reduction methods exist, among themmodal reduction [4] and Krylov subspace methods[2], and have been integrated into commercialsoftware [1, 15]. The linear mechanical contact problem constitutes a variational problem underunilateral constraints, and its discrete counterpart is mostly solved by the augmented Lagrangemethod [14]. While this method is very flexible with respect to various descriptions of thecontact condition, it may be difficult as a basis for model order reduction. The main issue withcontact problems is that the shape of the contact area is not known a priori, and correspondinglythe reduced equations need to account for this.In [3], a dynamic problem with a linear contact condition is discussed and both the displacementsand the Lagrange multipliers are reduced using methods such as the singular value decompositionor the non-negative matrix compression in order to find an online solution for a reduced saddle-point problem by means of the Lagrange multiplier method. In contrast, our approach reducesonly the primal displacements. In order to compute the projection matrix only the systemmatrices and the position vector of the applied force are required. This offline procedure isperformed only once. A major benefit is the fact that the trajectories of the full system arenot required. Moreover, the computed projection matrix can be re-used for a system with adifferent load that acts at the same position and in the same direction. After the reducedcontact problem is defined, we switch over from the space of reduced displacements to the spaceof the dual Lagrange multipliers, i.e., to the adjoint problem of the original variational problem.This leads to a Linear Complementary Programming (LCP) problem [8, 9, 19]. The complexityof this LCP problem scales with the number of constraints, which is usually small compared tothe number of displacements and thus, may be solved even in real time. The resulting update ofthe displacements is again performed in a reduced space. We think that the degree of reductionis optimal in the sense that the shape of the contact interface is preserved.For a further separation of the displacement degrees of freedom into interior and contact nodeswe additionally apply the method of Craig-Bampton [21]. It turns out that this separationprovides advantages in accuracy in some cases. The paper is organized as follows: In Section 2we recall the underlying equations of dynamic contact in linear elasticity. Section 3 describesthe solution of the contact problem by Linear Complementary Programming (LCP) both inthe static and the dynamic case. Section 4 describes the reduction of the problem by Krylovmethods and the extension for the Craig-Bampton method. The performance of the approach2s demonstrated on 2D geometries in Section 5. In this section we summarize the governing equations for dynamic contact in linear elasticity.We first define a generic model problem in semi-discretized form that applies to several differentfinite element discretizations and then provide some exemplary background.
We consider large-scale systems of differential-algebraic equations with unilateral constraints inthe form M ¨ q ( t ) + Kq ( t ) = f ( t ) + C T λ , (1) Cq ( t ) + b ≥ , λ ≥ , λ ∗ ( Cu ( t ) + b ) = 0 (2)where q ( t ) ∈ R N denotes the vector of nodal displacement variables, M ∈ R N × N the symmetricpositive definite mass matrix, K ∈ R N × N the symmetric positive semi-definite stiffnes matrix,and f ( t ) ∈ R N the load vector. Furthermore, C ∈ R m × N stands for the constraint matrix and b ∈ R m for a given initial clearance vector. The inequality signs in (2) are to be interpretedcomponent-wise and stand for the discretized non-penetration conditions of all variables in thecontact interface. Correspondingly, the vector of Lagrange multipliers λ ∈ R m enforces theconstraint in the dynamic equation (1). The multiplier λ , standing for the discretized contactpressure, is always positive and its inner product with the constraints satifies a complementaritycondition.In this paper, we will introduce a model order reduction method for the semi-discretized equa-tions (1) that projects the nodal variables q onto a much smaller space while explicitly preservingthe constraints and the Lagrange multipliers. For this purpose, the structure of the equationsand the properties of the matrices as stated above are an appropriate starting point. Frictionless,adhesive-free normal contact in combination with small deformation theory and a linear-elasticmaterial will lead to (1) if(i) the Lagrange multiplier method is used to enforce the non-penetration condition and(ii) the finite element method on matching meshes, i.e., node-to-node contact, is applied forthe discretization in space.Approaches like the augmented Lagrange method or the Nitsche method [25] lead to modifi-cations that will not be considered here. The same holds for additional friction effects andnode-to-segment contact.For a contact problem with k bodies that fits into the framework described so far, the massmatrix M will consist of k blocks on the diagonal that stem from the discretizations of theindividual bodies, and the same block-diagonal structure applies to the stiffness matrix K . Theconstraint matrix C will then typically exhibit a sparse structure where each row stands for3 node-to-node contact condition. We do not dive further into the details of various contactmodels and discretization schemes and refer instead to the standard references [17, 25]. But forthe sake of a self-contained presentation, we shortly sketch the setting of the classical obstacleproblem in the dynamic case. Assuming as before frictionless, adhesive-free normal contact and a linear-elastic material as wellas linear kinematics, we can express the dynamic contact problem in terms of the displacementfield u ( x , t ) ∈ R d and the contact pressure p ( x , t ) ∈ R where x ∈ Ω ⊂ R d , d = 2 or d = 3, isthe spatial variable and t ∈ [ t , T ] stands for the temporal variable. The boundary of the elasticbody is decomposed into ∂ Ω = Γ D ∪ Γ N ∪ Γ C where the latter represents the contact interface.Using the outward normal vector n ( x ) on the contact interface, the distance between the bodyand a given surface is described by the scalar gap function g ( x , u ) = (cid:0) ξ ( x ) − ( x + u ( x , t ) (cid:1) T n ( x ) , x ∈ Γ C . (3)Here, the point ξ ( x ) on the obstacle’s surface is obtained by projection of x in outward normaldirection. The contact problem in strong form is then given by the dynamic equations ρ u tt − div σ ( u ) = F in Ω × [ t , T ] (4)subject to standard boundary conditions σ ( u ) · n = τ on Γ N × [ t , T ] , u = 0 on Γ D × [ t , T ] (5)and furthermore subject to the contact conditions g ≥ , p ≥ , g · p = 0 on Γ C × [ t , T ] . (6)As initial data, u ( · , t ) = u and u t ( · , t ) = v (7)are prescribed with given functions u and v , respectively.In (4) and (5), the constant ρ denotes the mass density, F ( x , t ) ∈ R d the volume force, τ ( x , t ) ∈ R d the surface traction, and σ ( u ) ∈ R d × d the stress tensor given by σ ( u ) = E ν e ( u ) + νE (1 + ν )(1 − ν ) trace( e ( u )) I (8)with Young’s modulus E ≥
0, Poisson’s ratio − ≤ ν ≤ . e ( u ) = 12 ( ∇ u + ∇ u T ) ∈ R d × d . (9)In order to pass to a weak formulation of the obstacle problem, we introduce the function spaces V := (cid:110) v ∈ H (Ω) d : v = 0 on Γ D (cid:111) , (10) L := (cid:26) µ ∈ H / (Γ C ) (cid:48) : (cid:90) Γ C µw ds ≥ ∀ w ∈ H / (Γ C ) , w ≥ (cid:27) (11)4nd the abstract notation (cid:104) ρ u tt , v (cid:105) := (cid:90) Ω v T ρ u tt dx, a ( u , v ) := (cid:90) Ω σ ( u ) : e ( v ) dx, (cid:104) (cid:96), v (cid:105) := (cid:90) Ω v T F dx + (cid:90) Γ N v T τ ds. (12)The non-penetration condition in weak form with test function µ ∈ L is recast as0 ≤ (cid:90) Γ C gµds = (cid:90) Γ C g ( µ − p ) ds = − (cid:90) Γ C u T n ( µ − p ) ds + (cid:90) Γ C ( ξ − x ) T n ( µ − p ) ds. (13)The last two integrals give rise to the definitions of the bilinear form b c ( v , p ) := − (cid:90) Γ C v T n · p ds (14)on V × L and the linear form on
L(cid:104) m, p (cid:105) := (cid:90) Γ C ( ξ − x ) T n p ds. (15)Using these definitions, the weak form of the dynamic contact problem is stated as follows: Foreach t ∈ [ t , T ] find the displacement field u ( · , t ) ∈ V and the contact pressure p ( · , t ) ∈ L suchthat (cid:104) ρ u tt , v (cid:105) + a ( u , v ) = (cid:104) (cid:96), v (cid:105) + b c ( v , p ) for all v ∈ V ,b c ( u , µ − p ) + (cid:104) m, µ − p (cid:105) ≥ µ ∈ L . (16)The system (16) is discretized with respect to the spatial varibale by applying the standardGalerkin projection u ( x , t ) . = (cid:80) Ni =1 φ i ( x ) q i ( t ) = Φ ( x ) q ( t ) with basis functions φ i and nodalvariables q = ( q , . . . , q N ) to the displacement field. If the contact condition is simply expressedin terms of the distance between a node x j ∈ Γ C , i = 1 , . . . , m , and the obstacle, the integralover Γ C in the weak dynamic equation is replaced by a sum (cid:90) Γ C v T n · p ds . = m (cid:88) j =1 A j v ( x j ) T n ( x j ) · p ( x j , t ) , A j : area around node x j . (17)In the same way, the unilateral constraint is discretized. Putting finally λ ( t ) := ( p ( x , t ) , . . . , p ( x m , t ))as discrete pressure variable, the semi-discretized equations (1) and (2) follow in the usual way,with the matrix C ∈ R m × N being an indicator matrix for the nodes in the contact interfaceand b ∈ R m the offset of the contact. Note that only a few displacements are involved for thecontact condition and hence most of the entries of the matrix C are zero. A popular approach to solve the mechanical contact problem is the augmented Lagrange method[20, 14, 25], which is very flexible to different descriptions of the contact interface. However,in the context of model order reduction we advocate a formulation as linear complementarityprogramming (LCP) problem with corresponding solution methods. We treat first the staticcase and proceed then to the dynamic problem.5 .1 LCP for the Static Contact Problem
In the stationary case, we can recast the generic semi-discretized formulation (1) and (2) asminimization problem min q ∈ R N , λ ∈ R m + q T Kq − q T f − λ T ( Cq + b ) . (18)Note that the requirement λ ∈ R m + enforces the positivity of the contact pressure and, simulta-neously, the non-penetration condition. The corresponding KKT-conditions read Kq − f − C T λ = 0 , Cq + b ≥ , λ ≥ , λ ∗ ( Cq + b ) = 0 . (19)Next, assuming a positive-definite stiffness matrix, we eliminate the displacements q from theKKT-conditions (19) via q = K − ( f + C T λ ) (20)and insert this expression into the constraint equations. This results in CK − f + CK − C T λ + b ≥ , λ ≥ , λ ∗ ( CK − f + CK − C T λ + b ) = 0 . (21)Using the abbreviations A := CK − C T ∈ R m × m , B := CK − f + b ∈ R m , (22)the equations (21) read B + Aλ ≥ , λ ≥ , λ ∗ ( B + Aλ ) = 0 . (23)The LCP problem (23) may be solved by standard methods from constrained optimization,see below for more details. Note that tasks like the detection of active contact nodes are thustransferred to the LCP solver. In the transient case, the dynamic equations (1) must be discretized in time. As straightforwardand unconditionally stable method, we apply the implicit Euler scheme, i.e., the first orderbackward differentiation formula. For this purpose, the second order derivative is replaced bythe finite difference ¨ q ( t + h ) ≈ h (cid:0) q ( t + h ) − q ( t ) + q ( t − h ) (cid:1) (24)with time stepsize h . Note that the finite difference ˙ q ( t + h ) ≈ ( q ( t + h ) − q ( t )) /h for thevelocity is hidden in (24), and due to the second order time derivative, the implicit Euler leadshere to a two-step method. Since (24) possesses first order of accuracy only, one can apply otherintegration schemes instead, e.g. the generalized- α method, which is very common in structuralmechanics, cf. [24]. 6n this paper, however, we continue working with the implicit Euler method as our main interestfocuses on the model order reduction scheme. Moreover, dynamic contact involves frequentdiscontinuities that might impair the convergence of a higher order time integration method.Inserting (24) into (1) leads to M ( q ( t + h ) − q ( t ) + q ( t − h )) + h Kq ( t + h ) = h f ( t + h ) + h C T λ ( t + h ) (25)Assuming that the previous time steps are known, (25) may be resolved for q ( t + h ), q ( t + h ) = ( M + h K ) − ( h f ( t + h ) + h C T λ ( t + h ) + 2 M q ( t ) − M q ( t − h )) . (26)Now the same procedure as in the static case may be applied. Inserting (26) into the constraints(2) again leads to an LCP. With the definitions A := h C ( M + h K ) − C T , (27) B := C ( M + h K ) − ( h f ( t + h ) + 2 M u ( t ) − M u ( t − h )) + b (28)we again arrive at the LCP (23), which here has to be solved in each time step.We remark that the two-step time discretization requires initial values q ( t ) and q ( t + h )to start. E.g., if q and ˙ q as initial displacement and velocity are given, one can compute q ( t + h ) = q + h ˙ q by an explicit Euler step and then continue with the two-step formula(25). The Lagrange multiplier λ , on the other hand, does not require an initial value, andit is computed in each time step as an implicitly given function of q . In the terminology ofdifferential-algebraic equations, this means that q stands for the differential variables while thealgebraic variables λ possess no memory. If the unilateral constraints were replaced by equalityconstraints Cq + b = 0, the index of the resulting differential-algebraic equation would equal 3,which means that special care must be taken for the time integration [6, 16, 24].Another remark concerns the treatment of impact situations. In the presented time discretiza-tion, impacts are not specifically identified by means of computing an impact velocity and acorresponding coefficient of restitution. Instead, the LCP solution in each time step implicitlyenforces a plastic impact solution that might lead to a loss in kinetic energy. The core of our proposed solution algorithm consists of the efficient solution of the LCP problem(23). One way for solving LCPs is the usage of so-called NCP-functions φ NCP . An NCP-function φ NCP : R → R is characterized by the property φ NCP ( a, b ) = 0 ⇔ a ≥ , b ≥ , ab = 0 . (29)Two examples of NCP-functions are: φ NCP1 ( a, b ) = min( a, b ) , (30) φ NCP2 ( a, b ) = (cid:112) a + b − a − b. (31)7he first function (30) is not suitable for computational purposes because of its non-smoothness.The second NCP-function is also known as the Fischer-Burmeister-function [12]. Thus for solvingan LCP the following two formulations (32) and (33) are equivalent: B + Aλ ≥ , λ ≥ , λ ∗ ( B + Aλ ) = 0 . (32) ⇐⇒ φ NCP2 ( B + Aλ , λ ) = 0 . (33)Equation (33) may be solved by Newton’s method, but here we apply a much more efficientmethod in terms of Lemke’s algorithm (see [8], Sect. 4.4.5). Specifically, we use the Pythonimplementation of the algorithm [18]. In many applications, the number N of degrees of freedom of the discretized system (1) is largeand its numerical integration is not possible in real time. Model order reduction strategiesintroduce a reduced state v ∈ R n r with n r (cid:28) N , where v is defined by q = Qv , Q ∈ R N × n r . (34)We can interprete the map from q to v as a second Galerkin projection that acts on top of thefirst projection from the displacement field u to q . For a survey on model order reduction see,eg., [5]. One way to obtain the reduction matrix Q is to use modal reduction. Setting up the eigenvalueproblem in the unconstrained case ω M y = Ky (35)and taking the first n r eigenvectors, the matrix Q is then defined by Q = (cid:8) y , · · · , y n r (cid:9) . (36)A more recent technique are the Krylov subspace methods [2, 22]. The reduction is defined by Q = (cid:8) K − ω f , K − ω M K − ω f , · · · , ( K − ω M ) n r − K − ω f (cid:9) (37)The Krylov vectors can be computed by the Arnoldi algorithm, which delivers an orthonormalbasis of the subspace. Inserting the reduction (34), where Q is defined by (36) or (37), into thedifferential equation (1) and multiplying by Q T from the left, we obtain the reduced equationˆ M ¨ v + ˆ Kv = Q T (cid:0) f ( t ) + C T λ ( t ) (cid:1) (38)where ˆ M = Q T M Q , ˆ K = Q T KQ ∈ R n r × n r are the reduced matrices.8sing (38), moreover the reduced time discretization follows easily as v ( t + h ) = (cid:16) ˆ M + h ˆ K (cid:17) − (cid:16) h Q T f ( t + h ) + h Q T C T λ ( t + h ) + 2 ˆ M v ( t ) − ˆ M v ( t − h )) (cid:17) (39)We stress that (39) contains only operations with complexity O ( n r ) or O ( m ). Though Q ∈ R N × n r , the acting force f ( t + h ) has only a few entries and the Lagrange multiplier λ is still ofdimension m . Next, we introduce a modification of our approach which is inspired by the fact that often thedominant influence on the dynamics of the system stems from the nodes located on the contactinterface. In other words, we distinguish between the nodal variables on the contact interface,so-called master nodes and the ones outside of the contact interface, the so-called slave nodes.We adopt this technique from structural dynamics, see Craig and Bampton [21].Let q = ( q M , q S ) T ∈ R N the permuted solution vector such that we distinguish between masterand slave nodes. Then the corresponding permuted matrices and the permuted force vector aregiven by M = (cid:18) M MM M MS M SM M SS (cid:19) , K = (cid:18) K MM K MS K SM K SS (cid:19) , F = (cid:18) F M F S (cid:19) . (40)Due to the fact that the constraint matrix affects only the dynamics of the slave nodes it holds C = (cid:18) C M C S (cid:19) , C S = 0 . (41)The constraint matrix refers only to the master nodes. A crucial assumption is now that theinfluence acceleration term for the slave nodes in the dynamic equation can be neglected, andthe same for the force term. This yields a coupling equation between the master and slave nodevia K SM q M + K SS q S = 0 (42)Since we want to preserve the contact nodes and only reduce the slave nodes, we keep thestructure of the master nodes fixed, i.e, u M = 0 . Then due to (41) the dynamical system isgiven by M SS ¨ q S + K SS q S = F S (43)Instead for the full system, we compute the transformation matrix for the slave system (43)using the Arnoldi method, which yields Q S := (cid:8) K − SS F S , K − SS M SS K − SS F S , · · · , ( K − SS M SS ) n − K − SS F S (cid:9) . (44)Using (42), we obtain the complete transformation matrix that includes the relation betweenthe master and slave nodes, Q CB = (cid:18) I M − K − SS K SM Q S (cid:19) (45)9he transformation matrix Q CB with the partitioning of master and slave nodes in the fashionof the Craig-Bampton method provides an alternative for reducing the transient contact problem(39). Regardless whether we employ Q from the reduction (37) or Q CB , the computation of thetransformation can be done in an a priori offline phase and does not require solution trajectories,i.e., snapshots, of the full system. In this section we apply our approach to 2D problems with self-contact. The spatial discretizationuses standard bilinear shape functions on quadrilateral elements. The finite element methodand the reduction scheme are implemented by means of a custom Python script while Lemke’salgorithm is taken from an open source library [18].
We consider the unit square Ω = [0 , × [0 ,
1] under the assumption of plane stress. Theparameters in dimensionless form are ρ = 1 , E = 1000 , ν = 0 .
3. The left edge is fixed by meansof zero Dirichlet boundary conditions. We discretize the example by 1600 quadrilateral elements,which lead to N = 3386 DOF. An interior contact interface is predefined by a fixed number of m discretization points where each such point is represented by double nodes. This kind of datastructure allows us to distinguish the nodes placed on the contact interface from the remainingnodes. Note that the node-to-node non-penetration condition acts only in x -direction while theDOF in y -direction may move freely.An oscillating force acts on the right side of the domain which affects the contact interface byopening and closing the tear over time. The problem setup and the mesh are illustrated inFigure 1. Two sensors are placed for the data extraction in order to compare the displacementsof the full and the reduced model. Before starting the reduction method we eliminate the fixeddegrees of freedom due to the Dirichlet boundary condition. Overall, there are 4 m displacementvariables associated with the contact nodes due to the relation n c = 4 m First Example
In our first example the contact tear is located on the upper-right side of the domain. Thecontact is discretized by m = 12 points and there are n c = 48 contact variables respectivly. Weconsider an oscillating nodal force in horizontal direction (cf. Figure 1) with the load magitudegiven by f ( t ) = (cid:18) . . πt )0 (cid:19) , t ∈ [0 , . (46)We test our approach for three different numbers n r = 5 , ,
15 of the Krylov basis vectors.The resulting trajectories for the displacement in the first sensor node in the contact interface,Figure 2, show very good agreement of the reduced order model with the full order model.10igure 1: The depicted deformed state with an open contact tear occurs at the time point whenthe acting horizontal nodal load vector is pointing outward from the domain.As expected, the approximation quality increases with the dimension of the Krylov subspace.Nevertheless, the displacement period seems to be recovered even for small numbers of n r . Theresults for the second sensor node are also very satisfactory, as presented in Fig. 3.In our tests we also applied an outer force with the same magnitude as (46) but different anglegiven by f ( t ) = (cid:18) . . πt ) sin(0 . π )1 . . πt ) sin(0 . π ) (cid:19) . (47)In this case the y -displacement shows a stronger response. The comparison of the displacementsof the FOM and ROM is depicted in Figure 4 and 5 corresponding to the first sensor and to thesecond sensor, respectively. 11igure 2: The x - and y -displacement of the full order model (FOM) and the reduced order modelfor n r = 5 , ,
15 basis vectors. The displacement trajectory is evaluated in the first sensor nodelocated on the contact interface.Figure 3: The x - and y -displacements of the full order model (FOM) and the reduced ordermodel for n r = 5 , ,
15 basis vectors. The displacement trajectories are evaluated in the secondsensor node located on the contact interface. 12igure 4: x - and y -displacements for the FOM (full order model) and for the ROM (reducedorder model) with three different numbers of basis vectors, i.e., n r = 5 , , . First sensor nodeand load (47).Figure 5: x - and y -displacements for the FOM (full order model) and for the ROM (reducedorder model) with three different numbers of basis vectors, i.e., n r = 5 , , . Second sensornode and load (47) 13 .1.2
Second Example
Next, we consider another example where the contact interface is placed inside the domain, seeFig. 6. The corresponding displacements can be found in Fig. 7.Figure 6: Simulation scenario with interior tear. The depicted state is captured at the timeinstant when the horizontal nodal load is directed outside the domain, opening the interiorcontact tear.Figure 7: The x - and y -displacements of the FOM and ROM for n r = 10 , ,
30 basis vectorsrespectvely. The displacement trajectories are evaluated on the first sensor node located on theinner contact interface and correspond to the acting force (46) and Figure 6.14n contrast to the examples before, the double nodes appear in the interior of the domain. Theouter boundary and the surrounding inner area are represented by single nodes. The stiffnessof the surface and connectedness of the nodes produces an inner tear. Again, the reduced ordermodel performs well. At least n r = 20 basis vectors are required to obtain satisfactory results. We continue with a comparison of the reduction results with and without Craig-Bampton split-ting, which was discussed in Subsection 4.2. To guarantee a fair comparison of the two methods,we choose the same number of basis functions in both cases. In case of using no splitting, thebasis consists of only Krylov vectors, whereas in the splitting case, the basis consists of a fewKrylov vectors ( n k ) and additionally the contact nodes ( n c ), i.e., n r = n k + n c . The Figure 8shows the comparison of both reduction approaches with the full model, referring to the scenarioof Fig. 1. We can observe an increase of accuracy in case of using the splitting of Craig-Bampton.This holds also for the second sensor node and the displacements in y -direction.Figure 8: Upper figure: x -displacements of the FOM and the ROM, once without (Arnoldi) andonce with splitting (CB). In both cases, n r = 55 . In case of CB n r = n k + n c with n k = 7 and n c = 48 . Trajectory of first sensor node, which is placed on the contact interface. Lower figure:Lagrange multiplier in contact node 6.Furthermore, the Lagrange multipliers λ ( t ) exhibit an interesting behavior. We observe thatROM without splitting produces contact pressures that are quite different from those of theFOM, while the Craig-Bampton splitting yields values that agree well with the FOM. Bothapproaches guarantee positivity of the Lagrange multiplicator whenever the contact is activated.An explanation of this deviation with respect to the Lagrange multipliers in case of no splittingis as follows. The projection matrix Q in (34) alters the nodewise distribution of the pressure15n the contact variables considerably and involves spatial basis functions with large support.Splitting of variables, on the other hand, by means of Q CB in (45) preserves the contact variablesand thus the Lagrange multipliers. In Figure 8 the x -displacement in the sensor node 1 forboth reduction approaches and the FOM is displayed, along with the corresponding Lagrangemultiplier. The switching between contact, where the multiplier is positive, and positive gapwith vanishing multiplier can be nicely observed. In this paper we have presented a reduction method for dynamic contact problems in linearelasticity. The generic model elaborated in Section 2 does not cover all available methods butnevertheless is general enough to provide sufficient insight into the problem class. We planto extend the problem class in the future to situations like node-to-segment contact and non-matching meshes where the unilateral constraints might become nonlinear.The proposed model order reduction uses a Krylov subspace approach for the displacementvariables but preserves the Lagrange multipliers. This is intimately connected with the LCPproblem that we advocate for solving the contact problem in each time step. Obviously, thesavings in DOF are substantial as long as the number of contact constraints is small comparedto the number of displacement variables. If we consider solids and contact on outer or interiorboundaries, this holds true, but if we extend the method to shell problems where contact mightoccur all over the computational domain, the savings might be less substantial.It also turned out that the numerical approximation of the Lagrange multipliers, i.e., the contactpressure in the contact interface, is in general not accurate for the simple Krylov reductionscheme. By splitting of the nodes in the fashion of the Craig-Bampton method, a substantialimprovement can be achieved. Moreover, this modification leads also to better approximationsof the displacement variables. In any case, the reduction matrix is computed only once in anoffline procedure. There is no need to compute trajectories or snapshots of the full model toprovide data for the reduction scheme.Our future effort will focus on the development of methods for treating more complex geometriesand on the analysis of the proposed method class in terms of apriori and a posteriori errorestimates.
Acknowledgement : The authors would like to thank Meinhard Paffrath for his guidance onthe technical matters concerning the implementation. Furthermore, the authors want to thankChristoph Heinrich and his analysis group for valuable discussions.; 16 eferences
Applied Numerical Mathematics , 43:9–44, 2002.[3] M. Balajewiec, D. Amsallem, and C. Farhat. Projection-based model reduction for contactproblems.
Int. J. Numer. Meth. Engng , 106:pp. 644–663, 2015.[4] U. Baur, P. Benner, and L. Feng. Model order reduction for linear and nonlinear systems:a system-theoretic perspective. , MPIMD/14-07,2014.[5] U. Baur, P. Benner, and L. Feng. Model order reduction for linear and nonlinear systems:a system-theoretic perspective.
Archives of Computational Methods in Engineering , 2014.[6] K. E. Brenan, S. L. Campbell, and L. R. Petzold.
The Numerical Solution of Initial ValueProblems in Ordinary Differential-Algebraic Equations . SIAM, Philadelphia, 1996.[7] S. Chaturantabut and D.C. Sorensen. Discrete empirical interpolation for nonlinear modelreduction. 32:4316 – 4321, 01 2010.[8] R. W. Cottle, J. S. Pang, and Richard E. Stone.
The Linear Complementarity Problem .SIAM, 2009.[9] R.W. Cottle, J.S. Pang, and R.E. Stone. The linear complementarity problem.
San Diego,CA: Academic , 1992.[10] M. Eigner, T. Gilz, and R. Zafirov. Interdisziplin¨are produktentwicklung - modellbasiertessystems engineering.
PLM Portal , 2012.[11] C. Farhat, P. Avery, T. Chapman, and J. Cortial. Dimensional reduction of nonlinearfinite element dynamic models with finite rotations and energy-based mesh sampling andweighting for computational efficiency.
International Journal for Numerical Methods inEngineering , 98(9):625–662, 2014.[12] A. Fischer. A special newton-type optimization method.
Optimization , 24:pp. 269–284,1992.[13] A. Francavilla and O.C. Zienkiewicz. A note on numerical computation of elastic contactproblems.
International Journal for Numerical Methods in Engineering , 9:pp. 913–924,1975.[14] P. Gill, W. Murry, and M. Wright.
Practical Optimization
The Numerical Solution of Differential-AlgebraicEquations by Runge-Kutta Methods . Lecture Notes in Mathematics Vol. 1409. Springer,Heidelberg, 1989.[17] N. Kikuchi and J.T. Oden.
Contact Problems in Elasticity . SIAM, 1988.[18] A. Lamperski. Python: lemkelcp 0.1. https://pypi.org/project/lemkelcp/.[19] C. X. Li and S. L. Wu. A note on the unique solution of linear complementarity problem.
Cogent Mathematics , page DOI: 10.1080/23311835.2016.1271268, 2016.[20] M.J.D. Powell.
A method for nonlinear constraints in minimization problems . AcademicPress, in optimization ed. by r. fletcher edition, 1969.[21] D. J. Rixen. A dual craig–bampton method for dynamic substructuring.
Journal of Com-putational and applied mathematics , 168(1-2):383–391, 2004.[22] B. Salimbahrami and B. Lohmann. Order reduction of large scale second-order systemsusing krylov subspace methods.
Linear Algebra and its Applications , 415:385–405, 2005.[23] S. Salimbahrami and B. Lohmann. Order reduction of large scale second-order systemsusing krylov subspace methods.
Linear Algebra and its Applications , 415:385–405, 2006.[24] B. Simeon.
Computational flexible multibody dynamics . Springer, 2013.[25] P. Wriggers.