Explicit dynamics in impact simulation using a NURBS contact interface

In this paper, the impact problem and the subsequent wave propagation are considered. For the contact discretization an intermediate non‐uniform rational B‐spline (NURBS) layer is added between the contacting finite element bodies, which allows a smooth contact formulation and efficient element‐based integration. The impact event is ill‐posed and requires a regularization to avoid propagating stress oscillations. A nonlinear mesh‐dependent penalty regularization is used, where the stiffness of the penalty regularization increases upon mesh refinement. Explicit time integration methods are well suited for wave propagation problems, but are efficient only for diagonal mass matrices. Using a spectral element discretization in combination with a NURBS contact layer the bulk part of the mass matrix is diagonal.

used, where the geometry and its discretization coincide and the higher continuity delivers advantages for the enforcement of the contact constraints, see De Lorenzis et al. 8,9 However, the generation of an isogeometric discretization for complex 3D problems still has not reached the same level of automation as in the standard finite element (FE) analysis. Promising approaches are given in References 10,11 using triangular isogeometric discretizations.
According to Belytschko et al, 12 a hyperbolic system of partial differetial equations (PDEs) can be classified as a structural dynamics problem or a wave propagation problem. A structural dynamics problem is characterized by a frequency spectrum of the input, which is below the resolution limits of the mesh. Classical examples of a structural dynamics problem are seismic response or vibration problems. In a wave propagation problem, high frequency waves are part of the spectrum. The impact event is a classical example, where the sudden collision leads to high frequency waves in the material. While implicit methods are preferred for structural dynamics problems benefiting from a time step that is significantly larger than the critical time step, explicit methods are preferred for wave propagation problems due to the small time step, which is required to accurately track the high frequency waves in the body. Using explicit time integration for systems of ordinary differential equations (ODEs) with a diagonal mass matrix, no solution of a linear system is required and the explicit time integration becomes more efficient for wave propagation problems. Using a standard FE discretization, a lumping technique like the Hinton-Rock-Zienkiewicz 13 mass lumping is used to obtain a diagonal mass matrix. Alternatively, the spectral element discretization can be used, where the Lobatto integration points and FE nodes coincide. [14][15][16] This naturally provides a diagonal mass matrix. Using an isogeometric discretization along with explicit time integration, several techniques have been proposed like collocation Auricchio et al 17 or lumping procedures. [18][19][20] However, using higher-order isogeometric spatial discretizations a lumping procedure for the mass matrix gives only a limited accuracy, thus a modification of both the stiffness and mass matrices is required. Finally, for the wave propagation problem, the interaction between the spatial and time discretizations needs to be considered. Using the popular central difference method (CDM) with linear FEs for impact simulation results in apparently accurate results due to the compensation of spatial and time discretization errors, as shown in Otto et al. 3 Higher frequencies are only captured using smaller time steps, which leads to oscillations and increased numerical error.
The aim of this paper is to combine several methodologies to accurately simulate 2D/3D impact problems. The nonlinear and mesh-dependent penalty regularization for 1D impact problems, see Otto et al, 3 is extended to 2D and 3D impact problems. The coupled FE-NURBS approach, see Otto et al, 21 is extended to dynamic contact with a zero mass interface layer. The global system matrices, in particular the mass matrix is derived for this approach. The coupled approach uses an intermediate NURBS layer to compute the contact forces between the contacting bodies discretized by FEs (Figure 1). The bulk part is discretized by spectral elements, as a consequence, the advantages of a smooth isogeometric contact formulation are combined with efficient spectral element formulations that are simple to generate as well as very efficient for wave propagation problems. Mesh tying conditions are used to tie the NURBS layer to one of the bodies. The contact conditions are then enforced between the NURBS layer and the other body. Thus, the NURBS layer is the master side for both the mesh tying and the contact formulation. Smooth interelement boundaries and a smooth normal field can be provided. A complex 3D NURBS geometry is not needed and the tying conditions are computed only once and remain constant throughout the simulation. Using a well-established spectral element discretization, a slave body for tying contact mesh tying diagonal mass matrix is obtained for the bulk part of the discretization. The interpolation order of the spatial discretization can be easily increased and a higher order explicit time integration scheme, the classical Runge-Kutta (RK) method, is used. This paper is organized as follows: In Section 2, the contact conditions for frictionless contact and its variational formulation are summarized. Afterwards, the spatial FE and NURBS discretizations are discussed in Section 3. The contact discretization is described in Section 4 followed by the coupled FE-NURBS approach for impact simulation in Section 5. The explicit RK time integration for the simulation of wave propagation problems is given in Section 6 and the mesh-dependent penalty regularization in Section 7. Finally, the performance of the method is discussed for selected 2D and 3D examples in Section 8.

FRICTIONLESS CONTACT FORMULATION
A linear elastic material with Young's modulus E and Poisson's ratio is used. The contact problem is defined between a slave body and a master body for contact, Ω i with i = {s,m}. The deformation field of the slave and master bodies is described by where i is the mapping from the initial configuration X i to the current configuration x i with the displacements u i . Since small deformations are assumed, the strain is defined by The normal gap is given by where x s is a position on the slave contact zone and x m denotes the closest point to x s on the master side. The normal at x m is denoted by n m . For an elastic material, the frictionless contact problem can be written as a constrained energy minimization problem A sufficient condition for the equilibrium is given by the following equations where W kin and W int denote the kinetic and internal virtual work, respectively, W ext the external virtual work, W c the contact contribution to the overall virtual work, and u i the kinematically admissible variation of the displacement fields. The conditions in Equation ( 5b) are called the Karush-Kuhn-Tucker conditions, where the normal contact traction p N is an additional unknown in the system and is interpreted as a Lagrange multiplier. The virtual work of the normal contact traction is given by where Γ c denotes the potential contact surface on the slave body in the reference configuration. For detailed treatment, we refer to the monographs. 22,23 In this paper, the penalty method is used to regularize the contact conditions (5b). The normal contact traction p N is computed by The normal contact traction p N is a power function of the normal gap g N with degree q and factor a. The impenetrability condition g N ≥ 0 is enforced only approximately, since the contact traction increases with the penetration. The stiffer the function p N (g N ) is, the smaller is the penetration.

SPATIAL DISCRETIZATION
In the proposed spatially discretized model, a coupled FE-NURBS approach is used. 21 The NURBS layer is coupled to the FE mesh by mesh-tying conditions. The NURBS layer is the master side for mesh tying and contact conditions. The bulk part of the geometry is discretized by spectral FEs. In this section, first the FE discretization is discussed, followed by the description of the NURBS layer.

FE discretization
The spectral FE method was originally proposed by Patera 14 for fluid dynamic problems. As for the standard FE discretization, the domain is discretized by nonoverlapping elements. The shape functions of the spectral elements are computed by Lagrange interpolation at the nodes that coincide with the integration points of a Gauss-Lobatto quadrature scheme. Due to the identical positions of FE nodes and quadrature points an orthogonal relation for the shape functions is achieved: where k ( p ), k = 1, …, n are the shape functions of the spectral FE with n nodes and p denotes the integration point coordinates on the reference bi-unit domain. This orthogonality of the shape functions at the quadrature points results in a diagonal element mass matrix M e : The diagonal element mass matrices M e are assembled to a diagonal global mass matrix M. A Gauss-Lobatto quadrature with n integration points is accurate for polynomials up to degree 2n − 3, thus the diagonal mass matrix is the result of underintegration. Nevertheless, the total mass is conserved and the spectral convergence of the method is still ensured, see Maday et al. 24 Using an explicit time integration, a diagonal mass matrix M is advantageous in terms of computational efficiency.

Discretization with NURBS
The NURBS contact layer is represented by a curve for 2D problems or by a surface for 3D problems. In this section, the main terminology and notations are introduced. A more detailed description of NURBS can be found, for example, in Piegl and Tiller. 25 For mechanical problems, the first spatial discretization with NURBS was introduced by Hughes et al. 26 The integration of IGA in existing FE codes is described in Cottrell et al 27 and Borden et al. 28

NURBS basis functions
The NURBS basis functions are determined by the weights w j , j = 1, …, n, the degree p and the knot vector , with where the first knot is assumed to be 1 = 0 and the last knot to be n+p+1 = 1. The first and last p + 1 knots in are equal. The nonvanishing spans in the interior knot vector { p+1 , …, n+1 } can be interpreted as elements. In the case p+1 < p+2 < … < n < n+1 , the number of elements is n − p. One repeated interior knot reduces the number of elements by one. The NURBS basis functions are given by where the B-Spline basis functions are denoted by N i , i = 1, …, n. Efficient evaluation methods of N i are presented in Piegl and Tiller. 25 The local support of the NURBS shape functions include p + 1 elements or knot spans. This increases the bandwidth of the stiffness matrix if p > 1, but results in a  p−1 continuity of the basis functions.

NURBS curves and surfaces
A NURBS curve is given by where x i , i = 1, …, n are the control points. Since the first and last p + 1 knots in are assumed to be equal, the first and last control points are interpolated. The parametric coordinates of a NURBS surface are condensed in the vector A NURBS surface has a knot vector in each parametric direction d = 1, 2: where p d denotes the corresponding degree and n d the number of control points. The control points and weights are denoted by x kl and w kl , where k = 1, …, n 1 and l = 1, …, n 2 . Bivariate NURBS basis functions R kl are defined by a tensor product of the univariate basis functions in Equation (11): Finally, a NURBS surface is given by

Generation of the NURBS layer and the normal projection
In this work, an auxiliary NURBS layer is added to the contact problem of two FE bodies. Mesh-tying conditions are used to couple the NURBS layer with one of the FE bodies and contact conditions are applied between the NURBS layer and the second FE body. As presented in Otto et al, 21 the NURBS layer can be constructed in several ways. For simple geometries of the FE tying interface, as, for example, circles, arcs, lines, a NURBS representation can be found with a small number of elements. This NURBS representation is adapted to the FE mesh using the knot insertion technique. A matching parametrization of the NURBS is achieved, where each FE boundary in the tying zone is projected on exactly one NURBS element. Another possibility is generating the NURBS layer by interpolation of the FE nodes in the tying zone. The algorithm for NURBS interpolation is described in Piegl and Tiller. 25 First, the parameters i , i = 1, …, n ti are defined for which the NURBS S takes the values of the FE nodes x ti i . In this work, the chord length method is used to compute i and the knot vector is determined by averaging of i . 25 The control points are determined by solving a linear system of equations The interpolation is required in each parameter direction and is applicable to a wide range of geometries. It can be used in the manner of a black-box, but the number of elements in the NURBS layer is in general higher then using the knot insertion approach.
For the calculation of the normal gap between the FE mesh and the NURBS layer a normal projection is required. Newton's method is applied to determine the parameter of the projected point on the NURBS layer. Only for linear interpolations an analytical solution is feasible. A given point x s is projected onto the NURBS surface S( ) by minimizing the distance A sufficient condition for a minimum is d ′ ( ) = 0. This is equivalent to Here, S d denotes the derivative of S with respect to the coordinate d and Newton's method is applied to find the parameter vector . For a detailed description see the study of Piegl and Tiller. 25

CONTACT DISCRETIZATION
The normal gap between the slave FE domain and the NURBS layer is given by Here, x s i and R s i denote the FE nodal positions and the appropriate shape functions on the slave side. The control point positions and shape functions of the NURBS layer are x iga j and R iga j . The integration point s on the slave side is projected onto a point on the NURBS layer with the parameter and the normal vector n iga ( ). The variation of the normal gap Using the vector with the components the normal gap and its variation can be written in matrix form The vector contains the positions of the slave nodes and the master control points, its variation is denoted by x c . Only a limited number of FE and NURBS basis functions is not vanishing for specific parameters s and , which results in a sparse vector N.
The discretization of the virtual work of the normal contact traction Equation (6) is given by where the nodal normal gap g Ni and its variation are The tributary area A i of a node i is given by In this paper, the contact integrals are evaluated on Γ c (Equation (6)). The same notation is used for the discretized slave contact boundary in the reference configuration. The nodal contact pressure p Ni is defined by one of the following formulations: Gauss-point-to-surface (GPTS) formulation: In the GPTS formulation, see Fischer et al, 29,30 the contact is enforced at each predefined Gauss integration point located on the slave side. Only Gauss points with a negative gap are contributing to the contact force. The nodal contact pressure for the GPTS formulation is given by Mortar formulation: Using the mortar formulation the contact constraints are projected to the control points. 8,9 The nodal contact pressure is given by In the mortar formulation, the contact conditions are enforced by penalizing the nodal gap g Ni . The discrete contact pressure p h is given by In this paper, the contact constraints are enforced approximately by means of the penalty method. Using the mortar formulation, the nodal gap g Ni is penalized. In the GPTS formulation, contact constraints are enforced by the penalty function at each integration point.

Matrix form
The discretized virtual work of the normal contact traction (Equation (26)) is expressed in matrix form as: 21 The contact force f c is the product of the gap matrix G and the vector of the nodal contact pressures p N : The vector p N is assembled using the GPTS or the mortar formulations (Equations (29) or (30)). The gap matrix G is given by where The gap matrix G is split into submatrices D and M with

COUPLED FE-NURBS APPROACH FOR IMPACT SIMULATION
In the coupled FE-NURBS approach, 21 a NURBS layer is tied to one of the contacting bodies discretized by FEs. The NURBS layer is the master side for the mesh tying and contact conditions. Due to the smooth normal field of the NURBS layer, a robust and simple computation of the normal gap is ensured. Additionally, the NURBS layer provides smooth shape functions at element intersections, which leads to an accurate element-based contact discretization.

Mesh tying of the NURBS layer
In Otto et al, 21 two different methods are discussed for the tying of the NURBS layer: pointwise tying and mortar mesh tying. It is shown that the mortar tying is a robust method, giving accurate results for different NURBS layer parametrizations. The pointwise tying is computationally more efficient but only accurate for specific NURBS meshes with a matching parametrization or in cases when a node to node tying is achieved. Due to its robustness and general applicability, the mortar method is used for the tying of the NURBS layer to the FE mesh. The constraints resulting from the mortar mesh tying are given by where d ti are the FE node displacements at the tying zone and d iga are the NURBS control point displacements. The tying matrices D mt and M mt are given by Here, the FE shape functions at an integration point ti are condensed in the matrix N ti and the NURBS shape functions at the projected integration point are condensed in the matrix N iga with The number of FE nodes in the tying zone is given by n ti and the number of control points is given by n iga . The dimension dim quantifies whether it is a 2D or 3D problem. The Equation (37) may also be reformulated as Thus, the FE displacements d ti at the tying interface are calculated by means of the isogeometric displacements d iga .
A time-consuming part of the calculation of P is the inversion of D mt . Although the computation of P is performed only once in the simulation, the inversion for large matrices D mt is computationally not affordable and the storage due to the dense structure might exceed the computational resources. An alternative proposed in Wohlmuth 31 is to use bi-orthogonal dual basis functions on the slave side, which then renders the matrix D mt diagonal so that no inversion is needed to compute the operator P. The bi-orthogonality condition is typically applied for each slave element. Consequently, the computation of the whole bi-orthogonal basis becomes a problem with linear complexity. In this paper, the focus is placed on impact problems with a relatively small ratio between contact surface nodes and the total number of nodes, which allows to use standard FE shape functions. In this case, the matrix D mt is relatively small and its inversion can be carried out. In other situations, such as shell-like structures, where the modeling of the impacting wave is less relevant compared to the structural dynamics response (with much lower frequencies), implicit schemes might be more beneficial.

Coupled FE-NURBS formulation
In this work, the coupled FE-NURBS discretization, as introduced in Reference 21 for static problems, is applied to dynamic problems involving contact. Explicit time integration is used in this work. Thus, no solution of a system of nonlinear equations and its linearization is required. In this section, the dynamic equilibrium equation involving contact conditions and the mesh-tying constraints due to the NURBS layer are discussed. A linear elastic material model is assumed. The variational form of the dynamic equilibrium is given by The discrete kinetic and internal energies are denoted by W kin,h and W int,h . The discrete external energy due to Neumann boundary conditions or body forces is denoted by W ext,h . The contribution of the normal contact to the virtual work W c,h is given by Equation (26). The variation of the kinetic energy is expressed by where M is the mass matrix and d a block vector, containing the displacements of the NURBS layer d iga , the FE tying zone d ti and the remaining FE nodes d in : The acceleration vector is denoted byd. The variation of the internal and external energies is given by where f int,h and f ext,h are the internal and external force vectors. The contact contribution to the virtual work is given by Equation (32). Reformulating the matrix G and the nodal contact pressures p N with respect to the vector d, which virtually adds some zero entries and reorders the entries, results in the variational form of the dynamic equilibrium problem: where d is the kinematically admissible variation of the displacement vector d. According to Equation (40), the FE nodes at the tying interface are calculated by means of the NURBS layer displacements. This is represented by the constraint matrix C: where C = Applying Equation (46) to the variational form of the dynamic equilibrium Equation (45) results in where d ind is the kinematically admissible variation of the independent displacements d ind . The dynamic equilibrium equation with contact and mesh tying conditions follows from Equation (48): which is a second-order ODE. For the numerical solution the initial displacements and velocities are required as initial conditions. The mass matrix M is assembled according to the block structure of the displacement vector d in Equation (43). Additionally, the off-diagonal block entries of M are zero due to the spatial discretization with spectral FEs (Section 3.1). The NURBS layer does not have any mass, which results in the block entry M iga,iga = 0. Thus, M has the following structure: The matrix C T MC is expressed as The projection matrix P (Equation (40)) is generally nondiagonal, which results in a nondiagonal block matrix P T M ti,ti P. The size of the nondiagonal block matrix is (dim ⋅ n ti ) × (dim ⋅ n ti ). Using a NURBS layer on a narrow region on the boundary results in a relatively small nondiagonal submatrix. Obviously, there are mechanical problems, where the contact boundary is larger than the bulk part, for example, shell structures in contact. In these cases, the nondiagonal submatrix becomes larger than the diagonal submatrix of the mass matrix C T MC. In those cases, the advantage of using a bulk discretization with spectral elements is rather limited.

EXPLICIT TIME INTEGRATION
In this work, numerical methods are presented to model wave propagation effects after impact in linear elastic materials. Due to the impact event and its rough response, the wave spectrum contains high frequency waves. Implicit schemes are unconditionally stable and do not require a minimum time step for a stable time integration. Due to accuracy reasons, the time step of the implicit scheme has to be adapted to the highest frequency, to track the whole spectrum. Moreover, due to the contact conditions, a system of nonlinear equations has to be solved in the implicit method. The starting value for the numerical solution of this nonlinear system of equations is mostly the equilibrium of the last time step. Using a large time step in the implicit method means that the equilibrium state, which has to be computed in a recent time step, is probably far away from the last equilibrium state. This can result in convergence issues. The required small time step and the computational effort to solve a nonlinear system of equations in each time step makes implicit methods for nonlinear wave propagation problems very time consuming. Explicit methods are only conditionally stable and require a time step below a critical time step, which depends on the chosen explicit time integration scheme. The critical time step Δt c is determined by the highest frequency max in the wave spectrum: Δt c ∝ 1∕ max . This is required anyway to accurately track the whole wave spectrum. Using explicit schemes, no solution of a nonlinear system of equation is needed in each time step. Additionally, explicit time integration in combination with a diagonal mass matrix does not even require the solution of a system of linear equations in each time step, which makes the explicit time integration very efficient for wave propagation problems.
The dynamic semidiscretized equilibrium equation is a second-order ODE. A very popular explicit time integration scheme is the second-order CDM. As shown in Otto et al, 3 the popular combination of CDM with linear elements results in a compensation of spatial and temporal errors. Combinations of higher order methods in space and time are computationally more efficient.

Runge-Kutta time integration
As discussed in Otto et al, 3 where v is the velocity vector. The RK method applied to Equation (53) gives Here, Δt denotes the time step, a ij and b j denote constants defined by the well-known classical RK scheme of order k = 4. The coefficients a ij and b j are taken from table 1. The mass matrix M and the projection matrix P are constant throughout the simulation and the bulk part of the mass matrix C T MC is diagonal due to the spectral FE discretization (Equation (51)). To ensure the efficiency of the explicit time integration, the nondiagonal part of the mass matrix P T M ti,ti P is not inverted. Instead, it is prefactorized and a forward and backward substitution for the affected (dim ⋅ n ti ) degrees of freedom is carried out in each time step.

The critical time step
Explicit time integration methods require the time step to be smaller than a critical time step to guarantee the stable integration of the ODE. For linear wave propagation problems, the critical time step is not deformation dependent and is given for the classical RK method by where max is the maximum eigenvalue of the system For linear wave propagation, the stiffness and mass matrices K and M are constant. The dynamic equilibrium equation (49) is nonlinear. Thus, the maximum eigenvalue max is not constant throughout the simulation and is affected by the mesh density and the penalty formulation, which is used to enforce the contact conditions.

MESH-DEPENDENT PENALTY FORMULATION
In this work, the contact constraints are enforced by the penalty method. As shown in Otto et al, 3,21 a mesh-dependent penalty formulation is computationally advantageous. For dynamic problems, a nonlinear penalty regularization is proposed in Otto et al. 3 As shown in Figure 2, the penalty function is composed of a nonlinear part and a linear part. For the nonlinear part, a polynomial regularization p N (g N ) = a ⋅ g q N of degree q is used. This formulation has the advantage that the transition between no contact and contact is  q−1 smooth. The transition between the nonlinear and linear part is determined by the penetration g max . As soon as g max is reached, the penalty function becomes linear with the constant maximum stiffness k max .
The parameters of this penalty function are the penetration g max and the appropriate pressure p max = p N (g max ). The pressure p max is problem dependent and has to be estimated by the user. The penetration g max is chosen as where m is a user-defined constant, Δt c the critical time step without taking into account the contact conditions and v 0 the initial impact velocity. The factor a of the penalty function p N is determined from the equation p N (g max ) = p max and the maximum penalty stiffness is given by k max = p ′ N (g max ). The smaller m, the stiffer the formulation. A stiffer penalty formulation is also achieved by mesh refinement. The critical time step Δt c of the linear dynamic equilibrium problem is not deformation dependent, but choosing a finer mesh increases the maximum eigenvalue max and decreases the critical time step Δt c . Thus, also g max decreases for finer meshes and a stiffer regularization is achieved. In Reference 3, it is shown for 1D impact problems that for a constant factor m the critical time step of the penalty enforced impact problem remains constant upon mesh refinement. Thus a stiffer penalty formulation is achieved for finer meshes, but the critical time step of the impact problem remains constant when m is constant.
For 2D and 3D impact problems with penalty regularization, an estimation of the critical time step can be based on diagonalization of the penalty stiffness and the element eigenvalue inequalities, which is virtually equivalent to a 1D impact problem. 32 The bipenalty approach, where both inertia and stiffness terms are penalized to guarantee a constant critical time step upon the increase of the penalty stiffness, is valid for 1D impact problems. 33 Another possibility to guarantee a stable time integration is based on monitoring of energy balance or linear stability during the simulation and a possible restart with an adjusted time step.

Comparison of the GPTS and the mortar formulations
The static Hertz contact problem (Figure 3) is used to compare the GPTS and the mortar contact discretizations for the nonlinear penalty function p N = a ⋅ g 2 N . The NURBS layer is generated by interpolation and tied to the structured mesh region of the halfspace (Figure 3). In Temizer et al, 34 the enforcement of the contact conditions at each integration point resulted in a stiff formulation with convergence problems. In this section, the stiffness of the nonlinear penalty function is increased to compare to the results in the literature. The parameter p max = 592 of the penalty function is chosen slightly above the analytical maximum pressure of 591.43. The stiffness of the penalty function p N is controlled by the second parameter g max , which is chosen as −10 −5 or −10 −6 . Since p max is slightly above the maximum contact pressure, the parameter g max corresponds to the maximum penetration. An element-based contact discretization with eight Gauss-Legendre integration points per element is used. Figure 4 shows the contact stresses using the GPTS and the mortar formulations. The contact stresses are scaled by the analytical maximum contact stress p 0 and the contact width is scaled by the analytical maximum contact width a. Using g max = −10 −5 ( Figure 4A) the results of the GPTS and mortar methods are similar and convergence is achieved for both methods after seven iterations in Newton's method. Using g max = −10 −6 ( Figure 4B), the mortar method shows some more oscillations in the stresses, but converges in 10 iterations compared to 12 iterations needed for the GPTS method. For moderate penalty stiffness, the accuracy and convergence behavior of the mortar and GPTS methods are similar. Using higher penalty stiffness, the GPTS method needs more iterations but the stresses are free of artificial oscillations, as also reported by Kopačka 35 for a linear penalty formulation. Using an explicit time integration, no solution of a nonlinear system is required, thus the convergence behavior of the formulation is of minor importance. Using segmentation of the contact region, an accurate numerical integration can be achieved, where the oscillations in the contact stresses are reduced. 36,37 The segmentation is an accurate but computationally expensive technique, especially for 3D problems. In the FE-NURBS approach, the accuracy of the integration is increased by higher order quadrature or by an increased interpolation order of the NURBS layer. 21

2D cylinder-cylinder impact
In this section, a 2D plain strain impact problem of two collinear cylinders is considered and the GPTS discretization is used with eight Gauss-Legendre integration points per element. The FE mesh is shown in Figure 5A. The NURBS layer is tied to the lower cylinder in the structured mesh region and is constructed by interpolation of FE nodes. A cubic penalty regularization p N = a ⋅ g 3 N is used. The factor m = 15 is chosen for the penetration g max = m ⋅ (2v 0 )⋅Δt c . The penalty pressure at g max amounts to p max = 200. The time integration is done with the classical explicit RK method. The time step is chosen as Δt = 0.95Δt c , where Δt c is the critical time step of the dynamic equilibrium problem without the penalty term. According to Belytschko 12 a reduction factor of 0.95 is used. The wave propagation after impact is shown in Figure 5B. The critical time step Δt c,p for the nonlinear equation (49) is also computed during the simulation and the maximum deviation is 1−Δt c,p ∕Δt c = 2 ⋅ 10 −7 .
For stiffer regularizations m = 5 and 10, the ratio 1−Δt c,p ∕Δt c amounts to 2 ⋅ 10 −6 and 4 ⋅ 10 −7 , which shows that also for increased penalty stiffness the choice of the time step Δt = 0.95Δt c remains valid and ensures a stable time integration.

Comparison of linear and nonlinear penalty functions
In this section, the linear and nonlinear penalty regularizations are considered and compared to each other. For the 2D impact problem shown in Figure 5 no analytical solution exists. An approximate solution based on the conservation of momentum, where small deformations are assumed and no wave propagation is taken into account is given in Johnson 38 or Miyazaki et al 39 in the form of an ODE: Here, d and F c denote the indentation depth and the contact force. The contact width is given by a = √ Rd and the maximum contact stress for E = E 1 = E 2 and = 1 = 2 is given by  Figure 2, provides a smooth transition between no contact and contact state. In Figure 6 the mesh from Figure 5A is used along with a cubic penalty function p N = a ⋅ g 3 N . The penetration g max and the pressure p max are chosen as: g max = 15 ⋅ (2v 0 )⋅Δt c and p max = 200. For the linear penalty regularization, the penalty factor k max = p ′ N (g max ) is used. For both regularizations, the time step of Δt = 0.95Δt c is used. As shown in Figure 6A, the contact force using the linear penalty function is oscillating. Considering the contact stresses at one time step t = 0.4 ( Figure 6B), for both the linear and nonlinear formulations the stresses are similar.
The drop of the contact force at approximately t = 0.5 is induced by a stress wave causing an expansion of the material in the contact zone at that time.

Mesh-dependent formulation
Using a linear penalty regularization, the choice of the penalty factor is often heuristic. As shown in Otto et al. 21 for static contact problems, a mesh-dependent penalty formulation is more efficient. The error in the computations is not reduced by just increasing the penalty parameter for a constant mesh. After a certain threshold, the error even increases for higher penalty parameters. Using the nonlinear penalty regularization, the stiffness of the penalty function increases by mesh refinement. In Figure 7, the contact forces and the maximum contact stresses are illustrated for the problem in Figure 5. Three different meshes are used with the mesh structure corresponding to Figure Figure 7A). The contact stresses are given in Figure 7B. For finer meshes, the contact width is reduced and the contact stresses remain free of artificial oscillations.

Energy conservation
In Figure 8, internal E in , kinetic E kin and total E total energies are plotted for the problem in Figure 5. The total energy is constant throughout the simulation, which shows that energy is conserved in the given framework. Given the material parameters and the geometry dimensions in Figure 5, the total kinetic energy at the beginning of the simulation amounts to E kin = R 2 v 2 0 = 201.06, which coincides with the value of the total energy plotted in Figure 8.

Comparison of different orders of FE interpolation
As reported in De Lorenzis et al 8  when the degree of the FE interpolation and NURBS is increased. For all considered interpolation degrees, the contact stresses show similar results and the absolute values lie slightly above the quasi-analytical solution for the contact stresses.

3D cylinder-cylinder oblique impact
In this section, a 3D impact problem of crossed cylinders is considered. The meshes and parameters for the linear elastic material of the crossed cylinders are identical and shown in Figure 10. The 3D mesh is obtained by extrusion of the unstructured mesh and has 10 elements in the direction of the extrusion, as shown in the upper part of Figure 10. Quadratic FE and NURBS interpolations are used. The upper cylinder is located in the halfspace y ≥ 0 and the lower cylinder in the halfspace y ≤ 0. The NURBS layer is tied to the lower cylinder. The x − y cross section of the NURBS layer is an arc and is constructed by interpolation of FE nodes. The x coordinates of the interpolated FE nodes are lying in the range |x| ≤ 0.8R. In the z direction, the parametrization of the NURBS layer matches the FE mesh and extends in the coordinates |z| ≤ 0.9L. A nonlinear penalty regularization with p N = a ⋅ g 3 N is applied with g max = 5 ⋅ (2v 0 )⋅Δt c and p max = 250. In this example, the mortar and GPTS contact discretizations with the classical explicit RK time integration are used. The time step is chosen similar to the 2D example as Δt = 0.95Δt c . Due to the linear elastic material model, the critical time step Δt c is constant when the contact conditions and the appropriate penalty regularization are not taken into account.
The contact force F c and the maximum contact stress distribution at time t = 0.113 are shown in Figure 11A,B for the mortar and the GPTS methods. For the contact discretization, a Gauss-Legendre quadrature rule with six integration points in each direction is used. For the rather coarse mesh, the contact force F c and the contact stress are free of artificial oscillations. The contact stress has a symmetric distribution over the contact area, which corresponds to the symmetric setup. The GPTS and the mortar methods give virtually identical results. The GPTS formulation leads to slightly smoother stresses.
The deformations and the wave propagation until the time t = 0.113 as the maximum contact stress is reached are shown in Figure 12.

CONCLUSIONS
In this paper, two mechanical problems are considered: the impact problem and the subsequent wave propagation following the impact event.
For the contact discretization, a FE-NURBS approach is used. The auxiliary NURBS layer provides a robust and easy normal gap computation and an accurate element-based integration without the necessity to use full 3D NURBS geometries. An alternative to the element-based integration is segmentation, which is time consuming even for static 3D problems. A GPTS and a mortar contact discretizations are considered. Both methods are shown to provide accurate results for 2D and 3D impact problems. The numerical modeling of the impact event is intrinsically ill-posed due to the instantaneous change of velocities in the contact area, which leads to artificial oscillations in the numerical solution. In this work, a nonlinear penalty regularization is applied with a smooth transition from the no contact to the contact state. The nonlinear penalty function is defined by two mechanically motivated parameters, which make the nonlinear regularization mesh dependent. The stiffness increases upon mesh refinement. Compared to the linear penalty regularization, less oscillations appear in the contact forces. For finer meshes, the penalty stiffness is automatically increased without oscillations in the contact force or contact stresses.
The wave propagation problem requires a small time step for the time integration scheme, which makes implicit schemes not affordable. Explicit time integration schemes are efficient, provided the mass matrix is diagonal. Moreover, explicit time integration schemes are only stable for time steps below of a critical time step. For nonlinear problems, like the impact problem, the critical time step is not known a priori. In this work, a spectral FE discretization is used, thus the bulk part of the mass matrix is diagonal. Due to the coupled FE-NURBS discretization the part of the mass matrix corresponding to the contact region is not diagonal. This independent part of the mass matrix is prefactorized and only a forward and backward substitution is needed in each time step of the explicit time integration. As shown by the 2D and 3D examples, using the nonlinear mesh-dependent penalty regularization the critical time step is mostly influenced by the mesh density and not the penalty stiffness. Thus, it is sufficient for the linear elastic material model to compute the critical time step Δt c without taking into account the contact conditions. The simulations are run with the time step slightly below of Δt c .
So far this method is applied to linear elastic materials. The application to nonlinear elasticity with finite strain and friction is an interesting and from our point of view promising extension of the method.