Flexible multibody impact simulations based on the isogeometric analysis approach

Usually detailed impact simulations are based on isoparametric ﬁnite element models. For the inclusion in multibody dynamics simulation, e.g., in the framework of the ﬂoating frame of reference, a previous model reduction is necessary. A precise representation of the geometry is essential for modeling the dynamics of the impact. However, isoparametric ﬁnite elements involve the discretization of the geometry. This work tests isogeometric analysis (IGA) models as an alternative approach in the context of impact simulations in ﬂexible multibody systems. Therefore, the adaption of the ﬂexible multibody system procedure to include IGA models is detailed. The use of nonuniform rational basis splines (NURBS) allows the exact representation of the geometry. The degrees of freedom of the ﬂexible body are afterwards reduced to save computation time in the multibody simulation. To capture precise deformations and stresses in the area of contact as well as elastodynamic effects, a large number of global shape functions is required. As test examples, the impact of an elastic sphere on a rigid surface and the impact of a long elastic rod are simulated and compared to reference solutions.


Introduction
In machine dynamics impact problems often contain large nonlinear motions and small elastic deformations, if the flexible body is relatively stiff. This is the case for steel or aluminum bodies which allows the assumption of small linear deformations. Therefore, the floating frame of reference formulation can be applied to model impacts in flexible multibody systems. T The application of the floating frame of reference formulation requires global shape functions of the flexible bodies. One way to obtain the global shape functions is to model the flexible body with isoparametric finite elements. A major drawback of isoparametric elements is the discretization of the geometry. However, an exact representation of the geometry is essential for determining the occurring contact forces and stresses in an impact problem. The discretization error can be worse in a sliding contact. The contact bodies can get entangled due to an insufficient discretization as discussed in [2]. These artificial effects from isoparametric elements might yield a model behavior, which does not necessarily represent the reality.
This motivates the use of isogeometric elements [2]. Isogeometric elements have two advantages. Firstly, the local shape functions are nonlinear basis splines (B-splines) and allow the exact representation of the geometry without any discretization errors. Secondly, it is shown in [2] that high frequency modes of flexible bodies are represented more accurate with isogeometric elements compared to isoparametric elements. However, a downside of isogeometric elements compared to isoparametric elements is the increased computational cost. The evaluation of the nonlinear shape functions of isogeometric elements is more time consuming than the evaluation of the linear shape functions of isoparametric elements.
The literature proposes a variety of different contact algorithms for isogeometric elements [10,17,18]. There are two main types of contact methods. The first type are the node-to-surface algorithms. They include the node-to-surface algorithm of isoparametric elements. Several node-to-surface algorithms exist for isogeometric elements. One way is to check the exterior unique knots for contact [18]. Another possibility is to define a collocation method [10,18] which checks the collocation points for contact. The second type are integral-based descriptions of a contact, also referred to as knot-to-surface methods [10,17]. Using the Gauss integration, the contacts of individual Gauss points can be checked. This usually results in a higher number of evaluated points compared to the node-to-surface algorithms. Therefore, a knot-to-surface method is computationally more expensive but grants a higher accuracy than the node-to-surface algorithms. Since the focus in this work is on high accuracy of the geometry by using the IGA, a knot-to-surface method is used.
The use of the floating frame of reference formulation involves the description of deformations by global shape functions. Thereby, mostly the low frequencies are considered and thus the overall flexibility of the body is included. The straightforward approach involves the modal reduction. This method captures the overall flexibility and wave propagation inside the body in case of an impact. Whereas local deformations and stresses in the contact area cannot be recovered using a moderate number of global shape functions from modal reduction [19]. Another drawback of the modal reduction is that the solution of the impact does not converge to the exact solution due to missing local deformation. This behavior is observed in [16,20] when using the penalty method. Therefore, the Craig-Bampton method is used instead [3]. In addition to the normal modes, static shape functions are also included. Thus the deformation and stress in the contact area can be recovered very accurately. The cost of the Craig-Bampton method is the increased numerical stiffness due to the high frequency of the static shape functions and increased number of modes.
The aim of this work is to test isogeometric elements as an alternative approach to isoparametric elements in the context of the floating frame of reference approach. Therefore, it is shown how global shape functions can be determined from isogeometric finite element models, embedded in the floating frame of reference formulation and used, e.g., for impact simulations.
This work is organized in the following way. In Sect. 2 the concept of the floating frame of reference formulation is briefly summarized. The concept of the IGA is introduced in

Floating frame of reference formulation
A well-known approach for simulating flexible multibody systems is the floating frame of reference formulation [12,14]. Each flexible body of the system has a body-related moving reference frame K R . Within the inertial frame, the large nonlinear rigid body motion R r IR is defined, where the prescript 'R' denotes the description in the body reference frame. In this work, a Buckens frame [12] is utilized as floating frame. The assumed small elastic deformations of the bodies are linear and are described in the corresponding body reference frame. The position of an arbitrary point P on the flexible body with respect to the body reference frame is described by R c RP . The small linear elastic deformation can be reconstructed with the global shape functions and the n q elastic coordinates q e . The position R r IP of the point P within the inertial frame is expressed by in the reference coordinate system. The description of the point P in Eq. (2) can be seen in Fig. 1. The position, velocity, and acceleration of any point P can be expressed in terms of the position coordinates z I , velocity coordinates z II and accelerationsż II , which are defined as They consist of the rigid body deformation of the body reference frame with the translation R r IR and its velocity R v IR from the inertial frame to the body reference frame. Rotations are represented by rotation parameters β IR and the angular velocity R ω IR . The elastic motion is given by q e andq e . The equations of motion for a single flexible body are finally given by ⎡ ⎣ mE mc which are assembled by the standard input data (SID). The SID [12] is a well-known standard in providing elastic data in the context of the floating frame of reference formulation. In Eq. (4), the mass of the body is denoted by m, the center of mass relative to K R byc, the translational coordinate coupling matrix by C t , the rotational coordinate coupling matrix by C r , the mass matrix M e of the flexible body and the mass moment of inertia by I . The right-hand side vector h a is composed of the vector of surface forces h p , the discrete forces h d , the body forces h b , the inertial forces h ω and the internal forces Here K e represents the stiffness matrix and D e the damping matrix. Occurring contact forces are applied at the discrete forces h d .

Global shape functions from reduced isogeometric analysis models
A major question using the floating frame of reference formulation is how to obtain the global shape functions . This section presents the key idea of the IGA and the procedure to obtain the global shape functions from an IGA model. A detailed introduction to the IGA can be found, for instance, in [2]. A general way to determine the global shape functions is to generate a finite element model of the flexible body and then to identify the global shape functions form the finite element model using model reduction techniques. The generation of a two-dimensional isogeometric finite element model is briefly explained in the following.

Basis-splines
Two spaces are used in the IGA, the physical space and the parameter space. The shape functions of isogeometric elements are defined in the parameter space which can be seen on the left-hand side in Fig. 2. The parameter space is divided into elements and spanned by the knot vectors = ξ 1 ξ 2 . . . ξ n+p+1 and H = η 1 η 2 . . . η m+q+1 .
Thereby, p and q are the polynomial orders and n and m are the numbers of the basis functions N i,p and M j,q . The basis functions are based on B-splines and can be computed recursively by the Cox-de Boor algorithm. The shape functions in the ξ -direction are computed The calculation rule given by Eq. (7) and Eq. (8) is identical in η-direction. In practice, recursive functions are numerically inefficient. Therefore, a nonrecursive algorithm is used, which is suggested in [9].

NURBS
Besides the parameter space there is the physical space, which can be seen on the righthand side in Fig. 2. In the physical space the control net is defined. It is formed by control points P i,j which have a corresponding weight w i,j . To transform the parameter space into the physical space, the NURBS basis R p,q i,j is given by using above B-splines N i,p and M j,q . The NURBS basis and the control points then lead to the NURBS surface which gives the B-spline surface in the physical space. The degrees of freedom in the IGA correspond to the displacements u i,j of the control points whereby P 0 i,j represents the position of the undeformed control points and P i,j the position of the deformed ones. The deformation d of the NURBS surface can be computed with where the basis functions of the corresponding element are summarized in N .

Equations of motion for isogeometric elements
The procedure of setting up the equations of motion for isogeometric elements is almost identical to isoparametric elements. This is because the weak Galerkin method is also applied for isogeometric elements. Like for isoparametric models, the isogeometric model requires a preprocessing [2]. This includes the determination of the local mass and stiffness matrices of each isogeometric element. For relatively stiff materials linear elasticity can be assumed. The calculation rule is identical to isoparametric elements [1]. The local mass and stiffness matrices are given by where C is the material elasticity matrix and n e is the number of elements. The strain displacement matrix B is assembled by differentiating the NURBS basis R p,q i,j and using the Jacobian transformation [1]. The integration over each element e,j is achieved by the Gauss quadrature in the parameter space. The basis functions N and the strain displacement matrix B are evaluated at each Gauss point. As stated in [2], the same Gauss rule for a polynomial of the order p can be applied to a p-th order NURBS function. Therefore, the Gauss order p + 1 is chosen in ξ -direction and order q + 1 in η-direction. The global system matrices K e and M e are created by assembling the corresponding element matrices. The elements are connected by shared control points. The equations of motion of the isogeometric model in a body-attached nonrotated frame are then given by where the displacements of the control points are represented by u ∈ R n dof . The number of degrees of freedom is denoted by n dof and external forces by f e . It should be noted that Eq. (15) is only valid for a nonrotating body, otherwise a more complex inertial term would occur. Compared to isoparametric models, the system matrices K e and M e of isogeometric models are more dense since the number of shared control points is defined by (p + 1) × (q + 1). Therefore, more effort is required in the numerical solution. But since the system matrices are reduced before incorporating them into the floating frame of reference formulation, this disadvantage is not relevant in the simulation time.

Model order reduction of isogeometric elements
The most common approach for reducing the equations of motion (15) is the modal reduction. As shown in [16,20] this leads to inaccurate results in case of impact problems, since the local deformation in the contact region is not included in the reduced basis. An alternative approach is the Craig-Bampton method [3], which is a special case of the Component-Mode-Synthesis. The idea of the Craig-Bampton method is the combination of fixed-interface normal modes and constraint modes. The normal modes are determined by an eigenvalue problem and the constraint modes by unit displacements to each predefined interface control point. These predefined control points are located on the outside of the body in the contact area where later potential forces act. A detailed description of the Craig-Bampton method for isogeometric elements can be found in [8]. The procedure results in the shape function matrix ∈ R n dof ×nq . The resulting system consists of n q n dof degrees of freedom. The reduced mass and stiffness matrices are then given by where E is the identity matrix and ω i the natural frequencies.

Software implementation
The software solution used here to simulate flexible multibody systems is the MATLAB toolbox DYNMANTO [5], which is developed at the Institute of Mechanics and Ocean Engineering at the TUHH. The workflow to incorporate flexible bodies in multibody simulations is briefly summarized in the following. The geometry is modeled in a computer-aided design tool and imported in a finite element program such as ANSYS. The finite element data is then imported in the MATLAB toolbox RED [5]. RED reads the elastic data, applies a model order reduction technique and determines the SID. The flexible body is then defined in DYNMANTO using the SID to assemble the equations of motion (4).
This procedure changes only slightly for isogeometric elements. A geometry, which is based on B-splines, is designed with the in house MATLAB toolbox RIGA. Currently, only simple geometries can be modeled with RIGA. Importing any geometry from a CAD program is not possible at the moment. One problem is that standard CAD programs only map the outer surface of three-dimensional bodies in a two-dimensional parameter space. This is not sufficient for a full finite element model.
To enable compatibility with the workflow described above, the data structure of the ANSYS interface is adapted. Besides the global mass and stiffness matrix, the element table, the table of degrees of freedom and the coordinates of the control points are needed to enable the basic functionality. The use of the Craig-Bampton method additionally requires the specification of the identification number of the control points and the corresponding degrees of freedom of the interface control points.

Modeling of contacts using isogeometric elements
This section presents contact specific adjustments to the equations of motion (15) and aspects to increase the numerical efficiency. This includes the introduction of a penalty formulation for isogeometric elements.

Penalty formulation
The penalty formulation is chosen to model the contact in this work. The concept of the penalty formulation is the penalization of a penetration with a spring force. The spring constant c p is a tuning factor and needs to be chosen heuristically. Thereby, the penalty factor should be chosen large enough such that the results become independent of the chosen parameter [13]. Recapturing first the application of the penalty method for isoparametric elements, a set of exterior slave nodes is checked for contact. If a contact occurs, the node is loaded with the contact force f c = c p g n n, (18) where g n is the normal distance of the slave node to the master body, and n the normal vector.
In each time step of the simulation the body is checked for occurring contacts. This requires the deformed NURBS surface. Initially, the displacements u of the control points are recovered with using the global shape functions and the elastic coordinates q e . With Eq. (11) and Eq. (19), the position of the deformed control points R P i,j can be determined in the corresponding reference coordinate system of the flexible body. Finally, rotations of the body reference frame are considered using Tait-Bryan angles. The rotation parameters β IR = α β γ from Eq. (20) transforms a vector from the body reference frame to the inertial frame. To enable the check for occurring contacts, the deformed control points R P i,j are transformed to the inertial frame with The prescript 'I' denotes the description in the inertial frame.
In this work the contact is restricted to the contact of a body described by IGA with a rigid surface representing the master, see Fig. 3. A knot-to-surface based method using the Gauss integration [10,17,18] is chosen, where a set of n c slave elements is checked for contact. If the normal gap g n between the flexible body and the rigid surface is negative, the contact force I f c,j = e,j N c p g n n d e,j for j = 1...n c is applied with the penalty formulation to the slave body. After all slave elements have been checked for contact and the contact forces I f c,j have been determined, the resulting forces Flexible multibody impact simulations. . .
The forces f i,j,x and f i,j,y act on each individual control point in x-and y-direction, respectively. The forces I f i,j are described in the inertial frame, since the contact check is performed within the inertial frame. In order to determine the contact forces with respect to the flexible body, the forces I f i,j are transformed to the body reference frame with As mentioned before, the contact forces are assigned at the discrete forces h d in Eq. (4). The discrete forces h d consist of three components, including the components of translation h dt and rotation h dr of the rigid body, as well as the discrete forces h de affecting the elastic deformation. The discrete forces h d are then summarized in The contact resolution of isoparametric elements depends on the number of nodes in the contact area. In contrast to isogeometric elements the resolution depends primarily on the number of knots in the contact region. When the Gauss integration is used, the contact resolution additionally depends on the order of the B-splines. The resolution can be further increased by choosing a higher Gauss order. However, here the Gauss order p + 1 and q + 1 as defined in Sect. 3.3 is used. Thereby, p + 1 or q + 1 Gauss points are checked in each of the previously selected contact elements.

Increasing the numerical efficiency
One major drawback of the Craig-Bampton method is the numerical stiffness of the resulting equations of motion [11]. The frequencies of the normal modes are relatively low and the frequencies of the constraint modes relatively high. The high frequency Craig-Bampton modes are important to capture the local deformation of the contact region, however have no influence on the wave propagation. The wave propagation is determined by the normal Fig. 4 Impact of two identical spheres modes. One simple approach to increase the numerical performance is presented in [15] to efficiently damp the high frequency modes. This procedure has also been applied to reduced isoparametric element models [19] and is also used here for the reduced IGA model. By defining the Rayleigh damping coefficients the damping matrix can be determined. The frequencies f 1 and f 2 represent the frequency band of the normal modes. Thereby, frequencies between f 1 and f 2 are only minimally damped, which is necessary for capturing the wave propagation, where frequencies above f 2 are critically damped. The damping parameter κ is a tuning factor. Heuristic testing showed that the value κ = 0.001 does improve the numerical performance and has only a small influence on the elastic dynamic effects. The resulting damping matrix D e is then updated in the inner forces h e in Eq. (5).

Application Example I: impact of a sphere
As a first application example, the impact of two identical spheres is considered. The constellation is shown in Fig. 4. Both spheres have a radius of r = 10 mm. The chosen material is aluminum, therefore the Young's modulus is chosen as E = 7.28 × 10 10 Pa, the density as ρ = 2789 kg m 3 and the Poisson's ratio as ν = 0.33. Both spheres have an initial velocity of v 0 = 0.1 m s and gravity is ignored. In this section the proposed reduced isogeometric element model is compared to the analytical solution by Hertz [6]. Additionally, a comparison to a full nonlinear isoparametric model simulated in ANSYS is done.

Hertz impact
As a reference, the analytical solution for the impact of two spheres by Hertz is used [6]. For the impact of two spheres, they are modeled as mass points. During contact the contact force acts between the spheres, where δ denotes the displacement between the center of mass of the spheres. This corresponds to the overlap of the undeformed spheres, see Fig. 4. The equivalent Young's modulus E * and the equivalent radius r * are determined by The Young's modulus of the two spheres is denoted by E 1 and E 2 , the Poisson's ratio by ν 1 and ν 2 , the radius by r 1 and r 2 . Due to the symmetry of two identical spheres, the lower sphere, as shown is Fig. 4, is modeled by a rigid surface. An equivalent assumption can be found in [21] where the Hertz contact between two identical cylinders is modeled. The contact width a and the maximum pressure p 0 are given by Thus the Mise stress along the symmetry axis can be calculated as

Isogeometric model
To reduce the computational cost of the IGA impact simulation, the sphere is substituted by an axisymmetric semicircle. Thus a rotational symmetry model is used. This means that fewer interface control points are needed, the number of elastic coordinates n q is reduced and the number of points at which the contact is evaluated is reduced. However, the flexible body cannot perform a rigid body rotation, but only a rigid body translation along the symmetry axis. Due to symmetry, the lower sphere, as shown is Fig. 4, is modeled by a rigid surface. The order of the B-splines is given by p = 2 and q = 1. The knot vectors for a semicircle are  However, the definition of the order, knots and control points is only sufficient for the exact visualization of the geometry. To obtain an exact solution of the impact problem, the model must be refined, especially in the contact area [22]. Two common methods for refinement are the order elevation [7] and the insertion of knots [2]. The final model is depicted in Fig. 5a. The order is elevated to p = 4 and q = 3. The minimum edge length of an element in the contact zone is 50 µm. This results in 1 evaluation point per 10 µm, since the Gaussian order p + 1 is selected in the ξ -direction. The isogeometric model is reduced with the Craig-Bampton method. The number of normal modes is chosen as 6, and 10 interface control points on the lower side are used to determine the constraint modes. Therefore, the resulting number of elastic coordinates is n q = 26.
For comparison, the full IGA model is used. In order to include it in the same framework, the full IGA model is modally transformed using all modes. With n dof = 614 degrees of freedom this leads to n q = 614 elastic coordinates.

Isoparametric model
Beside the analytical solution and the proposed reduced isogeometric model, an isoparametric finite element model is used for comparison. Four-noded plane elements are utilized. The axisymmetric model is generated in ANSYS and visualized in Fig. 5b. The minimum edge length in the contact zone is 10 µm. As with the isogeometric model in Sect. 5.2, this results in 1 evaluation point per 10 µm. This isoparametric finite element model has n dof = 3588 degrees of freedom and due to the contact elements a full nonlinear analysis has to be performed.

Simulation results and analysis
Simulations of 0.1 ms are performed for analysis and comparison. The isoparametric model is solved with a fixed step size of 0.1 µs in ANSYS. Thereby, both the penalty formula-  Heuristic testing showed, that the numerical solution converged for the penalty factor c p = 5 × 10 18 N/m. This means that higher penalty factors yield the same results. An overview of both models can be found in Table 1. The simulations differ mainly in the short calculation time of the isogeometric model reduced with CMS compared to the isoparametric model. The isogeometric model which is modally transformed takes much longer than the isogeometric CMS model. The reason for this is the lower number of elastic coordinates in the CMS model. The differences in the results are discussed in the following. The solution by Hertz in Sect. 5.1 is used as a reference. In Fig. 6 the occurring contact force is displayed. The maximum absolute error of the isogeometric model is 1 N. The maximum absolute error of the isoparametric model is slightly lower, it is 0.7 N. There is no discernible difference between the CMS model and the full IGA model. The low frequency normal modes already map the overall elastic behavior, and the selected constraint modes are sufficient to describe the deformation in the contact area in this application example. However, the computation time of the CMS model is significantly lower compared to the full IGA model.
The relationship between contact force f c and displacement δ is shown in Fig. 7. The nonlinear relationship is well represented by all finite element models. Again, no difference can be seen between the IGA CMS model and the full IGA model.
The von Mises stress along the symmetry axis is given in Eq. (34) and visualized in Fig. 8. For the comparison, the time step was chosen at which the occurring stresses are maximum in the respective simulation. The maximum stress occurring in the flexible body is slightly better represented by the isoparametric model than by the isogeometric model. However, the  overall stress is better represented by the isogeometric model. Overall, the differences are very small and are of minor practical importance. The distribution of the von Mises stress in the contact area is shown in Fig. 9. It shows that the reduced IGA model reproduces the stresses predicted by the nonlinear isoparametric model very accurately.

Application Example II: impact of a long rod
Whilst the first application example exhibits quasistatic behavior, the second application example shows significant elastodynamic effects. For this second test case a long cylindrical aluminum rod with length of = 1 m and radius of r = 10 mm is used. At one end, the rod has a spherical end which impacts on a rigid surface. The wave speed in the rod can be computed, following, e.g., [4], with In the following it is investigated whether a reduced isogeometric model can represent the wave propagation, and this model is again composed to a full nonlinear isoparametric model in ANSYS.
T. Rückwald et al.  Testing showed again that the numerical solution of the reduced IGA model converged for the penalty factor of c p = 1 × 10 18 N/m and greater. An overview of the models is given in Table 2. As already observed in Sect. 5.4, the reduced isogeometric model again requires less computing time.
To monitor the wave propagation, the velocity in y-direction of a control point or a node is visualized in Fig. 11. For the sake of simplicity, no distinction will be made in the following between the terms "control point" and "node". Instead, the designation "point" is chosen. Both points are located on the symmetry axis at the top of the one meter long rod. The corresponding contact force of the impact is depicted in Fig. 12. The x-axis in Fig. 11 and Fig. 12 represents the integer multiples of the time it takes for the wave to travel one meter. The impact starts at the time stamp t 0 = 0 s. In the beginning, both points at the top end of the rod keep the initial velocity of the rod. The impact of the rod induces a wave, which reaches the monitored points at the time t = t 1 . As the wave has reached the end of the rod, the wave is reflected and moves back to the contact area. Between the times t = t 1 and t = t 2 the contact force remains almost constant. This shows that the rod, in contrast to the sphere, is not able to push itself off the rigid surface by decompressing the local deformation. Only the re-entry of the wave in the contact area at the time t = t 2 allows an end of  The temporal course of the velocity in Fig. 11 shows that both the isogeometric and isoparametric finite element models represent the wave propagation accurately. The wave speed from Eq. (37) is well mapped, which shows that the IGA is also capable of predicting elastodynamic behavior in a long rod well.

Application Example III: impact of a rotating body
In a third application example an impact problem with a large rotation is considered. As described in [14], the floating frame of reference formulation can model large rotations of flexible bodies. This application example is intended to show that this is also feasible with an IGA model. In the course of this, a rotating cylinder impacts on a rigid surface as depicted in Fig. 13.
The initial rotational velocity is ω 0 = 2π /s and the initial translational velocity in vertical direction is v 0 = 0.1 m/s. Due to the rotation, sliding friction is induced with the co- and the order of the B-splines is given by p = q = 2. The corresponding control points are and the weights are The order is elevated to p = q = 4 and the model is refined by inserting additional knots. Then, a model reduction is performed using the Craig-Bampton method. The number of normal modes is chosen as 10, and 26 interface control points on the lower side are used to determine the constraint modes. An overview of the final model is given in Table 3. The penalty factor c p = 1 × 10 17 N/m is determined in a series of simulation as described in Sect. 4.1.
As a reference the translational and rotational velocities after the impact can be determined with the principle of conservation of linear and angular momentum. An elastic impact is assumed and the moment of inertia of the cylinder about the z-axis is I = 1 2 mr 2 . With the defined initial velocities the velocities after the impact are determined by The calculated velocities can be used for validating an impact simulation. Figure 14 shows good agreement with the analytical solution in Eq. (44). The rotational velocity around the zaxis is visualized in Fig. 15. The excitation of the impact in tangential direction causes small vibrations in the cylinder. Nevertheless, the rotational velocity oscillates around the previously calculated value ω 1 in Eq. (45). Finally the position of the cylinder center of mass can be seen in Fig. 16. Due to the friction component of the contact force, the cylinder moves to the left after the impact.

Summary
Overall, it can be concluded that it is feasible to obtain the global shape functions for flexible multibody systems from isogeometric finite element models. By choosing the appropriate data interface, the global shape functions of the IGA can be smoothly included into a preexisting flexible multibody simulation toolbox. Due to the use of a model reduction of the IGA model a significant reduction of computation time is achieved compared to the nonreduced fully nonlinear isoparametric finite element model. The application of isogeometric elements for an impact simulation leads to accurate results in a quasistatic setup as presented for the sphere impact. The results of the isogeometric models agree very well with the isoparametric models and the analytical solution by Hertz. This can also be observed in the occurring contact forces and stresses in the flexible body. One reason for minimal differences in the results can be the reduction by the Craig-Bampton method. Additionally, the second application example of the impact of a long rod shows a good representation of an elastodynamic behavior when compared to an isoparametric reference model. A key feature of the floating frame of reference formulation, the handling of large rotations, is demonstrated in the third application example.
Funding Note Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.