Extension Effects in Compliant Joints and Pseudo-Rigid-Body Models

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 [6–9]. 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 [18–20].

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].

Shape deposition manufacturing (SDM) is a step-by-step method for depositing layers of desired materials in varying shapes to produce required parts [29–31]. 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.

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.

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, F_x and F_y, respectively. A moment M_z, 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, L₀, area A₀, and second moment of area I₀.

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.

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

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]$⁠.

where $k=I0/A0$⁠.

In this paper, the limits of f_bg 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.

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.

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 = L₀. 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 f_bg (i.e., very little elongation), the equations are solved without much difficulty and the solution process is rather robust. However, for large values of f_bg, only a good guess for the value of C can guide the solver to a feasible solution, especially for large combined loads.

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 f_bg = 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.

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 (f_bg = 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.

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 f_bg 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 f_bg 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.

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 f_bg, 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.

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.

• 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 f_bg. 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.

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 f_bg is presented in Fig. 7(a).

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 f_bg 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.

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 f_bg. The percentage error in the 1R and 3R models increases significantly with increase in f_bg, 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 f_bg 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 f_bg value increases. This means that the length of the extension spring increases with f_bg, which agrees with intuition. The value of $kθ$ remains constant, suggesting that the bending stiffness is unaffected by the change in f_bg, as expected. On the other hand, the value of k_ex 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 f_bg values require the use of PRBMs with extension springs, such as the RPR model, for analysis.

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 f_bg, as shown in Table 2.

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, f_bg, 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.

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.

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.