A CPRBM-based method for large-deflection analysis of contact-aided compliant mechanisms considering beam-to-beam contacts

Mechanism and Machine Theory 145 (2020) 103700

Contents lists available at ScienceDirect

Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmachtheory

Research paper

A CPRBM-based method for large-deflection analysis of contact-aided compliant mechanisms considering beam-to-beam contacts Mohui Jin a,b, Benliang Zhu c, Jiasi Mo a, Zhou Yang a,∗, Xianmin Zhang c, Larry L. Howell b a

College of Engineering, South China Agricultural University, Guangzhou, China Department of Mechanical Engineering, Brigham Young University, Provo, UT, USA c School of Mechanical and Automotive Engineering, South China University of Technology, Guangzhou, China b

a r t i c l e

i n f o

Article history: Received 30 July 2019 Revised 6 November 2019 Accepted 7 November 2019

Keywords: Compliant mechanisms Contact analysis Large deflection Chained pseudo-rigid-body model

a b s t r a c t Contact-aided compliant mechanisms (CCMs) utilize contact to achieve enhanced functionality. The contact phenomenon of CCMs increases the difficulties of their analysis and design, especially when they exhibit beam-to-beam contact. Considering the particularity of CCMs analysis, which is more about the mechanisms’ deformation, this paper presents a numerical method to analyze the large deflection and stress of the CCMs considering beam-to-beam contacts. Based on our previous work on beam-to-rigid contact, the large deformation of general beams in CCMs is modeled by using the chained pseudo-rigid-body model (CPRBM). An approximation based on the geometric information of CPRBM is proposed in this paper to rapidly determine the moving boundary curve for beam-to-beam contact constraints. The static equilibrium configuration of CCMs is solved by minimizing its potential energy function under the geometric constraints from the boundary curves of contacts. A formulation is also provided to evaluate the normal stress along the deformed beam based on the deformation of CPRBM’s torsional springs. Numerical examples and finite element analysis are used to verify the feasibility and accuracy of the proposed method. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Due to structural characteristics, the force and motion responses of compliant mechanisms [1] are usually smooth. As a particular kind of compliant mechanisms, contact-aided compliant mechanisms (CCMs) [2] utilize contact and large deflection to produce enhanced functionality, e.g., non-smooth path, motion and force-deflection characteristics [3]. This paper focuses on beam-based CCMs, in which the contact can be divided into two categories, i.e., beam-to-beam contact and beam-to-rigid contact. The beam-to-beam contact occurs when two beams of CCMs contact each other [4], while the beamto-rigid contact means that the CCMs interacts with rigid object [5] (Fig. 1). Thanks to their advantages, CCMs have been widely used to design compliant mechanisms with non-smooth motion [6], compliant joints [7–11], flapping wings [12–14], nonlinear stiffness actuators [15], surgical devices [16], and metamaterials [17]. ∗

Corresponding author. E-mail addresses: [email protected] (Z. Yang), [email protected] (L.L. Howell).

https://doi.org/10.1016/j.mechmachtheory.2019.103700 0094-114X/© 2019 Elsevier Ltd. All rights reserved.

2

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

initial deformed beam1

beam-to-beam contact

beam-to-rigid contact

rigid object beam2 Fig. 1. A beam-based CCM undergoing beam-to-beam and beam-to-rigid contacts.

However, it is hard to analyze and design CCMs because of their nonlinear deformations and the strong boundary nonlinearities caused by contact. In the field of computational contact mechanics, a lot of numerical methods have been developed to solve the contact problem, e.g., the matrix inversion method [18], variational inequality method [19–21] and finite element analysis (FEA) [22,23]. Paul and Hashemi [18] used the matrix inversion method to determine the normal contact pressure between railroad wheel and rail. The variational inequality method can determine the shape and size at contact by using well-developed optimization techniques. The two methods are based on the elastic half-space model so that the linear elasticity theory holds [24]. FEA has been the most popular method for large-deflection contact analysis, including static and dynamics problems [25,26]. FEA has also been used as an analysis tool for topology optimization [27] of CCMs [28–30]. Consequently, CCMs are usually analyzed and verified by the well-established FEA. Not only the large deformation of CCMs, but also the contact force/stress in contact zone can be obtained by the contact FEA. Although FEA is a powerful and general tool for analyzing the CCMs with given shape and dimension, efficient and focused analysis models of CCMs are still needed for the following reasons: (1) The number of analyses can be large in the synthesis process of CCMs, in which the shape and dimension of CCMs usually need to be changed or optimized. (2) FEA is general and applicable to a wide range of problems, but the more focused models can be better suited for being embed in shape or dimensional optimization programs. (3) Unlike general contact analysis, the analysis of CCMs is more about the contact’s impact on the nonlinear deformation and the required load of mechanisms rather than the contact force/stress in the contact zone. Consequently, a simple and easy-to-follow analysis method will be useful for the nonlinear analysis and design of CCMs, especially for the user without much knowledge of contact FEA. Alternatives to FEA are found in the literature, but they focus on beam-to-rigid contacts. For example, Cannon and Howell [7] modeled the required moment of a contact beam in a contact-aided revolute joint based on the relationship between the beam’s final curvature and moment. Moon [31] also used the Bernoulli-Euler beam theory to model the energy storage and required torque of a compliant contact-aided four-bar mechanism. Mehta et al. [32] proposed an analytical model for cellular structures by considering the contact effects. By assuming that only one contact point exists, Yin and Lee [33] proposed a numerical method based on elliptical integrals to analyze a large-deflection finger contacting an elliptic object. Lan and Lee [24] presented a contact model based on nonlinear constrained minimization to analyze deformation and contact forces of the compliant fingers with variable thickness. A six degree-of-freedom pseudo-rigid-body model (PRBM) was developed by Venkiteswaran et al. [34] to estimate the tip deflection and tip contact force of a 3D cantilever beam. Bilancia et al. [35] extended the chained beam-constraint-model [36] to the analysis of contact-aided cross-axis flexural pivot. Based on principal axes decomposition of compliance matrices, Chen et al. [37] proposed an approach for large-deflection analysis of spatial flexible rods and applied it to a point contact problem. Recently, we proposed a general method [38] for the nonlinear analysis of beam-based CCMs undergoing beam-to-rigid contacts. The chained PRBM (CPRBM) [39–42] was used to model the large deformation of general beams. The whole contact analysis problem was formulated as finding the minimum potential energy of the compliant system under the geometry constraints from beam-to-rigid contact. The beam-to-rigid contact constraint is to keep the joints of CPRBM on one side of a rigid boundary curve. Because beam-to-beam contact has to consider the interactions between two deformable bodies, where the boundary curves become movable, the large-deflection analysis of beam-to-beam contact is more complex than that of beam-to-rigid contact [4]. Motivated by the need to analyze and design CCMs with beam-to-beam contact, this paper extends the proposed CPRBM-based method to the CCMs undergoing beam-to-beam contacts. For the beam-to-beam contacts, the offset of beams’ thickness should be taken into account to form one CPRBM’s moving boundary curve which will be used to constrain the motion of the other beam’s CPRBM. However, since the configuration of CPRBM changes in every iteration, it will be time consuming to calculate the offset by using polynomial fitting and derivation. To solve the problem, an approximation based on the geometric information of CPRBM is proposed to rapidly determine the moving boundary curve for beam-to-beam contact constraint. The static equilibrium configuration of CCMs is solved by minimizing its potential energy function under the geometric constraints from the boundary curves of beam-to-beam and beam-to-rigid contacts. Based on the deformation of the torsional springs in the CPRBM, a formulation is also provided in this paper to evaluate the normal stress along the

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

3

Fig. 2. Modeling schematic of CPRBM: (a) straight beam and (b) curved beam.

deformed beam. Numerical studies and FEA comparisons are used to verify the feasibility and accuracy of the proposed method. One contribution of this paper is that we establish a technique for the beam-to-beam contact constraint and provide the formulation for normal stress evaluation. Based on that, we propose a simple and straightforward method to analyze the large deflection and stress of the CCMs considering beam-to-beam contacts. The rest of this paper is organized as follows. Section 2 describes the proposed method. Section 3 gives several numerical studies to verify the proposed method. The discussions and conclusions are presented in Sections 4 and 5. 2. A general method for analyzing beam-based CCMs 2.1. Problem formulation of contact analysis In this paper, we focus on the large-deflection analysis of general beams with frictionless contact and uniform cross section. A beam-based CCM that undergoes beam-to-beam and beam-to-rigid contacts is shown in Fig. 1. The flexible beams of the CCM are in general shapes, and will contact each other or rigid objects during their deformations. Based on the principle of minimum potential energy and our previous work, the static equilibrium problem of the CCMs can be formulated as

min sub ject to

PE c (θ ) ≤ 0 Z (θ ) = 0 cb (θ ) = 0 or ≤ 0

(1)

where PE is the potential energy function of the compliant system, c(θ ) is the vector that contains the geometric constraints from beam-to-beam or beam-to-rigid contact, Z(θ ) is the vector loop equation(s) of the mechanism, cb (θ ) is the vector that contains the constraint of other boundary conditions, and θ is the variable vector describing tangent angle of the links in CPRBM. By minimizing the potential energy function PE under these constraints, we can obtain the static equilibrium configuration of CCMs. 2.2. Modeling the large deformation of general beams 2.2.1. The CPRBM As shown in Fig. 2, a general beam is discretized into N small-length elements. Each element is modeled by a 1R PRBM, i.e., a torsional spring is placed at the middle point of the element. The position of each torsional spring can be determined approximately as

xsi =

xei−1 + xei , 2

ysi = f (xsi ),

i = 1, 2, · · · , N

(2)

where f(x) is the function describing the beam’s shape and xei is the coordinate of the ith element node. The stiffness constant ki of each torsional spring is approximately calculated as

ki =

EI Lei

i = 1, 2, · · · , N

(3)

4

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

where E is Young’s modulus, I is the moment of inertia, and Lei is the length of the ith element. Then, the general beam can be modeled as a series of rigid links connected end to end through these torsional springs, i.e., CPRBM. When the number of elements N is large enough, the resulting CPRBM can equivalently describe the nonlinear load-deflection characteristics of a beam with general geometry shapes and load boundary conditions. Not only the beams’ end position, but also their deformed shapes can be accurately described. 2.2.2. The potential energy of CPRBM For force-loaded cases, the potential energy function PE is the sum of torsional springs’ potential energy and the conservative external forces’ potential energy PEforce :

PE =

N  1 2 ki δ − P E f orce 2 i

(4)

i=1







where δ i is the deformation of the ith torsional spring and is calculated as δi = (θi − θi−1 ) − (θi − θi−1 ). θ i and θi are the final and initial tangent angles of the ith link, respectively. For the displacement-loaded cases, the potential energy function PE only contains the elastic potential energy of torsional springs.

PE =

N  1 2 ki δ 2 i

(5)

i=1

2.2.3. Position and stress evaluation The current position of each joint or endpoint of the beam in the deformed configuration is

xi =

i−1 

L j cos(θ j ),

j=0

yi =

i−1 

L j sin(θ j )

(6)

j=0

where Lj is the length of the jth rigid link. The normal stress of each element can be approximately evaluated by using the angular deformation of a corresponding torsional spring [43]. The maximum normal stress caused by bending in the ith element is

       Mit   ki δit   Et δi  σi =   =  2I  =  2Le  2I

(7)

i

where Mi is the moment of the ith torsional spring and t is the element thickness. 2.3. Geometric constraints of contacts 2.3.1. Geometric constraints of beam-to-beam contact The two general beams shown in Fig. 3(a) are taken as an example to illustrate the constraints of beam-to-beam contact. When the thickness of beams are considered, the contact pair is the lower edge of beam1 and the upper edge of beam2. In other words, the two edges should not intersect with each other during the mechanism’s deformation. Instead of constraining the two edges directly, we use a moving boundary curve gs (x) to constrain the neutral line of beam1, where gs (x) is an equidistant offset of the neutral line y(x) of beam2. The reason for this is that the deformed neutral line of beams are directly described by the current position of their corresponding CPRBM’s joints. The distance of the equidistant offset is (t1 + t2 )/2, where t1 and t2 are the thicknesses of beam1 and beam2, respectively. To create the moving boundary curve gs (x), we should determine the position of corresponding offset points bi for related joints p2i first. One approach is to calculate normal vector ni of the neutral line y(x) directly by using polynomial fitting and derivation. But this is time consuming, because the position of the related joints p2i change in every iteration. As shown in Fig. 3(b), we use an approximation to calculate current position of the offset points bi . The key of this approximation is to determine the tangent angle of the normal vector ni based on the geometric information of CPRBM, i.e., the joints’ current position p2i and the tangent angle θ i of links. At joint p2i , the normal vectors of the (i − 1 )th and ith links are denoted by v1 and v2 correspondingly. The angle varied from v1 to v2 is θi − θi−1 . Theoretically, the angle β varied from v1 to ni is

β=

γi−1 (θ − θi−1 ) γi−1 + γi i

(8)

where γi−1 and γ i are angular variation of the links’ corresponding curved segments in the neutral line y(x). However, γi−1 and γ i cannot be obtained directly from CPRBM. When the length of rigid links are small enough, γi−1 is considered to be equal to γ i . Approximately, the angle β can be calculated as

1 2

β ≈ (θi − θi−1 )

(9)

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

5

Fig. 3. Schematic of the contact constraints on the beams’ deformed configurations: (a) beam-to-beam contact, (b) schematic for determining the moving boundary curve, and (c) beam-to-rigid contact.

Then, the coordinates of the offset points bi are





1 π (t1 + t2 ) cos θi−1 ± + β 2 2   1 π ybi = y2i + (t1 + t2 ) sin θi−1 ± +β , 2 2 xbi = x2i +

i ∈ [m2L , m2R ]

(10)

where x2i and y2i are the coordinates of p2i . The sign of π /2 is determined by the angular variation from the vector of the (i − 1 )th link to v1 . To reduce the number of constraints, only the related joints that will potentially contact the other beam are considered for creating the moving boundary curve gs (x). The range of the related joints is [m2L , m2R ]. Based on the offset points bi , the moving boundary curve gs (x) can be determined using polynomial fitting. Then, the constraint of beam-to-beam contact is to keep the related joints p1j of beam1 at one side of gs (x). For the case shown in Fig. 3(a), the related joints p1j of beam1 should be above the moving boundary curve gs (x). Its beam-to-beam contact constraints are formulated as

c ( θ ) = gs ( x1 j ) − y1 j ,

j ∈ [m1L , m1R ]

(11)

where [m1L , m1R ] indicates the region of related joints p1j that will potentially contact the moving boundary curve gs (x). 2.3.2. Geometric constraints of beam-to-rigid contact For the beam-to-rigid contact as shown in Fig. 3(c), the contact object is considered to be rigid, whose boundary curve gm (x) remains unchanged during the mechanism’s deformation. That will make it easier to implement the constraint of beam-to-rigid contact. A fixed boundary curve gm (x) can be obtained based on the equidistant offset of the contact edge y(x). Then, the position of the related joints p1j in beam1 are constrained as Eq. (11). 2.4. Calculation of the reaction force at the input port for displacement-loaded cases The reaction force calculation technique proposed in our previous work is briefly described in this section. This technique is based on the principle of work and energy to calculate the reaction force at the input port of CCMs. Fig. 4 shows the nonlinear relationship between the input displacement din and reaction force Fre for a compliant system. To calculate Fre for the input displacement din , two input displacements d1 and d2 are selected on the left side and right side of din , respectively.

d1 = din (1 − ρ ),

d2 = din (1 + ρ )

(12)

The step between d1 /d2 and din is ρ din , where ρ is the step ratio. When ρ is small enough, the load-displacement relationship within [d1 , d2 ] can be regarded as linear. If the input displacement increases from d1 to d2 , the work increment W can be calculated as

W = Fre (d2 − d1 ) = Fre (2ρ din )

(13)

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

reaction force

6

F re ρd in ρd in d1

d in

d2

displacement

Fig. 4. Schematic for calculating the reaction force Fre of the displacement-loaded cases.

t=10

t=10

t=5

t=5

Fig. 5. Verification of the accuracy of calculated offset points bi .

W is transformed into potential energy of springs, which is equal to the potential energy difference between d1 and d2 :

W = PE = PE (d2 ) − PE (d1 )

(14)

The potential energy PE(d1 ) and PE(d2 ) can be obtained after determining the mechanism’s equilibrium configurations caused by d1 and d2 , respectively. Then, the reaction force Fre is calculated as

Fre =

P E ( d2 ) − P E ( d1 ) 2ρ din

(15)

3. Numerical studies In this section, the accuracy of the offset points calculation will first be verified. Then, four analysis examples will be used to illustrate the proposed method. The constrained minimization problems of all the examples are solved by using the fmincon function in MATLAB. 3.1. Verification of the offset points’ accuracy Since the position of offset points bi are calculated approximately, their accuracy should be verified first. Take the parabolic curve y = −0.06x2 + 10 (x ∈ [−10, 10]) as an example. Different numbers of CPRBM’s links are used in the discretization to determine the offset points bi for different distances of equidistant offset t. As shown in Fig. 5, the variable x of discretization varies sinusoidally within [−10, 10]. The angular variations γ i of the segments in the initial curve are different, which may reduce the accuracy of the calculated β . As a result, there are some deviations between the calculated offset points bi and the accurate bi , e.g., the case of N = 5 in Fig. 5. As the equidistant offset t increases, the deviations become larger. Fortunately, the calculated offset points bi are nearly on the accurate offset curves shown in Fig. 5. Since we determine the offset curve by using polynomial fitting, the fitting offset curve based on the calculated bi will be quite accurate. In other words, the deviations between the calculated bi and actual bi have little impact on the accuracy of the offset curve fitting. Moreover, when the number of CPRBM’s links increases, the deviations between the calculated bi and the actual bi can be also reduced, e.g., the case of N = 10 in Fig. 5. To illustrate the deviations between the calculated bi and the actual offset curve, the relative errors describing their distances along the y-axis are calculated for different offset distances t (t ∈ [1, 20]) and link numbers N (N ∈ [6, 30]). The

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

7

Fig. 6. Maximum relative errors between the calculated offset points bi and the offset curve for different offset distances t and link numbers N.

Fig. 7. A rotational hinge based on two semicircular beams.

maximum value among the relative errors of each case is plotted in Fig. 6. The relative errors are all below 0.04%. When N is less than 10, the relative error increases as the offset distances t increases. As the link numbers N increases, the offset distances t have little impact on the relative errors which are all below 0.01%. The result shows that the approximation in Eq. (9) can accurately calculate the angle β . This implies the accuracy of following offset points bi and the fitting offset curve. 3.2. A rotational flexible hinge with beam-to-beam contact A rotational hinge based on two semicircular beams is shown in Fig. 7. The two semicircular beams, beam1 and beam2, are symmetric with respect to the x-axis. The parameters of the two beams are: E = 2.3 × 109 Pa, thickness t = 2 mm, and width w = 8 mm. The shape of beam1 is defined as

(x − 26 )2 + (y − 26.25 )2 = 252 , mm

(16)

In the initial configuration, there is a gap between beam1 and beam2. The lower edge of beam1 will potentially contact the upper edge of beam2 during their deformations. The other two parts of the hinge, links AC and BD, are considered to be rigid. A pure moment M is applied to link BD. The two beams are modeled by CPRBM separately. Each beam is discretized uniformly by 40 elements. The potential energy and joint position of the two beams are formulated separately. The potential energy of the system is the summation of the two beams’ potential energy and the applied moment’s potential energy. The current positions of endpoints B and D are used to form the following vector loop equations:

Z = CA + AB − (CD + DB )

(17)

Since the rigid connection, the included angles between the last link of the two CPRBMs and link BD should stay constant during deformation. These are included in cb (θ ). The joints whose initial x coordinates are in the range of [19.75, 32.25] mm are considered as the related joints for beam-to-beam contact constraints. Unlike beam-to-rigid contact, the moving boundary curve gs (x) of beam-to-beam contact changes in every iteration of optimization. To implement beam-to-beam contact constraint, the formulations of the related joints’ position are included in the nonlinear constraint nonlcon for the function

8

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

Fig. 8. Comparison of the rotational hinge’s deformations obtained by CPRBM and nonlinear FEA. Link BD is not plotted.

Fig. 9. Nonlinear FEA results for the rotational hinge: (a) deformation and (b) maximum bending stress. Link BD is hidden.

fmincon. Based on the current position of the related joints in each iteration, the moving boundary curve gs (x) can be created to formulate the beam-to-beam contact constraints. In this case, the applied moment M is equal to 0.67 N · m. The deformed configurations of the two beams obtained by the proposed method are shown in Fig. 8, and compared with the nonlinear FEA results from ANSYS, whose corresponding graphic results are shown in Fig. 9. The figure shows that the two deformed CPRBMs are nearly on the two beams’ deformed configurations described by FEA. By comparing the results from the proposed method and FEA, the relative errors of points B and D’s displacements are 0.62% and 0.81%, respectively. The deformed configuration of the rotational hinge is accurately obtained by the proposed method. More details about the offset points bi and moving boundary curve gs (x) are also provided in Fig. 8. The number of the related joints is eight for each beam. The related joints in beam2 are utilized to determine the moving boundary curve gs (x). Thanks to the geometric constraints of beam-to-beam contact, the eight related joints in beam1 are kept above or on gs (x). Two related joints of beam1 are on gs (x), which indicates the contact region between the two beams. Based on the deformed configuration, the maximum normal stress caused by bending in every element is calculated. Then, the maximum normal stress along beam1 and beam2 can be obtained. Fig. 10 provides the stress along beam1 and beam2 with respect to their x coordinates, respectively. The calculated stress is compared with the combined stress and bending stress obtained

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

9

Fig. 10. Comparison of the rotational hinge’s stress obtained by CPRBM and nonlinear FEA. Left: the stress along beam1 with respect to its x coordinates. Right: the stress along beam2 with respect to its x coordinates.

beam-to-beam contact

A B

y

2

1

am

am

beam-to-rigid contact

g m( x )

be

be

rigid O

x

C

Fig. 11. The upper part of a compliant gripper undergoing two kinds of contacts.

from ANSYS. One can see that the calculated stresses from CPRBM match well with the bending stress of FEA. The calculated stresses are also close to the combined stress result that includes other kinds of stress, e.g., shear stress. 3.3. A compliant gripper undergoing two kinds of contacts A compliant gripper that is symmetric with respect to the x-axis is taken as the second analysis example. Due to symmetry, we only consider the upper part of this compliant gripper. As shown in Fig. 11, the analysis model consists of two flexible beams, i.e., beam1 and beam2. The parameters of the two beams are: E = 2.3 × 109 Pa, thickness t = 2 × 10−3 m, and width w = 8 × 10−3 m. beam1 (OA) is a straight beam defined by

√

√ 3 2 3 y= x+ t, 3 3



 √ t t 3 x∈ − , − 2 20 2

(18)

The motion of endpoint O of beam1 is constrained in its translational freedom along the y-axis and rotational freedom, while endpoint C of beam2 (BC) is fixed. The two beams are connected to each other by a solid block AB, which is considered to be rigid. The lower edge of beam1 will potentially contact the upper edge of beam2 during their deformations, i.e., beamto-beam contact. The shape of beam2 is a parabola defined by x = 14.64y2 + 0.05 (y ∈ [0, 0.05] m ). For the beam-to-rigid contact, the lower edge of beam2 will potentially contact a rigid and fixed object, whose boundary curve gm (x) for contact is a semi-circle defined as

gm ( x ) =



0.022 − (x − 0.07 )2

(19)

In CPRBM modeling, beam1 is discretized uniformly by 30 elements, while beam2 is discretized uniformly along the yaxis by 40 elements. The current positions of endpoints O, A, and B are used to form the following vector loop equations:

Z = CO + OA − (CB + BA )

(20)

In addition, the boundary constraints from the rigid connection between beam1 and beam2 should be considered too. For each beam, the angle between the last link of its CPRBM and link AB should stay constant during deformation. The joints of beam1 whose initial y coordinates are in the range of [0.03, 0.051] are considered as the related joints for beam-to-beam

10

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

(a)

(b) bi g s( x )

Fig. 12. Comparison of the compliant gripper’s deformations obtained by CPRBM and nonlinear FEA. (a) The whole mechanism and (b) the enlarged view of contact. Link AB is not plotted.

Fig. 13. Nonlinear FEA results for the rotational hinge: (a) deformation and (b) maximum bending stress. Link AB is hidden.

contact constraints, while the joints of beam2 whose initial y coordinates are in the range of [0.025, 0.05] are considered as the related joints for beam-to-beam contact constraints. For this case, the moving boundary curve gs (x) for beam-to-beam contact is created from beam1. The first two thirds of the joints in beam2 are considered to form the beam-to-rigid contact constraints. The mechanism is first loaded by a vertical force FA = 16 N at endpoint A. The deformed configurations of the two beams obtained by the proposed method are shown in Fig. 12(a), and compared with the nonlinear FEA results from ANSYS, whose corresponding graphic results are shown in Fig. 13. The two deformed CPRBMs in the static equilibrium configuration are nearly on the lines described by the FEA result. By comparing the results from CPRBM and FEA, the relative errors of points A and B’s displacements are 0.17% and 0.17%, respectively. The deformed configuration of this compliant gripper is obtained accurately by the proposed method. More details about the contacts are provided in the enlarged view shown in Fig. 12(b). The related joints in beam1 are utilized to determine the offset points bi and moving boundary curve gs (x) for the beam-to-beam contact constraints, in

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

11

Fig. 14. Comparison of the compliant gripper’s stress obtained by CPRBM and nonlinear FEA. Left: the stress along beam1 with respect to its y coordinates. Right: the stress along beam2 with respect to its y coordinates.

Fig. 15. Comparison of the reaction force obtained by CPRBM and nonlinear FEA.

which the related joints in beam2 are constrained to be below or on gs (x). The beam-to-beam contact region is indicated by the related joints in beam2 whose positions are on gs (x). For the beam-to-rigid contact, part of the related joints of beam2 in the deformed configuration are on the boundary curve gm (x). The figures show that the region of contacts indicated by the proposed method match well with the FEA results. Fig. 14 provides the stress along beam1 and beam2 with respect to their y coordinates, respectively. The calculated stress is compared with the combined stress and bending stress obtained from ANSYS. For this case, the combined stress is so close to the bending stress that the two set of data are overlapped in the figure. One can see that the calculated stresses obtained by the proposed method match well with the FEA results. Next, the mechanism is loaded by a vertical displacement y at endpoint A. Its direction is the same as shown in Fig. 11. By using the technique proposed in our previous work, we can obtain the reaction force for the input displacement y whose value varies from 1 mm to 40 mm. For this case, we set the step ratio to ρ = 0.2%. The calculated reaction forces are plotted in Fig. 15, and compared with the result from ANSYS. The relative errors between the two results are below 0.5%. Thus, the reaction force of the CCMs with beam-to-beam and beam-to-rigid contacts can be evaluated accurately by the technique. 3.4. A CCM with prestressed assembly To show the feasibility of analyzing CCMs with a prestressed assembly, a partially compliant mechanisms with beamto-beam contact will be analyzed in this section. The mechanism consists of one rigid crank ABC and two identical beams. The parameters of the two beams are: E = 2 × 1011 Pa, thickness t = 1 mm, and width w = 8 mm. Fig. 16(a) shows the initial state of the parts before assembly, including the rigid crank ABC and two straight beams, i.e., beam1 and beam2. The crank ABC is able to rotate around point A, whose parameters are OA = 10 mm, AD = 40 mm, and BD = DC = 20 mm. The lengths of beam1 and beam2 are 60 mm. In the assembly, one end of beam1 and beam2 are vertically fixed to the points (−t/2, 0 ) and (t/2, 0), while the other end of beam1 and beam2 are connected to the crank ABC by the rotational joints at points B and C, respectively. Since the maximum distance of OB or OC is less than the beams’ length, the two beams will be prestressed after assembly. To make the two beams contact each other, the maximum angle between BC and the end of beam1 (or beam2) is limited to 45◦ . Fig. 16(b) shows the equilibrium state of the mechanism after prestressed assembly, i.e., the two beams are symmetric with respect to the y-axis and α = 0◦ .

12

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

pinned

pinned

B

C

C

D beam1

beam2

beam2

beam1

D

45°

y

B

y

α

beam-to-beam contact A

A

f ix e d to

fix

ed

to O

x

O (a)

x

(b)

Fig. 16. A partially compliant mechanism with prestressed assembly: (a) initial state before prestressed assembly and (b) its equilibrium state after prestressed assembly.

C

C

B

B

A

A

O

O

C

C

A B A O O

B

Fig. 17. Four deformed configurations of the partially compliant mechanism with prestressed assembly.

This mechanism is displacement loaded at the crank ABC, and the input displacement is described by α . In CPRBM modeling, both beam1 and beam2 are discretized uniformly by 40 elements. To form the vector loop equations, the current endpoint of beam1 and beam2 should be in the positions of point B and C, respectively. Two additional boundary conditions are the constraints on the angle between BC and the last link of the two beams’ corresponding CPRBM. The joints of the two beams whose initial y coordinates are in the range of [0, 25] mm are considered as the related joints for the beam-to-beam contact constraints. The moving boundary curve gs (x) is formed based on the current configuration of beam2. For this case, the beam-to-beam contact constraints are formulated along the direction of x-axis, i.e., c (θ ) = x1 j − gs (y1 j ). By driving the crank ABC from α = 0◦ to α = 90◦ , the proposed method can obtain the final equilibrium configurations of the mechanism. Four of the deformed configurations are shown in Fig. 17. When α is less than 35◦ , the deformed configurations of the mechanism are similar to the case of α = 35◦ , in which large part of the two beams contact each other.

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

13

35

36

Fig. 18. The reaction moment and potential energy vs. the input displacement α of the crank ABC for the partially compliant mechanism with prestressed assembly.

B G

y beam-to-rigid contacts

beam1 beam2 beam-to-beam contact

A

C

α1 O

α2 F

D

x

Fig. 19. A two-DOF CCM based on five-bar mechanism and circular beams.

In the figure of α = 35◦ , most of the beam2’s offset points bi are on the neutral line of beam1. If α increases from 35◦ to 36◦ , the bistable characteristics make beam2 change abruptly to another configuration shown in Fig. 17. As a result, only a few parts of the two beams contact each other. The deformed configurations of the other cases are similar to the case of α = 36◦ . With the reaction force calculation technique, the reaction moment at the input can be determined for the deformed configurations. The relationship between the reaction moments and the input displacement α is shown in Fig. 18. In addition, it is easy for the proposed method to calculate the potential energy of the deformed configurations. Thus, the relationship between the potential energy and the input displacement α is also shown in the figure. The bistable characteristic of the mechanism is illustrated. The reaction moment and potential energy experience abrupt changes when α varies from 35◦ to 36◦ . 3.5. A two-DOF CCM The last example is to show how the proposed method can be used to model CCMs with more complex kinematics, e.g., the two-DOF CCM shown in Fig. 19. This CCM is a rigid-body five-bar mechanism coupled with two circular beams, i.e., beam1 and beam2. The parameters of the five-bar mechanism are OA=CD=50 mm and AB=BC=130.39 mm. The angular displacement of link OA and link CD are regarded as inputs of this five-bar mechanism. The initial tangent angles of links OA and CD are α1 = 120◦ and α2 = 60◦ , respectively. Beam1 is a circular beam coupled with link AB, whose center and radius are G(-25, 123) mm and R1 = 79.7 mm. The two ends of beam1 are connected to joint A and joint B, respectively. The center and radius of beam2 are F(70, 0) mm and R2 = 70 mm, whose two ends are connected to joint O and joint C, respectively. Other parameters of the two flexible beams are: E = 2 × 1011 Pa, thickness t = 1 mm, and width w = 8 mm. Both beam1 and beam2 are discretized uniformly by 40 elements for CPRBM modeling. The first link of the CPRBMs for beam1 and beam2 are connected to joint A and joint O, respectively. In addition to the variables of the two CPRBMs, the tangent angles of link AB and link BC should be considered as variables too. The objective function is the summation

14

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

B

B

C

C

A

A D

O

D

O

B

B

C C

A

A O

D

O

D

Fig. 20. Four deformed configurations of the two-DOF CCM. The dash lines represent the initial configuration of the mechanism, where the initial values of α 1 and α 2 are 120◦ and 60◦ , respectively.

of the two beams’ potential energy. The vector loop equations can be formulated according to the following geometric relationships.

OA + ABlink = OD + DC + CB = OA + ABbeam1 OD + DC = OCbeam2

(21)

For simplicity, three contact pairs are defined without considering the beam thickness, i.e., the beams are considered as lines. (1) The beam-to-beam contact between beam1 and beam2. The moving boundary curve is created from beam1 based on joints p1j , j ∈ [8, 26]. The joints of beam2, p2i i ∈ [11, 28], are subjected to the geometric constraints from the moving boundary curve. For this case, the beam-to-beam contact constraints are formulated along the direction of x-axis, i.e., c (θ ) = gs (y2i ) − x2i . The reasons for this will be explained at the end of this section. (2) The two beam-to-rigid contacts between beam1/beam2 and link BC. The last ten joints of the two beams are subjected to the geometric constraints from the moving link BC, which is defined by its tangent angle and the position of joint C. The two sets of related joints should not be above the straight line defined by link BC. The contact constraints can be formulated along the direction of either x-axis or y-axis. With given values of α 1 and α 2 , the static equilibrium configuration of the CCM can be obtained. Fig. 20 shows four deformed configurations of the two-DOF CCM, where α 1 is equal to 140◦ and α 2 varies from 90◦ to 130◦ . As link DC rotates anticlockwise, the deformation of the two beams increase and the contact zone of the beam-to-beam contact pair relocates along the beams. Meanwhile, the contact portion of the two beam-to-rigid contact pairs increase and deform into the shape of link BC. One can see that the three contact pairs’ impact on the two-DOF mechanism’s deformed configuration can be accurately evaluated by the proposed method. In the case of α1 = 140◦ and α2 = 130◦ , beam2 and link OA are overlapped because no contact constraint is applied to them. Since the relative movement between beam1 and beam2 is larger than that of other three examples, the constrained joints of beam-to-beam contact will possibly move out of the range of the moving boundary curve gs (x) in the direction of x-axis or y-axis, which will result in infeasible solutions. To avoid this, the beam-to-beam contact constraints should be formulated carefully for this case. If beam2 is regarded as the moving boundary curve, the constrained joints p1j close to joint A will possibly move out of the range of gs (x) in the direction of y-axis, e.g., the case of α1 = 140◦ and α2 = 90◦ . In addition, the constrained joints p1j close to joint B will possibly move out of the range of gs (x) in the direction of x-axis. If gs (x) is created from beam1, the constrained joints p2i close to joint C will possibly move out of the range of gs (x) in the direction of y-axis, while this problem can be avoided in the direction of x-axis.

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

15

Table 1 Comparison of the computation time of the proposed method and FEA. Methods

Rotational hinge

Compliant gripper

Contact FEA Contact FEA Non-contact FEA CPRBM CPRBM

4 s (80 elements/beam) 5 s (160 elements/beam) 4 s (80 elements/beam) 5 s (N1 = N2 = 40) 11.75 s (N1 = N2 = 60)

23 s (80 elements/beam) 40 s (160 elements/beam) 4 s (80 elements/beam) 4.32 s (N1 = 40, N2 = 30) 19.56 s (N1 = 60, N2 = 50)

4. Discussions All the numerical examples and FEA simulations in this paper are carried out on a computer with Intel Core i5 - 8500 (3.00 GHz) CPU, 8.00GB RAM, MATLAB R2018b, and ANSYS Workbench 19.1. Table 1 lists the computation time of the proposed method and FEA. The FEA results shown in Sections 3.2 and 3.3 are solved by discretizing each beam of the rotational hinge and compliant gripper with 80 beam elements, whose corresponding computational time are 4 s and 23 s, respectively. When the number of elements in FEA discretization is double, i.e., 160 elements per beam, the computational time for the two cases are 5 s and 40 s. For the rotational hinge case, the computation time of the proposed method is 5 s, which is larger than that of FEA. The reason for this is that the contact zone of the rotational hinge is so small that only a few contact constraints are activated in FEA due to its contact detection algorithm, whereas the proposed method considers all the contact constraints during optimization. Thus, the contact FEA simulation of the rotational hinge has the same efficiency as its FEA simulation without considering contact effects. For the case with more obvious contact effect, e.g., the compliant gripper, the contact effect will greatly increase the computation time of FEA. One can see that the contact FEA spends more time than the proposed method and non-contact FEA. When more elements are used for CPRBM modeling of the rotational hinge and compliant gripper, the number of variables and contact constraints will increase, which will increase the computation time of the proposed method. For example, the computation time of the two cases will increase to 11.75 s and 19.56 s if the CPRBM discretizations are set to the data shown in the last row of Table 1. The computation time of the proposed method can be reduced by discretely selecting the related joints for contact constraints, i.e., without considering all the joints in the contact zone. To show the material property’s impact on the accuracy of the proposed method, three different materials (ABS, titanium, and steel) are used for the rotational hinge and compliant gripper cases. For the three materials, the relative errors of the rotational hinge are all below 0.8%, while the relative errors of the compliant gripper are all below 0.17%. It should be pointed out that the beam-to-beam contact analysis in ANSYS is not as easy as that of surface-to-surface or beam-to-surface contacts. One of the key settings is that the contact radius of the beams must be set to half of the beam’s thickness because the contact beam surface in ANSYS is assumed to be the surface of an equivalent cylinder. The users have to carefully set the real constant of beams in ANSYS APDL or use additional command code in Workbench. For the CCM with prestressed assembly, the analysis procedure of FEA becomes more complicated. One feasible method is to use multiple load steps and two auxiliary moments. The key modeling process is as follows: (1) Model the mechanism without prestressed assembly and let beam1 and beam2 be in the vertical direction. Define the revolute joints at points B and C and let the other end of the two beams be freed. (2) In the first load step, the rotational freedoms of the two beam’s freed end are fixed, and their positions are constrained to (−t/2, 0 ) and (t/2, 0), respectively. Meanwhile, α is set to zero. To avoid buckling at the two beams and control the deformed configuration, two auxiliary moments act on the joints B and C, respectively. (3) In the second load step, set the auxiliary moments to zero and keep α = 0. The beams will deform to the static equilibrium configuration. (4) In the third load step, add the rotational displacement on crank ABC. The whole process in ANSYS requires careful modeling, or it will be easy to suffer from convergency problems. Compared with the above analysis process in ANSYS, the proposed method can solve the problem using the same procedures as the cases without prestressed assembly. The proposed method regards the crank as a rigid body directly, whereas the contact analysis of ANSYS often suffers from convergence problems when the stiffness property of a body is set to rigid. It regards the contact analysis of CCMs as finding the static equilibrium configuration of CPRBM under the geometric constraints from contacts’ boundary curves. The contact constraints are all formulated based on the geometric information of contact pairs, and no artificial spring or coefficient is used in the contact. The whole process including the modeling of CPRBM is straightforward and easy-to-follow. The method can be easily implemented by the user without much knowledge of non-linear analysis and contact FEA. In addition to the examples shown in this paper, the method can be applied to the problems that consider contacts between three or more planar beams by utilizing multiple beam-to-beam contact pairs. Furthermore, the framework of the proposed method can be extended to the case of contact beams with variable cross-section. One limitation of the proposed method is that the contact constrained joints should move within the range of boundary curves. Otherwise, the fmincon function cannot obtain the feasible solution. 5. Conclusions Based on our previous work on beam-to-rigid contact, a numerical method is proposed in this paper to analyze the large deflection and stress of the CCMs undergoing beam-to-beam and beam-to-rigid contacts. The proposed approximation can

16

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

rapidly and accurately determine the moving boundary curve for beam-to-beam contact constraints. The static equilibrium configuration of CCMs can be accurately obtained by minimizing the potential energy function of CPRBM under the geometric constraints from contacts. The normal stress along the deformed beam can be evaluated easily according to the torsional springs’ deformation in CPRBM. With an accuracy comparable to FEA, the proposed method is more straightforward in analyzing CCMs with large deflection, partial compliance, or prestressed assembly. It can be an alternative to FEA, especially in the conceptual design phase of CCMs and for the user without much knowledge of contact FEA. Declaration of Competing Interest We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted. Acknowledgments The authors gratefully acknowledge the financial support from the Young Innovative Talents Program of Universities in Guangdong (grant nos. 2018KQNCX018 and 2018KQNCX021), the CSC Scholarship, and the National Natural Science Foundation of China (grant no. 51675189). References [1] L.L. Howell, Compliant Mechanisms, John Wiley & Sons, Inc, New York, 2001. [2] N.D. Mankame, G.K. Ananthasuresh, Contact aided compliant mechanisms: concept and preliminaries, in: ASME 2002 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers, Montreal, Canada, 2002, pp. 109–121. [3] N.D. Mankame, G.K. Ananthasuresh, Topology optimization for synthesis of contact-aided compliant mechanisms using regularized contact modeling, Comput. Struct. 82 (15–16) (2004) 1267–1290. [4] P. Litewka, Finite Element Analysis of Beam-to-beam Contact, Springer Science & Business Media, Berlin Heidelberg, 2010. [5] P. Kumar, A. Saxena, R.A. Sauer, Computational synthesis of large deformation compliant mechanisms undergoing self and mutual contact, J. Mech. Des. 141 (1) (2019) 12302. [6] N.D. Mankame, G.K. Ananthasuresh, A novel compliant mechanism for converting reciprocating translation into enclosing curved paths, J. Mech. Des. 126 (4) (2004) 667–672. [7] J.R. Cannon, L.L. Howell, A compliant contact-aided revolute joint, Mech. Mach. Theory 40 (11) (2005) 1273–1293. [8] P.A. Halverson, L.L. Howell, A.E. Bowden, A flexure-based bi-axial contact-aided compliant mechanism for spinal arthroplasty, in: ASME 2008 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers, Brooklyn, USA, 2008, pp. 405–416. [9] J.R. Montierth, R.H. Todd, L.L. Howell, Analysis of elliptical rolling contact joints in compression, J. Mech. Des. 133 (3) (2011) 31001. [10] T.G. Nelson, R.J. Lang, S.P. Magleby, L.L. Howell, Curved-folding-inspired deployable compliant rolling-contact element (D-CORE), Mech. Mach. Theory 96 (2016) 225–238. [11] T.G. Nelson, J.L. Herder, Developable compliant-aided rolling-contact mechanisms, Mech. Mach. Theory 126 (2018) 225–242. [12] Y. Tummala, A. Wissa, M. Frecker, J.E. Hubbard, Design and optimization of a contact-aided compliant mechanism for passive bending, J. Mech. Robot. 6 (3) (2014) 31013. [13] J.P. Calogero, M.I. Frecker, Z. Hasnain, J.E. Hubbard Jr., A dynamic spar numerical model for passive shape change, Smart Mater. Struct. 25 (10) (2016) 104006. [14] J. Calogero, M. Frecker, Z. Hasnain, J.E. Hubbard, Tuning of a rigid-body dynamics model of a flapping wing structure with compliant joints, J. Mech. Robot. 10 (1) (2018) 11007. [15] Z. Song, S. Lan, J.S. Dai, A new mechanical design method of compliant actuators with non-linear stiffness with predefined deflection-torque profiles, Mech. Mach. Theory 133 (2019) 164–178. [16] K.W. Eastwood, P. Francis, H. Azimian, A. Swarup, T. Looi, J.M. Drake, H.E. Naguib, Design of a contact-aided compliant notched-tube joint for surgical manipulation in confined workspaces, J. Mech. Robot. 10 (1) (2018) 15001. [17] L.A. Shaw, S. Chizari, M. Dotson, Y. Song, J.B. Hopkins, Compliant rolling-contact architected materials for shape reconfigurability, Nat. Commun. 9 (1) (2018) 4594. [18] B. Paul, J. Hashemi, Contact pressures on closely conforming elastic bodies, J. Appl. Mech. 48 (3) (1981) 543–548. [19] G. Fichera, Problemi elastostatici con vincoli unilaterali: il problema di signorini con ambigue condizioni al contorno, Mem. Accad. Naz. Lincei 8 (7) (1964) 91. [20] G. Duvaut, J.L. Lions, Les inéquations en mécanique et en physique, Dunod, Paris, 1972. [21] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, PA, 1988. [22] P. Wriggers, Computational Contact Mechanics, John Wiley & Sons, Chichester, 2002. [23] D. Sun, C. Liu, H. Hu, Dynamic computation of 2d segment-to-segment frictionless contact for a flexible multibody system subject to large deformation, Mech Mach Theory 140 (2019) 350–376. [24] C.-C. Lan, K.-M. Lee, An analytical contact model for design of compliant fingers, J. Mech. Des. 130 (1) (2008) 11008. [25] Q. Wang, Q. Tian, H. Hu, Dynamic simulation of frictional contacts of thin beams during large overall motions via absolute nodal coordinate formulation, Nonlinear Dyn. 77 (4) (2014) 1411–1425. [26] Q. Wang, Q. Tian, H. Hu, Dynamic simulation of frictional multi-zone contacts of thin beams, Nonlinear Dyn. 83 (4) (2016) 1919–1937. [27] B. Zhu, X. Zhang, H. Zhang, J. Liang, H. Zang, H. Li, R. Wang, Design of compliant mechanisms using continuum topology optimization: a review, Mech. Mach. Theory 143 (2020) 103622. [28] B.V.S.N. Reddy, S.V. Naik, A. Saxena, Systematic synthesis of large displacement contact-aided monolithic compliant mechanisms, J. Mech. Des. 134 (1) (2012) 11007. [29] A. Saxena, A contact-aided compliant displacement-delimited gripper manipulator, J. Mech. Robot. 5 (4) (2013) 41005. [30] P. Kumar, R.A. Sauer, A. Saxena, Synthesis of C0 path-generating contact-aided compliant mechanisms using the material mask overlay method, J. Mech. Des. 138 (6) (2016) 62301. [31] Y.-M. Moon, Bio-mimetic design of finger mechanism with contact aided compliant mechanism, Mech. Mach. Theory 42 (5) (2007) 600–611. [32] V. Mehta, M. Frecker, G.A. Lesieutre, Stress relief in contact-aided compliant cellular mechanisms, J. Mech. Des. 131 (9) (2009) 91009. [33] X. Yin, K.-M. Lee, Modeling and analysis of grasping dynamics for high speed transfer of live birds, in: Second IFAC Conference on Mechatronic Systems, Berkeley, CA, 2002.

M. Jin, B. Zhu and J. Mo et al. / Mechanism and Machine Theory 145 (2020) 103700

17

[34] V.K. Venkiteswaran, J. Sikorski, S. Misra, Shape and contact force estimation of continuum manipulators using pseudo rigid body models, Mech. Mach. Theory 139 (2019) 34–45. [35] P. Bilancia, G. Berselli, S. Magleby, L. Howell, On the modeling of a contact-aided cross-axis flexural pivot, Mech. Mach. Theory 143 (2020) 103618. [36] G. Chen, F. Ma, G. Hao, W. Zhu, Modeling large deflections of initially curved beams in compliant mechanisms using chained beam constraint model, J. Mech. Robot. 11 (1) (2019) 11002. [37] G. Chen, Z. Zhang, H. Wang, A general approach to the large deflection problems of spatial flexible rods using principal axes decomposition of compliance matrices, J. Mech. Robot. 10 (3) (2018) 31012. [38] M. Jin, Z. Yang, C. Ynchausti, B. Zhu, x Zhang, L.L. Howell, Large deflection analysis of general beams in contact-aided compliant mechanisms using chained pseudo-rigid-body model, J. Mech. Robot. (2020), doi:10.1115/1.4045425. [39] L. Saggere, S. Kota, Synthesis of planar, compliant four-bar mechanisms for compliant-segment motion generation, J. Mech. Des. 123 (4) (2001) 535–541. [40] J. Pauly, A. Midha, Pseudo-rigid-body model chain algorithm, part 1: introduction and concept development, in: ASME 2006 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers, 2006, pp. 173–181. [41] J. Pauly, A. Midha, Pseudo-rigid-body model chain algorithm, part 2: equivalent representations for combined load boundary conditions, in: ASME 2006 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers, 2006, pp. 183–190. [42] R.P. Chase Jr, R.H. Todd, L.L. Howell, S.P. Magleby, A 3-d chain algorithm with pseudo-rigid-body model elements, Mech. Based Des. Struct. Mach. 39 (1) (2011) 142–156. [43] L.L. Howell, A. Midha, A method for the design of compliant mechanisms with small-length flexural pivots, J. Mech. Des. 116 (1) (1994) 280–290.