Compliant members come in a variety of shapes and sizes. While thin beam flexures are commonly used in this field, they can be replaced by soft members with lower aspect ratio. This paper looks to study the behavior of such elements by analyzing them from the view of beam theory for 2D cases. A modified version of the Timoshenko beam theory is presented which incorporates extension and Poisson's effects. The utility and validity of the new approach are demonstrated by comparing against Euler–Bernoulli beam theory, Timoshenko beam theory, and finite-element analysis (FEA). The results from this are then used to study the performance of pseudo-rigid-body models (PRBMs) for the analysis of low aspect ratio soft compliant joints for 2D quasi-static applications. A parallel-guiding mechanism comprised of similar compliant elements is analyzed using the new results to validate the contribution of this work.

## Introduction

The analysis of beam elements has many applications in mechanical engineering. Beam theory methods have evolved over the years in various forms. While the Euler–Bernoulli beam equations are still used for a multitude of static applications, more recent forms of beam theory include the variational-asymptotic method for finite-elements (VABS) [1], and geometrically exact beam theory [2]. FEA is an another widely used method for calculating deformations of the mechanical elements. In this paper, an approach will be presented to determine the deflection of soft joints in planar mechanisms from the view of beam theory. This method will incorporate shear and extension effects into the analysis of compliant mechanisms and will provide a more accurate theoretical formulation than those previously used for analysis in this field.

The groundwork for introducing shear and rotational inertia into beam equations was laid down by Timoshenko [3]. Over the years, many researchers have used and embellished it. Chaisomphob et al. [4] developed a finite displacement theory accounting for axial effects, and Goto et al. [5] derived elliptic integral solutions for the same. However, they did not consider the change in the cross-sectional area due to Poisson's effect. There have also been many papers dealing with the shear correction factor, and analyzing its accuracy for various cases [69]. In this paper, the results by Cowper [10] will be used since these were developed for static applications.

Compliant mechanisms have proved to be extremely useful for a variety of applications, from flapping wing micro-air vehicles [11,12] to shape morphing applications [13] and micromechanisms [14]. In many cases, the flexible members in compliant mechanisms exhibit beam-like behavior. Such compliant elements can be analyzed using different methods [15], such as the beam theory approach [16], the beam constraint model [17], or PRBMs [1820].

The PRBM approach has a number of advantages. It transforms compliant elements into their rigid-body counterparts, and allows the user to take advantage of the well-established equations used in rigid-body kinematics. It is a good numerical approximation requiring less computation when compared to some of the other approaches such as topology optimization [21]. While some of the basic groundwork for this approach was laid down by Howell and Midha [22,23], many researchers have worked on this approach over the past two decades and expanded the knowledge base significantly. Extension springs were used to incorporate general loading cases and also initially curved pinned–pinned beams [24,25]. More elements were added for special cases such as beams with points of inflection [26], and other improvements have been made to increase the accuracy [27,28].

### Extension Effects in Compliant Mechanisms.

Shape deposition manufacturing (SDM) is a step-by-step method for depositing layers of desired materials in varying shapes to produce required parts [2931]. Polymeric materials are usually used in the deposition process. The method allows multiple materials to be used and also allows embedding of other parts such as sensors and actuators. The process eliminates the need for seams, fasteners, or other forms of assembly, i.e., the assembly that occurs in situ during the deposition of material.

Two compliant parallel-guiding mechanisms fabricated by SDM are shown in Fig. 1(a). The mechanism on the left is fabricated using urethane plastic of flexural modulus 2.413 GPa. The one on the right uses urethane rubber for the beams with a 100% modulus of 2.068 MPa. As shown in Fig. 1(b), both the mechanisms have similar performance under the action of a horizontal load. Therefore, the same performance goals were achieved by the second design with a mechanism half the size of the first one. However, the soft beams of low aspect ratio have significant extension effects, which are not incorporated into the conventional methods of analysis in compliant mechanisms.

Fig. 1
Fig. 1
Close modal

To tackle this problem, this paper will focus on using the beam theory approach for the analysis of soft compliant elements. The formulation derived here will aim to provide a set of beam equations that can be readily solved while capturing most of the significant effects. Timoshenko beam theory assumes no sectional change. In this paper, it is assumed that the only sectional change is the Poisson's effect which is more significant than the effect of bending. The study will be restricted to two-dimensional quasi-static cases, but this is still applicable for a large variety of compliant mechanisms. The joints or beam elements are expected to undergo a change in length of less than 10%, but this can lead to significant error if one uses conventional beam theory or PRBMs for analysis. The second part of this paper will look at the effect of introducing extension elements in PRBMs, and study how they can reduce the error in predicting the behavior of these kinds of compliant elements.

In Sec. 2, the beam theory formulation will be derived. Some specific loading cases will be examined to illustrate the consequence of shear, extension, and Poisson's effects on the results. Conventional PRBMs will then be compared to the beam theory results, and the resulting increase in error will be studied. A set of three PRBMs with extension spring optimized for beam geometry will be presented, showing a significant decrease in error, followed by a discussion of the salient points of the results.

## The Modified Timoshenko Beam Formulation

In this section, the equations for the new beam theory approach will be derived. Most conventional compliant beams have a length to width ratio of 20:1 or higher. Here, compliant members of relatively low aspect ratio (between 20:1 and 3:1) will be considered and treated as beam elements. The change in cross section of the beam due to elongation of the beam axis will be incorporated into the formulation.

Consider the initially straight cantilever beam of uniform cross section shown in Fig. 2 with general loads at its free end. The force at the free end is divided into its horizontal and vertical components, Fx and Fy, respectively. A moment Mz, also acts at the tip of the beam. The inset shows a segment of the beam at a distance $s∈[0,L]$ along the length of the beam. It is subjected to an axial force P(s), a shear force V(s), and a bending moment M(s). The beam element is rotated by an angle θ to the x-axis, and the cross section is rotated by a angle ψ to the y-axis. It is assumed that the cross section remains planar and is not distorted by the loads acting on the beam. The initial dimensions of the beam are length, L0, area A0, and second moment of area I0.

Fig. 2
Fig. 2
Close modal
Upon loading, the length and cross section of the beam are subject to change. At any point along the length of the beam, let A(s) and I(s) be the area and second moment of area of the beam, respectively. For a Timoshenko beam
$M(s)=EI(s)ψ′(s)$
(1)
$ψ(s)=θ(s)+V(s)κGA(s)$
(2)
where E is the modulus of elasticity of the material, κ is the Timoshenko shear coefficient, and G is the shear modulus. Based on force and moment equilibrium equations, the internal loads can be calculated as
$M(s)=Mz+Fy(a−x(s))−Fx(b−y(s))$
(3)
$V(s)=−Fy cos θ(s)+Fx sin θ(s)$
(4)
$P(s)=Fx cos θ(s)+Fy sin θ(s)$
(5)

where a and b are the x and y coordinates of the tip of the beam. The axial extension is represented using λ, the extension ratio. It is defined as $λ=1+e$, where e is the engineering strain.

The Poisson's effect is represented in a nonlinear form shown in the equation given below
$ww0=(ll0)−ν=λ−ν$
(6)

where w is a dimension orthogonal to the length l, and ν is the Poisson's ratio of the material. The subscript 0 represents the initial dimension.

It follows that the area and second moment of area of the cross section are calculated as
$I(s)=I0 λ(s)−4ν$
(7)
$A(s)=A0 λ(s)−2ν$
(8)
Differentiating Eq. (1) with s yields
$M′(s)=EI′(s)ψ′(s)+EI(s)ψ″(s)$
(9)
and similarly differentiating Eq. (2) gives
$ψ′(s)=−V(s)A′(s)κGA(s)2+V′(s)κGA(s)+θ′(s)$
(10)
$ψ″(s)=2V(s)A′(s)2κGA(s)3−2A′(s)V′(s)κGA(s)2 −V(s)A″(s)κGA(s)2+V″(s)κGA(s)+θ″(s)$
(11)
In order to obtain a relationship for $M′(s)$, Eq. (3) can be differentiated to obtain
$M′(s)=−Fy cos θ(s)+Fx sin θ(s)$
(12)
where we can substitute the geometric relationships
$dxds=cos θ(s); dyds=sin θ(s)$
Equations (7) and (8) can be substituted into Eq. (11). The new equation for $ψ″(s)$ along with Eq. (12) can be substituted into Eq. (9) to obtain Eq. (13). V(s) and P(s) can be substituted from Eqs. (4) and (5) to obtain the expanded form in terms of the forces Fx and Fy
$V(s)=1A0GκEI0λ(s)−4ν−2[A0Gκλ(s)2θ″(s)−2ν{2ν+1}λ(s)2νλ′(s)2V(s)+2νλ(s)2ν+1λ″(s)V(s)−4A0Gκνλ(s)θ′(s)λ′(s)+λ(s)2ν+2{θ″(s)P(s)−θ′(s)2V(s)}]$
(13)
In the axial direction, for any segment of the beam, Hooke's law gives $σ(s)=Eε(s)$ and the definition of axial stress, $σ(s)=P(s)/A(s)$ (σ is the axial stress and ϵ is the axial logarithmic strain). Also, $ε=log λ$. These relationships, along with Eq. (5) can be substituted into Hooke's law to give
$Fx cos θ(s)+Fy sin θ(s)=EA0 λ(s)−2ν log λ(s)$
(14)
Two boundary conditions can be established for $θ(s)$. At the fixed end of the beam, the cross section cannot rotate, and therefore, $ψ(0)=0$. Also, at s = L, the bending moment $M(s)=Mz$, since the moment arms of Fx and Fy are zero. Evaluating Eq. (2) at s = 0 and Eq. (10) at s = L gives
$θ(0)=−V(0) λ(0)2νκGA(0)$
(15)
$Mz λ(L)4νEI0=λ(L)2νV′(L)κA0G+θ′(L)+2νV(L)λ(L)2νλ′(L)κA0Gλ(L)$
(16)

Equations (13) and (14) with the boundary conditions Eqs. (15) and (16) form a system of differential-algebraic equations (DAEs). The set of DAEs is solved numerically to obtain interpolating functions for $λ(s)$ and $θ(s)$ in the domain $[0,L0]$.

For a uniform beam, it is possible to convert the above equations into a dimensionless form. The nondimensional form [18] is important for dividing the problem into two; a problem of load design and a problem of geometric design. For convenience, let us define nondimensionalized loads as
$αx=FxL02EI0; αy=FyL02EI0; β=MzL0EI0$
In addition, another parameter which we will call the beam geometry factor, fbg, is defined as
$fbg=I0A0L02$
The Beam Geometry Factor: The parameter fbg is very significant for the remainder of the discussion in this paper. It relates the extension of the beam element along the neutral axis to the deformation due to bending. A large value of fbg means that the elongation and shear effects in the beam cannot be ignored under general loading, whereas a small value of fbg suggests that the beam deformation is mainly determined by the bending moment equation (classical Euler beam theory). fbg can also be expressed as the square of the slenderness ratio
$fbg=(kL0)2$
(17)

where $k=I0/A0$.

In this paper, the limits of fbg are set to values of $10−5$ and $10−2$. For a beam with a circular cross section, this would mean a ratio between the diameter and the length varying between $1:80$ and $1:2.5$. This should cover most of the types of beams used in various applications.

To obtain the tip deflection in terms of the coordinates x, y, and the angle ψ, the following relationships are used. Equation (20) is derived by setting s = L0 in Eq. (2) and substituting expressions for the other variables
$xtip=∫0L0λ(s)cos θ(s)ds$
(18)
$ytip=∫0L0λ(s)sin θ(s)ds$
(19)
$ψtip=θ(L0)−1κ[2fbg(1+ν) λ(L0)2ν×(αy cos θ(L0)−αx sin θ(L0))]$
(20)

### Solution Process.

The equations derived for the modified Timoshenko beam form a set of DAEs in a boundary value problem (BVP) format. This makes the solution process nontrivial. The equations are too complex to consider closed form or analytical solutions, which are difficult to obtain even for Euler–Bernoulli beams with combined loads. Numerical solutions may be more practical for a problem such as this.

Notice that Eq. (14) is an algebraic relationship between $λ(s)$ and $θ(s)$. However, the solution for $λ(s)$ in this equation is in the form of a LambertW function [32], also called the “product-log” function. While it is possible to substitute this solution of $λ(s)$ into Eq. (13) and convert it into a system of ordinary differential equations (ODEs), it would still contain a set of terms with the product-log function and make the numerical solution process more cumbersome. It may be possible to develop a simpler relationship between $λ(s)$ and $θ(s)$ in Eq. (14) using a series expansion of the terms, and a few preliminary attempts were made to tackle this problem. However, they did not provide any improvements to the solution process.

The numerical solution for DAEs in the form of initial value problems (IVPs) is available in many commercially available solvers. This is possible because the boundary conditions only need to be satisfied at the first point, and therefore, the differential equations can be integrated forward along the direction of the independent variable (in this case, the length of the beam s). However, for BVPs, the solution process usually follows a “shooting method,” where a “guess” is made for the solution at the initial position [33]. For a highly complex system, it can be very difficult to arrive at the right guess. For a problem as complex as the one discussed here, there is a very little literature on solving BVP DAEs. A few numerical solvers have been built for DAEs [34,35]. Ascher and Spiteri developed the code named COLDAE [36] which was transferred to the open source language R by Soetaert et al. [37] in the package BVPSolve. However, the highly coupled relationship between the variables λ and θ made it difficult to solve these equations using this package.

The solver used in this paper is Wolfram Mathematica's NDSolve, which is only capable of solving DAEs in the IVP format. To work around this, the idea behind the shooting method is used. A fifth equation for the initial value of $θ′(s)$ is proposed and the value assumed to be a variable, C
$θ′(0)=C$
(21)

This equation is used as a boundary condition at the initial point, s = 0. The DAE system comprising of Eqs. (13)(15) and (21) is solved as a family of solutions with the variable C. Equation (16) is then solved to determine of value of C that satisfies the constraint at s = L0. This gives the solution of $λ(s)$ and $θ(s)$ over the range $s∈[0,L0]$.

For the solution process, it is seen that for small values of fbg (i.e., very little elongation), the equations are solved without much difficulty and the solution process is rather robust. However, for large values of fbg, only a good guess for the value of C can guide the solver to a feasible solution, especially for large combined loads.

### Beam Theory Results.

The applicability and validity of the modified beam theory formulation can be demonstrated by comparing it against conventional beam theory methods and FEA. For long slender beams, the results are similar to those obtained from Euler–Bernoulli equations because shear and extension effects are negligible compared to bending. However, for aspect ratios of 10:1 and lower, these effects cannot be neglected. The equations for Euler beam model and Timoshenko beam model are given in the Appendix.

The FEA was performed in abaqus software using B21 beam elements with the geometric nonlinearities allowed, and each beam is divided into 100 elements. The material was set to be linear isotropic with unit elastic modulus and Poissons ratio of 0.4, and the length of the beam was unity. The loads on the beam were specified using αx, αy, and β, as defined in Sec. 2. The beam is of rectangular cross section and Poisson's ratio ν = 0.4. This gives κ = 0.85 as defined by Cowper [10]. The beam theory formulation presented in this paper will be called modified Timoshenko beam theory (MTBT) for the remainder of the paper.

Figure 3(a) shows the tip deflection of a beam with fbg = 0.005 under the action of combined loads. αy = 1 and β = 0.5 are held constant while x is varied in the range [−1,7]. For high values of αx, the extension effects are noticeable. It can also be seen that the error, as shown in Fig. 3(b), is much lower than the Euler beam or Timoshenko beam, especially as αx increases. In this case, the length of the beam increases by 3–5% across the range of loads.

Fig. 3
Fig. 3
Close modal

The results of similar comparisons under compressive loads are shown in Fig. 4. A rectangular beam of unit length with width, w = 0.3 and thickness, t = 0.4 (fbg = 0.0075) was subjected to αx = −1, β = 0, and αy varying in the range [0,0.5]. The plot at the top shows the error in prediction compared to FEA for calculating the x-coordinate of the beam tip, while the bottom one shows the results for the y-coordinate. The MTBT works very well to approximate the FEA results for the x-coordinate, and does much better than both the Euler beam and Timoshenko beam for the y-coordinate.

Fig. 4
Fig. 4
Close modal

Figure 5 shows the extension in a beam of varying cross section under the action of a transverse load (αy) as calculated using MTBT. The extension is presented for five different beams with values of fbg varying from 0.0001 to 0.01. Four loading cases were used: $αy∈$ [0.6,1.2,1.6,2]. The results demonstrate the relationship between fbg and the amount of extension in the beam, even for a transverse load. While it is negligible for long beams, the extension in short joints can be anywhere between 5% and 10%, which cannot be ignored.

Fig. 5
Fig. 5
Close modal

These results indicate significant extension or compression effects in soft joints of low aspect ratio. It stands to reason that traditional PRBMs would be insufficient for accurately predicting the behavior of such compliant elements. In the following section, the MTBT solution is used to determine the values of the pseudo-rigid-body (PRB) parameters. The beam theory results are reproduced for four different values of fbg, for a large range of loading cases. These results are then used to determine the optimal values of the PRB parameters for three PRBMs with extension springs. The performance of the new models is compared against traditional PRBMs in literature to demonstrate the effect of extension in soft joints.

## Comparing PRBMs

PRBMs provide a convenient and intuitive method for the analysis of compliant mechanisms. While most of the pioneering work was done by Howell and Midha [22,23,38], there have been quite a few advances in this area, and many others have adapted the idea to be used for different situations with varying accuracy. In this paper, the following three PRBMs have been studied in detail. The values of the parameters and the schematics of the models are given in Table 1.

Table 1

PRBMs being compared

PRBMSchematicPRB parameters
1R$γ=0.85Kθ=2.25EIL$
2R$γ0=0.1 γ1=0.54 γ2=0.36Kθ1=3.40EIL Kθ2=1.58EIL$
3R$γ0=0.1 γ1=0.35 γ2=0.4 γ3=0.15Kθ1=3.51EIL Kθ2=2.99EIL Kθ3=2.58EIL$
PRBMSchematicPRB parameters
1R$γ=0.85Kθ=2.25EIL$
2R$γ0=0.1 γ1=0.54 γ2=0.36Kθ1=3.40EIL Kθ2=1.58EIL$
3R$γ0=0.1 γ1=0.35 γ2=0.4 γ3=0.15Kθ1=3.51EIL Kθ2=2.99EIL Kθ3=2.58EIL$
• 1R model: This model was defined in Refs. [23] and [38]. It consists of one torsion spring of stiffness $Kθ$ attached to two rigid segments of length γ and $1−γ$. The values of the PRB parameters vary with the loading condition, and it is only defined for beams with end loads and no moments.

• 2R model: This model has two revolute joints with torsion springs and three rigid segments. It was defined in Ref. [39]. It also has parameter values that vary with load.

• 3R model: It was defined in Ref. [18]. This model consists of four rigid segments joined by torsion springs on revolute joints. Unlike the other two, it is load-independent.

Since this paper deals with the errors that arise in PRBMs due to axial extension, three PRBMs with extension springs will be considered as possible solutions to the problem. These models will have optimized PRB parameters for different values of fbg. The process of determining the optimum values of the PRBMs is similar to the discussions in previous work done by the authors [19,40]. The errors are calculated based on the deflection of the beam tip in terms of its position and orientation. The PRBMs are listed along with their optimized values in Table 2.

Table 2

Optimized PRBMs for varying fbg under general loads

Optimized values
PRBMSchematicPRB parametersfbg = 0.0001fbg = 0.001fbg = 0.005fbg = 0.01
PRγ0.3010.2970.2750.251
$kθ$1.4801.4881.5401.592
kex3.1772.2401.8871.378
RPγ0.3020.2980.2810.259
$kθ$1.4741.4821.5121.553
kex0.0860.8710.8700.955
RPRγ0.2050.2010.1870.169
$kθ$2.0352.0322.0312.024
kex0.2240.6380.9780.997

Optimized values
PRBMSchematicPRB parametersfbg = 0.0001fbg = 0.001fbg = 0.005fbg = 0.01
PRγ0.3010.2970.2750.251
$kθ$1.4801.4881.5401.592
kex3.1772.2401.8871.378
RPγ0.3020.2980.2810.259
$kθ$1.4741.4821.5121.553
kex0.0860.8710.8700.955
RPRγ0.2050.2010.1870.169
$kθ$2.0352.0322.0312.024
kex0.2240.6380.9780.997

Note: Values of $kθ$ and kex are dimensionless, and should be multiplied by factors EI/L and EA/L, respectively.

In order to compare the accuracy of the PRBMs, two loading cases were studied. At first, a qualitative assessment of the beam deflection was performed under the action of vertical forces on the beam tip. Vertical loads contribute to bending and axial elongation. The loads applied were in the range $αy∈[0,2]$. The second case looked at combined horizontal and vertical forces and a moment. These loads were in the range $αx∈[−1,1], αy∈[−1,1]$ and $β∈[−0.5,0.5]$ and were used to study beams with more complex beam curves. Of special importance are beams with points of inflection, where the bending moment equals zero.

Figure 6(a) shows the tip deflection as αy goes from 0 to 2. It has previously been shown that the 3R model is very close to the classical Euler beam theory and is more accurate than the 2R and 1R models [39]. However, it is noticeable that when extension effects start to weigh in, the error in the PRB results gradually increases. On the other hand, the optimized PRBMs with extension springs demonstrate negligible error, as is evident in Fig. 6(b). The change in error with increase in fbg is presented in Fig. 7(a).

Fig. 6
Fig. 6
Close modal
Fig. 7
Fig. 7
Close modal

The change in error in calculating tip deflection using the three PRBMs for combined loads is given in Fig. 7(b). Table 2 lists the optimized PRB values when fbg varies from 0.0001 to 0.01. A large set of loading points was necessary for the optimization in this case to ensure that the error is calculated over a population of data that is representative of the beam behavior under such loads.

### Discussion.

Figures 6(a) and 6(b) clearly show that the models with the extension springs are more accurate in predicting the tip deflection for large fbg. The percentage error in the 1R and 3R models increases significantly with increase in fbg, as is evident from Fig. 7(a). The same figure also shows that the RPR and RP models are better suited for softer beams (For the PRB models, R and P refer to revolute and prismatic joints, respectively). This leads us to conclude that the 3R model has very high accuracy for long beams, whereas the thicker beams are better represented by models with extension springs. For the combined loading case, it is again noticeable from Fig. 7(b) that the PR model is not suitable for such beams as its error remains the same. However, the RPR and RP models show continuous decrease in error as fbg increases. The low error of the 3R model for long beams can be attributed to the higher degrees-of-freedom of the 3R model as opposed to the RP and PR models.

There are a couple of interesting points to be inferred from Table 2. It can be noted that for the RPR model, which is the most accurate, the value of γ decreases as the fbg value increases. This means that the length of the extension spring increases with fbg, which agrees with intuition. The value of $kθ$ remains constant, suggesting that the bending stiffness is unaffected by the change in fbg, as expected. On the other hand, the value of kex increases, demonstrating that extension is affected by beam geometry.

Generally speaking, there is a need to consider the error induced due to beam elongation when using PRBMs. While the 3R model is highly accurate for thin long beams, softer beams with larger fbg values require the use of PRBMs with extension springs, such as the RPR model, for analysis.

## Discussion

The soft joints analyzed in this work use a constant value for the modulus of elasticity, E, in the beam theory approximation as well as the PRBMs. This assumption is reasonable since many elastomeric materials are used in elongation as high as 400%, whereas the values experienced by these joints are well below 50%. Such material characteristics can possibly be incorporated by converting material models for hyperelastic materials [41,42] into simple stress–strain forms and using them in these equations.

The discussion is limited to static analysis since this is widely used in the design of compliant mechanisms. Dynamic analysis of compliant mechanisms using PRBMs can also be achieved by including inertial effects. The solution process for the MTBT approach detailed here is only capable of dealing with beams with a single inflection point. However, it may be possible to use prior knowledge of the number of inflection points to change the solution, as shown by Zhang and Chen [43].

The current formulation of the MTBT makes the solution process of the BVP DAE system rather complex. A good guess for the slope variable, C, is needed to obtain the solution when large combined loads act on the beam. The relationship between λ and θ is a key bottleneck here. It would be beneficial if the BVP DAE system can be simplified into an ODE, or a better solver is devised.

The FEA used for comparison captures the shear in the beam since the elements are constructed using Timoshenko beam theory. It also takes into account the extension in the beam as well as the Poisson's effect, both of which are also captured by the modified Timoshenko beam presented Sec. 2. The PRBM is capable of reproducing the effect of linear expansion due to the extension spring. However, it does not directly capture the resulting change in cross section by Poisson's effect. This is embedded within the values of the PRB parameters, which change based on the value of fbg, as shown in Table 2.

## Conclusions and Future Work

In the first part of this paper, a beam theory approach was proposed for the analysis of compliant members with a low ratio of length to width. The proposed theory incorporates extension and Poisson's effects into Timoshenko beam theory, and has been proven to work for 2D quasi-static applications. The solution process is complicated by the relationship between the bending, shear, and extension effects, and a shooting method with a numerical ODE solver is used for calculations. The results from this approach are then highlighted for a couple of different cases. The development of the beam theory approximation (MTBT), along with the final equations and method of solution, is one of the major contributions of this work.

The rest of the paper makes use of these results to study of accuracy of PRBM, and measures to be taken for analyzing short beams of soft materials. The beam geometry factor, fbg, is used to categorize beams as short and long. The extension effects in short beams are studied and compared to results from the conventional PRBMs. It is found that for short beams, extension springs in PRBMs lead to a significant decrease in error during analysis. This is also validated using the example of a parallel-guiding mechanism. The study of the effect of extension (in relation to the aspect ratio of the joints) on the error in the PRBMs, while illustrating the difference in the PRB parameters as a function of the change in shape, is a second contribution of the work presented here.

These concepts can be further extended by studying the properties of similar soft materials and incorporating those effects into PRBMs. The beam theory approximation presented in this paper can be further validated by experiments and testing, which can be rather difficult when dealing with elastomeric materials such as the one shown in Fig. 1(a). It is possible to extend the use of PRBMs to the design of compliant mechanisms, and also incorporate nonlinear materials and members into the process. This will open up the field of compliant mechanisms to more fabrication processes. By utilizing a wide range of materials, compliant mechanisms may have a gateway into many more applications on different length scales.

## Acknowledgment

The authors would like to extend their thanks to the discussion with Prof. Wenbin Yu from Purdue University on the beam theory. This material is based upon the work supported by the Air Force Office of Scientific Research under Contract No. AFOSR FA9550-12-1-0070 and the National Science Foundation under Grant Nos. CMMI-1144022 and CMMI-1161841.

### Appendix

The equations for Euler and Timoshenko beam models referred to in Sec. 2.2 are given below. These equations were used to calculate the results used for comparison in Figs. 3(a), 3(b), and 4.

#### Euler–Bernoulli Beam

$ODE→ 1L2d2θds2=1EI(Fx cos θ−Fy sin θ)BC→ θ(0)=0 θ′(L)=MzLEIxtip=∫0L cos θds ytip=∫0L sin θds$

#### Timoshenko Beam

$ODE→ (1+Fx cos θ+Fy sin θκGA)d2θds2+Fy cos θ−Fx sin θκGAdθ2ds=1EI(Fx cos θ−Fy sin θ)BC→ θ(0)=FyκGA θ′(L)=MzLEI(1+Fx cos θ+Fy sin θκGA)xtip=∫0L cos θds ytip=∫0L sin θds$

## References

1.
Popescu
,
B.
, and
Hodges
,
D. H.
,
2000
, “
On Asymptotically Correct Timoshenko-Like Anisotropic Beam Theory
,”
Int. J. Solids Struct.
,
37
(
3
), pp.
535
558
.
2.
Yu
,
W.
, and
Blair
,
M.
,
2012
, “
GEBT: A General-Purpose Nonlinear Analysis Tool for Composite Beams
,”
Compos. Struct.
,
94
(
9
), pp.
2677
2689
.
3.
Timoshenko
,
S.
,
1921
, “
On the Correction for Shear of the Differential Equation for Transverse Vibrations of Prismatic Bars
,”
Philos. Mag. Ser. 6
,
41
(
245
), pp.
744
746
.
4.
Chaisomphob
,
T.
,
Nishino
,
F.
,
Hasegawa
,
A.
, and
Alygamalaly
,
A.-S.
,
1986
, “
An Elastic Finite Displacement Analysis of Plane Beams With and Without Shear Deformation
,”
Doboku Gakkai Ronbunshu
,
1986
(
368
), pp.
169
177
.
5.
Yoshiaki
,
G.
,
Tomoo
,
Y.
, and
Makoto
,
O.
,
1990
, “
Elliptic Integral Solutions of Plane Elastica With Axial and Shear Deformations
,”
Int. J. Solids Struct.
,
26
(
4
), pp.
375
390
.
6.
Stephen
,
N. G.
,
1980
, “
,”
ASME J. Appl. Mech.
,
47
(
1
), pp.
121
127
.
7.
Kaneko
,
T.
,
1975
, “
On Timoshenko's Correction for Shear in Vibrating Beams
,”
J. Phys. D: Appl. Phys.
,
8
(
16
), pp.
1927
1936
.
8.
Jensen
,
J. J.
,
1983
, “
On the Shear Coefficient in Timoshenko's Beam Theory
,”
J. Sound Vib.
,
87
(
4
), pp.
621
635
.
9.
Hutchinson
,
J. R.
,
2000
, “
Shear Coefficients for Timoshenko Beam Theory
,”
ASME J. Appl. Mech.
,
68
(
1
), pp.
87
92
.
10.
Cowper
,
G. R.
,
1966
, “
The Shear Coefficient in Timoshenkos Beam Theory
,”
ASME J. Appl. Mech.
,
33
(
2
), pp.
335
340
.
11.
Bejgerowski
,
W.
,
Gerdes
,
J. W.
,
Gupta
,
S. K.
,
Bruck
,
H. A.
, and
Wilkerson
,
S.
,
2010
, “
Design and Fabrication of a Multi-Material Compliant Flapping Wing Drive Mechanism for Miniature Air Vehicles
,”
ASME
Paper No. DETC2010-28519.
12.
Venkiteswaran
,
V. K.
, and
Su
,
H.
,
2014
, “
Optimization of Mechanism Design of Flapping Wing MAV
,”
55th AIAA/ASMe/ASCE/AHS/SC Structures, Structural Dynamics, and Materials Conference
, American Institute of Aeronautics and Astronautics, Reston, VA, Paper No.
AIAA
2014-0573.
13.
Lu
,
K.-J.
, and
Kota
,
S.
,
2003
, “
Design of Compliant Mechanisms for Morphing Structural Shapes
,”
J. Intell. Mater. Syst. Struct.
,
14
(
6
), pp.
379
391
.
14.
Zhou
,
L.
,
Marras
,
A. E.
,
Su
,
H.-J.
, and
Castro
,
C. E.
,
2014
, “
DNA Origami Compliant Nanostructures With Tunable Mechanical Properties
,”
ACS Nano
,
8
(
1
), pp.
27
34
.
15.
Howell
,
L. L.
,
Magleby
,
S. P.
, and
Olsen
,
B. M.
,
2013
,
Handbook of Compliant Mechanisms
,
Wiley
,
New York
.
16.
Howell
,
L. L.
,
2001
,
Compliant Mechanisms
,
Wiley
,
New York
.
17.
Awtar
,
S.
,
Slocum
,
A. H.
, and
Sevincer
,
E.
,
2006
, “
Characteristics of Beam-Based Flexure Modules
,”
ASME J. Mech. Des.
,
129
(
6
), pp.
625
639
.
18.
Su
,
H.-J.
,
2009
, “
A Pseudorigid-Body 3r Model for Determining Large Deflection of Cantilever Beams Subject to Tip Loads
,”
ASME J. Mech. Rob.
,
1
(
2
), p.
021008
.
19.
Venkiteswaran
,
V. K.
, and
Su
,
H.-J.
,
2015
, “
A Parameter Optimization Framework for Determining the Pseudo-Rigid-Body Model of Cantilever-Beams
,”
Precis. Eng.
,
40
, pp.
46
54
.
20.
Venkiteswaran
,
V. K.
, and
Su
,
H.-J.
,
2015
, “
Effect of Beam Geometry on the Accuracy of Pseudo-Rigid-Body Models
,”
ASME
Paper No. DETC2015-47267.
21.
Frecker
,
M. I.
,
Ananthasuresh
,
G. K.
,
Nishiwaki
,
S.
,
Kikuchi
,
N.
, and
Kota
,
S.
,
1997
, “
Topological Synthesis of Compliant Mechanisms Using Multi-Criteria Optimization
,”
ASME J. Mech. Des.
,
119
(
2
), pp.
238
245
.
22.
Howell
,
L. L.
, and
Midha
,
A.
,
1994
, “
A Method for the Design of Compliant Mechanisms With Small-Length Flexural Pivots
,”
ASME J. Mech. Des.
,
116
(
1
), pp.
280
290
.
23.
Howell
,
L. L.
, and
Midha
,
A.
,
1995
, “
Parametric Deflection Approximations for End-Loaded, Large-Deflection Beams in Compliant Mechanisms
,”
ASME J. Mech. Des.
,
117
(
1
), pp.
156
165
.
24.
Saxena
,
A.
, and
Kramer
,
S. N.
,
1998
, “
A Simple and Accurate Method for Determining Large Deflections in Compliant Mechanisms Subjected to End Forces and Moments
,”
ASME J. Mech. Des.
,
120
(
3
), pp.
392
400
.
25.
Edwards
,
B. T.
,
Jensen
,
B. D.
, and
Howell
,
L. L.
,
1999
, “
A Pseudo-Rigid-Body Model for Initially-Curved Pinned-Pinned Segments Used in Compliant Mechanisms
,”
ASME J. Mech. Des.
,
123
(
3
), pp.
464
468
.
26.
Kimball
,
C.
, and
Tsai
,
L.-W.
,
2002
, “
Modeling of Flexural Beams Subjected to Arbitrary End Loads
,”
ASME J. Mech. Des.
,
124
(
2
), pp.
223
235
.
27.
,
M. H.
,
2001
, “
Variable Parametric Pseudo-Rigid-Body Model for Large-Deflection Beams With End Loads
,”
Int. J. Non-Linear Mech.
,
36
(
7
), pp.
1123
1133
.
28.
Chen
,
G.
,
Xiong
,
B.
, and
Huang
,
X.
,
2011
, “
Finding the Optimal Characteristic Parameters for 3r Pseudo-Rigid-Body Model Using an Improved Particle Swarm Optimizer
,”
Precis. Eng.
,
35
(
3
), pp.
505
511
.
29.
Cham
,
J. G.
,
Bailey
,
S. A.
,
Clark
,
J. E.
,
Full
,
R. J.
, and
Cutkosky
,
M. R.
,
2002
, “
Fast and Robust: Hexapedal Robots Via Shape Deposition Manufacturing
,”
Int. J. Rob. Res.
,
21
(
10–11
), pp.
869
882
.
30.
Dollar
,
A.
, and
Howe
,
R.
,
2006
, “
A Robust Compliant Grasper Via Shape Deposition Manufacturing
,”
IEEE/ASME Trans. Mechatronics
,
11
(
2
), pp.
154
161
.
31.
Bailey
,
S. A.
,
Cham
,
J. A.
,
Cutkosky
,
M. R.
, and
Full
,
R. J.
,
2000
, “
Biomimetic Robotic Mechanisms Via Shape Deposition Manufacturing
,”
Robotics Research—International Symposium
, Vol.
9
, pp.
403
410
.
32.
Corless
,
R. M.
,
Gonnet
,
G. H.
,
Hare
,
D. E. G.
,
Jeffrey
,
D. J.
, and
Knuth
,
D. E.
,
1996
, “
On the LambertW Function
,”
,
5
(
1
), pp.
329
359
.
33.
Ascher
,
U. M.
, and
Petzold
,
L. R.
,
1998
,
Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations
,
SIAM
,
.
34.
Petzold
,
L. R.
,
1982
, “
A Description of DASSL: A Differential/Algebraic System Solver
,” IMACS World Congress, pp.
430
432
.
35.
Tolsma
,
J.
, and
Barton
,
P. I.
,
2000
, “
DAEPACK: An Open Modeling Environment for Legacy Models
,”
Ind. Eng. Chem. Res.
,
39
(
6
), pp.
1826
1839
.
36.
Ascher
,
U.
, and
Spiteri
,
R.
,
1994
, “
Collocation Software for Boundary Value Differential-Algebraic Equations
,”
SIAM J. Sci. Comput.
,
15
(
4
), pp.
938
952
.
37.
Soetaert
,
K.
,
Cash
,
J.
, and
Mazzia
,
F.
,
2012
,
Solving Differential Equations in R
,
Springer Berlin Heidelberg
,
Berlin, Heidelberg
.
38.
Howell
,
L. L.
,
Midha
,
A.
, and
Norton
,
T. W.
,
1996
, “
Evaluation of Equivalent Spring Stiffness for Use in a Pseudo-Rigid-Body Model of Large-Deflection Compliant Mechanisms
,”
ASME J. Mech. Des.
,
118
(
1
), pp.
126
131
.
39.
Yu
,
Y.-Q.
,
Feng
,
Z.-L.
, and
Xu
,
Q.-P.
,
2012
, “
A Pseudo-Rigid-Body 2r Model of Flexural Beam in Compliant Mechanisms
,”
Mech. Mach. Theory
,
55
, pp.
18
33
.
40.
Kalpathy Venkiteswaran
,
V.
, and
Su
,
H.-J.
,
2016
, “
A 3-Spring Pseudo-Rigid-Body Model for Soft Joints With Significant Elongation Effects
,”
ASME J. Mech. Rob.
,
8
(
6
), p.
061001
.
41.
Arruda
,
E. M.
, and
Boyce
,
M. C.
,
1993
, “
A Three-Dimensional Constitutive Model for the Large Stretch Behavior of Rubber Elastic Materials
,”
J. Mech. Phys. Solids
,
41
(
2
), pp.
389
412
.
42.
Yeoh
,
O. H.
,
1993
, “
Some Forms of the Strain Energy Function for Rubber
,”
Rubber Chem. Technol.
,
66
(
5
), pp.
754
771
.
43.
Zhang
,
A.
, and
Chen
,
G.
,
2013
, “
A Comprehensive Elliptic Integral Solution to the Large Deflection Problems of Thin Beams in Compliant Mechanisms
,”
ASME J. Mech. Rob.
,
5
(
2
), p.
021006
.