DAS-2D: a concept design tool for compliant mechanisms

Compliant mechanisms utilize the deformation of the elastic members to achieve the desired motion. Currently, design and analysis of compliant mechanisms rely on several commercial dynamics and finite element simulation tools. However, these tools do not implement the most recently developed theories in compliant mechanism research. In this article, we present DAS-2D (Design, Analysis and Synthesis), a conceptual design tool which integrates the recently developed pseudo-rigid-body models and kinetostatic analysis/synthesis theories for compliant mechanisms. Coded in Matlab, the software features a kinematic solver for general rigid-body mechanisms, a kinetostatic solver for compliant mechanisms and a fully interactive graphical user interface. The implementation details of all modules of the program are presented and demonstrated with four different case studies. This tool can be beneficial to classroom teaching as well as engineering practices in design of compliant mechanisms.


Introduction
Rigid-body mechanisms (Uicker et al., 2003) transfer motion of rigid links connected by kinematic joints. Likewise, compliant mechanisms (Howell, 2001) perform the same tasks via the deflections of the flexible members. It is essential during the design process to perform kinematic and static analysis on mechanisms to understand if displacements, velocities, and accelerations are acceptable for the desired function and to determine the required loads to accomplish the predetermined function. In the past two decades, three most important methods have been developed to analyze elastic deformation of compliant mechanisms: pseudo-rigid-body (PRB) model approach, continuum mechanics and finite element approach.
The PRB model approach substitutes compliant links with a series of (typically two to four) rigid segments joined with torsion or linear springs. Howell and Midha (Howell, 1991;Howell and Midha, 1995) had developed the earliest PRB 1R model which replaces elastic members with two rigid-body links connected with a torsion spring. "R" and "P" represent a revolute and a prismatic joint for a PRB model. PRB 1R model also has two parameters γ and K , characteristic ra-dius factor and equivalent spring stiffness respectively, which needs to be tuned for different loading conditions. Saxena and Kramer (1998) improved the 1R model by switching to two linear springs to limit the change of the characteristic radius factor. Su (2009) proposed a PRB-3R model that is load independent within a load range. Saggere and Kota (1999) developed a PRB-FSM (Finite-Segment) model in which the beam is divided into (n + 1) ≥ 3 segments joined by n torsional springs of stiffness EI / l.
In terms of computer-aided design of mechanisms, there has been a number of static solvers that have been developed throughout the years and KINMAC and STATMAC by Paul (1979) are maybe the first examples of a kinematic and static analysis tool. SAM (Rankers and Schrama, 2002) is a commercial software capable of kinematic and static analysis limited to the rigid-body linkages. More recently, CoMeT (Petri, 2002;Culpepper and Kim, 2004) was developed for static analysis of compliant mechanisms which is not limited to the planer mechanisms. However, direct stiffness method was employed in the software which is only accurate for small deflections of elastic members limiting the capabilities of the software. Geometrically Exact Beam Theory (GEBT) (Yu and Blair, 2012) is a FEM tool for static analysis of Published by Copernicus Publications. 136 O. A. Turkkan and H.-J. Su: DAS-2D: a concept design tool for compliant mechanisms structures but it is not possible to analyze mechanisms that contain moving joints. Similarly, SPACAR (Jonker and Meijaard, 1990) is a software package for dynamic modeling and analysis of multibody flexible systems with finite element analysis.
WorkingModel 2D is a commercial software for multibody dynamic simulation of planar physical systems. It has been often used in design and simulation of planar mechanisms. Adams (MSC Software) is an another commercial software capable of dynamic simulation of planar mechanisms. MATLAB's SimMechanics toolbox is a powerful multi-body dynamics solver however the design process with block diagrams is not entirely intuitive and the toolbox lacks of proper interactivity between the user and the simulation. Yue et al. (2012) have developed virtual reality user interfaces that communicate with SimMechanics solver for interactive design of planar mechanisms. Ch Mechanism Toolkit (Cheng and Trang, 2005) enables users to create and perform kinematic analysis on rigid body mechanisms utilizing the Ch language. However, for each individual mechanism, a different class file must be prepared and kinematic equations must be manually inputted.
Given numerous advances in compliant mechanism theories in the past two decades, there is no design software that is dedicated to design of compliant mechanisms. In this paper, we aim to fill this gap by developing a conceptual design tool that integrates kinetostatic analysis and synthesis theory for design of compliant mechanisms. Kinematics of rigid-body mechanisms will also be a side product of this software since they are considered as a special subset of compliant mechanisms. The software was developed in MATLAB in order to take advantage of the built-in functions, such as nonlinear equation solvers and optimization routines. Also, the software tool has a fully interactive graphical user interface to aid users during the design process.
First section of the paper explains the theory behind the kinematic and static analysis modules of the software. Then, a detailed explanation of the implementation of these theories for the different analysis modules are presented with flowcharts and algorithms. An overview of the graphical user interface is given in the next section. Finally, four case studies are illustrated to demonstrate and verify different analysis modules.

Theory and mathematical formulation
It is well known that analysis of compliant mechanisms involves simultaneously solving a set of kinematic equations coupled with static force equilibrium equations which are called kinetostatic (kinematic and static) equations. In this section, we present the mathematical formulation of kinetostatic equations.  Figure 1. Two of the three available independent loops that are required to analyze the mechanism are shown.

Kinematic vector loop closure equations
Complex number method is one of the most commonly used methods in kinematic analysis of rigid mechanisms: where i = √ −1 is the complex number unit. In general, the vector loop equation for a mechanism with l independent closed loops can be written in the complex form: where coefficient C ki can be −1, 0 or 1. For instance, the mechanism shown in Fig. 1 has l = 2 independent loops. The two vector loop equations are: Each complex equation will produce equations in x and y directions, while producing four equations at total, for the mechanism shown in Fig. 1. If one link is driven, these four equations must be solved simultaneously to determine the final configuration.

The pseudo-rigid-body model for compliant links
Compliant mechanisms have at least one compliant link that can be deformed under external forces. Pseudo-rigid-body models (PRB) have been widely used in the analysis and synthesis of compliant mechanisms. In this approach, n rigid segments are joined with torsion and/or linear springs.
where each row represents a different segment in the flexible beam. Each row consists of length of the segment γ , the axial stiffness of the spring k ex and the magnitude of the torsion spring k θn . This representation of compliant beams is adopted in the software and each individual flexible beam is represented with a predefined but customizable PRB matrix. Different PRB models have varying accuracies and thus, the accuracy of static analysis greatly depends on the PRB models used in the analysis. The accuracy of pseudo-rigidbody models can be illustrated with a basic example of a single beam loaded from one end and fixed to wall at the other end. The beam is discretized using PRB-3R and PRB-FSM (11 segments) methods and compared with Bernoulli-Euler large deflection equation and GEBT. The first test case is a large force and moment acting in the same direction and the next case is an opposing force and moment. Table 1 shows good agreement between PRB models and the analytical model. The error with the analytical model comes from the discretization of the compliant beam (number of segments) and the superiority of the PRB-FSM model over the PRB-3R model comes from the number of rigid-link segments used in discretization of the elastic beam.
Typically, PRB models with higher number of segments yield more accurate results. PRB-FSM model with more than ten segments will be always reasonably accurate for any compliant mechanism regardless of the loading condition. However, employing PRB-FSM model can increase the degrees of freedom of a compliant mechanism beyond a limit where an analysis is not possible within an acceptable time frame. There exists some PRB models (Su, 2009; that are made of three or four segments and are very accurate for a range of loading conditions. These models produce faster results without compromising the accuracy, although special care must be taken not to exceed the optimized loading limit of the PRB model. Sometimes, extension or shortening of the compliant members is necessary in order to drive a compliant mechanism and therefore, PRB models with linear springs must be employed. PRB-FSM model with linear springs will be again accurate but slow. Venkiteswaran and Su (2014) constructed a library of PRB models with linear springs that had been optimized for short beams, but the accuracy of the models with more than three segments are not guaranteed for long slender beams.

Kinetostatic analysis of compliant mechanisms
Kinetostatic analysis is defined as determining the deflection of a compliant mechanism upon a given load or vice versa. Several methods exist for derivation of kinetostatic equations (coupled kinematic and static equilibrium). Direct stiffness approach (Melosh, 1963;Hughes, 2012;Ranzi et al., 2004) can be used to perform linear static analysis. This method is currently being used in static analysis of rigid-body (SAM, Rankers and Schrama, 2002) and compliant mechanisms (CoMeT, Culpepper and Kim, 2004). Although it has advantage of simple formulation, the stiffness approach is based on linear models and it is not accurate for large deflections.
On the other hand, energy methods (Reddy, 2002) can be used for nonlinear static analysis. Virtual work principle was employed to calculate deformation of parallel mechanisms (Xiangzhou et al., 2007;Jesús Cervantes-Sánchez et al., 2012). Optimization of potential energy has been used for a wide variety of applications from molecular mechanics (Pillardy et al., 2001) to human walking (Zarrugh et al., 1974). In this article, we implement the energy based approach for formulation of the kinetostatic equations. The total work done on a mechanism is the sum of work done by external forces (forces and moments) and internal forces (springs), Total work done on the mechanism by spring forces is the negative of the work done on the springs: where n and m are the number of linear springs and the number of torsion springs, respectively. The first sum represents the total work done by linear springs, whereas the second sum represents the work done by torsion springs. Total work done on the system by external forces (along a path − → r ) and external moments is calculated as: where l and s are the number of external forces and the number of external moments, respectively. The first sum represents the total work done by external forces, whereas the second sum represents the work done by external moments. Potential energy of the system is the negative of the work function; i.e., E = −W = −U −V . The main principle of the virtual work leads to δE = δE δq 1 δq 1 + . . . + δE δq n δq n = −δW = 0.
Equation (9) states that a stationary potential energy is the necessary and sufficient condition for equilibrium of a conservative system. This means that equilibrium position (stationary) point can be obtained from minimization of the potential energy. This optimization must be subjected to the kinematic constraint equations defined in Eq. (2). Finally, the optimization problem is mathematically formulated as: where (k, ψ) equals to (k e , x) and (k θ , θ ) for linear and torsion springs, respectively. Equation (10) is very similar to the minimization of potential energy (MinPE) method by Aten et al. (2011). MinPe method handles the external loads by adding force equilibrium equations at the links where loads are applied. The force balance equations may bring additional optimization variables and implementation difficulties compared to the Eq. (10). However, the MinPE method has the advantage of having the joint loading information without any further post-processing after minimization.
If the direction of a force is fixed during minimization at Eq. (10), the force integral will reduce to F x (x −x 0 )+F y (y − y 0 ) where F x and F y are the force magnitudes in the x and y directions, x and y are the distance to the origin during the optimization, and x 0 and y 0 are the initial distance to the origin. Since −F x x 0 −F y y 0 is constant during minimization, the accuracy of the optimization will not be affected if the optimization is done in one or more steps. If follower forces that change direction with the applied link exists, the force integral will be path dependent. The magnitude of all loads must be incrementally increased (Algorithm 1) in order to have better accuracy. The accuracy of the force integral will then be parallel with the number of increments. Similar to the non-follower forces, the accuracy of the moment integral is independent of the number of increments.

Implementation
In this section, we describe the implementation details of DAS (Design, Analysis and Synthesis) 2D, a computer-aided design tool for planar compliant mechanisms. The compliant mechanism design software is implemented in MATLAB using object oriented programming (Fig. 3) while taking advantage of the built-in nonlinear equation solver, optimizer and plotting tools. Also, graphical user interfaces have been implemented to facilitate the design process.

Graphical representation of mechanisms
Graph theory is often employed in representing mechanisms (Dobrjanskyj and Freudenstein, 1967;Chen and Lin, 1998;Shi and McPhee, 2000). This representation is widely used in classification and type synthesis of rigid-body mechanisms. The two different graphical representations of the same slider-crank mechanism are showed in Fig. 4. At the top, Freudenstein and Maki (1979) give an example of the common practice in which vertices and edges represent the links and the kinematic joints, respectively. Since rigidly connected vertices are represented with a single vertex, this representation is not well suited for implementation. Also, only binary joints can be handled with this representation   without breaking multiple joints into binary joints. A new representation is adapted over the traditional method and in this new representation, vertices represent the nodes (or the joints) and the edges denote the links between the joints. This new representation is parallel with how individual links are 5 stored in the software.

Kinematic analysis
The main challenge of kinematic analysis is finding the independent loops for formulating the kinematic constraint equations. In the graph theory, the topology of any mechanisms 10 can be represented with an adjacency matrix and the number of independent kinematic loops can be calculated with the Euler's formula. Once the adjacency matrix and the number of independent kinematic loops are determined, back edges of the graph can be found by a depth first search algorithm 15 (Algorithm 2). The back edges are then converted to base cycles by adding the parent nodes until the start and end nodes are the same.
Algorithm 2 Depth first search for finding the back edges of a graph. After all independent kinematic loops are found, nonlinear kinematic equations of the mechanism can be derived. The 20 kinematic equation for the ith link can be written as: www.mech-sci.net/7/1/2016/ Mech. Sci., 7, 1-14, 2016 without breaking multiple joints into binary joints. A new representation is adapted over the traditional method and in this new representation, vertices represent the nodes (or the joints) and the edges denote the links between the joints. This new representation is parallel with how individual links are stored in the software.

Kinematic analysis
The main challenge of kinematic analysis is finding the independent loops for formulating the kinematic constraint equations. In the graph theory, the topology of any mechanisms can be represented with an adjacency matrix and the number of independent kinematic loops can be calculated with the Euler's formula. Once the adjacency matrix and the number of independent kinematic loops are determined, back edges of the graph can be found by a depth first search algorithm (Algorithm 2). The back edges are then converted to base cycles by adding the parent nodes until the start and end nodes are the same.
Algorithm 2 Depth first search for finding the back edges of a graph. After all independent kinematic loops are found, nonlinear kinematic equations of the mechanism can be derived. The kinematic equation for the ith link can be written as:  where Z i0 and λ i are constant design parameters for the link and f is a function for which input is the unknown parameter(s) needed to define the link. As an example, for a binary link with two pin joints, the term Z i0 is zero, λ i = r i is length of the link and f = (cos θ i , sin θ i ) T where θ i is the link angle. Following this notation, for a mechanism with n links and l independent loops, the kinematic constraint equations can be systematically derived as n vector equations or 2n scalar equations: where coefficient C ki can be −1, 0 or 1. By applying Eq. (13), kinematic equations will be formed and they can be solved simultaneously using a nonlinear equation solver with the number of inputs equal or less than the degrees of freedom of the mechanism. The interactive kinematic analysis capabilities of the software was previously demonstrated at Turkkan and Su (2014).

Static (force) analysis
Static force analysis can be described as determining the relationship between the external loading and the deflection of the compliant mechanism. Figure 5 illustrates the schematic view of the static analysis problem and the flowchart of the static solver. The module starts with converting compliant links into a series of rigid-body links using a pseudo-rigidbody model as described earlier. Then, kinematic equations are derived based on the independent kinematic loops. A breadth-first search algorithm is used to determine the deflection of the point where the external force is applied. Finally, the optimization in Eq. (10) and Algorithm 1 are employed in the last optimization step to determine the resulting configuration of the mechanism.

Distance analysis
Contrary to the static force analysis, the distance analysis defined as determining the required load(s) that will result in a prescribed deformation. As shown in Fig. 6, the desired deformation can be translation of a node or rotation of a link. The output of this module will be the magnitudes of the unknown loads.
After designing a mechanism, the desired motion of a node or a link and unknown loads on the mechanism are defined. Then, the magnitudes of unknown loads are obtained by an optimization process which minimizes the difference between current deformation and the desired deformation.

Mechanical advantage analysis
Mechanical advantage is defined as the ratio of the output force / moment over the input force / moment, e.g. M o /M i for the case in Fig. 7. It is one of the most important design criterion in compliant mechanisms. Mechanical advantage analysis is to find the balancing load for a given input load. For a mechanism, one load is selected as the input force and another load is selected as the output force. The input force is assumed to be an unit load. Then, an optimization process shown in Fig. 7 is initialized. At each optimization step, magnitude of the output force is guessed and resulting configuration is calculated via force analysis. The displacement of the nodes are minimized to zero in order to ensure that the mechanism stays at the initial design configuration.

Load and energy curve plotting
The energy/torque vs. the input load curve gives the designer a visualized way to evaluate the quality of a compliant mechanism. Algorithm 3 illustrates the procedure for determining energy and torque values. The algorithm starts with creating a equally spaced displacement list between the initial and the target displacement (link rotation or a slider displacement). At each iteration, the corresponding link angle or slider position is fixed and the potential energy of the system is minimized while subjected to kinematic constraints of the mechanism. The total potential energy is calculated and finally, torque is calculated by taking derivative of the potential energy stored in springs with a central finite difference method.

Graphical user interfaces
To expedite the design process, a graphical user interface has been developed with MATLAB and interfaces for the some of the design and analysis modules are shown in Fig. 8. Panel (a) shows a step in the design process. The PRB model of a compliant link can be specified using a dedicated user interface shown in panel (b). Panel (c) shows the kinematic analysis. Users can interact with a mechanism via mouse operations (dragging, clicks) which trigger the kinematic analysis class. Similarly, a mechanism can be animated within a range of input values using this module. Panels (d), (e) and (f) display the load analysis, distance analysis and energy analysis interfaces, respectively.

Case studies
In this section, four case studies are presented to demonstrate the capability and the validity of the kinetostatic analysis software described in the previous sections.

Parallelogram flexure mechanism
To test the static solver, we give an example of a parallelogram flexure mechanism which consists of two identical compliant beams (length = 100 mm, in-plane width = 1 mm, out-plane depth = 20 mm, E = 72 GPa), as shown in Fig. 9. Normalized (dimensionless) loads exerting on the rigid link are defined as: Static analysis is performed on the mechanism for a large range of load magnitude (f x , f y , m z ∈ [1, 10]). Two different PRB models are employed and the two analysis results are compared with Beam Constraint Model (BCM) Awtar et al. (2006) analytical model and the Abaqus FEM software. Figure 9 shows the software output when normalized loads are equal to 10 and it can be seen that quite large y displacement (≈ 28 mm) of the rigid link is present. Therefore, this range of normalized loads are satisfactory in investigating accuracy of the software for small and large displacements.

Compound multibeam parallelogram mechanisms
Compound multibeam parallelogram mechanisms (CMPM) are composed of multiple parallelogram flexure mechanisms. Li (2015, 2016) developed an analytical model for nonlinear deformation of CMPMs. The study found out that the configuration at Fig. 11 has nonlinear displacement characteristics at the primary direction due to load-stiffening effect. Since lengths of the compliant members must increase during deformation, only PRB models with linear springs can be employed during analysis. The geometric properties of a compliant member are as follows: in-plane width = 1 mm, out-plane depth = 10 mm and length = 40 mm. The elastic modulus and Poisson's ratio of the elastic members are 2.3 GPa and 0.37, respectively. The force is applied in the primary direction at the rigid link of length 30 mm. Three different PRB models (PRB-RPR, PRB-2(PR)R, PRB-FSM with 10 springs; Table A1) are employed in the analysis and the results are compared with analytical model Figure 11. A typical compound multibeam parallelogram mechanism loaded at the primary direction.
at Fig. 12. The analysis requires approximately 0.15, 0.2 and 0.4 s on a computer with i7 processor at 3.4 GHz for the PRB-RPR, PRB-2(RPR)R and PRB-FSM, respectively. All PRB models were able to capture load stiffening effect in the primary translation direction, although PRB-RPR model fails to capture the deflection accurately. The analytical model is very accurate until very large forces, whereas PRB-FSM model follows the Abaqus results closely in all load levels. PRB-2(PR)R model is not as accurate as the PRB-FSM model for smaller load levels, but converges to Abaqus results at higher load levels.  (Hao and Li, 2015) and ABAQUS software. f y is normalized as shown in Eq. (14).

Distance analysis
We considered a fully compliant bistable mechanism (Wilcox and Howell, 2005) shown in Fig. 13 to test the distance analysis module. The long beams are approximated with PRB models that contain three linear springs and two torsion springs. The linear spring is utilized to capture the linear elongation of the long beams. The short compliant beams are represented with the PRB-3R model. The cross sections of the short and long elastic beams are comparable in size but the elastic modulus of the short beam is approximately 30 times larger resulting in short beams being much stiffer than the longer beams. The length of the center rigid link is 30 mm. The target vertical displacement of the mechanism is set to y ∈ [1, 10] mm and three different PRB models (PRB-RPR, PRB-2(PR)P, PRB-2(RP)R; Table A1) were chosen to represent the long compliant beam. Figure 13 also shows the analysis results compared with the Abaqus simulation result. The newly developed model, PRB-2(RP)R, is very close to the Abaqus result. Although the other two models come closer to the Abaqus result for large deflections, they fail to capture the general trend in the Abaqus result.

Load and energy curve plotting
A compliant four bar bistable mechanism (Fig. 14) with a flexible link was designed and 90 • clockwise rotation of the crank link was specified as the target motion. The lengths of the links are 15, 37.0842 and 43.294 mm for the crank, coupler and compliant links, respectively. The width and depth of the compliant beam is 1.5 and 5 mm, respectively. Steel (E = 207 GPa) was chosen as the design material of the compliant beam. Figure 14 shows the comparison of the output of the program compared with the Adams software. The PRB-FSM (10 torsion springs) model was used to represent the elastic link. The difference between the torques for the two softwares is around 6 % for the negative peak and 10 % for the positive peak. Since the location of the stable and unstable points is essentially a kinematics problem, (un)stable points happen at same angle values in the both cases. Also, it was found that the torque values are highly sensitive to the mesh of the elastic link in the Adams. Although not shown, if the same PRB model is implemented in Adams, both tools give the exact same results.

Conclusions
This paper presents DAS-2D, a kinetostatic solver for conceptual design of planar compliant mechanisms. The recently developed theories in compliant mechanism research such as pseudo-rigid-body models have been implemented. Static force analysis is performed by minimizing the total potential energy of the system. Considering rigid-body mechanisms as a special case of compliant mechanisms, kinematic analysis routines based on vector loop closure equations have been developed. Implementation details of the different modules are presented and they are demonstrated with four representative case studies. Future work includes the incorporation of kinematic and static synthesis of complaint mechanisms and modeling compliant members with analytical models that can provide strain energy calculation for the Eq. (10) such as Beam Constraint Model (Awtar and Sen, 2010) and Chained Beam Constraint Model (Ma and Chen, 2015). Also, the PRB models optimized by Venkiteswaran and  were found to be inadequate for slender beams and a new PRB model with linear springs will be optimized to be employed in the software for analysis of slender beams.   Table A1. List of the PRB Matrix values for the PRB Models used in the paper. First row of the PRB matrix of a PRB-FSM model with n torsion springs is 0.5/n ∞ ∞ , the last row (n + 1th) is 1/n n ∞ and the remaining rows are 1/n n ∞ . If extension effect is desired, k ex will be 2n, 2n and n for the first, the last and the remaining rows of the PRB Matrix of the PRB-FSM model with n torsion springs, respectively. The parameters of the PRB-RPR and PRB-2(PR)P are obtained from Venkiteswaran and Su (2014)