Derivation of a superelement with deformable interfaces – applied to model flexure joint

Design and optimization, as well as real time control, of flexure mechanisms require efficient but accurate models. The flexures can be modelled using beam elements and the frame parts can be modelled using superelements. Such a superelement efficiently models arbitrarily shaped bodies by few coordinates, using models obtained by model order reduction. The interfaces between the frame parts and the flexures often experience considerable deformation which affects the stiffness. To define the interface deformation in a reduced order model, this paper derives a multipoint constraint formulation, which relates the nodes on the deformable interface surface of a finite element model to a few coordinates. The multipoint constraints are imposed using a combination of the Lagrange multiplier method and master–slave elimination for efficient model order reduction. The resulting reduced order models are used in the generalized-strain multi-node superelement (GMS) that was defined in (Dwarshuis et al. in Multibody Syst. Dyn. 56(4):367–399, 2022). The interface deformations can be coupled to the cross-sectional deformation of higher order beam elements (i.e. beam elements of which the deformation of the cross-sections is explicitly taken into account). This paper applies this technique to model flexure joints, where the flexures are modelled with beam elements, and the frame components and critical connections using the GMS. This approach gives generally over 94% accurate stiffness, compared to nonlinear finite element models. The errors were often more than 50% lower than errors of models which only contain beam elements.


Introduction
Design and optimization, as well as real-time high bandwidth control, of flexure mechanisms require efficient but accurate models.The often long and slender flexures can be modelled using beam models.Sophisticated beam elements [2,3] for the modelling of flexures have been derived and implemented in the generalized strain formulation [4].The frame parts in flexure mechanisms, which often have a complex shape, can be modelled efficiently using superelements.A superelement linearly describes deformation of arbitrarily shaped parts with only a few coordinates.In [1] a superelement has been formulated in the generalized strain formulation, referred to as the generalized-strain multi-node superelement (GMS).It has been applied to model the frame parts of flexure mechanisms, showing efficient modelling with relatively good accuracy.However, the interfaces of the GMS were defined to be rigid, whereas in reality the interfaces between frame parts and the flexures often experience considerable deformation.For better accuracy, the deformation of the interfaces of the frame parts that are modelled using the superelement should be taken into account.
The stiffness and inertia properties of superelements are generally obtained using model order reduction methods.These methods reduce the number of degrees of freedom in the finite element model of a component by describing the deformation using a limited number of deformation modes.Overviews of the different methods can be found in [5][6][7].The surfaces of the component to which other components are connected are called interface surfaces.The deformation of an interface surface in the finite element model is described by the displacements of all the nodes on the surface.Most conventional model order reduction techniques (e.g. the techniques proposed by Hurty [8] and Craig and Bampton [9]) take all these displacements into account in the reduced model.However, this results in large reduced order models if there are a lot of nodes on the interface surface.
Interface reduction methods can be used to reduce these models further.They describe the deformations of the interface surfaces using a limited number of modes.Overviews of the different methods can be found in [10,11].One method is to perform interface reduction on the assembled reduced model [12].One disadvantage is that this assembled model may have many degrees of freedom, such that this reduction can take a lot of computation time.Moreover, as the reduction is applied after assembly, modifying one of the components requires recomputing the reduction.One solution is to apply interface reduction on subsets of the full assembled reduced model [13].Another solution is to reduce the interfaces of the individual components before assembly [14], which involves two challenges.First, it requires reduced order models of the components with a sufficiently accurate description of the deformation around the interfaces, as, for example, discussed in [15,16].Second, the compatibility between the interfaces of two connected components should be enforced, as, for example, discussed in [17].Sometimes the compatibility cannot be simply enforced, for example, in a dissipative interface [18,19].
Considering an interface to be completely rigid allows describing the displacement of the interface by the displacement and rotation of a single master node [20], which is often called condensation node.The existence of the condensation node in this approach makes the resulting reduced order models suitable for use in multibody analysis (see, e.g.[21]), because the positions of the condensation nodes of two connected components can be coupled in a geometrically nonlinear analysis.This approach was also applied in the GMS [1].
The 'prior basis function method' [10,22] extends this method by adding a linear combination of deformation fields to the rigid interfaces in order to describe interface surface deformation.This approach is used in the current paper.These deformation fields are hereinafter called interface deformation fields.The generalized coordinates that describe the amount of this deformation, together with the six coordinates of the condensation node, are called the condensation coordinates.
The compatibility between two components can be imposed by choosing the same interface deformation fields for both components and relating the corresponding condensation coordinates.This approach has been applied on a two-node superelement [23].However, apart from in this paper, interface deformation in geometrically nonlinear multibody simulations has rarely been investigated in literature to the best of the author's knowledge.
The relation between the nodes on the interface surface and the condensation coordinates is called a multipoint constraint.References [21,[24][25][26][27][28][29] show formulations for a multipoint constraint without interface deformation fields.Two types are considered: the 'rigid multipoint constraint' and the 'interpolation multipoint constraint'.Both types are extended in this paper in order to apply them for deformable interfaces.The results of the interpolation multipoint constraint are compared with the results of existing literature.The rigid multipoint constraint will be referred to as 'exact multipoint constraint' because the term 'rigid' is confusing in this context.
• Exact multipoint constraint (EMPC).The interface surface displaces and deforms exactly as prescribed by the displacements of the condensation node and interface deformation fields.• Interpolation multipoint constraint (IMPC).The interface surface is completely free to deform.The condensation coordinates follow the average motion of the surface.
The multipoint constraints can be imposed using the penalty function method, the Lagrange multiplier method (also referred to as the 'dual formulation') and master-slave elimination (also referred to as the 'primal formulation') [24,25].
• A disadvantage of the penalty function method is that it requires the selection of a suitable penalty factor.This selection is nontrivial, compromising between accuracy and computational stability [25,30].• A disadvantage of the Lagrange multiplier method is that it increases the number of unknowns, whereas master-slave elimination decreases the number of unknowns, both in proportion to the number of constraint equations.This is not a big issue for the IMPC as the number of constraint equations of this constraint is much lower than the number of degrees of freedom in the finite element model.However, for the EMPC, the number of constraint equations can be much higher as it scales with the number of nodes on the interface surface.• A disadvantage of master-slave elimination is that it requires a suitable selection of a set of dependent coordinates to avoid singularity in the equations.For the EMPC, all nodes on the interface surface are dependent.However, the selection is nontrivial in case of the IMPC.For the IMPC the selection can be based on physical insights [24,28], but this is shown only for multipoint constraints without interface deformation fields.The selection can also be avoided by computing the null-space of the constraint relations [21].However, this may require a lot of computation time and has a negative effect on the sparsity of the constrained finite element matrices.
To a large extend, all these disadvantages can be avoided by using master-slave elimination to impose the EMPCs and the Lagrange multiplier method to impose the IMPCs, which will be detailed in this paper.
In higher-order beam theories, the deformation of the cross-section of the beam is taken into account [31][32][33][34].This deformation can be interpreted as interface deformation.The cross-sectional deformation is especially important in thin-walled beams with an open profile [35][36][37][38].The resulting beam elements have extra degrees of freedom at both nodes that define the amount of cross-sectional deformation.This deformation can be coupled to the deformation of the interface surface of a superelement if the deformation fields are equal.This method was applied in the linear structural analysis of frames [39,40] and for concept modelling of vehicles [41,42].
This paper shows how interface deformation can be defined in the geometrically nonlinear superelement (GMS) of [1].To obtain a reduced order model for the GMS, a multipoint constraint for deformable interfaces is derived and imposed efficiently using a combination of Lagrange multipliers and master-slave elimination.The interface deformation of the GMS can be defined consistently with that of the cross-sectional deformation of connected higher order beam elements.This is applied in order to analyse flexure joints using a combination of GMSs and beam elements.
These models are more efficient to evaluate than using nonlinear finite element models, because the models contain far less elements.Although the superelements may be derived from complex finite element models, which implies that the model order reduction may take considerable computation time, this time will be significantly less than the time to evaluate a finite element model of the full flexure joint for the following reasons: • The model order reduction of the frame parts only requires linear finite element calculations, where the finite element model of the full flexure joint with large deformation requires nonlinear calculations.• Only part of the flexure joint has to be modelled using finite element models.
• In a dynamic simulation, a full finite element model should be evaluated at each time step, where the model order reduction has to be applied only once.
The cost of this increase in efficiency is that the models will be slightly less accurate.However, the models can be significantly more accurate than the models which consider the frame parts to be completely rigid, as validated in this paper.
Section 2 shows how a reduced order finite element model with deformation of interfaces can be obtained using multipoint constraints.Section 3 summarizes the GMS formulation and explains how the reduced order model can be used in this superelement.Section 4 briefly introduces higher order beam theory, to show how the cross-sectional deformation is related to interface deformation of the superelement.The formulation is validated with five examples in Sect. 5.

Reduced finite element model
This section defines how the Craig-Bampton reduced order model including interface deformation is obtained.Section 2.1 defines the multipoint constraints and Sect.2.2 derives the constraint equilibrium equation from which the reduced order model is obtained.

Multipoint constraints
Figure 1 shows a finite-element model (henceforth FE-model), with three interfaces.The nodes of the model will be referred to as FE-nodes.The frame i is the global reference frame of the FE-model.The interface surfaces are coloured dark grey and to each interface a condensation node is attached, visualized by a frame.The vector r i,k g defines the undeformed position of an FE-node g (lower index) with respect to condensation node k (second upper index), defined in the orientation of frame i (first upper index).The vector u i,i k defines the global displacement of node k, and ψ i,i k defines its rotation, which are assumed to be small.For a rigid interface the displacement of each FE-node, g, on the interface surface can be predicted linearly based on the displacement of the condensation node.This predicted displacement ûi,i g can be expressed as The tilde defines the skew-symmetric matrix which is related to the cross-product, such that for two arbitrary 3 × 1 vectors a and b, the following relations hold: Deformation of the interface surface can be added to this predicted displacement using interface deformation fields ω r k,k g .These fields are a user-defined function of the position on the interface surface.The fields are multiplied by coordinates, collected in the vector α k .The resulting predicted displacement can be written as where R i k is the rotation matrix that defines the rotation of node k with respect to the global frame, and This is summarized in Fig. 2. The effect of the displacement of the condensation node and the deformation fields is combined to obtain the predicted displacements.These predicted displacements can be used to define two types of multipoint constraints.The first type, the exact multipoint constraint (EMPC), implies that the displacements on the interface surface should equal the predicted displacement, resulting in three constraint equations for each FE-node on the interface surface: The second type of multipoint constraint is the interpolation multipoint constraint (IMPC).It is defined as follows: the interface surface is free to deform and the condensation coordinates are defined such that the squared error between the expected displacements and the real displacements of all FE-nodes on the interface surface is minimized.The squared error is defined in a cost function as where w g is a weighting factor for the FE-node, which should be chosen proportional to the part of the interface surface that the node represents.The constraints that minimize this cost 3) and then enforcing the derivative to be zero: (2.7) The IMPC will typically underestimate the stiffness and is the most logical choice if the interface is connected to a component which is more flexible around the interface.On the other hand, the EMPC will typically overestimate the stiffness and is the most logical choice if the interface is connected to a component which is much stiffer around its interface.This difference between the two types of multipoint constraints was also noted in [24].
Equations for the IMPC without interface deformation fields are also given in [29,43].These results are derived by assuming a relation between the forces on the FE-nodes on the interface surface and the forces and moments on the condensation node.The result is a quite long expression which is slightly different from Eq. (2.7).Another simplified relation for the IMPC is given in [24][25][26][27][28].In these papers, the condensation node is placed in the centre of the interface surface and the translational displacement and rotation of the condensation node is written as It can be shown that the translational displacement of the condensation node for both these methods corresponds to the constraint in Eq. (2.7), but the resulting rotation is different.To evaluate this difference, Fig. 3 shows two interface surfaces, both with four FE-nodes that have equal weighting.Condensation node k is placed in the centre of the interface surfaces.The four FE-nodes are displaced in the x-direction, according to a rigid rotation β around the z-axis.The resulting rotations of condensation node k are given in Table 1.The IMPC of Eq. (2.7) as used in this paper gives the expected rotation β around the z-axis for both shapes, whereas the equations used in the other papers give results that depend on the width b.

Unconstrained static equation
The equations of an unconstrained FE-model can be written as where M i FEM and K i FEM are the mass and stiffness matrix of the FE-model; u i,i FEM , üi FEM and F i FEM are the displacements, accelerations and forces of the FE-nodes respectively, all expressed in the orientation of the global frame i of the FE-model.Besides the displacements of the FE-nodes, there are also condensation coordinates.The displacements of these coordinates will become the boundary displacements in the reduced method and are defined in a vector p i,i bnd as The unconstrained static equation in terms of all displacements can be written as where F i bnd are the loads on the boundary coordinates.The forces applied to the FE-nodes, F i FEM , are zero in the derivation of the Craig-Bampton reduced model.

Constrained static equation
The model can contain EMPCs and IMPCs.The equations for all IMPCs in the model, as defined in Eq. (2.7), are combined to (2.12) These constraints are applied to Eq. (2.11) using the Lagrange multiplier method, ⎡ ⎣ 0 0 where λ are the Lagrange multipliers.
The EMPCs will be solved using master-slave elimination.The equations of all EMPCs are combined to where the terms in B A come from the constraint equations as defined in Eq. (2.5).The vector u i,i f contains the displacements of the FE-nodes that are not on the interface surface of an EMPC, so B B is just a Boolean matrix that relates these displacements in u i,i FEM to u i,i f .Equation (2.14) can be used to write a relation in terms of all displacements and Lagrange multipliers, Applying this coordinate transformation to Eq. (2.13) gives

.16)
Computing the matrix-products results in the constrained static equation (2.17)

Craig-Bampton reduction
A Craig-Bampton reduced model contains boundary modes and internal modes.The boundary modes are related to the stiffness of the boundary-displacements p i,i bnd and can be obtained by Guyan reduction.The internal displacements can be expressed in terms of the boundary displacements using the last two rows of Eq. (2.17): This result can be substituted into Eq.(2.14) by which the boundary modes i bnd are obtained that relate the boundary displacements to the displacements of the FE-nodes, (2.19) The internal Craig-Bampton modes are the natural modes of the component where the boundary coordinates are fixed.These are obtained by solving the constrained eigenvalue problem of the inner part of the stiffness and mass matrix, namely Only the internal modes in the desired frequency range have to be selected.The internal modes i int for all FE-displacements are obtained using Eq.(2.14) by noting that the displacements p i bnd are zero for the internal modes, where φ i,desired are the modes φ i in the desired frequency range and η int is the vector with the generalized coordinates of the internal modes.
All the Craig-Bampton modes (Eqs.(2.19) and (2.21)) are combined into Using these modes, the reduced stiffness and mass matrices in the orientation of the global frame of the FE-model can be written as All . (2.24)

Implementation in the superelement formulation
This section summarizes the coordinates by which the configuration of the GMS [1] is defined and shows how the reduced model of the previous section can be used for this superelement.
Figure 4 shows a GMS with four interface nodes.It is defined with respect to global frame O.Note that this frame is different from the global frame in the FE-model that was used to obtain the reduced order model.The global position of a node k is defined by vector r O,O k (indices are defined in a similar way to the undeformed positions r i,k g of the FE-model in Sect.2.1).The rotation matrix R O k defines the orientation of node k (lower index) with respect to global frame O (upper index).It depends on three independent parameters.The exact parameterization is not relevant to this overview and is therefore not detailed.The six independent parameters that define the global position and orientation of node k are stored in vector q O,O k .The configuration of the GMS is fully defined by the absolute nodal coordinates of all interface nodes, {q O,O IF 1 , . . ., q O,O IF N }, in combination with the generalized coordinates of the internal deformations, q int .
The undeformed position and orientation of the GMS is defined by its element frame j .The coordinates of the element frame are dependent coordinates, i.e. they do not appear in the equation of motion.The user can define six relations that define the position and orientation of the frame as functions of the absolute nodal coordinates and internal coordinates, The simplest option is to define the frame in one of the interface nodes, but other options are also possible as detailed in [1].
After the position of the element frame is obtained, the local coordinates of each interface node can be obtained from its global coordinates as (3. 2) The elastic displacement of an interface node can then be obtained by subtracting the undeformed position from the local coordinates.The elastic displacement of the internal deformation modes equals the generalized coordinates, p j,j k = p j,j k q j,j k . (3.3) These displacements define the deformation of a GMS and should be related to the displacements of the reduced model, p i,i  All , which are defined in Sect.2.1.The displacements of the interface nodes, p j,j k , are the displacements of the condensation nodes in the reduced model.The displacements of the internal deformation modes in the reduced model, η int , are part of the internal deformation modes in the GMS.The displacements related to interface deformation fields, p α , do not explicitly appear in the displacements of the GMS.
However, the values of the warping coordinates do not depend on the position of the element frame, therefore they can be treated as internal displacements.The vector with all displacements can be written as These are the same displacements as in the reduced model in Eq. (2.23).However, they are defined in the orientation of a different frame.The displacements are related by a rotation matrix, in which [R i j ] consist of 2N IF times the 3 × 3 rotation matrix R i j and an identity matrix corresponding to the length of vector q int .The reduced stiffness and mass matrix can be expressed in the orientation of the element frame by applying this rotation to the matrices of Eq. (2.24), These matrices, in combination with the local positions of the condensation nodes in undeformed configuration are the required input to define a GMS.The displacements, p j All , are a result of the multibody simulation.These displacements can be applied to the constrained FE-model to obtain strain and stress results.
To ensure that the interface surfaces of two connected bodies in the multibody simulation match, the condensation nodes of their reduced order models should be defined in the same position with respect to the surface.Furthermore, the interface deformation fields of both reduced order models should be defined equivalent to ensure that the interface deformations of both components match.

Summary of higher-order beam elements
In conventional beam theory, the cross-section is assumed to be undeformed.Therefore, the global position of each point on the cross-section can be obtained by the coordinates of its elastic line.The position, r O,O g of a node g at the cross-section at side p can be written as (see Fig. 5) In the higher-order beam elements, deformation of the cross-section is typically added using multiple deformation fields which are functions of the position in the cross-section ω r p,p g , multiplied by axial coordinates α (s).The values of these coordinates at node p are denoted by α p .These coordinates are extra degrees of freedom at this node, which can be coupled to a connected beam element.The position at the cross-section on interface p can be defined by adding this deformation to Eq. (4.1), This formulation of the positions is similar to the expected displacement on the interface surface of a superelement, as defined in Eq. (2.3).These elements can be connected by coupling the node of the beam element to the interface node of the GMS.In order to enforce the compatibility of the interface surface deformation, the same interface deformation fields ω r p,p g should be chosen for both interfaces.In this paper we will only use one deformation field for the beam-elements, namely the axial warping caused by torsion.The corresponding warping fields for a thin rectangular cross-section and an I-profile are given in Fig. 6.

Validation
The GMS is validated using the multibody software SPACAR [44,45] with traditional beam elements [2] and beam elements in which warping due to torsion is included [38], referred Fig. 5 Higher-order beam element, to show how the global position of a node g on an interface can be expressed in terms of the beams coordinates Fig. 6 Warping fields of a thin rectangular cross-section and an I-shaped cross-section; the latter is an approximation based on thin-walled beam theory to as warping beam elements.In this section the GMS does not allow interface deformation when it is connected to a traditional beam element The GMS connected to a warping beam element contains interface deformation according to the connected element.The finite element models to obtain the stiffness properties of the GMSs are assembled using the PDE-toolbox in MATLAB.
First, the connection of two I-profiles gives an example of thin-walled beams where the modelling of cross-sectional deformation is essential to obtain accurate stiffness results.Secondly, the GMS is used to study the clamping of a flexure.Next, examples of a folded flexure and a cartwheel joint show how the GMS can be used to accurately model the connection between multiple flexures.
Finally, the GMS is used in a beam-problem in order to perform a convergence analysis.

Connection of I-profiles
Torsional warping is especially important in thin-walled beams with an open profile [35].
Figure 7 shows a horizontal I-profile that is connected to a vertical I-profile.Using only beam elements, the warping of the horizontal profile at this connection can only be considered either completely constrained or completely free.The use of a GMS allows a more precise analysis.The profiles are made of steel (Young's modulus 200 GPa, Poisson ratio 0.3) and a torsional moment of 1 000 Nm is applied to the horizontal profile.The vertical profile is exactly constrained: the global x-and y-displacement of both ends are constrained, and at the upper node the axial displacement and torsional rotation are constrained.The construction is modelled in six different ways, listed in Table 2.All the used beam elements have a length of 50 mm.Using more beam elements does not significantly contribute to the accuracy.All finite element meshes are all generated using ANSYS with quadratic tetrahedrons with a size of 5 mm.The 15 flexible modes of the GMS are defined using the free-free-option described in Sect.3.4.2 of [1], and the interface constraints are imposed using IMPCs.A finite element model of the whole structure is used as a reference.
Table 2 shows the resulting rotation angle at the position of the applied moment and the maximum stress, based on which the following observations are made: • Modelling the structure with a single GMS (case 'b') gives almost the same results as the linear finite element model (case 'a'), because the GMS is based on the same finite element model.• The stiffness cannot be computed accurately by modelling the horizontal I-profile using only beam elements, where the warping at the connection is either completely constrained or completely released (case 'e' and 'f' respectively).• Using the GMS in combination with warping beam elements (case 'd') gives more than 96% accuracy in stress and stiffness, which is over three times more accurate than the result obtained with the GMS with traditional beam elements (case 'c').This indicates that connecting the torsional warping of beam elements to the deformation of the interfaces of the GMS can increase the accuracy significantly.• The stiffness of the GMS with 10 warping beams (case 'd') is slightly too high.This is caused by the fact that the part of the horizontal I-profile modelled with the GMS is relatively small compared to the size of the cross-section.Some deformation is present close to this connection that is not modelled with the warping beam element.

The clamping of a flexure
To investigate the warping behaviour of the interface between the GMS and a warping beam element the clamping of a flexure is modelled, see Fig. 8.The flexure is made of steel (Young's modulus 200 GPa, Poisson ratio 0.3), has a length of 12 mm, a width of 10 mm and a thickness of 0.5 mm.A torsional moment of 0.25 Nm is applied at the unclamped side.The first 10 mm of the unclamped side of the flexure is modelled with five beam elements.The clamped side is modelled in six different ways, see Fig. 8:   The GMSs and finite element models are modelled by quadratic tetrahedrons with a mesh-size of 0.125 mm for the part of the flexure and 0.5 mm for the part of the block.The left side of the GMS is always imposed with an EMPC, the side that is connected to a beam element is imposed using an EMPC or an IMPC. Figure 9 shows the resulting rotation and maximum stress based on which the following observations can be made: • Constraining the torsional warping has a significant influence on the torsional stiffness, as the rotation in case 'a' is about 20% higher than the other cases.• There is only a small difference of about 2% between the rotation of the cases without a block (case 'b', 'c' and 'g') and the cases with the block (case 'd', 'e' and 'h').This indicates that the warping at such an interface can be considered to be fully constrained.The fillet (case 'f' and 'i') does add about 8% stiffness.• Comparing case 'b' with 'c' and case 'd' with 'e' indicates that modelling part of the flexure with a GMS gives about the same stiffness as the beam element.• Comparing the rotation of case 'c' with 'g', case 'e' with 'g' and case 'f' with 'i' shows that the rotation of the finite element model is consistently about 0.1 degrees higher than the rotation obtained using the GMS and beams.This is mainly because the finite element model is slightly more compliant at the side where the moment is applied.After compensating for this, the resulting stiffness obtained using the superelement with warping beams is more than 97% accurate.• The distribution of the stress around the clamp is very similar in all cases, except for case 'a'.The stress distribution of the finite element models at the side where the moment is applied does not correspond to the stress of the beam-elements, because the multipoint constraint used in the finite element model causes some deformation which does not correspond to the beam model.• The stress in the sharp corners of the block without fillet (case 'd', 'e' and 'h') becomes theoretically infinite high, therefore the maximum stress obtained by the GMS differs significantly from the stress of the finite element model.For the other two shapes, the maximum stress of the GMS is closer to that of the finite element model: the errors with an EMPC are lower than 3%, the errors with an IMPC are lower than 13%.

Folded flexure
A folded flexure (see Fig. 10) has a high stiffness in the vertical translational direction and is compliant in the other five directions.Contrary to what its name suggests, a folded flexure is not necessarily manufactured by folding a flat strip.The folded flexure in Fig. 10 is made of steel (Young's modulus 200 GPa, Poisson ratio 0.3), the fold is modelled using a GMS and both sides by warping beam elements.The GMS is modelled with quadratic tetrahedrons with a mesh-size of 0.18 mm, resulting in models that consist of about 23 000-162 000 elements.The interface constraints are imposed with EMPCs.The fold is modelled with a radius or a thickening, as shown in Fig. 11.The radius is also approximated using six very short beam elements.The thickening could not be modelled using beam elements.A finite element model of the full folded flexure is used as a reference.
Figure 11 shows the compliance in two directions.The compliance in the support direction of the flexure (the z-direction) increases significantly with the radius.This increase in compliance can be obtained with about 80% accuracy by the GMS, and also using six warping beams to approximate the fold.An applied moment around the x-axis causes torsion of the clamped side of the folded flexure.This is because the torsion of the two sides of the folded flexure interact through the warping around the fold.This effect is affected by the size of the thickening, mainly because the thickening increases the resistance against warping.The effect is modelled with 95% accuracy using the GMS.

Cartwheel joint
Figure 12 shows a cartwheel joint.The flexure is made of steel (Young's modulus 200 GPa, Poisson ratio 0.3).The warping of each of the four flexures interact with each other at the connection.This part is modelled using a GMS, built from quadratic tetrahedrons with a mesh-size of 0.18 mm, resulting in a model with about 56 000 elements.The interface constraints of the GMS are imposed with EMPCs.The result with only beam elements is obtained by assuming completely constrained warping of the four flexures in the connection.A finite element model of the full flexure is used as a reference.
Figure 13 shows the rotational compliance around the global z-axis as a function of the rotation.After some rotation, this compliance depends on the torsional stiffness of the flexures and therefore also depends on the modelling of the warping around the connection.After a 20-degree rotation around the y-direction, the resulting error for the GMS is about 50% smaller than the error obtained by only using warping beam elements.

Beam with a tip-load
Figure 14 shows a beam-structure that is modelled using multiple GSMs.The beam has length of 200 mm, a width of 10 mm and a height of 5 mm.It is made of steel (Young's modulus 200 GPa, Poisson ratio 0.3).The beam is fixed at the base.At the tip a force is Fig. 13 Rotational compliance around the z-axis of the cartwheel Fig. 14 Beam modelled using multiple GMSs. Figure also shows the cross-section with local coordinate-system.For visualization, the mesh is coarser than the mesh used to obtain the results applied of 1 000 N in the z-direciton and 500 N in the y-direction.The GMSs are modelled by quadratic tetrahedrons with a mesh-size of 1.5 mm.Each GMS has two interfaces which are either modelled without interface deformation fields or with deformation fields up to the second order in the local coordinates of the face, as shown in Table 3.The deformation of the cross section at the both ends of the beam is fixed.The element frames of the GMSs are attached to the interface node that is closest to the base.
Figure 15 shows the resulting vertical position of the tip, from which three observations can be drawn: • If no deformation fields are used, using EMPCs results in an error which increases with the number of elements.This is because all interfaces between the GMSs are completely rigid in this case.Including deformation fields gives much better results.• If no deformation fields are used, using IMPCs gives only a small error for a large number of GMSs.The beam is only slightly too compliant, mainly because the deformation at both ends of the beam cannot be constrained in this case.• If interface deformation fields are included and a sufficient number of beams is used, the GMS gives comparable results to normal beam elements.
Figure 16 shows how the error in the tip-position converges as function of the number GMSs that is used to model the beam.This is because the amount of deformation per GMS  The second plot shows the time to obtain the reduced order model for one GMS.This time is almost inversely proportional to the number of GMSs as the number of elements per GMS is also inversely proportional to the number of GMSs, which is about 25 000 elements if one GMS is used and about 2 500 elements if 10 GMSs are used.
The computation time of the multibody simulation increases proportional to the number or GMSs.Including deformation fields significantly increases the computation time as this results in extra degrees of freedom.The current example shows a much higher computation time for the model order reduction than for the multibody-simulation.However, note that this can be different in other cases as these computation times are for example highly dependent on the number of elements per GMS and the number of bodies in the multibodysystem, respectively.

Conclusions
Superelements compute the small deformation of arbitrarily shaped components efficiently, using the results of model order reduction techniques.A multipoint constraint has been derived that can be used to obtain a reduced order model with deformable interfaces.The formulation gives more consistent results than other multipoint constraint formulations for rigid interfaces in literature.The multipoint constraint is imposed using a combination of the Lagrange multiplier method and master-slave elimination to allow for efficient model order reduction.
The resulting reduced order models with deformable interfaces are used in the GMS, a superelement in the generalized strain formulation.The interface deformation can be defined consistent to the deformation of the cross-section of higher-order beam elements.In this way, structures with slender parts can be modelled efficiently and accurately using beam elements in combination with GMSs.
This paper combined the GMS with beam elements in which the axial warping due to torsion is included.The GMS gives accurate results if the deformation in the GMS is small and sufficient interface deformation fields are used.The stiffness of a frame consisting of two I-profiles was modelled 96% accurately, and the maximum stress over 98% accurately.The GMS was also applied to model the critical parts of a single clamped flexure, a folded flexure and a cartwheel joint.In these models the errors in the stiffness were below 6%, typically at least twice as accurate as the stiffness modelled using only beam elements.

Fig. 1
Fig.1Reducing a finite element model of an I-profile connection.Frame i is the global reference frame of the finite element model, the other frames are condensation nodes; g is an FE-node on the interface surface that is related to condensation node k.For visualization, the mesh is coarser than the mesh used to obtain the results in Sect.5.1

Fig. 2
Fig.2Overview of the multipoint constraints, for an example of a one-dimensional interface with seven nodes.The predicted displacements are obtained based on the displacement of the condensation node and the deformation fields.These predicted displacements are used to obtain an EMPC or IMPC

Fig. 4
Fig. 4 Coordinates in a GMS with global frame O, element frame j and four frames that indicate the interface nodes

Fig. 7
Fig. 7 Connection of two I-profiles modelled by a GMS and 10 beam elements.The colours in the GMS show the von Mises stress.Dimensions are given in mm.Displacements are magnified by a factor of 10. (Color figure online)

Fig. 8 von
Fig. 8 von Mises stress in the clamping of a flexure subjected to a torsional moment, modelled in nine different ways.(Color figure online)

Fig. 9
Fig. 9 Rotation and maximum von Mises stress of the clamped flexure

Fig. 10
Fig.10 Folded flexure.For visualization, the mesh is coarser than the mesh used to obtain the results, and the sides show only three warping beam elements where 20 warping beam elements per side were used to obtain the results

Fig. 15
Fig. 15 Vertical position of the tip of the beam as function of the number of GMS.Note that the x-axis has a logarithmic scale

Fig. 16 Fig. 17
Fig.16 Results of the beam-problem where the constraints are imposed using IMPCs, as a function of the number of GMSs.The first plot shows the error-convergence of the tip-position (16 GMSs are used to obtain a reference).The second plot shows the required computation time to obtain the reduced order models that are used by the GMS.The grey line in this plot shows the time to assemble the full stiffness matrix[K i FEM ].The other two lines show the time to compute the reduced stiffness matrix,[K i  All ].The third plot shows the required computation time of the nonlinear multibody simulation of the beam-model that consist of multiple GMSs Fig.17 Remaining error after each iteration step during the solution of the beam-problem that is modelled with 10 GMSs.The total force is applied in four steps

Table 1
Rotations of condensation node k for the displaced interfaces of Fig.3

Table 3
The 12 interface deformation fields that are used in the beam-problem, ω = ω x ω y ω z