A multinode superelement in the generalized strain formulation

Design and optimization of flexure mechanisms and real time high bandwidth control of flexure based mechanisms require efficient but accurate models. The flexures can be modeled using sophisticated beam elements that are implemented in the generalized strain formulation. However, complex shaped frame parts of the flexure mechanisms could not be modeled in this formulation. The generalized strain formulation for flexible multibody analysis defines the configuration of elements using a combination of absolute nodal coordinates and deformation modes. This paper defines a multinode superelement in this formulation, i.e., an element having its properties derived from a reduced linear finite element model. This is accomplished by defining a local element frame with the coordinates depending on the absolute nodal coordinates. The linear elastic deformation is defined with respect to this frame, where rotational displacements are defined using the off-diagonal terms of local rotation matrices. The element frame can be defined in multiple ways; the most accurate results are obtained if the resulting elastic rotations are as small as possible. The inertia is defined in two different ways: the so-called “full approach” gives more accurate results than the so-called “corotational approach” but requires a special term that is not available from standard finite element models. Simulations show that (flexure based) mechanisms can be modeled accurately using smart combinations of superelements and beam elements.


Introduction
Design and optimization of flexure mechanisms and real-time high-bandwidth control of flexure-based mechanisms require efficient, but accurate models. The flexures can be modeled using beam elements [1][2][3]. Sophisticated beam elements [4,5] for the modeling of flexures have been derived and implemented in the generalized strain formulation [6]. However, the frame parts of flexure mechanisms can have complex shapes (see, for example, the spherical joint [2] in Fig. 11) and can therefore not be modeled with beam elements.
Efficient modeling of arbitrarily shaped parts requires reduced order models. Reduced order models are finite element models, the (linear) deformation of which is reduced to a few generalized coordinates [7,8] (also referred to as component mode synthesis [9]). By using these models, arbitrarily shaped bodies can be defined by a single element with few degrees of freedom. In several formulations, such an arbitrary shaped element is called a superelement [10][11][12][13].
A two-node superelement is derived in the generalized strain formulation [12,13]. However, many frame parts are connected to more than two other components such that these frame parts cannot be modeled using the two-node superelement. This paper introduces a superelement with an arbitrarily number of interfaces, which will be referred to as Generalized-strain Multinode Superelement (GMS). Figure 1 gives an overview of the geometrically nonlinear multibody formulations and nonlinear finite element formulations in order to show the relation of the GMS with respect to existing formulations, as detailed below. All formulations define the configuration with respect to a global reference frame (also known as inertial frame). The formulations are categorized based on the type of coordinates that are used as degrees of freedom. In this paper the term "degrees of freedom" refers to the unknown coordinates that appear in the equation of motion. Other overviews of the different formulations can be found in [14][15][16][17]. In contrast to these papers, the current overview does not distinguish between multibody analysis (i.e., modeling physical components as a whole) and finite element analysis (i.e., partitioning physical components in multiple standard elements) as most formulations can be used for both analyses. The terminology in the literature about these two analyses is slightly different from each other, in this paper the following terms are used: "element" is used for modeling parts, "node" for the connections between the elements, and "element frame" is the local frame that defines the position and orientation of an element. "Absolute coordinates" are coordinates with respect to the global frame, while "local coordinates" are coordinates defined with respect to the element frame.
Two categories of formulation can be distinguished based on the degrees of freedom that are used to define the large motion of an element (the two columns in Fig. 1). One category uses the absolute nodal coordinates of the elements, i.e., the position (and orientation) of the nodes with respect to the global frame. The formulations in the other category use an element frame for each element to define its large motion. In these formulations the element frame is typically referred to as the floating frame and the formulations are referred to as the floating frame formulation, see, e.g., [18]. An advantage of using absolute nodal coordinates is that constraints can be easily applied: elements can be connected to each other by sharing nodes, and the displacement of some nodes can be prescribed. In other words, applying constraints eliminates degrees of freedom. In the floating frame formulations, the constraint equations are generally nonlinear relations between the coordinates of the element frames and the deformation coordinates of the elements. These equations are generally solved using the Lagrange multiplier method, increasing the total number of unknowns in the equation of motion. An advantage of the floating frame formulation is that small elastic deformation of Fig. 1 Overview of formulations for finite element and multibody analysis, categorized based on the coordinates that are used as degrees of freedom the element can be described linearly, relative to the element frame. This facilitates the use of arbitrarily shaped reduced-order models in the formulation, see, e.g., [19][20][21]. In [22], a slightly modified floating frame formulation was proposed, based on the model order reduction described in [23], which results in a constant mass matrix at the expense of additional deformation modes. The use of reduced order models is also possible for all formulations that use an element frame, as indicated in Fig. 1. Because of this potential to model arbitrary shaped bodies, the floating frame formulation is most often used in multibody simulations if the displacements due to elastic deformation of the physical components are small. The absolute nodal coordinate formulations are the preferred method for (nonlinear) finite element simulations, the finite element models generally contain many (standard) elements, making it important to keep the number of degrees of freedom as low as possible.
The second division of the categories in Fig. 1 defines whether an element based on absolute nodal coordinates uses an element frame (often called corotational frame in these formulations) to describe the rigid rotation and to define the elastic deformation relative to this frame. In contrast to the floating frame formulation, the coordinates of these frames are not necessarily degrees of freedom (i.e., they do not appear as unknowns in the equation of motion) but their coordinates are implicitly defined as functions of the degrees of freedom.
The method that applies this approach is the corotational formulation, see, e.g., [24][25][26]. Because the elastic deformation can be described linearly to this frame, reduced-order models can be used in the corotational formulation, see, e.g., [11,27,28]. The inertial frame formulation (see, e.g., [29]) does not use element frames to distinguishes the rigid motion from the flexible motion. Therefore, the nonlinear Green-Lagrange strain definition is used, which is valid under large rigid motions. An "absolute nodal coordinate formulation" that is developed by Shabana [30] is a formulation in terms of absolute nodal coordinates that defines the orientations of the nodes using slopes. In this formulation corotational elements [30,31], as well as inertial frame elements [32], can be developed.
The two rows in Fig. 1 define whether the degrees of freedom contain generalized coordinates related to deformation modes in order to define the stiffness. In case of the floating frame formulation, the inclusion of these deformation modes is the only way to include flexibility of the elements. In the absolute nodal coordinate formulations, the deformation of elements is already implicitly defined by the displacements of the nodes.
The generalized strain formulation [6,33,34] (also referred to as the natural modes' approach [35][36][37]) defines the deformation of an element using deformation modes. The generalized coordinates associated to these modes are expressed as analytical functions of the absolute nodal coordinates. These generalized coordinates are called "generalized deformations" (in other literature sources also referred to as "generalized strains" although they are related to displacements instead of strain). The generalized deformations remain constant under rigid body motion. Using proper definitions, the deformation modes can be given a physical meaning like the elongation of a beam element. The constitutive law is expressed in terms of the deformation modes. Because the generalized deformations are independent of rigid motions, the resulting constitutive equations are linear or relatively simple nonlinear equations, in contrast to the inertial frame formulation. Rigid elements can be modeled by applying constraints on all the deformation modes. Also part of the deformation modes can be constraint to keep only the most important flexibility. This is an advantage compared to the inertial frame and corotational frame formulations which only allow the modeling of flexible elements. The inertia forces of the element are defined using the absolute nodal coordinates.
A challenge in the generalized strain formulation is the definition of suitable deformation modes. For many default elements deformation modes are defined, like trusses, beams [4,38], hinges [38], and wheels [39]. Also a two-node superelement [12,13] was formulated, based on the deformation modes of beam elements. However, a superelement with more than two interface nodes has not been derived. Therefore, arbitrarily shaped bodies that are connected to more than two other parts cannot be easily modeled in the generalized strain formulation.
This paper presents a multinode superelement by introducing an implicit element frame, using the relations derived for the corotational superelement in [27]. The coordinates of the element frame are not part of the degrees of freedom, but the coordinates can be obtained from the degrees of freedom with a Newton-Raphson iteration. Deformation modes are defined using the local coordinates of the nodes to make the superelement applicable in the generalized strain formulation. Section 2 summarizes the generalized strain formulation and introduces the notation used throughout the paper. Section 3 formulates the superelement. An expression for the deformation modes can be chosen by the user, Sect. 4 shows three general ways to define these modes. The superelement is validated with examples in Sect. 5. Fig. 2 Generalized strain formulation of a two-dimensional beam element, L 0 is the undeformed length, ε 1 defines the elongation, and ε 2 , ε 3 define the bending

Summary of the generalized strain formulation
This section presents the generalized strain formulation, using the two-dimensional beam element in Fig. 2 as an example. Detailed derivations for this specific element will not be given as the purpose of this section is only to give an impression of the generalized strain formulation in general. The details for many element types can be found in the literature cited in the introduction. (2.1)

Notation for coordinates
The lower index "All" will be used to define all coordinates that define the configuration of an element. These coordinates are also denoted by x, In the three-dimensional case, x and q O,O All have a slightly different meaning: in q O,O All the orientations are expressed as the finite rotations around the x, y, and z-axes where the orientations in x are expressed using Euler parameters [40]. Note that due to the nonvectorial nature of rotations, the vector q O,O All does not exist in three dimensions. However, only the virtual change of this vector is used in the derivations.

Stiffness in terms of deformation modes
The generalized strain formulation defines the deformation of the element using deformation modes. The generalized coordinates associated to these deformation modes are called generalized deformations and denoted by ε. In case of the two-dimensional beam, three deformation modes can be defined of which one defines the elongation and two define the bending (see Fig. 2). The generalized deformations are explicit functions of the nodal coordinates, 3) The generalized force associated to the deformation modes is denoted by σ and is related to the generalized deformations by the constitutive law. If this relation is linear, it can be written as σ = Sε, (2.4) where S is a constant stiffness matrix.

Inertia
The inertia is modeled based on the coordinates x such that for the resulting forces f on all nodes we can write where M (x) is the mass matrix and h (x,ẋ) contains the convective inertia terms. The dot and double dots on x define the first and second derivative with respect to time, respectively.

Equation of motion
Using the principle of virtual work, we obtain where f a is the force that is applied on the nodes, or the reaction force in case of prescribed coordinates, and (. . . ) T defines the transpose; δ denotes the virtual change of a variable. This is the equation of motion of one element in terms of the absolute nodal coordinates in combination with the generalized deformations. One way to solve this is by substituting δε = D ,x δx, in which D ,x defines the derivatives of the generalized deformations which can be obtained analytically for each element, In this case the generalized deformations are only used implicitly. The equation of motion can also be defined for a set of degrees of freedom that include (part of) the generalized deformations. In this way (part of) the generalized deformations of the element can be constraint which allows the modeling of (partly) rigid bodies. This is detailed in [41,42] and Appendix A of [39].

Derivation of the superelement
This section derives the GMS. The configuration of the GMS is defined by the absolute coordinates of the interface nodes and by generalized coordinates of any internal modes. Together these are the configuration coordinates. An element frame j defines the rigid body motion. The coordinates of the element frame are not part of the configuration coordinates. They do not appear in the equation of motion, but they can be determined for a given set of absolute configuration coordinates. Once the position of the element frame is computed, the other element-dependent functions and matrices can be derived, e.g., the generalized deformations (which are a function of the local nodal coordinates) and the mass matrix.
Derives a relation between the virtual change of the element frame and local configuration coordinates in terms of absolute configuration coordinates  Table 1 shows an overview of the steps in derivation of the GMS in this section. Section 3.1 relates the local configuration coordinates (i.e., the coordinates with respect to the element frame) to the absolute configuration coordinates and the position of the elementframe. Section 3.2 defines the displacements in terms of these local coordinates. In Sect. 3.3 these displacements are used to define the generalized deformations and the position of the element frame. Sections 3.4 and 3.5 present some relations between the virtual change of the different coordinate types. The stiffness and inertia terms are derived in Sects. 3.6 and 3.7, respectively.

Local configuration coordinates in terms of absolute configuration coordinates
is the virtual change of the orientation of frame k and the tilde defines the skew-symmetric matrix: The virtual change of the absolute position and orientation of node k can be expressed in terms of the local coordinates and the change of the element frame. This is derived in Appendix A.1. Combining Eqs. (4.4) and (4.5) gives with the definitions: Note that terms between square brackets can have a slightly different meaning than the same term without brackets in this paper. Combining Eq. (3.3) for all interface nodes gives where the subscript "IF" refers to the coordinates of all interface nodes, and: Henceforth a single upper index, as used in ˆ j rig , defines the frame in which the variable is expressed, unless the index is explicitly specified differently.
In Fig. 3, the deformation of the right most part of the GMS is typically not affected by the positions of the interface nodes. This deformation can be described using internal modes with generalized coordinates q int . The generalized coordinates of internal modes are not expressed in the orientation of a specific frame such that their values in local and absolute coordinates are equal. These coordinates can be added to the vector with all interface coordinates in Eq. (3.5), where the subscript "All" refers to the configuration coordinates, and: This equation relates the absolute configuration coordinates to their local coordinates through the coordinates of the element frame. Equation (3.7) can also be rewritten to express the local coordinates in terms of absolute coordinates: (3.9)

Displacements in terms of local configuration coordinates
The displacement of interface node k, expressed in the orientation of element frame j , is denoted by p j,j k , see Fig. 3. It is composed of displacements and rotations, where r j,j k is the undeformed position of k with respect to the element frame and ψ j,j k defines the orientation of k. In the undeformed configuration, the local orientations of the interface nodes are defined to be zero. The rotation ψ j,j k will be specified by means of the local rotation matrix. A rotation matrix can be defined as the matrix exponential of the skewsymmetric matrix of the rotation vector [43]: where n j,j k is the unit rotation axis and φ j k the magnitude of the rotation k with respect to element frame j . By assuming that the elastic rotation of node k is small, the rotation matrix can be approximated by the first-order Taylor expansion (3.12) Based on this approximation, the rotation ψ j,j k will be implicitly defined using the offdiagonal terms of this local rotation matrix, The virtual change of this rotation can be related to the virtual change of the local coordinates of node k, see Appendix A.2. Based on Eq. (A.8), the virtual change of the displacement p j,j k can be expressed as (3.14) The matrix H j k is defined in the appendix and equals the identity matrix for zero rotation of the node k. This equation can be combined for all interface nodes: (3.15) The displacements of the internal modes are defined by their corresponding coordinates q int . Combining this relation with Eq. (3.15) gives an expression for the virtual change of all displacements in terms of the change of the local configuration coordinates:

Generalized deformations and the position of the element frame
This section relates the generalized deformations ε to the displacements that are derived in Sect. 3.2. This will also result in a relation for the coordinates of the element frame. The displacement vector p j,j All describes the elastic deformation. However, it also describes the six rigid body motions as it includes the displacements of all interface nodes. Therefore, p j,j All can be linearly related to six rigid body motions in combination with the elastic deformations that are described by ε: where η rig are the six coordinates of the six rigid body motions and the constant matrix The matrix with deformation modes, j flex , can be defined by the user. Section 4 describes three general methods to define these deformation modes. In the remaining part of this section we will assume that this matrix is known. Also the matrix j rig0 is known at the start of the simulation as it can be computed based on the local positions of the interfaces of the undeformed element. This means that also the matrices V  (3.19) where N All is the number of configuration coordinates, N IF is the number of interface nodes, and N int is the number of internal modes. The rigid body motion of the element is described by the coordinates of its element frame. It can therefore not also be described by the rigid modes as this will result in a singular system. This means that η rig should be zero which defines six constraints on p j All : (3.20) Based on these six constraints, we can find the position and orientation of the element frame for a given set of absolute configuration coordinates q O,O All . However, an explicit relation does not exist in general, so that it has to be solved based on a Newton-Raphson iteration. Substituting Eq. (3.16) in the virtual change of the constraint in Eq. (3.20) gives We want to find the position of the element frame for a given set of absolute coordinates, i.e., δq O,O All = 0. Therefore, Eq. (3.9) can be used to obtain Using this equation, we can update the position of the element frame using the following Newton-Raphson procedure: The hat on q O,O j emphasizes that this vector fundamentally does not exist. As noted in Sect. 2.1, the rotation in this vector is parameterized by finite rotations. This does not work for large rotations in three dimensions. However, the orientation can be defined by a rotation matrix or Euler parameters, and the vector can be updated in an equivalent way.
Once the position of the element frame is found, the displacements p j,j All can be obtained using Eq. (3.10) after which the generalized deformation is obtained from Eq. (3.18), (3.24)

Virtual change of the element frame and local configuration coordinates as function of absolute configuration coordinates
Section 3.3 defined the position of the element frame based on the absolute coordinates using the six constraints in Eq. (3.20). Once the element frame is in the right position and orientation, these constraints can also be used to write the virtual change of the element frame as a function of the virtual change of the absolute coordinates. The constraints imply that also their virtual change should stay zero as defined in Eq. (3.21). By substituting Eq. (3.9) into this relation, we obtain The deformation of the GMS is assumed to be small, for this small deformation, the matrices H j and j rig will typically only change slightly. This indicates that the term V j rig H j j rig is close to the identity matrix, which implies that it is invertible. Therefore, the equation can be rewritten to relate the virtual change of the element frame to the virtual change of the absolute configuration coordinates: (3.26) A physical interpretation of the 6 × 6N IF matrix Z j is that it defines the rigid body motion as a function of an arbitrary motion expressed in the local frame. By substituting Eq. (3.26) into Eq. (3.9), we obtain the change of the local configuration coordinates as a function of the absolute configuration coordinates: A physical interpretation of T j is that it removes the rigid body motion from an arbitrary motion, leaving the flexible motion of the coordinates. More elaborate geometric interpretations of the matrices j rig , Z j , and T j are given in [44].

First derivative of the generalized coordinates
This Note that δq O,O All contains the virtual change of finite rotations for each interface node. Rotations in three dimension should be specified using a parameterization like Euler angles or Euler parameters. In this paper Euler parameters are used. The orientation of node k with respect to frame O will be parameterized by λ O k . The absolute coordinates with this parameterization are given by x: The virtual change of finite rotations can be related to the virtual change of the parameterization (see, e.g., [45]) resulting in an equation like where G is a function that is in case of Euler parameters, (3.33)

Stiffness matrix
The GMS uses the stiffness and mass matrix of a linear finite element model that is reduced using Craig-Bampton modes [7] (i.e., boundary and internal modes). Note that the vector with displacements, p j,j All , indeed consists of these boundary displacements and internal modes. The result of the reduced model, expressed in the orientation of element frame j is like where M j All is the constant reduced mass matrix, K j All the constant reduced stiffness matrix, and f j All defines applied forces on the modes. The potential energy of this model is By substituting Eq. (3.17) with η rig = 0, the stiffness matrix in terms of the deformation modes, S, can be obtained, which can be used in the equation of motion, Eq. (2.7):

Inertia terms
The inertia terms can be derived in different ways. This paper presents two approaches: the "corotational inertia" defines the energy based on the corotated mass matrix and derives the inertia vector using Lagrange's equation. The "full inertia" derives the global acceleration of the material points in the body and obtains the inertia terms by integrating these accelerations over the volume. Both approaches neglect the higher-order terms in deformation and result in the same mass matrix, but a different convective inertia. The approaches are consistent to two of the approaches described in [10].

Corotational inertia
The kinetic energy of the reduced linearized finite element model is Based on Lagrange's equation, the total inertia forces, H, can be defined as function of the kinetic energy. These inertia forces should equal the inertia forces defined in Eq. (2.5): Substituting Eq. (3.38) and rewriting gives the convective inertia as The full expression is derived in Appendix A.3. This result is similar to the result obtained in the superelement of [13].

Full inertia
For an initial undeformed configuration, the global velocity or virtual change of an arbitrary point s in the GMS can be written in terms of its absolute nodal coordinates: where the 3 × N All matrix j All,s defines the mode-shapes evaluated at position s. Note that the mode-shapes contain a linear combination of all the Craig-Bampton boundary modes such that rigid body motion is also included.
Based on the principle of virtual work, the total inertia force is implicitly defined by the volume integral where ρ is the material density. The resulting total inertial is derived in Appendix A.4, resulting in Eq. (A.33): The first term in this expression is the global mass matrix times the acceleration, see Eq. (3.38). The convective inertia vector according to the full approach is  30), where it can be seen that it involves an integral that cannot be computed from the finite element matrices that are commonly available in a linear finite element analysis. The consistent derivation of this integral requires evaluating a specific term for each element in the finite element model, but the integral can also be estimated by using a lumped mass approximation.

Comparison
The second terms of the convective inertias of both approaches as defined in Eqs. However, the remaining terms in both approaches are different. This difference exists because the corotational approach uses the energy in the discretized, reduced form to obtain the inertia forces, where the full approach derives the inertia forces from the continuum and applies the model order reduction afterwards. Figure 4 visualizes this. In both approaches the total energy is conserved. However, the corotational approach implicitly assumes that the inertia forces can be written in terms of the reduced mass matrix M j All . The full approach shows that the exact evaluation of the inertia forces requires the term N j All which cannot be expressed in terms of the reduced mass matrix. Sections 5.1 and 5.2 of this paper further evaluate the differences. In [46] (Sect. 5.3) a more elaborate derivation of the inertia terms is given, which also shows that the inertia terms cannot be written in terms of the finite element mass matrix.

Local interface displacements
The first option directly relates the generalized deformations, ε, to the local displacements of all interface nodes except one. If we, for example, exclude the first interface node, the matrix with flexible modes becomes This choice causes the orientation of the frame to be the orientation of the remaining interface node. Therefore, a logical choice in combination with these deformation modes is to place the element frame in the remaining interface node. A disadvantage of this option is that the result will depend on the interface node chosen. An advantage is that the position of the element frame does not have to be found by the Newton-Raphson iteration outlined in Sect. 3.3 as its coordinates are simply the coordinates of the related interface node. This also simplifies some of the other relations that have been defined in Sect. 3.

Natural modes of the free body
For this option the natural modes are extracted from the reduced model with free motion by solving the eigenvalue problem The vector j free−rig defines the first six eigenmodes, which are rigid body modes with zero eigenvalues. Note that this matrix is not necessarily exactly the same as j rig0 , but is spans the same space. The modes j free−flex are defined to be j flex . An advantage of this choice is that the mode-shapes define the natural modes which means that the modes with the higher natural frequencies can be constraint in the GMS. Another advantage is that the stiffness matrix as derived in Eq. (3.36) becomes diagonal, which simplifies the evaluation of the stiffness equation (Eq. (2.4)), where the vector ω defines the eigen frequencies, and the matrix diag ω 2 is part of the diagonal matrix . The element frame can be positioned anywhere, and its position will not affect the results. A classic choice is to place it at the center of mass.

Frame attached to a material point
This option gives the position of the element frame a physical meaning, it is attached to a material point. This option was also used in [27]. The frame is chosen to be in the center of mass. Based on the finite element model, the displacement and rotation of the material point that is located at the position of the element frame, p j,j FFR , can be expressed as a function of the displacements p j,j All . This displacement should be zero if the element frame is attached to this material point:  The matrix j flex should be defined in such way that the remaining part of Eq. (4.5) holds. The natural modes of the free motion, j free−flex , will be used to define this flexible modes. However, in order to make sure that the top right term of Eq. (4.5) holds, the effect of rigid body motion should be subtracted from these modes. Note that the matrix T j as defined in Eq. (3.27) subtracts the rigid body motion from an arbitrary motion. The flexible modes will therefore be defined using this matrix in undeformed configuration: (4.6) By using the relation V j FFR j rig0 = 1, it can be shown that this indeed satisfies the rightupper part of Eq. (4.5): The stiffness matrix of this method equals the stiffness matrix of the method in Sect. 4.2. This can be shown by the fact that the stiffness matrix multiplied by rigid body modes equals zero, K j All j rig0 = 0, and, using Eq. (4.6), This means that this method shares the advantages with the method in Sect. 4.2: the mode shapes are related to the natural modes and the stiffness matrix is diagonal. However, the mass matrix of this method is different from the mass matrix in Sect. 4.2.

Validation
The GMS is validated using the multibody software SPACAR [38,41]. A rigid rotating beam demonstrates the importance of the convective inertia terms. The application of the GMS in dynamic simulation is shown by a slider-crank case. A static cantilever beam shows the effect of different definitions of the element frame. In these first three examples, the mass and stiffness of the GMS are obtained using a finite element model of beam elements. These examples are validated using the beam elements defined in [4]. Deformation due to shear is neglected in these examples. In the fourth and fifth case, the GMS is used in a spherical flexure joint and a misaligned cross-flexure to show the usefulness in flexure based mechanisms. The flexures are modeled with the beam element described in [47].

Rigid rotating beam
This section shows an example to evaluate the convective inertia terms. Figure 6 shows a beam that rigidly rotates around point A, which moves with a constant velocity v in the xdirection. The beam is modeled using a single GMS with either one or two interface nodes. The first interface node is positioned in point A, the second optionally in point B or C. The element frame is fixed to the material point in the center of the beam (see Sect. 4.3). Table 2 shows the required force F on the beam when it is horizontally oriented (i.e., the position shown in figure). The results obtained with the corotational inertia and the full inertia correspond to the centrifugal force on a rotating beam, i.e., mLω 2 /2. In absolute nodal coordinates based finite element simulations, the convective inertia is generally neglected. For some element-types this gives exact results, for other elements this only results in a small error if the elements are small. However, in this rotating beam example, using no convective inertia (i.e., using only the mass matrix times the accelerations) does only give correct results if the center of mass is exactly in the center of both interface nodes. This illustrates the importance of the convective inertia term. Table 3 shows inertia terms to compare the corotational inertia with the full inertia. These are the terms in two dimensions, so the inertia forces on each interface point contains three terms: for the translational x-direction, the y-direction, and the rotational direction around the z-axis. The total inertia force H in the x-direction always equals −mLω 2 /2 for both approaches as also given in Table 2. The corotational inertia results in a rotational term which dependents on the overall velocity v. This is a nonphysical result, in the first place because the overall velocity of a mechanism should not affect the inertia forces. Secondly, because a rigid rotating component should only experience centrifugal inertia forces. In this case with a rigid beam there is no effect as the extra bending moments at both interface nodes cancel each other. However, in a flexible beam element, the extra bending moments will affect the bending deformation of the element as shown in Sect. 5.2. The full inertia only results in inertia in the x-direction and is independent of the overall velocity v.

Two-dimensional slider-crank
This example evaluates the accuracy of the GMS in a dynamic simulation. A twodimensional slider-crank problem that was also analyzed in [27,42,48] is shown in Fig. 7. Its physical properties are given in Table 4. The rigid crank is initially horizontally oriented  to the right, and rotates with a constant angular velocity of 150 rad/s. The flexible connector between the crank and slider is initially undeformed and has an initial velocity corresponding to the velocity of the crank (i.e., its initial velocity is a clockwise rotation of 75 rad/s around the slider). The mass of the slider is half of the mass of the connector. Figure 8 shows the midpoint deflection of the connector perpendicular to the undeformed connector divided by the length of the connector. The connector is modeled in seven different ways: a. 10 serial connected beam elements, this serves as reference-case; b. 2 serial connected beam elements; c. 2 GMSs with full inertia, the result is identical to the case where no convective inertia is modeled; The following observations can be made: • The GMS with full inertia gives almost the same results as the GMS without convective inertia (plotted by a single line). This indicates (together with the results in the rotating beam problem) that the convective inertia terms can be neglected in beam-like components; • The errors of the GSM with full inertia and without convective inertia (case c) with respect to the reference case, are in the same order as the errors obtained by beam elements and the corotational superelement (cases b and f). • The GMS with corotational inertia gives a significant different deflection compared to the other results. Using four elements instead of two, the results are much closer to that of the other simulations. This indicates that the corotational inertia method gives small inaccuracies if the elements are large.

Static equilibrium of a cantilever beam
This example evaluates the influence of the frame position on the accuracy and computation time. A hollow circular cantilever beam is subjected to a vertical tip force. The length of the beam is 1 m, the outer radius of the cross-section 0.01 m and the wall thickness 0.001 m. The Young's modulus is 70 GPa. Figure 9 shows configurations obtained by three different methods: the GMS, a beam element [4], and the corotational superelement of [27]. In all three cases, the beam is modeled using three serial connected elements, a reference is obtained by ten serial connected beam elements. The modes of the GMS are defined using the option described in Sect. 4.3: "frame attached to a material point" and the frame is located in the center of the element. For a force of 10 000N, the result for the GMS did not converge due to the large deformation in the left most element. In general, all three methods give the same results. Small differences are visible for the larger deformations because this results in large deformation per element. The difference in the results between the GMS and the corotational superelement are caused by the matrix H j that was introduced in Eq. (3.14).
In the derivation of the corotational superelement, this matrix was neglected by assuming small deformations. One of the advantages of the GMS is that the position of the element frame can be defined in different ways, see Sect. 4. Figure 10 shows results to study the effect of different positions. In the first three frame-options, the frame is placed in an interface point and the modes are chosen by the option of Sect. 4.1 ("local interface displacements"); in the fourth frame-option, the method of Sect. 4.3 is used ("frame attached to a material point").
The results indicate that positioning the frame at the center of the element gives the most accurate results. The most important reason is that the elastic rotational displacements are the smallest in this case. Defining an extra interface node in the center of the elements increases the number of degrees of freedom and significantly increases the computation time, especially when many elements are used.
Placing the element frame in the center without defining an extra interface node requires to update the frame in each loadstep using the Newton-Raphson procedure defined in Sect. 3.3, but this only slightly increases the computation time. Figure 10 gives the total computation time that was required for these Newton-Raphson updates. These times are approximately the same as the time difference between the total computation time of case d and the total computation time of the cases a and b. On average, 3.8 iteration steps were required to find the coordinates of an element frame. Figure 10 shows the total number of iteration steps in one simulation, where one step took on average 7.3 · 10 −5 s. Figure 11 shows the serial stacked spherical joint that was introduced in [2]. The most important dimensions of this flexure joint are given in Table 5. The flexure joint consist of six folded flexures which are connected to three frame parts. These folded flexures are placed in such way that lines through the folds coincide in the center of the joint. In this way,   Figure 12 shows the support-stiffness in the vertical direction (z-direction). The results indicate that compliance of the frame parts is significant with respect to the total compliance and that this compliance can be modeled accurately using the GMS. Figure 13 shows eigenfrequencies for the case that the base and end effector are fixed to the ground at their triangular-shaped face. The first three eigenfrequencies are rotations of the ring. These eigenfrequencies are related to a low stiffness and therefore almost not influenced by the flexibility of the connecting parts. The other three eigenfrequencies are influenced by the flexibility of the connecting parts which is modeled accurately using the GMS. Figure 14 shows an overconstrained cross-flexure with a misalignment in the overconstraint direction that was described in [49,50]. The two flexures are each modeled using eight beam elements with torsional warping (as defined in [47]) to model the thin part and one beam element to model the thick part which is used for attachment to the frame parts. The flexures are made of steel (Young's modulus 200 GPa, Poisson ratio 0.3), have a thickness   Figure 15 shows the first natural frequency (in which the shuttle rotates around the indicated rotation axis) as function of the misalignment. The results show that the compliance of the shuttle has significant influence on the result and this effect can be modeled with the GMS. Also the shear deformation in the beam elements has a significant effect. The exclusion of shear has a similar effect as a rigid shuttle on the first natural frequency.

Conclusions
A superelement has been presented which can be used to model arbitrarily shaped parts with multiple interface nodes in the generalized strain formulation. The deformation is defined linearly with respect to a local frame, where rotational displacements are defined using the off-diagonal terms of local rotation matrices. The coordinates of the frame are not part of the degrees of freedom, but can be obtained by a Newton-Raphson iteration, as a function of the degrees of freedom. This frame can be defined in multiple ways. Simulations show that this definition of the frame can have significant influence on the results. More accurate results are obtained if the elastic rotations with respect to the element frame are small. Two methods are presented to define the inertia: simulations show that the "full approach" gives more accurate results than the "corotational approach," however, the full approach includes terms that cannot be derived from a standard reduced finite element model. The paper shows that complex components with slender parts can be modeled accurately using a proper combination superelements and beam elements.

Appendix A: Derivations
This appendix shows some derivations for the formulation presented in Sect. 3.

A.1 Relating the virtual change of absolute and local coordinates of an interface node
This section shows how the virtual change of the absolute position and orientation of an interface node can be expressed in terms of the virtual change of its local coordinates and the virtual change of the coordinates of the element frame.
An absolute position r O,O k of node k can be defined by means of the coordinates of the element frame (see Fig. 3): where R O j defines the absolute orientation of element frame j . The virtual change of a rotation matrix can be expressed as where δθ O,O j defines the virtual change in finite rotations of frame j with respect to the global frame O. The tilde defines the skew-symmetric matrix of a vector which is related to the cross-product such that, for two arbitrary 3 × 1 vectors a and b, the following relations hold:ã The where n x , n y and n z are unit-vectors in the x, y, and z-direction, respectively. The virtual change of the first term can be expressed using the virtual change of a rotation matrix as For undeformed elements, the matrix R j k equals the identity matrix which means that the matrix H j k also equals identity for undeformed elements.

A.3 Derivation of the corotational inertia
This section derives the full expression for the corotational convective inertia as given in Sect. 3.7.1. The corotational convective inertia vector, as defined in Eq. (3.40) is where ρ is a function, namely (A.20) The terms for the individual interface rotations can be combined for the full matrix into a single equation such that the second term in δ [B]ẋ becomes

A.4 Derivation of full inertia
This section derives the full inertia, which was implicitly defined in Sect. 3