Direct Kinetostatic Analysis of a Gripper with Curved Flexures

Micro-electro-mechanical-systems (MEMS) extensively employed planar mechanisms with elastic curved beams. However, using a curved circular beam as a flexure hinge, in most cases, needs a more sophisticated kinetostatic model than the conventional planar flexures. An elastic curved beam generally allows its outer sections to experience full plane mobility with three degrees of freedom, making complex non-linear models necessary to predict their behavior. This paper describes the direct kinetostatic analysis of a planar gripper with an elastic curved beam is described and then solved by calculating the tangent stiffness matrix in closed form. Two simplified models and different contributions to derive their tangent stiffness matrices are considered. Then, the Newton–Raphson iterative method solves the non-linear direct kinetostatic problem. The technique, which appears particularly useful for real-time applications, is finally applied to a case study consisting of a four-bar linkage gripper with elastic curved beam joints that can be used in real-time grasping operations at the microscale.


Introduction
Compliant mechanisms have been adopted in several applications for centuries, because of several well-known advantages, such as the absence of sliding friction, backlash, and significant wear, with the advantage of requiring minimum effort assembly. The first complete and systematic description of their characteristics appeared in the literature about thirty years ago, around 1994, when Midha and Howell introduced a classification [1] that has been used until today. Another significant step forward in developing the compliant mechanisms was taken a couple of years later, when the Pseudo-Rigid-Body equivalent Model (PRBM) was introduced to evaluate the elasticity of compliant mechanisms with significant deflection capabilities [2]. More materials concerning the compliant mechanisms were presented in 2001 [3].
The actual configuration of a compliant mechanism depends not only on the applied forces and torques but also on its geometric characteristics, with kinematic and mechanical coupling and non-linearity problems that especially arise in case of large deformations.
A comprehensive survey on many recent techniques for modeling the kinetostatic and dynamic behavior of flexure-based compliant mechanisms has been recently presented [26].
Kinetostatic models of complex plane compliant mechanisms have been developed in both micro and macro scale devices by using a wide variety of different linear and non-linear methods, such as, for only representative example, those based on: Loop-closure equations and the static equilibrium conditions for multi-loops compliant mechanisms [27]; chained beam constraint model and geometric parameter optimization, specially conceived for translational motion [28]; a combination of a beam constraint model, load equilibrium conditions, and geometric compatibility equations, specially conceived for 3-PPR compliant parallel mechanisms [29]; a kinetostatic modeling approach that integrates the screw theory with the energy method, with consequent avoidance of the problem of finding a solution to the equilibrium equations of nodal forces and the possibility of taking into account the parasitic deformations in space [30]; an extension of the chained beam constraint model specially revisited to analyze flexible beams of effective variable length [31]; a mathematical formulation of the compliance matrix method, combined with the inverse kinematic, specially introduced for modeling the flexure-based parallel compliant mechanisms with multiple actuation forces [32]; the adoption of three representations of multiple segments 2D beam models, namely, beam constraint model, linear Euler-Bernoulli beam, and PRBM [33]; the adoption of a new two-colored digraph representation of planar flexure-based compliant mechanisms for the automatic generation of the kinetostatic equations [34]; the introduction of virtual flexure hinges, link-flexure incidence matrices, and path matrices to generate automatically the formulation of the kinetostatic equations [35].
A more specific contribution has been dedicated in 2013 [36] to the solution of the problem of inverse kinetostatic analysis of a compliant four-bar linkage with flexible circular joints and pseudo-rigid bodies. This problem was attacked by extensively applying the theory of curved beams to the flexible parts, which gave rise to the closed-form symbolic expression of the compliance matrix, and by applying the static balance equations to both the elastic and pseudo-rigid parts.
The present investigation explores the opposite problem of direct kinetostatic analysis of a planar gripper with circular flexures. Two possible models based on the static balance of flexures are provided. The first linear model considers the static equilibrium in the undeformed configuration, while the second considers the balance in the deformed configuration. Both models simplify the fully non-linear model by exploiting a constant stiffness matrix of the undeformed curved beam element. These models allow us to find the tangent stiffness matrix in closed form as the sum of different contributions. Furthermore, dividing the tangent stiffness matrix into its contributions allows for evaluating each term's importance and setting strategies to speed up convergence. Furthermore, through a validation process of the fully non-linear model results performed on a case study, it will be possible to ascertain how the simplifications still provide accurate values in almost the entire mechanism's range of motion. As known, the tangent stiffness matrix is the heart of an iterative solving method. It is the basis of many implicit integrators widely used for the study of flexible mechanisms, such as the generalized α-method [37] or the HHT-method [38].
The main target of this article is to create two simplified models: • Solving the problem of the direct kinetostatic analysis of planar grippers with curved beams; • Being reliable in terms of motion accuracy and actuation forces; • Being computationally efficient to extend the formulation for real-time applications.
Any Finite Element Analysis (FEA) or Multibody Dynamics Simulation (MBDS) package is very reliable for solving any general problem in kinetostatic analysis numerically. Despite this, the availability of a ready-to-use independent algorithm to solve the direct kinetostatic problem gives rise to the possibility of implementing it in any real-time applications. Nevertheless, the MBDS Adams software has been used herein for validation purposes.
The paper is divided into the following sections. Section 2 gives the fundamentals of the curved beam model. Section 3 outlines the kinetostatic analysis. Two simplified linear and partial non-linear models are developed, and their tangent stiffness matrices are obtained in closed form. Section 4 includes a detailed case study description. Section 5 compares and validates the two models and gives essential insights into convergence and computational burden. Finally, Section 6 gives the concluding remarks.

The Adopted Curved Beam Model
Flexures employed in this context are curved beams. It has been demonstrated that curved beams can provide large rotations while maintaining small errors in terms of displacements of its center [39], as it is typical for classic revolute pairs. This feature is important to guarantee finite rotations of the end-effector in monolithic structures such as MEMS-based grippers. Furthermore, a linear model is capable of faithfully reproducing the displacements and in-plane rotation of the curved beam tip up to rotations of approximately ±20 • . This feature has the considerable advantage of using a constant stiffness matrix, as will be recalled below.
In the following, the curved beam compliance matrix, and its inverse stiffness matrix, will be recalled from [36]. Let us consider a curved beam with a circular profile of radius r f and beam characteristic angle θ f , as displayed in Figure 1. First, let us consider the general- containing the displacementξ f and the rotation angle φ f of the end section due to the deformation. Then, introducing the generalized wrench array w f = [F T f , M f ] T containing the force vector and the torque applied to the end section, the compliance matrix C f , derived in [36], follows from and depends only on the geometric and structural parameters of the curved beam, i.e., . . .
where E is Young's modulus and I is the area moment of inertia, assumed both constant for the circular profile. In the following sections, the inverse of the compliance matrix C f , i.e., the stiffness matrix K f will be employed to write the kinetostatic equations of planar mechanisms with curved beams. Furthermore, the stiffness matrix will be expressed in its locale frameŜ f attached to the end-section of the curved beam in the undeformed configuration, as shown in Figure 1.

Kinetostatic Analysis
Hereafter, all vectors denoted with the hat will refer to the undeformed configuration, while the same vectors will indicate the deformed configuration without the hat. Position vectors of local frames as well as rotation matrices of frames describing the orientation of bodies in the undeformed configuration are constant.
A curved beam links two components, as it happens for the two bodies displayed in Figure 2. First, consider the undeformed system composed of two rigid bodies, identified by the reference framesŜ i andŜ j , and by the flexuref . The vectorsr i andr j denote the positions of the body reference frame origins with respect to the fixed reference frame Σ. In contrast, the vectorsŝ i f andŝ j f , respectively, indicate the distance vectors going from the body-reference frame origins to the attachment points of the curved beam to the bodies. In the undeformed configuration, the position vectorp f going from the attachment point on body i to that on body j is obtained through the following expression Then, consider a generic configuration in which the bodies undergo finite displacements and rotations, and the flexure is deformed. For the assumption of rigid bodies, it follows thatŝ ≡s j f where the superscript denotes the reference frame in which the vector is expressed. In the previous expressions,s i f ands j f have been introduced to simplify the notation. Then, the following closure equation stands, where A i and A j are the rotation matrices mapping S i and S j into Σ and p f is the distance vector between the two flexure extremities in the deformed configuration.
Let us introduce the deformation vector x f containing the deformations of the flexure due to the displaced configuration described in Section 2. As known from the continuum mechanics, this vector can be represented using either the material or the spatial description of motion. In the following, only the material description is implemented. Therefore, expressing x f in the frameŜ i f of Figure 2, it follows where p f has been pulled back to the undeformed configuration as required in the material description of motion. Then, as recalled in Section 2, the circular flexure model requires x f to be expressed in the frameŜ f j instead ofŜ i f , therefore whereÂ f j is the constant rotation matrix mappingŜ f j intoŜ j . The rotation angle φ f due to the flexure deformation reads where θ andθ, respectively, are the rotation angles of the bodies in the spatial and material configurations and ∆θ, ∆θ denote the corresponding relative rotation angles. The generalized displacement array of the curved beam f , already introduced in Section 2, becomes

Jacobian of the Deformation Vector
Since the direct kinetostatic analysis will be solved using an iterative procedure, the variation of the generalized displacement array must be calculated. Using Equation (5), the variation δx f is whereĀ = ∂A/∂θ while δr, δθ are the variations of the body coordinates in the deformed configuration. The variation δx f is evaluated in the reference frame Σ but can be easily expressed inŜ f j remembering thatÂ f j andÂ j in Equation (6) are constant, i.e., δx Considering the angles, the variation of Equation (7), leads to Finally, the variations can be combined to form the variation of the generalized deformation vector δψ f . The latter satisfies the following expression where δq ij = [δq T i δq T j ] T is a 6-dimensional vector containing the variations of the body coordinates of bodies i and j andJ f is the (3 × 6) Jacobian matrix, defined as In derivingJ f , the property A T iĀ i =1 has been employed, being a particular skew-symmetric matrix used to define the cross-product in the planar case.

Kinetostatic Equations
The kinetostatic equations of the system require the static equilibrium of a curved beam. Consider the layout of Figure 3 showing the static balance of a curved beam connecting the bodies i and j. The deformation of the beam yields force and moment applied on the section S f j that must be equilibrated at sectionŜ i f . A first simplified model, hereafter referred to as the linear model, performs the balance in the undeformed configuration and leads to the following expressions where A f is the rotation matrix mappingŜ f j toŜ i f and d f is the distance vector between the two sections, respectively, defined as In Section 2, the stiffness matrix of the curved beam has been derived considering the undeformed configuration, meaning that the stiffness model is linear and cannot capture the geometrical non-linearity coming from the change of configuration during the beam deformation. Despite this, the balance in the deformed configuration can be modified including the tip displacement due to deformation. Referring to Figure 3, the tip displacement can be included in deriving the moment M i f at the first sectionŜ i f , therefore modifying the previous Equation (15) into where F f j and M f j are referred to the deformed configuration. It is noteworthy that this partial non-linear model is not the geometrically exact fully non-linear model of the curved beam since the forces and moment at sectionŜ f j are still obtained using a linear stiffness model for the curved beam. In the following, either the linear or the partial non-linear model will be included to derive the kinetostatic equations of a planar mechanism. Let us consider the body i in its deformed configuration, as displayed in Figure 4. The flexures have been removed and replaced with their reaction forces and torques where the minus signs come from Newton's third law. The static balance of body i requires that the following system be satisfied where F i and M i are the external force and torque applied to body i, respectively. From the balance Equations (15) and (17), the forces and moments coming from the flexures have been expressed in the undeformed configuration and are now turned into the deformed one. From Figure 4, it can be found that A 1i F i2 , hence it follows that Considering the frame invariance of the scalar equation of moments, the final system reads Then, denoting with and with R i the residual vector of body i, the final kinetostatic model is where w i = [F T i , M i ] T is the generalized force vector or wrench acting on body i, K 1 and K 2 are the stiffness matrices of the curved beams connected to the body, and T 2 is the matrix defined in Equation (17). If the latter is replaced by the matrix T 2 of the linear model in Equation (15), the kinetostatic model turns into The kinetostatic equations of a complex multibody system are derived by assembling the residual vectors of all rigid bodies. The final system can be cast in the form where R indicates the global residual vector, w is the global wrench of external forces and torques, and q is the global vector of generalized coordinates. The elastostatics model offers two types of analyses: the inverse and the direct kinetostatic analysis. The deformed configuration is the input, and the global wrench is the output in the inverse analysis, as has been already described in [36]. The solution of the inverse kinetostatic analysis is straightforward and does not require an iterative procedure. In the direct kinetostatic analysis, the forces and moments applied to the system are known, while the final configuration of the deformed mechanism is sought. This highly non-linear problem can be solved using an iterative procedure such as the Newton-Raphson method described in Algorithm 1. ← given threshold to terminate iterations 2: k = 1 ← iteration number 3: q (k) =q ← set the solution guess value 4: procedure NEWTON-RAPHSON ITERATIVE METHOD(q (k) ,q,w) 5: R(w, q (k) ) ← calculate the residuals as in Equation (24) 6: K T (q (k) ) ← calculate the tangent stiffness matrix as in Equation (26) 7: 10: exit procedure 11: else 12: k ← k + 1 13: goto step 4 14: end if 15: end procedure

Tangent Stiffness Matrix Determination
Suppose that the external forces and moments are fixed in space. Then, considering the residual of the partial non-linear model of Equation (23), the tangent stiffness matrix is where If the residual of the linear model of Equation (24) is used instead, the tangent stiffness matrix turns into The Appendix A reports the expressions for the terms of K Ti .

Case Study: The Four-Bar Linkage
In this section, a compliant gripper with curved beams, shown in Figure 5a, is studied. The monolithic structure of the mechanism can be reduced to two in-parallel four-bar linkages, as revealed in Figure 5b. For symmetry along the vertical axis, only half a mechanism will be analyzed, i.e., the right side of Figure 5b. The layout of Figure 6 has been plotted following the notation presented in Section 3. Body 1 is the frame, here considered fixed, while bodies 2, 3, and 4 are moving rigid bodies. Considering grippers with CSFH joints, one link is a part of the monolithic structure enclosed between two circular beams.
In an MEMS device, the entire monolithic structure can deform. However, the assumption of rigid links coupled to flexible circular beams has been verified using FE models. It is fully justified since the links are at least one order of magnitude stiffer than the circular beam flexures.  First, vectorsp f can be calculated knowing the undeformed configuration, i.e.,q, thereforê Similarly, the vectors p f of the flexures in the deformed configuration, i.e., q, are Since body 1 is the frame, r 1 =r 1 and A 1 =Â 1 . Furthermore, if the frame S 1 of body 1 is coincident with Σ, it follows that r 1 = 0 and A 1 = 1. Here, these matrices are written for the sake of completeness.
Notice that q could be one of the iterative solutions q (k) employed in the Newton-Raphson algorithm. The deformation vectors x f in the material description and expressed in the local frames of the undeformed flexures are The angular deformation φ f is obtained as The expressions (31) and (32) allows for determining the flexure generalized deformations The transformation matrices T f of Equation (17) are defined as where A f and d f can be found for each curved beam using Equation (16). Expressions similar to T f , not reported for brevity, can be written to determine T f of Equation (15). Then, the matrices N of Equation (22) are Setting the external wrenches w i applied at the mass centers G i of the rigid bodies, the four residual vectors R i , i = 1, . . . , 4, are The residuals are employed to form a system of 12 non-linear kinetostatic equations, i.e., that must be solved using an iterative procedure. To calculate the tangent stiffness matrix necessary to apply the Newton-Raphson algorithm, let us define the JacobiansJ The JacobiansJ can be mapped using Boolean matrices to the final dimension of the system, i.e., Then, following the expression (A2) of K I Ti reported in the Appendix, it yields The final expression for the first part of the tangent stiffness matrix is The expressions for K I I Ti can be obtained starting from z i in Equation (A7), i.e., where the matrices G of Equation (A4) are defined as Therefore, K I I T becomes Finally, the third part of the tangent stiffness matrix K I I I T is derived through the where F is obtained taking the force vector from the flexure wrench w The block-matricesJ f x , or J (Ŝ f j ) f x , are derived taking only the first two rows, i.e., the block of x f , from the corresponding Jacobian matrices. Using the equations of (A12), the matrices V f = V f i V f j can be obtained and therefore, the matrix K I I I T becomes The final expression for the tangent stiffness matrix is K T = K I T + K I I T + K I I I T . If the DOFs of the first body, i.e., the fixed frame, are removed by imposing fixed boundary conditions, the final form of K T will be a (9 × 9) matrix with the following pattern

Numerical Application
Referring to Figure 6, body 1 is fixed while body 2 is actuated through a vertical force applied at its center of mass. Finally, the end-effector is attached to body 3. Following the layout of Figure 6, all geometric and structural parameters necessary for direct kinetostatic analysis of the case study are reported in Table 1.

Comparison and Validation
First, let us consider the linear flexure model expressed through Equation (24) of Section 3. Considering an initial actuation force F 2y = −60 (µN), the latter is increased until the x-coordinate of the end-effector becomes zero, i.e., the gripper clamp is completely closed. The gripper deforms as displayed on the left side of Figure 7, wherein the limit cases and the undeformed configuration are reported. Then, let us consider the partial non-linear flexure model expressed through the Equation (23) of Section 3. Performing the same simulations, the workspace becomes that of the right side of Figure 7. The values of F 2y for which the clamp is closed, respectively, are F 2y = 63 (µN) for the linear model and F 2y = 77 (µN) for the partial non-linear model. The two models have been compared to a model implemented using the commercial multibody software MSC Adams © . The model includes rigid links and flexible circular beams. The latter are combined to the links' endpoints through fixed connections. The model has been developed to be comparable with the Matlab model. The only difference pertains to the flexures since the curved beams have been modeled using a two-dimensional geometrical non-linear representation for beam-like structures. Compared to the previous Matlab models, this representation is fully non-linear, as the stiffness matrix of each flex-ure is updated during deformation. The actuation force has been applied through step functions, and two simulations have been carried out to achieve convergency; one for the forward movement of closing and the other for the backward movement of opening. Figure 8 shows the two simulations and the trajectories accomplished by the end-effector point. As for Figure 7, the deformed configuration has also been plotted to understand the deformation process better. The three models have been compared in Figure 9 in terms of the end-effector trajectory, also referred to as the mechanism's workspace (left subplot) and actuation forces (right subplot).
Two materials have been modeled for the curved beams: silicon with Young modulus E y = 100 (GPa) and nylon with Young modulus E y = 3.84 (GPa). Considering the same rectangular cross-section, whose dimensions are reported in Table 1, the silicon bending stiffness is EI = 26.0417 (µNµm 2 ) while the nylon bending stiffness is EI = 1 (µNµm 2 ). The two materials have different properties in terms of elasticity. Silicon has a brittle behavior, while nylon has an anisotropic hyperelastic or visco-hyperelastic behavior. The two bending stiffness values should be seen as extreme cases to test the proposed method. With this premise in mind, the two materials will be assumed to have both isotropic linear behavior, while the different degrees of non-linearity of the models will only concern the geometric stiffness.
The first row of plots in Figure 9 pertains to the silicon while the second one is the nylon. First, let us consider silicon. Observing the top-left subplot, the three arc-shaped trajectories of the end-effector reveal relevant differences only in the final part of the path. The influence of the fully non-linear flexures becomes more evident in the opening movement, where the trajectories become more distant. Compared to the linear model, the top-right subplot reveals that the flexure non-linearities introduce a stiffening effect, and the actuation force required to produce the same displacement grows. It can be observed that the partial non-linear model is stiffer than the fully non-linear model in the forward closing movement and softer in the opening movement. The three models have equal stiffness only at the undeformed configuration where the actuation force is zero. Then, let us consider nylon. Observing the bottom-left subplot, the end-effector trajectory follows a trend similar to the previous case. The bottom-right plot reveals differences in force range, as could be expected considering the lower bending stiffness of the nylon. Now, the non-linearities are more pronounced, and three inflection points appear for the fully non-linear flexures, which are totally absent in the remaining cases. Despite this, the plot is similar to the simplified cases.
Excluding the limit points of the workspace during the opening movement, the simplified models provide excellent results. Added to this is that the proper workspace is usually limited by other constraints such as the electrical interfaces or the maximum stress in the material, thus making the three models closer than they might appear in Figure 9.
For example, the Conjugate Surfaces Flexure Hinges (CSFH) employed in MEMS micro-grippers have a rotation range limited to ±20 • to prevent the silicon from breaking [39,40]. This range is displayed in the opening-closing movement of Figure 9. However, this range is further limited to about ±2 • by other phenomena coming from the electrostatic actuation such as sticking-friction anomalies, the pull-in or the impossibility to generating high actuation forces [41].

Shape Optimization
The simplified models have been employed to perform the cross-section optimization of the curved beams, as shown in Figure 10. These plots can be used in various ways; for example, knowing the maximum actuation force the electrical interface can produce, the section parameters can be chosen to cope with this value. Another example could be related to the choice of the section parameters based on the maximum allowable stress of the material or its fatigue limit. Likewise, the optimization could affect other structural parameters or mechanism lengths.

Tangent Stiffness Matrix
Since the tangent stiffness matrix is the key element of the direct kinetostatic analysis, in the following, a detailed analysis of the tangent stiffness matrix and its role in the Newton-Raphson algorithm convergence is detailed. As already described in Section 4, after imposing fixed boundary conditions on body 1, the tangent stiffness matrix turns into a 9 × 9 symmetric matrix. Referring to the partial non-linear model, K T has the expression reported in Table 2 in the undeformed configuration. Let us consider the mechanism in the final deformed configuration obtained by applying F 2y = 50 (µN). In this case, the expression of K T is reported in Table 3, while the percentage difference between the deformed and undeformed case is provided in Table 4. It can be observed that relevant differences appear during the deformation process, especially in the off-diagonal components. Furthermore, the matrix K T is no longer symmetrical in the deformed configuration.  Since K T is composed of three terms, it is legitimate to ask what the contribution of each term is. As reported in Table 5, the first term K I T is the closest to the final expression of K T . The second term K I I T is reported in Table 6. Remembering the expression for z i in Equation (A7) and the form of K I I T in Equation (A8), observing Table 6, it can be found that at the equilibrium, i.e., if and only if the residuals are zero, the following expression stands For the case study, only body 2, and therefore z 2 , has components different from zero. Finally, the third term K I I I T is reported in Table 7. It can be noticed that only the components of the inner moments, i.e., due to the flexures, are different from zero. This feature comes from the particular form of V f i , or V f j , in Equation (A12).  To better focus on the importance of the tangent stiffness matrix in achieving solution convergence, let us keep the tangent stiffness matrix constant and equal to that obtained in the undeformed configuration for the entire simulation, i.e., K T = K T (q). The results of Figure 11 reveal that the number of iterations necessary to achieve convergence grows exponentially as the input load increases and no longer converges beyond F 2y = 76 (µN). The same simulation has been repeated by updating the tangent stiffness matrix following four strategies based on: • The complete matrix K T starting each simulation from the undeformed solution; • The first term K I T starting each simulation from the undeformed solution; • The complete matrix K T starting each simulation from the previous converged solution; • The first term K I T starting each simulation from the previous converged solution. It can be seen that the approaches using the previous converged solution reduce the number of iterations. Similarly, using the complete matrix instead of the first-term approximated matrix results in fewer iterations. The results are displayed in Figure 12.
It is interesting to observe how these trends translate into a computational burden. The reduced number of iterations provided by the strategies based on the previous converged solution translates directly into savings in computation time. The two strategies based on the first term of the tangent stiffness matrix lead to a higher number of iterations but, simultaneously, require a smaller number of variables to be determined and save on calculation times, as shown in Figure 13. The peaks observed in Figure 13 come from memory allocation and other inner processes of Matlab during the first computation. On the other hand, it can be observed from Figure 12 that a higher number of iterations do not correspond to these CPU times. The results have been obtained using an HP workstation equipped with an Intel Xeon CPU @3.20GHz with 32 GB of RAM.  This section is concluded by giving some insights into the computational time obtained using Adams. In Table 8, the CPU times obtained in the opening-closing movement for the Adams fully non-linear model, the Matlab linear model, and the Matlab partial nonlinear model are compared. Adams simulations have been performed by disabling the graphic display. The comparison has been carried out for both silicon and nylon. As can be observed, the CPU time of the simplified methods is from 30 to 500 times faster than Adams. It is noteworthy that the simulation time for nylon is ten times faster than silicon for the Adams fully non-linear model.

Conclusions and Discussion
The tangent stiffness matrix has been used as a conceptual base to solve the direct kinetostatic problem of planar grippers with curved beams. Two models have been presented to cope with flexure deformations. The first linear model considers the flexure equilibrium in the initial undeformed configuration, while the second partial non-linear model considers the equilibrium in the deformed configuration. Both methods do not include a fully non-linear geometric description of the curved beam flexure whose stiffness matrix is kept constant and equal to that obtained in the undeformed state. The tangent stiffness matrix has been divided into sub-parts to facilitate both the theoretical treatment and the numerical implementation. The linear model led to two sub-parts, while the partial non-linear model introduced a further third sub-part. Both models were tested and compared with a fully non-linear model obtained using the commercial software MSC Adams. The results proved to be in good agreement on most of the mechanism's workspace, except for the extreme areas wherein the geometric non-linear effects become relevant. The same case study was used to show the method's potential; for example, in conducting a shape optimization of the flexure cross-section. Finally, the importance of each term of the tangent stiffness matrix in the convergence process was detailed in terms of the number of iterations required to achieve convergence and computational load.
From what has been outlined, the proposed method offers various advantages: 1.
The results of Figure 13 suggest a possible extension to real-time applications of micro and nano-grippers. It is known that the control often requires simplified models to be executed quickly by the control unit. Often these models are obtained by linearizing the equilibrium equations around one or more operating points. Using models with reduced complexity would allow more efficient control strategies such as control in the operating space, inverse dynamics control, pre-calculated torque control. Furthermore, the closed form helps creating more efficient reduced order models [42][43][44].

2.
The tangent stiffness matrix is obtained in closed form. This feature prevents the use of numeric differentiation, making the convergence process of the direct kinetostatic solution more robust. Furthermore, splitting the expression of K T allowed for identifying its most basic terms and calibrating the compromise between the number of iterations and calculation time. The calculation times are considerably reduced by using only the first term of the tangent stiffness matrix and recalculating it at each iteration of the Newton-Raphson algorithm described in Algorithm 1.

3.
The tangent stiffness matrix can be employed to develop a dynamics model to study vibrations. The tangent stiffness matrix is the core of implicit time integration methods primarily employed in flexible multibody dynamics [45]. Shape optimization takes further advantage of the closed form of K T opening scenarios to gradient-based constrained optimization problems based on the kinetostatic analysis.

4.
Both the two simplified models employ curved beams modeled by a constant stiffness matrix. Despite this, the curved beams guarantee finite displacements/rotations in the mechanism, allowing for the expansion of the reachable workspace. The model remains reliable for most of the mechanism's workspace. The results are accurate in the functional area except for the limit zones of the workspace in which physical constraints usually prevent motion. When the maximum rotations of the curved beams exceed ±20 • the constant stiffness hypothesis can no longer capture the geometric nonlinearities, and the results deviate from the actual case.

5.
Although the proposed method is valid only for planar cases, it can be extended to other compliant mechanisms with constant stiffness flexures without changing the mathematical background. where V f is a (3×6) matrix while the 6-dimensional vector