Kinematically started efficient position analysis of deformed compliant mechanisms utilizing data of standard joints

Topology optimization of a flexure-based mechanism requires the properties of the mechanism in several deformed configurations. This paper presents a fast and accurate method to compute these configurations. It is generally applicable on mechanisms with complex standard flexure joints. First kinematic equations of the mechanism are derived by allowing the mechanism to move only in the directions for which it is designed. Secondly the configurations of the joints are approximated based on the rotations of the elements by which the joints are modeled. These orientations are obtained by a parameterization based on a priori knowledge of standard flexure joints. Finally, the resulting approximation is used as initial guess to obtain the configuration accurately, after which relevant properties like stiffness can be derived. For a manipulator with three complex joints the computation time was reduced up to a factor of 65 compared to a conventional method. When for optimization purposes an approximation is acceptable, the computation time can be reduced by a factor of 600, using a linear description of the deformation that remains in the first part of the method. © 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


a b s t r a c t
Topology optimization of a flexure-based mechanism requires the properties of the mechanism in several deformed configurations. This paper presents a fast and accurate method to compute these configurations. It is generally applicable on mechanisms with complex standard flexure joints. First kinematic equations of the mechanism are derived by allowing the mechanism to move only in the directions for which it is designed. Secondly the configurations of the joints are approximated based on the rotations of the elements by which the joints are modeled. These orientations are obtained by a parameterization based on a priori knowledge of standard flexure joints. Finally, the resulting approximation is used as initial guess to obtain the configuration accurately, after which relevant properties like stiffness can be derived. For a manipulator with three complex joints the computation time was reduced up to a factor of 65 compared to a conventional method. When for optimization purposes an approximation is acceptable, the computation time can be reduced by a factor of 600, using a linear description of the deformation that remains in the first part of the method.

Introduction
Flexure-based mechanisms have become increasingly complex, especially for the cases where a large motion should be realized. These mechanisms typically consist of multiple flexure joints, and each flexure joint can consist of more than thirty leaf springs, see Fig. 1 for an example [1] . Moreover when large deformation is considered in the leaf springs, each leaf spring should be modeled by using multiple elements. In the end many elements are required to obtain an appropriate model of flexure based mechanisms.
Design optimization is common in flexure based design [1][2][3][4][5] as the performance of these mechanisms highly depends on the configuration of the mechanism and the geometry of the flexure joints. These optimizations require an analysis of the mechanical properties of the mechanism, e.g. the support stiffness, the maximum stress, the buckling load and the eigen frequencies of the mechanism. These properties should be evaluated over the full range of motion. Therefore a number of relevant deformed configurations is chosen for which these properties are evaluated in each optimization iteration step. Computing these relevant configurations currently requires the solution of nonlinear equations by iterative methods, the high eigen frequencies in order to increase the allowable time step in dynamic simulations. In this paper a similar technique will be used to improve static computations, i.e. a technique that first solves a large part of the motion and estimates the remaining part by a linear approximation.
All in all, there is no proper reduction technique to obtain a model for flexure based mechanisms that is accurate for large deformation and that does not require a priori obtained knowledge of the complete mechanism.

Approach
This paper presents a combination of two methods. In the first place a method to kinematically obtain the configuration of a flexure-based mechanism by only allowing the motion for which the mechanism is designed. This configuration is used as a starting point to obtain the real deformed configuration based on static equilibrium. In this way the need for a cumbersome iterative procedure to obtain the deformed configuration is avoided. This method will be referred to as the Kinematically Started Deformation method (KSD-method).
The motion for which the mechanism is designed is henceforth referred to as intended motion/deformation . In other papers this is referred to as degrees of freedom, despite the fact that the mechanisms can also move slightly in the other directions. This other motion (i.e. the motion in the support-directions) will be referred to as unintended motion/deformation . In order to distinguish intended and unintended deformation, it appears to be useful to consider flexure based mechanisms as an assembly of flexure joints and very stiff links. The distinction between intended and unintended deformation will also be made in each flexure joint.
In the KSD-method, each flexure joint is modeled by a small finite element model. The method requires the internal configurations of the deformed flexure joints in the mechanism to be obtained. This requires a time-consuming iterative procedure in general. However, this step can be applied much more efficiently, using the fact that flexure joints are typically designed according to standard assemblies, e.g. the cross flexure [ 23 , 24 ] the cart wheel [ 24 , 25 ], the butterfly hinge [ 25 , 26 ], the infinity hinge [2] or the 3x-Infinity hinge in Fig. 1 [1] . Other standard assemblies can be found in [27] . Joints based on these standard assemblies will be referred to as standard flexure joints .
Therefore it is considered to a feasible solution to setup a database that contains information about these standard flexure joints in order to speed up the static calculation of all mechanisms containing these standard flexure joints. In other words, we propose a data-driven technique. In the literature review above, it was stated that using a data-driven reduction technique requires much a priory computation time. However, the proposed data-driven technique only requires data of default parts, which means that the databases does not have to be updated before each design optimization. Another disadvantage of data-driven techniques mentioned above is that the accuracy is sensitive to variations in the dimensions of the model. This disadvantage is resolved by the second method in this paper, which is based on a priory obtained parameters that are valid throughout the full range of commonly encountered dimensions of the flexure joints.
The second method presented in this paper obtains the largely deformed configurations of standard flexure joints, based on a limited number of parameters that are not sensitive to changes in the dimensions of that joint. It appears that the orientations of the elements are a good choice for these parameters as the orientations are independent of most of the dimensions of the flexure joint and provide enough information to obtain a configuration that is close to the deformed configuration. A body that is described by this technique will be referred to as an Element Orientation based Body (EOB).
Once the internal configuration in the EOB is obtained, its stiffness matrix can be obtained and this stiffness matrix can be reduced by model order reduction. In this paper, the term element is reserved to refer to a modeling part of which the deformation can be described by analytic relations, e.g. beam elements or plate elements. A body is a modeling part that is based on a reduced finite element model. So each flexure joint is a body that is described by multiple elements.
The main idea of the KSD-method is to prevent the need for a large number of iterations to obtain the deformed configuration of a mechanism , where the main idea of the EOB is to prevent the need for an iterative approach to find the internal configuration of the joints . With this combination the required simulation time is reduced significantly.
The proposed method shares various concepts with the methods in the literature review. The KSD-method first solves for a main part of the deformation, i.e. the intended deformation. This is similar to the PRBMs, the semi-analytic-deflectionmethod [12] and the GMP-technique. In contrast to PRBMs, flexure joints are not simply modeled as ideal hinges, but parasitic motion is taken into account, and the KSD-method does offer the opportunity to approximate the unintended deformation. In contrast to GMP, the KSD-method separates the large and small contribution to the deformation based on physics, i.e. in intended and unintended motion. Another difference is that GMP is used in dynamics and the KSD-method solves static equations. In contrast to the semi-analytic-deflection-method, the KSD-method can be applied to mechanisms with flexure joints that are arbitrarily complex.
The EOB computes the internal configurations of joints based on a priory simulations, similar to the data-driven reduction techniques. The main novelty of the EOB is the type of data that is used, i.e. the element orientations. Therefore the data is insensitive to changes of most of the dimensions. Section 2 explains the KSD-method and Section 3 gives details about the EOB. The KSD-method with EOBs is validated in Section 4 by the analysis of a single cross flexure, a fourbar-mechanism with four cross flexure joints and a manipulator with three 3x-infinity joints. The paper ends with the most important conclusions.
Please cite this article as: K. Dwarshuis

Static deformation starting from a kinematic approximation
The KSD-method computes the static deformed configuration of a mechanism in five steps. Details about the steps are given in the following subsections. In the first two steps the configuration is estimated purely based on kinematic relations. In these steps, only the intended deformation is described and a database with data of the standard flexure joints is used. Based on these two steps, stiffness matrices of the joints can be obtained. In steps 3 and 4 these stiffness matrices are used to update the estimation of the configuration by static equilibrium. In these steps, the deformation in the unintended directions is described linearly.
Step 5 computes the configuration based on static equilibrium, taking the configuration after step 4 as an initial guess. Because this initial guess is already quite good, only a few iteration steps are required in this fifth step. Two conditions should be satisfied in order to perform the KSD-method. In the first place, the system should be kinematically determinate, as explained in Section 2.2 . Secondly, data of all the joints in the mechanism should be present in a database.
The steps are schematically visualized in Fig. 2 for a two-dimensional fourbar mechanism with four cross flexures. Each cross flexure is modeled by a finite element model with three flexible beam elements per leaf spring and four rigid elements to link the flexible elements. The configuration is described by the positions and rotations of the nodes. The nodes at which the joint is connected to other parts are referred to as interface nodes and the nodes inside the joints are internal nodes . The displacements related to these nodes are interface displacements and internal displacements respectively. The term displacement refers to the combination of rotations and translational displacements in this paper.

Database with information of standard flexure joints
This subsection explains what kind of data of a joint is stored in a database and how this information can be obtained. The joint is modeled by a finite element model. The position and orientation of one of the interface nodes of this model is fixed. The displacement of the other interface node is prescribed in the intended directions, on a finite number of values over the range of motion. The resulting configurations of the flexure joint for each prescribed value are obtained based on static equilibrium, where the forces and moments in the unintended directions are zero. These conditions result in a unique configuration of the flexure joint for each prescribed value of the intended direction. Fig. 3 shows the finite element model of the cross flexure that is used as an example in this section. It is fixed at interface node A. A cross flexure is meant to allow rotation, so the rotation of interface node B is the intended direction. This intended direction can for example be prescribed on 10 intervals in the range between −30 °and 30 °, assuming the reaction forces in the x-direction and y-direction on interface node B to be zero.
The resulting configurations are used to derive two kinds of functions that are stored in the database for use in step 1 and 2 of the KSD-method. The functions for step 1 are kinematic relations between both interface nodes of the joint. These are formulated as constraints for the unintended deformations. In case of the cross flexure, these functions can be chosen to be the horizontal and vertical positions of interface node B with respect to interface node A as functions of the intended  deformation θ , so the functions x ( θ ) and y ( θ ). These functions can be derived based on the simulation data, for example by a least square polynomial fit. The second kind of functions that are stored in the database are the rotations of the elements as functions of the intended deformation. These functions are also obtained based on the same simulation data. Section 3.1 gives more details about the way to describe element orientations.
Note that static equilibrium is used to obtain the necessary data for the database, but that the resulting functions in the database are purely kinematic. One of the advantages of using kinematic relations is that they do not depend significantly on most of the dimensions of the joint, e.g. the kinematic relations of the cross flexure joint do not depend on the width and thickness of the flexures.
However, some design parameters may affect the functions in the database. For the cross flexure, the overall length of the flexure (i.e. the distance between both interface nodes) will affect the constraint functions x ( θ ) and y ( θ ). However, this can be easily avoided by storing dimensionless functions in the database: x (θ ) = x/L and ȳ (θ ) = y/L , where L is the length of the flexure. But also the angle α between both flexures (as indicated in Fig. 3 ) influences the results. This means that the functions in the database explicitly depend on α: x ( θ , α) , ȳ ( θ , α) So, although we wanted to obtain data that is insensitive to the dimensions of the joint, still some dimensions may have to be considered explicitly.

Step 1 -estimate the interface displacements by kinematics
In step 1 the positions and orientations of the interface nodes of the mechanism are approximated by a kinematic approach based on two conditions. In the first place the displacement of the end-effector in the actuated directions, q end , is prescribed. The second condition is that the joints only allow motion in the intended direction. This is done by using the constraint equations for the standard flexures in the database. In this way the flexure joints do not have to be modeled as simple ideal pivot joints, but parasitic motion is also taken into account.
These two conditions are sufficient to obtain the configuration based on kinematic equations. However, this can only be solved uniquely if the mechanism is kinematically determinate (in other literature also referred to as a mechanism that is not underconstrained). This means that the KSD-method does not work for mechanisms that are kinematically indeterminate. On the other hand, if the mechanism is statically indeterminate (overconstrained) the number of unknown displacements is less than the number of equations. This means that some of the constraint equations are redundant. By removing these redundant equations the interface displacements can still be found. This system of nonlinear kinematic equations is solved by a Newton-Raphson iteration. Although this is an iterative process it can be computed very fast with respect to the computation of the solution of the full finite element model of the mechanism.

Step 2 -estimate the internal displacements by kinematics and obtain stiffness matrices of the joints
In step 2 the internal displacements are obtained from the interface displacements of step 1. These displacements will then be used to obtain the stiffness matrices of the joints. This can be done by using the Element Orientation based Body (EOB) description which is explained in Section 3 . The result is a relation for each joint like: where F nodes and q nodes are the forces on the nodes and the displacements of the nodes of the joint, respectively. K nodes is the stiffness matrix of the joint. The term F nodes is a virtual force that appears because this stiffness relation is linearized around a deformed configuration and not in the undeformed configuration.
Note that the internal configuration of a joint can also be obtained without using the EOB. The most straightforward method to find the internal configuration is by solving the equilibrium equation of the finite element model of the joint, using the displacements of the interface nodes of step 1 as boundary conditions. Once the internal configuration is found, also the stiffness relation of Eq. (2.2) can be found. However, solving the equilibrium equation will require an iterative process that will take a lot of computation time with respect to the solution of the EOB. Therefore in this paper EOBs are used to obtain the stiffness matrices of the joints.

Step 3 -update the interface displacements based on static equilibrium
Step 3 updates the interface displacements by using the stiffness matrices from step 2. From this step on the unintended deformation is not constrained anymore and force-displacement relations are used, in contrast to step 1 and 2 where kinematic relations were used. This means that also load on constrained directions will have influence on the result of this step. Also the compliance of the stiff connecting links (e.g. the three links in the fourbar mechanism) can be modeled in this step.
In step 3A the stiffness matrices of the joints that specifies the stiffness between the interface nodes are obtained by using the boundary modes of the Craig-Bampton method [28] . The extra force term F , has to be taken into account during this reduction. Eq. (2.2) is split in interface coordinates of the joint (normally called boundary coordinates in the Craig-Bampton method, indicated by b ) and internal coordinates of the joint (indicated by i ): By assuming no applied forces on the internal nodes, so F i = 0 , this equation can be rewritten to: where: K T is the Craig-Bampton reduced stiffness matrix of the joint. F T is a virtual force that is applied on the interface nodes of the body.
In step 3B the interface displacements are computed by using these Craig-Bampton reduced stiffness matrices. The interface displacements of step 1 are used as the starting point for this step. Because these displacements are close to the displacements that are obtained in this step, only a few iteration steps are required to solve this static problem.

Step 4 -estimate the internal displacements based on static equilibrium
In step 4, for each joint the internal displacements are updated by using the interface displacements of step 3 and the joint stiffness matrices of the joints of step 2. An equation for the internal displacements can be derived from Eq. (2.3) : Step 5 -update all displacements based on static equilibrium In step 5, the interface and internal displacements are updated simultaneously by solving for static equilibrium by an iterative approach. The interface displacements of step 3 and the internal displacements of step 4 are used as a starting point. Because the starting configuration is close to the deformed configuration, only a few iteration steps are required. The system of equations that is solved in this step is exactly the same as the system of equations that is solved in the conventional approach, i.e. obtaining the deformed configuration by an iterative process starting from the undeformed configuration. This means that step 5 will give the same result as the conventional method where it is generally faster as it will involve fewer iteration steps.

Reduced KSD-method and full KSD-method
In the kinematic iterative process in step 1 and in the static equilibrium iterative process in step 3 a mechanism is considered in which each joint is considered as a single body and where only the interface displacements were obtained. This mechanism will be referred to as the reduced mechanism . In contrast, in step 5 the full mechanism is considered. In the full mechanism each element in the finite element models of the joints is considered as a separate modeling part and the interface displacements and internal displacements are computed simultaneously. The result after step 4 is a good approximation of the displacements. Using this as final result will be referred to as KSD-reduced . Using all 5 steps will be referred to as KSD-full .
Step 3 may seem unnecessary for KSD-full, as step 3 updates the displacements in unintended directions which are generally small. However, using step 3 to update the displacements of the interface nodes gives the advantage of performing some more iteration steps with the reduced mechanism instead of the full mechanism, which is much faster. Moreover, the inaccuracy of the obtained interface displacements in step 1 can be significant, such that a static equilibrium iteration on the full mechanism in step 5 can become instable.

Element orientation based approach to obtain the internal configuration in standard joints
This section shows the approach to obtain the internal configuration and the stiffness relations of an EOB, which is used to apply step 2 of the KSD-method. The position and orientation of the interface nodes should be known on beforehand. So the EOB is valuable in combination with the KSD-method, as the interface displacements of the joints are computed in step 1 of the KSD-method. Based on the interface positions and rotations and a database, the element orientations are obtained ( Section 3.1 ). Using these orientations a configuration is derived that is close to the deformed configuration, referred to as the near configuration ( Section 3.2 ). In this near configuration the stiffness matrix of the body is derived ( Section 3.3 ).

Obtain the element orientations
The orientations of the elements are described by functions of the intended deformation that are stored in the database.
These functions may depend on some of the design parameters. In case of the cross flexure, the angle α (see Fig. 3 ) appears to have a significant influence on the rotation of the elements, as noted in Section 2.1 . The intended deformation is already known after step 1 of the KSD-method, the element orientations can therefore be obtained by simply evaluating the functions that are stored in the database: ˆ β k = f k ( Intend ed d e f ormation, Design P aram ) , where ˆ β k are the parameters that are used to describe the rotation of element k . The meaning of the hat on β is explained below. If the intended deformation is only in one plane, for example the deformation of a 2D or 3D cross flexure, the rotation of each element can be described by the rotation about one axis. If the elements rotate about multiple axis other techniques are required to describe the orientation. In that case ˆ β k can be for example a vector with Euler parameters or Euler angles.
Instead of using a database with functions to obtain the element orientations, we could also define functions that define other properties of the joint by which the internal configuration and stiffness matrix of the joint could be obtained. One other option is to store some relevant terms of the stiffness matrix. However, these terms will depend on the material properties of the joint and are very sensitive to the thickness and width of all the individual flexures. That means that the database should provide different functions for all the relevant designs of the flexure joint. The configuration of the joint is less sensitive to these design-parameter variations, therefore storing kinematic functions is preferred.
Another option is to store functions for the displacements (rotational and translational) of all the internal nodes in the body. This would immediately describe the internal configuration and provide enough information to obtain the stiffness matrix. However, the translational displacements are more sensitive to variations in the dimensions of the joint than the element rotations. All in all, the element orientations are probably the most dimension independent parameters that provide enough information to obtain the stiffness matrices of the elements and the internal configuration of the body as the remaining part of this section will show. Fig. 4 (a) shows the near configuration. The displacements of the nodes in the near configuration are defined with respect to one of the interface nodes, in this case with respect to node A. The orientations of the elements in the near configuration are equal to the orientations that are described by ˆ β. The orientations of two elements that are connected to each other are not necessarily the same and therefore the orientation of the intervening node cannot be uniquely defined. Therefore each element is described by its own nodes in the near configuration, these nodes will be referred to as markers . So the near configuration of the joint in Fig. 4 (a) is described by 24 markers as it consist of 12 elements that all have two markers.

Obtain the near configuration
To obtain the translational displacements of the markers, the condition is set that the elements are undeformed in the near configuration.
where ˆ u k,q and ˆ u k,p are the translational displacement of marker q and p respectively on element k . The hat indicates that the related variable describes a displacement from the undeformed configuration to the near configuration, these displacements will be referred to as rigid displacements . The vector d k, pq is the distance from marker p to marker q in the undeformed configuration. The matrix ˆ R k is the rotation matrix that defines the rotational part of the rigid displacement is only a function of ˆ β k . Note that in case of one or more closed loops of elements, there will be points at which the markers of two coupled elements are not exactly at the same location because they are related to node A by a different chain of elements. This is also the case in the cross flexure, therefore the locations of both markers at node B are not exactly the same. For clarity this effect is exaggerated in Fig. 4 (a).

Obtain stiffness matrix in the near configuration
This section describes how the stiffness matrix can be obtained in the near configuration by which a better approximation of the deformed configuration can be obtained. Because the orientation of each element in the near configuration should be close to the orientation of the deformed configuration, the displacements between both configurations can be described linear. Note that there can be deviation in the locations of both configurations, but large translations do not make the force-displacement relation nonlinear. The force-displacement relation for an element k is where F k is the vector with the forces and force moments on the markers of element k expressed in absolute coordinates. So this vector consist of six numbers for each marker of a three dimensional element k . K k is the global stiffness matrix of element k which can be obtained by rotating the local stiffness matrix, using the rotation matrix ˆ R k of the element. The vector q k contains the translational and rotational displacements on the markers of element k . The bar indicates that these are displacements from the near configuration to the deformed configuration, which will be referred to as the flexible displacements . Note that there is no force required for the rigid displacements, this means that the force in Eq. (3.3) is the total force that is required to have element k in its deformed configuration. The relations for all elements can be combined in one equation:  5) where N is the number of elements by which the joint is modeled. The total displacement is the sum of the rigid and the flexible displacement. This is shown for the two markers of one element in Fig. 5 . So, the flexible displacement is the total displacement minus the rigid displacement: Note that, in case of large rotations about multiple axis, the rotational parts in these displacement vectors cannot be added directly due to the non-vectorial nature of rotations. Therefore the rotations in q all and q all should be defined in such way that their difference equals the finite rotations in q all . In this paper only large planar deformation is considered, such that the rotation in q all and q all can be described by the rotation about the axis perpendicular to this deformation plane. Substituting Eq. (3.6) in Eq. (3.4) gives: in which the term: is constant. This term compensates for the fact that the rigid motions in the first term on the right hand side of the equation are computed linear. This term can be computed directly as the rigid displacements q all are the displacements from the undeformed to the near configuration that was obtained in Section 3.2 .
In order to obtain the deformed configuration, the constraints between the elements should be applied. These constraints are applied by defining a Boolean matrix L such that: where q nodes is the matrix with the total displacements of the nodes in the joint. The forces on the nodes can then be related like: This means that the relation between the forces and the displacements on the nodes can be written like: F nodes = K nodes q nodes + F nodes , K nodes ≡ L T K all L. (3.11) This is the static equilibrium equation, which is valid around the near configuration. The equation is identical to eq. (2.2) that should be obtained in step 2 of the KSD-method.

Numerical validation
The EOB is validated by an analysis of a 3D cross flexure ( Section 4.2 ) and the KSD-method with EOBs is validated by a fourbar mechanism with four cross flexure joints ( Section 4.3 ). Section 4.1 describes the cross flexure model that is used in these sections. The KSD-method is further analyzed by a spatial manipulator with more complex joints in Section 4.4 .  The used stiffness matrices of the flexible 3D beam elements are derived in Appendix A . The performance of the KSDmethod is compared to that of a conventional method that is described in Appendix B . Both, the KSD-method and the conventional method are implemented in matlab and both methods solve the same finite element models of the mechanisms.

Cross flexure model
The finite element model of the cross flexure that was used in the simulations is shown in Fig. 6 . The relevant properties of the flexure are given in Table 1 . Both connector blocks are modeled by rigid elements and the interface nodes are exactly in the center of these connector blocks. The leaf springs were modeled with 1 to 5 serial connected equally sized elements. The intended rotation is described by the angle between node A and B around the z -axis. It is denoted by θ and referred to as the (intended) rotation. For each of the models with 1 to 5 elements per leaf spring, simulations were run where the intended rotation was varied from 0 to 75 °in steps of 5 °and the displacements of all nodes were obtained. These displacements are used as simulation data to obtain the two required types of functions.
In order to perform step 1 of the KSD-method, five constraints for the unintended deformation are required. Four of these constraints can be defined without using data of a priory simulations (the chosen constraints are different from the constraints for the 2-dimensional cross flexure, introduced in Section 2.1 ): • The displacement in the z -direction of node B with respect to node A should be zero. • The rotation around the x -axis of node B with resprect to node A should be zero • The rotation around the y -axis of node B with resprect to node A should be zero • The rotation around the z -axis of node A and B with respect to the line through both nodes should be equal in size and opposite, as shown in Fig. 6 .
The fifth constraint is the length of the cross flexure (the distance between both interface nodes) as a function of the intended deformation. This length was parameterized as a function of the intended rotation by a sixth order polynomial least squares fit to the simulation data. The odd terms in this polynomial are zero. The obtained relation is found to be almost independent of the cross flexure dimensions, except for the angles at which the leaf springs are positioned with respect to the local x -axis, which is denoted by α in Fig. 6 where D is the distance between both interface nodes, and D 0 is this distance in undeformed configuration.
The rotations of the elements around the z -axis were fitted by third order polynomials as a function of θ . So for the rotation of each element k a relation is obtained like: The global stiffness matrices of the beam elements as described in Appendix A depend on the orientations of the nodes, therefore also the orientation of the nodes are parameterized by a third order polynomial, similar to Eq. (4.2) .

Results for single cross flexure
The accuracy of the results of the linearized equilibrium in the EOB, Eq. (3.11) was analyzed, which is representative for the accuracy of KSD-reduced. Fig. 7 shows the force moment that has to be applied on the cross flexure as a function of the intended rotation. The force moment is normalized by dividing it by the linear approximation [23] , given by: where M is the applied moment on node B around the z -axis and I is sum of the second moment of area for all the three leaf springs:

(4.4)
An analytic reference for the single cross flexure, which is based on the assumptions of infinite axial stiffness of the flexures and no shear, was adapted from [23] . Fig. 7 (a) shows that in the default case, the results of the EOB are exactly similar to the results of the finite element model of the joint: The line of the FEM result in this figure is exactly behind the line of the result of the EOB. Also the stiffness in the unintended directions of the EOB was found to match the result of the finite element method perfectly.
The applicability of the EOB is found to be bounded by two kinds of limitations. The first kind of limitation is that the difference in stiffness between the inner flexure and the two outer flexures may not vary too much. Fig. 7 (b) shows the effect of wider outer leaf springs. In the default case, w 2 = 0 . 02 m , the line of the reference case is hidden behind the line of the EOB. If the width of the outer flexures is four times as large, the standardized moment is a little lower for both cases, but this effect for the EOB is smaller than for the reference case. So there is some variation between the EOB and the reference, which is caused by the fact that the parameterizations of the rotations ˆ β are not accurate for a large variation in the ratio between these stiffness of the inner and the outer flexures. The result indicates that it should be carefully considered for which ranges in the dimensions the parametrization holds.
The second kind of limitation is on the amount of unintended deformation. The results become inaccurate if this deformation becomes so large that the assumption of linear unintended deformation is not valid anymore. Fig. 8   intended direction for large intended deformation. Therefore the applied force is scaled by dividing it by the stiffness in this direction: where F z is the applied force on node B in the z -direction. The z -displacement is normalized by dividing it by the length of the leaf springs, L . The force moment that is required for the intended deformation is scaled by dividing it by the linear stiffness: Note that this scaled force moment is different from the normalized force moment that was introduced in Eq. (4.3) . The results in Fig. 8 indicate that the deformation in the z -direction is approximated quite accurately with the EOB, but that it shows almost no effect of the disturbance force on the applied force moment, where this effect can be quite significant. Fig. 9 shows the fourbar mechanism. The cross flexures are modeled with the properties given in Table 1  flexures are assumed to be infinitely stiff. The desired motion of the fourbar is described by the displacement of the center of the upper link in the x -direction. The mechanism is actuated by a force on this point in the x -direction. Fig. 10 shows the normalized applied force as a function of the normalized displacement for the full mechanism and for the reduced mechanism. The displacement, d , is scaled by dividing it by the height H of the mechanism, the force is normalized by dividing it by the linear approximation, given by:

Fourbar-mechanism
where F is the force that is applied on the center of the upper link. The result shows that KSD-reduced gives accurate results.
The computational efficiency of the KSD-method is examined by computing the deformed configuration of the fourbar mechanism where the endpoint is displaced by d = 0 . 3 m . The computation stops if the norm of the vector with resulting forces and force moments on the nodes is smaller than 5 · 10 −7 , in which the forces are expressed in N and the moments in Nm. Fig. 11 (a) shows the number of iteration steps that are required for the conventional approach, KSD-reduced and KSD-full. In each iteration step the independent coordinates are updated one. Inside each iteration step, another iterative algorithm updates a set of dependent coordinates. For the full approach, the number of iterations is split in the number of iterations on the reduced mechanism for step 3 of the KSD-method (which is of course identical to the number of iterations for KSDreduced) and the number of iterations on the full mechanism for step 5. The conventional approach requires more iterations than the KSD-method, especially in case of a large number of elements per leaf spring. This is because of the higher number of independent parameters that describe the configuration, increasing the nonlinearity of the system of equations. Fig. 11 (b) shows the average time per iteration. The solution time per iteration applied on the reduced mechanism is almost independently of the number of elements per leaf spring. In the case of the fourbar each iteration applied on the reduced mechanism was solved in approximately 30 ms. The time per iteration on the full mechanism is higher. The average time per iteration on the full mechanism in step 5 of the full approach is shorter than the time per iteration of the conventional approach. It is reduced by a factor of about 2 for the cases with 4 or 5 elements per leaf spring. The reason for this is that the displacements in these iterations are much smaller, and therefore it takes shorter to update the dependent coordinates in each iteration.
In order to study the effect of forces in the unintended directions, an external moment around the z-axis is applied on the end-effector, i.e. the center of the top-bar. Fig. 12 (a) shows the resulting rotation of the upper link for the full model and for the reduced model. This error is result of the linearization of the unintended motion of the cross flexures. Even without external moment, there is a difference of about 20% in the rotation. This unintended motion is present as the cross flexures are not excited to a pure moment. Fig. 12 (b) shows the required number of iterations and the total simulation time for the conventional method and step 5 of KSD-full. It indicates that for both methods the required number of iterations increases with an increasing force in the unintended directions. The KSD-method is about 5 times faster over the full range of the external moment in unintended direction. It was observed that in simulations with a force out of plane (e.g. a force in the z-direction on the center of the top link) the computation of the KSD-method does not converge in step 5.  Overall KSD-reduced is about 20 times faster than the conventional method for 4 or 5 elements per leaf spring, respectively. KSD-full is about 5 times faster than the conventional method for 4 or 5 elements per leaf spring. The main reason for this increase in efficiency is the reduced number of required iterations.

Spatial three-link manipulator
The KSD-method is applied to the manipulator with three 3x-infinity joints and three links, as shown in Fig. 13 . The first joint can rotate around the global z -axis and the other two rotate initially around the global y -axis.  These constraints are expressed as a fourth order polynomial of the intended rotation which is a least squares approximation of the simulation data. The rotations of all the elements and nodes are also approximated by a fourth order polynomial in the intended rotation, in order to perform step 2 of the KSD-method. The manipulator has three intended degrees of freedom, therefore three parameters are required to define the position of the end-effector. Two different parameter sets are analyzed to describe the position: • Joint-case: the rotational motion of each of the three joints; • Location-case: the location of the end-point of the third link in x, y , and z -direction.
Both cases are essentially different in terms of the reaction forces. The joint-case assumes three reaction moments around the joints, the location-case assumes three linear forces on the end-point of the third link. This means that the motion in the unintended directions is different. Actually there is no unintended motion in the joint-case as all joints are actuated by a pure force moment. Therefore the KSD-method is expected to perform better in the this case.
Simulations are run in which the rotation of all three joints is simultaneously varied in the range from minus 45 to 45 °. Fig. 13 (b) shows some of the resulting configurations. In a first run the joint-case was performed. In a second run, the location-case was performed, using the locations for the end-point of the third link obtained from the joint-case. Fig. 15 shows the required number of iterations and the total simulation time. It shows that the KSD-method indeed requires more iterations for the location-case. This holds surprisingly also for the conventional method. Fig. 15 (d) shows that using KSDfull instead of the conventional method reduces the simulation time up to a factor of 65. KSD-reduced reduces the total simulation time by a factor up to 600. The joint-case has almost 100% accuracy (except from some numerical and interpolation errors). The error on the estimated reaction force of the location-case is less than 2% and the error in displacements less than 0.3%.

Conclusions
In this paper we presented a Kinematically Started Deformation method (KSD-method) to obtain static equilibrium of flexure mechanisms efficiently, by exploiting prior knowledge on the intended degrees of freedom. The relatively large deformation in the intended degrees of freedom is kinematically approximated as a first step in the iteration process.
The flexure joints in the mechanism were modeled with Element Orientation based Bodies (EOBs). The configuration of an EOB is approximated based on a parameterization of the rotations of the elements and this approximation is refined using static equilibrium. Using element orientations makes the parameterization suitable for joints with a large range of dimensions and therefore applicable in design optimizations.
The KSD-method with EOBs provides a way to use a database with information of standard joints. As a consequence, the application of the method is limited to mechanisms that consist of joints of which information exist in the database. This required information consists of relations between the interface points during intended deformation and the rotations of the elements as function of the intended deformation. Another limitation is that the mechanism should be kinematically determined in order to perform the first step of the KSD-method, but this is typical the case for flexure based mechanisms. A last limitation is that unintended deformation (deformation in supporting directions) should be low. This is not a large limitation as the purpose of a flexure based designs are related to high support stiffness. Although in this paper the method  3, 2020;23:13 ] was only applied to mechanisms that were modeled by beam-elements, the method is not fundamentally limited to this modeling approach. The efficiency of the KSD-method is compared to that of a conventional method. For a fourbar mechanism the required computation time was decreased up to a factor of 5 to obtain the accurate deformed configuration. For a manipulator with three 3x-infinity joints the computation time reduced up to a factor of 65. The computed configuration of the KSD-method is exactly the same as the configuration found by the conventional method as the KSD-method solves in the end exactly the same equations. However, the KSD-method can also be used to obtain the approximated deformed configuration in which the deformation in the unintended direction of the flexure joint is assumed to be linear. The time reduction to find this approximation with respect to the conventional method is a factor up to 20 for the fourbar-mechanism and a factor up to 600 for the manipulator. The proposed method potentially saves orders of magnitude of valuable time during the optimization of flexure mechanisms in the conceptual phase.

Acknowledgement
This work is part of the research program HTSM 2017 with project number 16210, which is partly financed by the Netherlands Organisation for Scientific Research ( NWO ).

Appendix A. Stiffness beam elements
In this appendix the stiffness matrix of a beam element is derived that is used in this paper. The used beam formulation and its derivation are relatively similar to the formulation described in [29] .

A.1. Element configuration
A.1.1. Nodal coordinates Fig. 16 shows a beam element in its initial configuration and its deformed configuration. The configuration of the beam is defined by twelve independent nodal coordinates: where r p and r q define the locations of node p and q respectively with respect to the global fixed frame. The rotation matrices that define the orientations of the nodes, R p and R q , are parameterized by λ p and λ q respectively. The initial orientation of the element is given by [ n x n y n z ] . The unit vectors at the nodes of the elements, as visualized in Fig. 16 can therefore computed by: n p y = R p λ p n y , n q y = R q λ q n y , n p z = R p λ p n z , n q z = R q λ q n z .

.2. Deformation coordinates
The deformation of the beam is described by six mode shapes. These modes are specified by deformation coordinates, ɛ , that are an explicit function of the nodal coordinates: The chosen functions for these deformation coordinates are: 2 L 0 asin n p z · n q y − asin n p y · n q z ( torsion ) ε 3 = −L 0 asin n p z · n l ε 4 = L 0 asin n q z · n l ε 5 = L 0 asin n p y · n l ε 6 = −L 0 asin n q where: The resulting mode shapes are visualized in Fig. 17 . The difference between these deformation functions and the deformation functions introduced in [29] , are the arcsine-terms in the second till the sixth function. The arcsine-terms are included to make the deformation functions linear dependent on the angle between the unit vectors. These arcsine-terms were also introduced in the 3D beam formulation in [30] and correspond to most of the 2D corotational beam formulations like [ 31 , 32 ].

A.1.3. Relation between change of nodal and deformation coordinates
The change of displacements of the nodes is given by: where δφ p and δφ q are the virtual rotations of node p and node q respectively. Note that the orientations are defined different in the nodal coordinates. However, there exists relations between these virtual rotations and the change of δλ p and δλ q . These relations depend on the choice of parametrization of the rotations and is therefore not discussed in this appendix. The relation between the change of the deformation and the virtual displacements is given by a matrix D : In order to derive an expression for D , based on the definition of the deformation coordinates in Eq. (A.4) , we note that the change of a rotation matrix can be expressed by the virtual rotations like:   Here, E is the modulus of elasticity and G is the shear modulus. A is the area of the cross section, I p is the polar moment of area of the cross section and I ȳ and I z are the second moments of area of the cross section with respect to the principal y and z -axis respectively. k x is the torsion correction factor according to Saint-Venant's theory. The shear factors are given by 16) in which k ȳ and k z are shear correction coefficients [33] . The derivation of this stiffness matrix is based on the exact solution of the equilibrium equations of the Timoshenko beam for the case that only load is applied on the nodes of the beam element. For further details the reader is referred to [29] .

A.2.2. Global stiffness matrix
The forces at the nodes, F p and F q and the force moments T p and T q are combined in one vector with nodal forces: These forces are in the same direction as the nodal displacements, q , of the element. Therefore, the relation between the nodal forces, F , and the elastic load, σ, can be derived by the principle of virtual work using Eqs. (A.7) : where the displacements, q , are assumed to be small. K is the global stiffness matrix of the beam element that was used in the simulations in this paper.

Appendix B. Conventional method
To analyze the performance, the KSD-method is compared to a conventional method which solves the full finite element model of the mechanisms. As this conventional method uses exactly the same model that was used in the KSD-method, the final results of both methods are exactly the same, but the required simulation time may be different.
The theoretical background of this conventional method is the same as the theoretical background that has been implemented in spacar [34] , as this software has often been used to model flexure based mechanisms [1][2][3][4][5] . It is particularly suited to obtain the relevant mechanism properties that are required in a design optimization, like maximum stress and eigen frequencies. The configuration is described by a set of independent coordinates and a set of dependent coordinates that are kinematically related to the independent coordinates. This theoretical background of spacar is implemented matlab , similar to the KSD-method, in order to have a fair comparison between both methods.
The number of iterations in the conventional methods is automatically determined as shown in Fig. 18 . The total displacement of the end-effector, d tot is first tried to be solved in one step. One iteration is performed, which means that the independent coordinates are updated based on the linearized equilibrium equations and the dependent coordinates are updated based on the underlying kinematic relations. After an iteration an error, ɛ , is computed based on the resulting forces on the nodes. The next step is decided based on this error, which has three options: • The error is larger than a threshold, ɛ max : in this case the iterative process will probably not converge and the step size in terms of displacement of the end-effector, d step , is reduced by a factor of two. • The error is smaller than ɛ max , but above a certain accuracy-threshold, ɛ acc : the equilibrium equations are linearized again and a new iteration is performed to increase the accuracy. • The error is smaller than ɛ min : the computation for the current step, d step , is completed, therefore either the full process is completed, i.e. the total displacement d tot has been reached, or a new step on the displacement of the end-effector is applied.
Please cite this article as: