# Extending Timoshenko Beam Theory for Large Deflections in Compliant Mechanisms

Ke Wu, Gang Zheng, Guimin Chen

► To cite this version:

Ke Wu, Gang Zheng, Guimin Chen. Extending Timoshenko Beam Theory for Large Deflec- tions in Compliant Mechanisms. Journal of Mechanisms and Robotics, 2023, 15 (6), pp.061012. 10.1115/1.4056501 . hal-03917098

HAL Id: hal-03917098

https://inria.hal.science/hal-03917098

Submitted on 3 Jan 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Extending Timoshenko Beam Theory for Large Deflections in Compliant Mechanisms

**Ke Wu**

Univ. Lille, Inria, CNRS
Centrale Lille, UMR 9189 CRISTAL
F-59000 Lille, France
Email: ke.wu@inria.fr

**Gang Zheng***

Univ. Lille, Inria, CNRS
Centrale Lille, UMR 9189 CRISTAL
F-59000 Lille, France
Email: gang.zheng@inria.fr

**Guimin Chen**

State Key Laboratory for Manufacturing Systems Engineering and Shaanxi Key Laboratory of Intelligent Robots
Xi'an Jiaotong University, Xi'an, Shaanxi, 710049, China
Email: guimin.chen@gmail.com

Compliant Mechanisms (CMs) have presented its inher-
ently advantageous properties due to the fact that CMs uti-
lize elastic deformation of the elementary flexible members
to transfer motion, force and energy. Previously, the clas-
sic Euler-Bernoulli beam theory is the most used theory in
terms of modeling large beam deflections in CMs. However,
it has some assumptions that may decrease the modeling ac-
curacy, such as ignoring the shear strain and the axial strain
of cross sections. In this paper, to take into account the shear
and axial strains, we adopt Timoshenko beam theory along
with some modifications to consider the axial elongation. To
simplify the complexity of the proposed governing boundary
value problem (BVP), we transform the BVP into an explicit
formulation and use weighted residual methods to numeri-
cally approximate the solution. We first focus on the single-
beam deflection of a straight beam and an initially curved
beam (ICB) using Euler Bernoulli beam theory, Timoshenko
beam theory and solid mechanics to analyze the contribu-
tions of the influences of shear and axial strains in beam de-
flections. Then, we prove the feasibility of the proposed mod-
eling strategy via mechanism synthesis for a bi-stable mech-
anism and an ICB-based parallelogram mechanism. Finally,
the deduction of the mathematical model and the numerical
results are provided along with brief analysis on the mechan-
ical performances of the studied CMs.

# 1 Introduction

## 1.1 Compliant mechanisms

Compliant mechanisms (CMs) have been serving as a
promising alternative for current mechanical applications

due to their inherent advantages [1]. Conventional rigid-
body mechanisms utilize the cooperation of kinematic pairs
and rigid links to provide the transfer of motion, force and
energy. CMs, as a unique concept of mechanical systems,
can realize the same function merely through the deforma-
tion of the built-in elementary flexible and compliant mem-
bers [1]. This powerful concept has presented many desired
characteristics over rigid ones: compliance due to the flexi-
bility of the used elastic material [1]; monolithic design that
results in increased motion precision with no concerns about
errors from assembly [2], simplified manufacturing process
[1] [2] [3], reduced weight [1] [3], cost reduction [1] [2] [3]
and less maintenance [1] [3]. These mentioned desired char-
acteristics have made CMs a gradually mature component
in many mechanical applications, such as high-accuracy po-
sitioning motion stages [4], compliant kinematic joints [5],
bi/tri-stable mechanisms [6] and other more complex appli-
cations that integrate the former three [7].

## 1.2 Commonly used modeling theories and the corresponding methods for modeling beams in CMs

As mentioned in the previous section, the built-in flexi-
ble members, which are normally straight beams and initially
curved beams (ICBs) [8], play an important role in realizing
the functions of CMs. Therefore, accurately and efficiently
modeling beams will be a common focus either for perfor-
mance predictions or model-based optimization [9]. There
are two major aspects: first, the fundamental modeling the-
ory should be properly chosen; second, and the correspond-
ing analytical or numerical methods that can implement the
mentioned theory are logically needed.

1. Solid mechanics: solid mechanics is supposed to be the most general theory to model any deformable object, and it is often implemented via FEM [10]. This is

<footnote>*Corresponding author: Gang Zheng University of Lille, Inria, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000 Lille, France, Tel: +33 3 59 57 79 53, Fax: +33 3 59 57 78 50.</footnote>
due to the fact that general (irregular) deformable bodies normally don't have a global governing equation. A deformable body needs to be discretized into a finite number of fine mesh elements, and Euler-Lagrange method is then applied to these elements, resulting in a high-dimensional differential equation to govern the overall deformation of the studied object.

2. **Euler-Bernoulli beam theory:** in the current literature, Euler-Bernoulli beam theory [11] is the most used theory to govern the deflection of a slender beam. For example, Howell used elliptical integrals to handle Euler-Bernoulli beam equation [1], and then Chen proposed a comprehensive elliptical solution for large beam-deflection problems for CMs [12]; Pseudo rigid body model, termed as PRBM, was also used to model Euler-Bernoulli beam where the deflection of a slender beam is physically represented by that of a rigid mechanism [1]; Awtar linearized Euler-Bernoulli beam equation for the cases where slender beams are deflected within an intermediate range, resulting in the commonly used Beam Constrained Model (BCM) [13]; chained BCM [14] and chained PRBM [15] are both based on Euler-Bernoulli beam theory to take care of large deflections of slender beams; in FEM-based commercial software, Euler-Bernoulli beam differential equations can be directly solved by FEM. Recently, weighted residual methods have also been proposed to solve the boundary value problems that are formulated by Euler Bernoulli beam theory [8] [16].

3. **Timoshenko beam theory:** in the state-of-the-art literature of CMs, there are researchers who have completed some impressive work regarding Timoshenko beam theory [11]. For the intermediate range of deflection, Guimin extended Awtar's work (BCM) by taking into account the shear strain of cross sections via Timoshenko beam theory [17], called Timoshenko beam constrained model (TBCM). For large beam-deflection problems, chained-TBCM is logically a promising option but no one has ever touched it to the best of our knowledge in the current literature. Besides, FEM is also a good tool to implement the governing BVPs behind Timoshenko beam theory for modeling large deflection of beams.

As just stated above, solid mechanics is used to analyze the general deformation of the studied object so the bending, axial elongation, shear deformation of the flexible beam and other possible deformation in beam deflections could be all considered. However, FEM normally develops a high-dimensional model for the studied object, which is logically tedious to solve. This explains why solid-mechanics-based FEM is more used as a verification method (accurate) than an optimization tool (time-consuming) in the field of CMs [8].

To increase the efficiency of beam-modeling process, geometrically nonlinear Euler Bernoulli beam theory is commonly adopted where just one ordinary differential equation can govern the large planar deflection of a slender beam [16]. It's straightforward and simple since it only considers the

<figure><img src="image_2.png" /><figcaption>Fig. 1: Small (linear) deflection of Timoshenko beam compared to Euler-Bernoulli beam subjected to beam-end loading</figcaption></figure>

major contribution of beam deflection: bending. Besides, in some cases, the axial elongation and shear deformation of the studied beam cannot be neglected especially when designing high-precision CMs or bi-stable mechanisms [17]. For example, compliant parallelograms work under some stretching payload [4], and the forward critical force needs to be precisely characterized in bi-stable mechanisms [17].

In previous studies, linear Timoshenko beam theory is normally used for analyzing the small deflection of stubby straight beams since the shear strain of cross sections make rather important contributions in these cases. In this paper, we utilize geometrically nonlinear Timoshenko beam theory to handle the bending and shear strains of cross sections in large deflections of straight beams and ICBs. To take the axial elongation into account, we take one step further to add the linearly elastic assumption along the beam axis to the Timoshenko beam equations.

In this paper, we use Galerkin method to solve the governing equations derived from the proposed modified geometrically nonlinear Timoshenko beam theory for beam deflections [8]. Via this modeling strategy, we will be able to synthesize CMs as well as conduct model-based topology optimization for CMs more efficiently.

## 1.3 Structure of the paper

In the following part, we first recall the linearized Timoshenko beam theory and present the deduction of the modified geometrically nonlinear Timoshenko beam theory where large beam-deflection problems can be formulated as a non-linear BVP. Then, we focus on the single-beam deflection of a straight beam and a circularly curved beam, via Euler Bernoulli beam theory, Timoshenko beam theory and solid mechanics to analyze the contributions of shear and axial strains in beam deflections. Finally, we proceed to model 2 representative CMs: a straight-beam-based bi-stable mechanism and an ICB-based compliant parallelogram to the prove the feasibility of the proposed modeling strategy.

# 2 Modeling of large deflection via geometrically nonlinear Timoshenko beam theory

## 2.1 Recall of linear Timoshenko beam theory

First of all, we would like to recall the linear Timoshenko beam theory, the formulation of which is presented below with its corresponding graphical demonstration shown
in Fig. 1:

$$
\begin{align*}
\frac{d^3 y}{dx^3} &= \frac{d^2 \theta}{dx^2} = \frac{V(x)}{EI} = \frac{\frac{dM(x)}{dx}}{EI} \\
\frac{d^2 y}{dx^2} &= \frac{d \theta}{dx} = \frac{M(x)}{EI} \\
\frac{dy}{dx} &= \theta(x) = \int_0^x \frac{M(x)}{EI} dx - \frac{k_t V(x)}{GA} = \int_0^x \frac{M(x)}{EI} dx - \frac{k_t \frac{dM(x)}{dx}}{GA} \tag{1}
\end{align*}
$$

with

$$
M(x) = M_o + F_x b + F_y (a - x) \quad (2)
$$

where E denotes the Young's modulus; I denotes the second
moment of inertia of the cross section ($I = \frac{1}{12}h^3w$); A denotes the area of the cross section ($A = wh$); G denotes the
shear modulus ($G = \frac{E}{2(1+v)}$ for isotropic materials); $k_t$ denotes Timoshenko shear coefficient [18] ($k_t = \frac{12+11v}{10+10v}$) and
v is the Poisson's ratio of the material; $\frac{d^2y}{dx^2}$ denotes the linearized curvature of the beam; $V(x)$ denotes the shear force at the
cross section; $M(x)$ denotes the moment exerted at point x due to external forces and moments, which in this case is
the beam-end loading: $F_x$, $F_y$ and $M_o$ (see Fig. 1); L, h and
w denote the length, depth and width of the beam; a and b denote the coordinates of the beam end. Note that this linear
theory is only applicable to small deflections, equivalently
when $\frac{dy}{dx}$ is small, without considering geometric nonlinearity. Therefore, Eq. (2) can be simplified:

$$
M(x) = M_o + F_y(L-x) \quad (3)
$$

where $a \approx L$ and $b \approx 0$ in small deflections. Besides, $M_o$
should be carefully handled since a relatively larger $M_o$ can
easily break the assumption made for the linearized curva-
ture as well as small deflections. Although the effect of axial
force $F_x$ in bending is ignored as presented in Eq. (3), it
plays an important role in stress stiffening even in small de-
flections. Here, we proceed with $M_o$ being 0, which results
in a very common beam-deflection case, beam-end concen-
trated loading: $M(x) = F_y(L-x)$. Rearranging this equation
and Eq. (1), we will be able to arrive at the following:

$$
\begin{align}
\frac{dy}{dx} &= \theta(x) = \frac{k_t F_y}{GA} + \frac{F_y L}{EI} x - \frac{F_y}{2EI} x^2 \tag{4} \\
y(x) &= \frac{k_t F_y}{GA} x + \frac{F_y L}{2EI} x^2 - \frac{F_y}{6EI} x^3 \nonumber
\end{align}
$$

**Remark 1.** In Eq. (1), $\frac{M(x)}{EI}$ is the resulting curvature due to bending, which is also the major contribution to beam deflections. Different from Euler Bernoulli beam theory, it is assumed that the cross section of the beam is not normal to the beam axis any more since the shear strain of the cross section is considered. $\frac{k_t V(x)}{GA}$ denotes the extra rotation of the

<figure><img src="image_3.png" /><figcaption>Fig. 2: Large deflection of Timoshenko beam compared to Euler-Bernoulli beam under beam-end loading</figcaption></figure>

cross section due to the local shear force, which is taken into
account by linear Timoshenko beam theory. Thus, ignoring
the shear strain of the cross section, we can easily derive the
governing equations of linear Euler Bernoulli beam theory:

$$
\begin{align}
\frac{d^2 y}{dx^2} &= \frac{d \theta}{dx} = \frac{M(x)}{EI} = \frac{F_y(L-x)}{EI} \nonumber \\
\frac{dy}{dx} &= \theta(x) = \int_0^x \frac{M(x)}{EI} dx = \theta(x) = \frac{F_y L}{EI} x - \frac{F_y}{2EI} x^2 \tag{5} \\
y(x) &= \int_0^x \frac{F_y L}{EI} x - \frac{F_y}{2EI} x^2 dx = \frac{F_y L}{2EI} x^2 - \frac{F_y}{6EI} x^3 \nonumber
\end{align}
$$

As mentioned before, both linear Timoshenko beam the-
ory and linear Euler Bernoulli beam theory are valid only
if analyzing small deflections. Precisely, since the curva-
ture κ(x) = d²y/dx² [1 + (dy/dx)²]⁻³/₂, then by assuming the Y-
axis deflection of the beam is small, i.e. dy/dx ≈ 0, we have:

κ(x) = \frac{d^2 y}{dx^2} Similarly, x is chosen as the independent vari-
able due to the fact that the deformed beam axis still ap-
proximately coincides with x axis in small-deflection cases
(see Fig. 1), yielding Eqs. (1)(4) for Timoshenko beam and
Eq. (5) for Euler Bernoulli beam. It is clear that the men-
tioned previous work above (linear Timoshenko beam the-
ory and linear Euler Bernoulli beam theory) is all based on
the small-deflection assumption where the linearized formu-
lation of curvature is only valid for small deflections.

2.2 Geometrically nonlinear Timoshenko beam theory
with axial linearly elastic assumption

Here, we would like to start with geometrically nonlin-
ear Timoshenko beam theory. Dropping the assumption for
small deflections stated in Section 2.1, the formulation needs
to be transferred from the global coordinate system (see Fig.
1) to the local body frame (see Fig. 2) where the independent
variable is as well redefined from x to s. Then, we can still
consider the rotation angles (θ<sub>b</sub>(s) and θ<sub>s</sub>(s)) of each cross
section caused by bending and shear strains respectively. For
the curvature of beams due to bending, we have the same lin-
early elastic assumption as Euler Bernoulli beam equation:

$$
\frac{d\theta_b}{ds} = \frac{M(s)}{EI(s)} + \frac{1}{R(s)} \quad (6)
$$
where

$$M(s) = M_o + F_y(a - x(s)) + F_x(b - y(s)) \quad (7)$$

$I(s)$ implies the second moment of inertia varies along the beam axis and $R(s)$ implies the initial curvature varies along the beam axis as well. In terms of rotation angle caused by shear strain at each cross section, we have: $\theta_s(s) = \frac{k_t}{GA(s)}(F_y \cos\theta(s) + F_x \sin\theta(s))$ where $A(s)$ implies the cross sections vary along the beam axis. Therefore, the total rotation angle at each cross section is described as follows:

$$\theta(s) = \theta_b(s) + \theta_s(s) \quad (8)$$

If we take one step further to consider the linearly elastic elongation of each cross section, we have

$$\varepsilon_a(s) = \frac{F_y \sin\theta(s) - F_x \cos\theta(s)}{EA(s)} \quad (9)$$

For each cross section, its infinitesimal length can be reformulated as $(1 + \varepsilon_a(s))ds$. Therefore, we can arrive at the following by considering axial strain in deflected beams:

$$\frac{dx}{ds}(s) = (1 + \varepsilon_a(s)) \cos\theta(s); \quad \frac{dy}{ds}(s) = (1 + \varepsilon_a(s)) \sin\theta(s) \quad (10)$$

Differentiating Eq. (6) considering the relationships of Eqs. (7)(10), we can arrive at:

$$\begin{aligned} \frac{d^2\theta_b}{ds^2} &= \frac{\frac{dM(s)}{ds}EI(s) - M(s)E\frac{dI(s)}{ds}}{E^2I(s)^2} + \frac{\frac{dR(s)}{ds}}{R(s)^2} \\ &= \frac{\frac{dM(s)}{ds}}{EI} - \frac{M(s)}{EI(s)^2}\frac{dI(s)}{ds} + \frac{\frac{dR(s)}{ds}}{R(s)^2} \end{aligned} \quad (11)$$

$$\begin{aligned} \frac{d^2\theta_b}{ds^2} &= \frac{\frac{dM(s)}{ds}}{EI} - \frac{M(s)}{EI(s)^2}\frac{dI(s)}{ds} + \frac{\frac{dR(s)}{ds}}{R(s)^2} \\ &= -\frac{F_y}{EI}\frac{dx}{ds}(s) - \frac{F_x}{EI}\frac{dy}{ds}(s) + \frac{\frac{dR(s)}{ds}}{R(s)^2} \end{aligned} \quad (12)$$

Besides, there are 2 boundary conditions defined from the physical conditions: $\theta_b(0) = 0$ and $\frac{d\theta_b}{ds}(L) = \frac{M_o}{EI} + \frac{1}{R(L)}$. Rearranging Eqs. (8)(10)(12), we will be able to obtain the final formulation of the desired boundary value problem:

$$\begin{aligned} \text{D.E.} \quad \frac{d^2\theta_b}{ds^2} &= -\frac{F_y}{EI}\frac{dx}{ds}(s) - \frac{F_x}{EI}\frac{dy}{ds}(s) + \frac{\frac{dR(s)}{ds}}{R(s)^2} \\ \theta_s(s) &= \frac{k_t}{GA}(F_y \cos\theta(s) + F_x \sin\theta(s)) \end{aligned} \quad (13)$$

$$\theta(s) = \theta_b(s) + \theta_s(s)$$

$$\text{B.C.} \quad \theta(0) = 0, \quad \text{and} \quad \frac{d\theta_b}{ds}(L) = \frac{M_o}{EI} + \frac{1}{R(L)}$$

BVP (13) governs the deflection of a straight beam or an ICB under beam-end loading. As is stated above, we can also consider some other engineering scenarios in CMs, such as varying cross sections, varying initial curvature, distributed loading along the beam body and so on. It is proved that all the mentioned engineering cases can be formulated like Eq. (7) by rearranging its right-hand terms into the corresponding total moment $M(s)$ and defining $I(s)$, $A(s)$, $R(s)$ and other physical parameters as functions of $s$ [16]. In BVP (13), the boundary conditions define a beam is fixed at one end and is exerted by external forces and moment at another free end. Besides, other different beam-end conditions shown in Fig. 3 yield different modifications on BVP (13) (mainly the BCs) so different cases need to be studied individually. For char-

<figure><img src="image_4.png" /><figcaption>Fig. 3: Different types of beam ends</figcaption></figure>

acterization of the deflected beam shape, we derive the following via Eqs. (10):

$$x(s) = \int_{0}^{s} \frac{dx}{ds}(\xi)d\xi; \quad y(s) = \int_{0}^{s} \frac{dy}{ds}(\xi)d\xi \quad (14)$$

As shown in BVP (13) and Eq. (14), the proposed constitutive beam theory, which is termed as modified geometrically nonlinear Timoshenko beam theory incorporating axial elongation, takes care of bending, shear and axial strains under the linear elasticity framework, which theoretically improves the accuracy of the model compared to the classic Euler Bernoulli beam theory. In the following section, we will present the advantages through numerical modeling over other modeling strategies. For calculating the stored energy of the studied beam, we have the following formulation:

$$E_p = E_{pt} + E_{pa} \quad (15)$$

with $E_{pt} = \int_0^L \frac{1}{2}EI(\frac{d\theta}{ds})^2 ds$ and $E_{pa} = \int_0^L \frac{1}{2}EA(\frac{F_y \sin\theta(s) - F_x \cos\theta(s)}{EA})^2 ds$ where $E_{pt}$ denotes the stored energy due to bending and shear and $E_{pa}$ denotes the stored energy from axial strain. Note that the shear strain energy is not individually stated above since the shear influence has been added into the curvature change $\frac{d\theta}{ds}$ according to BVP (13).

### 3 Analysis of beam-deflection modeling and its numerical validation

In this section, we aim to present the enhancement of the model through several cases of beam-deflection numerical modeling. We first nondimensionalize the governing equation to discover its mathematical essence as well as dig into
its physical characteristics for better understanding of how different types of strains contribute to the final deflection. Note that the nondimensionalized BVPs or ODEs are mathematically and physically equivalent to the original ones, and the nondimensionalization is essentially just a linear transformation of ODEs or BVPs. Then, we logically modify the governing equation to reach a more straightforward (explicit) approximate solution to the numerical multi-solution trap of the original complex BVP (13). Finally, several numerical tests are conducted to numerically prove the improvement of the proposed modeling strategy compared to other models.

## 3.1 Prerequisites

First, we proceed to nondimensionalize the Timoshenko beam governing equation (13). Therefore, the dependent variable $\theta$ and the independent variable $s$ are nondimensionalized via:

$$ 
\begin{aligned}
\hat{\theta} &= \frac{\theta}{1} = \frac{\theta_b + \theta_s}{1}; \quad \theta = 1\hat{\theta}; \quad \theta_b = 1\hat{\theta}_b; \quad \theta_s = 1\hat{\theta}_s \\
\hat{s} &= \frac{s}{L}; \quad s = L\hat{s}
\end{aligned}
\qquad (16) $$

Note that $\theta$, $\theta_b$ and $\theta_s$ are already non-dimensional but we still mark $\hat{\theta}$, $\hat{\theta}_b$ and $\hat{\theta}_s$ as the nondimensionalized ones for universal notation in this paper. Similarly, we have the same case for $\varepsilon_a(s)$ in Eq. (9) and Eq. (10):

$$ \frac{d\hat{\varepsilon}}{d\hat{s}}(\hat{s}) = (1 + \hat{\varepsilon}_a(\hat{s})) \cos \hat{\theta}(\hat{s}); \quad \frac{d\hat{\varepsilon}}{d\hat{s}}(\hat{s}) = (1 + \hat{\varepsilon}_a(\hat{s})) \sin \hat{\theta}(\hat{s}) \qquad (17) $$

where $\hat{\varepsilon}_a(\hat{s}) = \frac{f_y \sin \hat{\theta}(\hat{s}) - f_x \cos \hat{\theta}(\hat{s})}{12} \eta^2$, with $\eta$ being the slenderness ratio defined as $\eta = \frac{h}{L}$. Then, we can nondimensionalize (13) using (16) accordingly:

$$ 
\begin{aligned}
\text{D.E.} \quad \frac{d^2 \hat{\theta}_b}{d\hat{s}^2} &= -f_y \frac{d\hat{\varepsilon}}{d\hat{s}}(\hat{s}) - f_x \frac{d\hat{\varepsilon}}{d\hat{s}}(\hat{s}) + \frac{\frac{dr(\hat{s})}{d\hat{s}}}{r(\hat{s})^2} \\
\hat{\theta}_s(\hat{s}) &= \frac{(12+11v)}{60} \eta^2 (f_y \cos \hat{\theta}(\hat{s}) + f_x \sin \hat{\theta}(\hat{s})) \quad (18) \\
\hat{\theta}(\hat{s}) &= \hat{\theta}_s(\hat{s}) + \hat{\theta}_b(\hat{s}) \\
\text{B.C.} \quad \hat{\theta}(0) &= 0, \quad \text{and} \quad \frac{d\hat{\theta}_b}{d\hat{s}}(1) = m_o + \frac{1}{r(1)}
\end{aligned}
$$

where $f_x = \frac{F_y L^2}{EI}$; $f_y = \frac{F_x L^2}{EI}$; $m_o = \frac{M_o L}{EI}$; $r(\hat{s}) = \frac{R(\hat{s})}{L}$; $\eta = \frac{h}{L}$; $\hat{s} \in [0, 1]$; $s \in [0, L]$. In terms of characterizing the nondimensionalized beam shape, we can use the following nondimensionalized formulation according to Eq. (17): $\hat{x}(\hat{s}) = \int_0^\hat{s} \frac{d\hat{s}}{d\hat{s}} d\hat{\xi}$ and $\hat{y}(\hat{s}) = \int_0^\hat{s} \frac{d\hat{y}}{d\hat{s}} d\hat{\xi}$ where $\hat{x}(\hat{s})$ and $\hat{y}(\hat{s})$ are the nondimensionalized coordinates along the beam axis. Essentially, the nondimensionalized BVP (18) governs all large deflection scenarios where different materials, different beam geometries and different loading are all considered. Under this framework of nondimensionalization, all physical parameters can be represented by nondimensionalized values. This

framework provides both feasible and logical prerequisites for the following sections.

**Remark 2.** As long as the isotropic material used is settled, BVP (18) gradually converges to the following if $\eta \to 0$:

$$ 
\begin{aligned}
\text{D.E.} \quad \frac{d^2 \hat{\theta}}{d\hat{s}^2} &= -(f_y \cos \hat{\theta}(\hat{s}) + f_x \sin \hat{\theta}(\hat{s})) + \frac{\frac{dr(\hat{s})}{d\hat{s}}}{r(\hat{s})^2} \\
\text{B.C.} \quad \hat{\theta}(0) &= 0 \\
\frac{d\hat{\theta}}{d\hat{s}}(1) &= m_o + \frac{1}{r(1)}
\end{aligned}
\qquad (19) $$

with the nondimensionalized formulation of characterizing the beam shape:

$$ \hat{x}(\hat{s}) = \int_{0}^{\hat{s}} \cos \hat{\theta}(\hat{\xi}) d\hat{\xi}; \quad \hat{y}(\hat{s}) = \int_{0}^{\hat{s}} \sin \hat{\theta}(\hat{\xi}) d\hat{\xi} \qquad (20) $$

The above BVP (19) with Eq. (20) is exactly the nondimensionalized governing equation of Euler Bernoulli beam theory [19]. We can easily draw the conclusion that higher $\frac{h}{L}$ causes larger shear influences in beam deflection, which agrees with the statement: shear influences can not be ignored when analyzing the deflection of stubby beams, so Timoshenko beam theory is commonly used in this case. Similarly, higher $\frac{h}{L}$ also results in larger influences of axial elongation as well according to Eq. (17) and Eq. (18).

**Remark 3.** Logically, if we ignore the shear in the studied beam, BVP (18) can be deduced as follows:

$$ 
\begin{aligned}
\text{D.E.} \quad \frac{d^2 \hat{\theta}_b}{d\hat{s}^2} &= -f_y \frac{d\hat{x}}{d\hat{s}}(\hat{s}) - f_x \frac{d\hat{y}}{d\hat{s}}(\hat{s}) + \frac{\frac{dr(\hat{s})}{d\hat{s}}}{r(\hat{s})^2} \\
\text{B.C.} \quad \hat{\theta}(0) &= 0, \quad \text{and} \quad \frac{d\hat{\theta}}{d\hat{s}}(1) = m_o + \frac{1}{r(1)}
\end{aligned}
\qquad (21) $$

BVP (21) along with Eq. (17) is termed as modified geometrically nonlinear Euler Bernoulli beam theory incorporating axial elongation.

## 3.2 Numerical modeling through Galerkin method

As stated in BVP (18) and BVP (19), the governing equations of beam deflection are essentially a highly nonlinear BVP. Compared to the simple BVP (19) (one nonlinear ODE subjected to 2 constraints), BVP (18) is a system of ODEs subjected to 2 constraints, which is more difficult to handle. Directly numerically approximating the solution of BVP (18) may result in multi solutions due to its highly nonlinearity. Therefore, we need a rather simple and straightforward strategy to solve the problem by making most of the physical information both in and behind BVP (18).

In the cases of large deflection, the dominant deformation of the beam lies in bending whereas the shear strain
along the beam is the secondary and way less important con-
tributor. Therefore, we can simplify the governing equa-
tion in BVP (18) by assuming $\theta_b \approx \theta$, equivalently $\hat{\theta}_b \approx \hat{\theta}$.
On the contrary, we closely look at BVP (18) again in the
small-deflection cases. Obviously, $\cos\theta_b \approx \cos\theta \approx 1$ and
$\sin\theta_b \approx \sin\theta \approx 0$ when the magnitude of $\theta$ is small. There-
fore, we can still assume $\theta_b \approx \theta$ in small-deflection cases.
Then, we can rewrite Eqs. (17):

$$
\begin{align}
\frac{d\hat{x}}{ds}(\hat{s}) &\approx \frac{d\hat{x}^*}{ds}(\hat{s}) = (1 + \hat{\epsilon}_a^*(\hat{s})) \cos \hat{\theta}_b(\hat{s}) \nonumber \\
\frac{d\hat{y}}{ds}(\hat{s}) &\approx \frac{d\hat{y}^*}{ds}(\hat{s}) = (1 + \hat{\epsilon}_a^*(\hat{s})) \sin \hat{\theta}_b(\hat{s}) \tag{22}
\end{align}
$$

where $\hat{\epsilon}_a^*(\hat{s}) = \frac{f_y \sin \hat{\theta}_b(\hat{s}) - f_x \cos \hat{\theta}_b(\hat{s})}{12} \eta^2$. Now, we can logically rewrite the governing BVP into the following:

$$
\begin{align}
\text{D.E.} \quad \frac{d^2 \hat{\theta}_b}{ds^2} &= -f_y \frac{d\hat{x}^*}{ds}(\hat{s}) - f_x \frac{d\hat{y}^*}{ds}(\hat{s}) + \frac{-\frac{dr(s)}{ds}}{r(s)^2} \nonumber \\
\hat{\theta}_s(\hat{s}) &= \frac{(12+11\nu)}{60} \eta^2 (f_y \cos \hat{\theta}_b(\hat{s}) + f_x \sin \hat{\theta}_b(\hat{s})) \tag{23}
\end{align}
$$

$$
\text{B.C.} \quad \hat{\theta}_b(0) = 0, \text{ and} \quad \frac{d\hat{\theta}_b}{ds}(1) = m_o + \frac{1}{r(1)}
$$

Essentially speaking, we bypass the need of implicitly ap-
proximating θs in BVP (18) to explicitly approximate θs
through BVP (23), which hugely decreases the complexity
of approximation process. As proposed in [8], we have up
to nearly 10 methods to treat this BVP. In this paper, we use
Galerkin method to handle BVP (23), and a brief recap is
presented here.

According to Weierstrass’s first theorem, a given differ-
entiable function for a random closed interval can be always
approximated by a polynomial of a proper order. Therefore,
we then define a function ˆθ<sub>b</sub>(s) which is an n<sup>th</sup>-order poly-
nomial, to approximate the solution of the studied BVP (23):

ˆθ<sub>b</sub>(s) = Σ<sub>i=0</sub><sup>n</sup> c<sub>i</sub>s<sup>i</sup> for i = 0, 1, ..., n where there are n + 1 un-
known coefficients to solve. Next, we can define an approx-
imating residual function of BVP (23) utilizing the polyno-
mial solution ˆθ<sub>b</sub>(s) as follows:

$$
R_e(\hat{s}) = \frac{d^2 \hat{\theta}_b}{ds^2} + f_y \frac{d\hat{x}^*}{ds}(\hat{s}) + f_x \frac{d\hat{y}^*}{ds}(\hat{s}) + \frac{\frac{dr(s)}{ds}}{r(s)^2} \quad (24)
$$

The essence of weighted residual method is to force the
residual Eq. (24) to zero in some average sense over the
domain: ∫10 WᵢRₑ(s)ds = 0. Here, we choose Wᵢ taken from
the derivative of the approximating function ̂θ₉(s) according
to Galerkin method: Wᵢ(s) = ∂θ₉(s) / ∂cᵢ = sᵢ for i = 1, 2, 3...n + 1.
Besides, ̂θ₉(s) also needs to satisfy the boundary conditions
defined in BVP (23): ̂θ₉(0) = 0 and dθ₉/ds(1) = m₀ + 1/r(1). So
far, we have altogether n + 3 equations to solve n + 1 poly-
nomial coefficients cᵢ (i = 0, 1, 2, 3...n) via Newton-Raphson

method. Like all numerical methods, the initial guess is
of great importance in terms of convergence using Newton-
Raphson method. Here, we provide 2 ways to initialize the
unknowns.

1. **Direct way:** we randomly initialize the values of $c_i$, which yields a random shape of $\hat{\theta}_b(\hat{s})$. The direct way is suitable for the scenario where the problem tends to have only one single solution, i.e., it contains only global minima, such as general large deflection of beams [16].

2. **Indirect way:** we choose a globally initial guess for $\hat{\theta}_b(\hat{s})$, such as trigonometric functions, which yields the corresponding $c_i$ for the $\hat{\theta}_b(\hat{s})$. The indirect way goes for the scenario where we reasonably initialize the initial guess (possibly near the final solution) for the final solution via using the existing information, such as nonlinear buckling of beams where there are normally multisolutions for a given load [20].

As long as we have the approximate solution $\hat{\theta}_b(\hat{s})$ available, we can explicitly obtain the approximation of $\hat{\theta}_s(s)$ through the closed-form formulation below according to BVP (23):

$$
\hat{\theta}_s(\hat{s}) = \frac{(12 + 11\nu)}{60} \eta^2 (f_y \cos \hat{\theta}_b(\hat{s}) + f_x \sin \hat{\theta}_b(\hat{s})) \quad (25)
$$

Finally, we can reach the modified approximate solution via:

$$
\hat{\theta}(\hat{s}) = \hat{\theta}_b(\hat{s}) + \hat{\theta}_s(\hat{s}) \quad (26)
$$

Finally, replacing $\hat{\theta}(\hat{s})$ with $\hat{\theta}(\hat{s})$, we can just follow Eqs. (17)(26) to characterize the deflected beam shape.

3.3 Deflection analysis on shear strain and axial elonga-
tion

In this section, 2 cases (straight beams and circularly
curved beams) are studied via different physical models. In
Table 1, we have provided different models implemented by
different numerical or analytical methods. These models can
be classified into 3 major groups: Euler Bernoulli beam mod-
els, Timoshenko beam models and solid mechanics. All the
models are used to compare the results of the proposed mod-
eling strategy (NTBA-GM), BVP (23) where solid mechan-
ics serves as a verification theory. Then, we will be able to
figure out the importance of shear and axial strains in beam
deflections and why we need to take care of them under some
specific modeling scenarios. As stated in Section 3.1, the
slenderness ratio η serves as the key contributor of shear and
axial strains to beam deflections. Therefore, we would like
to both validate and reveal how η affects the deflection of
beams via plotting the η-ŷ relationship. The loading condi-
tion is schematically demonstrated in Fig. 4.

Besides, under this nondimensionalization framework,
the analytical solutions of both linear Timoshenko and Eu-
ler Bernoulli beam theories Eq. (4) and Eq. (5) need to be
nondimensionalized for more general comparisons as well.
<table><caption>Table 1: Acronyms of different beam models</caption><thead><tr><th></th><th>Description</th><th>Governing equations</th><th>Model considerations</th><th>Refs</th></tr></thead><tbody><tr><td>EB-AS</td><td>Linear Euler Bernoulli beam theory<br>(Analytical solution)</td><td>Eq. (29)<br>Eq. (30)</td><td>Bending strain<br>within small-deflection range</td><td>[21]</td></tr><tr><td>NEB-GM</td><td>Geometrically nonlinear Euler Bernoulli beam theory<br>(Galerkin method)</td><td>BVP (19)</td><td>Bending strain<br>within large-deflection range</td><td>[16]</td></tr><tr><td>EBA-FEM</td><td>Linear Euler Bernoulli beam theory<br>incorporating axial elongation<br>(Finite element method)</td><td>See<br>reference</td><td>Bending and axial strains<br>within large-deflection range</td><td>[22]</td></tr><tr><td>NEBA-GM</td><td>Geometrically nonlinear Euler Bernoulli beam theory<br>incorporating axial elongation<br>(Galerkin method)</td><td>BVP (21)</td><td>Bending and axial strains<br>within large-deflection range</td><td>N/A</td></tr><tr><td>TB-AS</td><td>Linear Timoshenko beam theory<br>(Analytical solution)</td><td>Eq. (27)<br>Eq. (28)</td><td>Bending and shear strains<br>within small-deflection range</td><td>[23]</td></tr><tr><td>TBA-FEM</td><td>Linear Timoshenko beam theory<br>incorporating axial elongation<br>(Finite element method)</td><td>See<br>reference</td><td>Bending, shear and axial strains<br>within large-deflection range</td><td>[22]</td></tr><tr><td>NTBA-GM</td><td>Geometrically nonlinear Timoshenko beam theory<br>incorporating axial elongation<br>(Galerkin method)</td><td>BVP (23)</td><td>Bending, shear and axial strains<br>within large-deflection range</td><td>N/A</td></tr><tr><td>SD-FEM</td><td>Solid mechanics<br>(FEM)</td><td>See<br>reference</td><td>All sorts of deformation<br>within large-deflection range</td><td>[22]</td></tr></tbody></table>

Apart from SD-FEM, the constitutive equations of all beam theories are under the framework of linear elastic assumption.

Using the same strategy presented before for nondimension- alization, Eq. (4) can be rewritten as the following:

$$ \begin{aligned} \frac{d\hat{y}}{d\hat{x}}(\hat{x}) &= \hat{\theta}(\hat{x}) = \frac{12+11\nu}{60}\eta^2 f_y + f_y\hat{x} - \frac{1}{2}f_y\hat{x}^2 \\ \hat{y}(\hat{x}) &= \frac{12+11\nu}{60}\eta^2 f_y\hat{x} + \frac{1}{2}f_y\hat{x}^2 - \frac{1}{6}f_y\hat{x}^3 \end{aligned} \quad (27) $$

To evaluate the beam-end rotation and vertical coordinate, we set $\hat{x} = 1$:

$$ \begin{aligned} \frac{d\hat{y}}{d\hat{x}}(1) &= \hat{\theta}(1) = f_y\left(\frac{12+11\nu}{60}\eta^2 + \frac{1}{2}\right) \\ \hat{y}(1) &= f_y\left(\frac{12+11\nu}{60}\eta^2 + \frac{1}{3}\right) \end{aligned} \quad (28) $$

which is termed as TB-AS. Likewise, the analytical solution of linear Euler Bernoulli beam theory, Eq. (5) can be arranged below under nondimensionalization:

$$ \frac{d\hat{y}}{d\hat{x}} = \hat{\theta}(\hat{x}) = f_y\hat{x} - \frac{1}{2}f_y\hat{x}^2, \quad \text{and} \quad \hat{y}(\hat{x}) = \frac{1}{2}f_y\hat{x}^2 - \frac{1}{6}f_y\hat{x}^3 \quad (29) $$

If $\hat{x} = 1$, we can end up with the beam-end rotation and ver-

tical coordinate via Eq. (29):

$$ \frac{d\hat{y}}{d\hat{x}}(1) = \hat{\theta}(1) = f_y - \frac{1}{2}f_y = \frac{1}{2}f_y, \quad \text{and} \quad \hat{y}(1) = \frac{1}{2}f_y - \frac{1}{6}f_y = \frac{1}{3}f_y \quad (30) $$

which is termed as EB-AS.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Studied loading conditions</figcaption></figure>

## 1. Deflection analysis of a straight beam

Here, we assume $v = 0.33$, $f_x = 0$, $f_y = 0.2$, $m_o = 0$, $\frac{1}{r} = 0$, $\frac{w}{L} = 0.1$ and $\eta = \frac{h}{L} = \frac{iL}{L} = i$ ($i = 0.1, 0.2, 0.3...1$).

As shown in Fig. 4a, a straight beam is fixed at one end and subjected to a transverse load at another free end
<figure><img src="image_6.png" /><figcaption>Fig. 5: Deflections of straight beams subjected to pure transverse forces</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6: Deflections of circularly curved beams subjected to pure vertical forces</figcaption></figure>

where $\hat{y}(1)$ is evaluated. Since the deflection is evaluated within a small range, EB-AS and TB-AS are used for comparisons as well. In Fig. 5, it can be noticed that all results of Timoshenko theory (TBA-FEM, NTBA-GM and TB-AS) present the impact of shear strain in the deflection of straight beams compared to the results of Euler Bernoulli beam theory (EBA-FEM, NEB-GM, NEBA-GM and EB-AS). As $\eta$ increases, the impact of shear strain increases as well, which perfectly agrees with the statement that shear cannot be neglected in stubby beams. The results of TBA-FEM, NTBA-GM and TB-AS capture the same trend as that of SD-FEM, and the errors are mainly due to the common assumption of beam theories: plane stress assumption. Generally speaking, if $w$ gets smaller, the results of SD-FEM will converge to that of TBA-FEM and NTBA-GM since the strict plane stress assumption is gradually satisfied in SD-FEM. However, the axial elongation doesn't play a major role in the studied case of a straight beam since the transverse load $f_y$ is very much normal to the beam axis.

## 2. Deflection analysis of a circularly curved beam

Similarly, we assume $v = 0.33$, $f_x = 0$, $f_y = 0.2$, $m_o = 0$, $r = \frac{2}{\pi}$, $\frac{w}{L} = 0.1$ and $\eta = \frac{h}{L} = \frac{iL}{L} = i$ ($i = 0.1, 0.2, 0.3...1$).

As shown in Fig. 4b, a circularly curved beam is fixed at one end and subjected to a vertical load at another free end where $\hat{y}(1)$ is evaluated.

There are 3 groups of results numerically obtained from different models: NEB-GM; EBA-FEM and NEBA-GM; TBA-FEM, NTBA-GM and SD-FEM presented in Fig. 6. Only considering bending of the beam (NEB-GM) results in easily observable errors to SD-FEM in the deflections of circularly curved beams. As is mentioned, EBA-FEM considers axial strain whereas NEB-GM does not, resulting in axial stretching caused by the vertical load $f_y$. As $\eta$ increases, $\hat{y}(1)$ increases as well, which numerically validates the theoretical conclusion drawn from BVP (18): larger $\eta$ causes larger impact of axial elongation. The closest results to SD-FEM are logically the ones of the models considering bending, shear and axial strains (TBA-FEM and NTBA-GM), which agrees with results in Fig. 6. In particular, the results of TBA-FEM, NEBA-GM and SD-FEM almost cover each other because the model of SD-FEM more or less meets the plane stress assumption. Logically, if the ratio $\frac{w}{L}$ get larger, the differences among TBA-FEM, NEBA-GM and SD-FEM will get larger. This is because the plane stress assumption will not be valid anymore and tend to transform into the plane strain assumption if $\frac{w}{L}$ reaches a very large value.

## 4 Mechanism synthesis: bi-stable compliant mechanisms and compliant parallelograms

In this section, we aim to model 2 types of CMs: bi-stable compliant mechanisms and compliant parallelograms under the proposed modeling framework BVP (13) where bending, shear and axial strains are all considered.

To verify the feasibility of the proposed modeling strategy in modeling straight beams and ICBs, a straight-beam-based compliant mechanism and an ICB-based compliant parallelograms are chosen to be synthesized as examples here. Following the basic modeling principles of solid mechanics, 3 types of governing equations need to be developed, which are constitutive equations of flexible members, force-equilibrium equations and geometric-relationship equations. Numerical results are provided with corresponding analysis both on the characteristics of the studied CMs and the choices of modeling strategies.

### 4.1 Bi-stable compliant mechanisms

Bi-stable compliant mechanisms experience nonlinear buckling and nonlinear post-buckling when operating, which is difficult to model. In this section, a typical type of lumped bi-stable compliant mechanisms is modeled where the built-in flexible members are modeled BVP (13). This means, besides bending, the axial and shear strains are both considered. In terms of the numerical implementation, the explicitly numerical strategy for nondimensional modeling to handle BVP (23) is used here for dimensional modeling.

As shown in Fig. 7, the bi-stable mechanism is com-
<figure><img src="image_8.png" /><figcaption>Fig. 7: Schematic diagram of bi-stable compliant mechanisms</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8: Detailed force diagram of the classic compliant bi-stable mechanism</figcaption></figure>

posed of 2 rigid links, 1 rigid motion stage and 4 flexible beams. State 1 refers to its fabricated shape where the mechanism is symmetrically fabricated and fixed at the two sides. The detailed geometric parameters are all provided in Fig. 7. The mechanism has 3 states where State 1 and State 3 are the stable equilibrium states whereas State 2 is the unstable equilibrium state. Subjected to the external force $F_o$, the mechanism works from State 1 through State 2 finally arriving at State 3. Here, we are interested in the relationship between $F_o$ and the displacement $\Delta Y$ of the reference point R. Since the structure of the bi-stable mechanism is symmetrical, we can simplify the modeling by working on the half of it as shown in Fig. 8.

### 4.1.1 Static modeling

Normally, we start with the **constitutive equations** of each flexible member first. According to BVP (13), we have the following governing equation for Beam #1:

$$ 
\begin{aligned}
\text{D.E.} \quad \frac{d^2\theta_{b1}}{ds^2} &= -\frac{F_{y1}}{EI} \frac{dx_1}{ds}(s) - \frac{F_{x1}}{EI} \frac{dy_1}{ds}(s) \\
\theta_{s1}(s) &= \frac{k_f}{GA} (F_{y1} \cos\theta_1(s) + F_{x1} \sin\theta_1(s)) \\
\theta_1(s) &= \theta_{b1}(s) + \theta_{s1}(s)
\end{aligned}
\quad (31) $$

$$ \text{B.C.} \quad \theta_1(0) = \theta_o \quad \text{and} \quad \frac{d\theta_{b1}}{ds}(L_1) = \frac{M_{o1}}{EI} $$

For Beam #2, we have a similar BVP as follows:

<figure><img src="image_10.png" /><figcaption>Fig. 9: Force-displacement relationship of the studied bi-stable mechanism</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10: Elastic energy stored-displacement relationship of the studied bi-stable mechanism</figcaption></figure>

$$ 
\begin{aligned}
\text{D.E.} \quad \frac{d^2\theta_{b2}}{ds^2} &= -\frac{F_y}{EI} \frac{dx_2}{ds}(s) - \frac{F_x}{EI} \frac{dy_2}{ds}(s) \\
\theta_{s2}(s) &= \frac{k_f}{GA} (F_y \cos\theta_2(s) + F_x \sin\theta_2(s)) \\
\theta_2(s) &= \theta_{b2}(s) + \theta_{s2}(s)
\end{aligned}
\quad (32) $$

Then, we can characterize the deflected shapes of Beam #1 and Beam #2 in their local coordinate systems respectively according to Eq. (14): $x_1(s) = \int_0^s \frac{dx_1}{d\xi}(\xi)d\xi$, $y_1(s) = \int_0^s \frac{dy_1}{d\xi}(\xi)d\xi$, $x_2(s) = \int_0^s \frac{dx_2}{d\xi}(\xi)d\xi$ and $y_2(s) = \int_0^s \frac{dy_2}{d\xi}(\xi)d\xi$. Besides, we also have the following **force-equilibrium relationships**: $F_x = F_{x1}$, $F_y = F_{y1}$, $F_o = 2F_y$ and $M_{o1} = M_o + F_y(L_r \cos\theta_1(L_1) + x_2(L_2)) + F_x(L_r \sin\theta_1(L_1) + y_2(L_2))$ and **the geometric relationships**: $(L_1 + L_r + L_2) \cos\theta_o = x_1(L_1) + L_r \cos\theta_1(L_1) + x_2(L_2)$ and $\Delta Y = (L_1 + L_r + L_2) \sin\theta_o - (y_1(L_1) + L_r \sin\theta_1(L_1) + y_2(L_2))$ where $\Delta Y$ denotes the displacement of the reference point R. Finally, to calculate the energy stored in this mechanical system along its workspace, we just need to follow Eq. (15).

### 4.1.2 Numerical results

To verify the modeling methodology, we use solid-mechanics-based FEM as the reference. First, the geometric
and material parameters chosen for the studied straight beam is shown in the following: $E = 200 \times 10^9$ Pa; $v = 0.30$; $w = 0.015$ m; $h = 0.005$ m; $L_1 = 0.25$ m; $L_r = 0.25$ m; $L_2 = 0.25$ m; $\theta_o = 0.1266$ rad; $k_t = \frac{12+11v}{10+10v}$; $G = \frac{E}{2(1+v)}$; To analyze the mechanical performance of bi-stable mechanisms, the relationship between the force and the displacement of point of force is of great importance. Therefore, $F_o - \Delta Y$ is plotted in Fig.9 for further analysis. Besides, we also present the relationship between the elastic energy $E_p$ stored in this system and the displacement $\Delta Y$ as well (see Fig.10).

### 4.1.3 Results analysis

In Fig. 9 and Fig. 10, many interesting insights behind the plotting can be discovered. Clearly, the results in these 2 figures both verify the existence of 2 stable equilibrium states (State 1 and State 3) and 1 unstable equilibrium state (State 2). In Fig. 9, 3 equilibrium states lie in the line $F_o = 0$ that implies at these 3 states the mechanism can stay equilibrium without external inputs. States 1 and 3 have positive stiffness while State 2 has negative stiffness, which means States 1 and 3 are stable whereas State 2 is not. According to the minimum total potential energy principle, the potential energy of States 1 and 3 reaches the local minimum as well as stay in a stable equilibrium state (Fig. 10). On the contrary, State 2 stays at the unstable state since its potential energy reaches the local maximum. Note that since $\theta_o$ is small in this numerical example, reference point R will experience a large axial force change within a small displacement. Therefore, the force-displacement curve near the State 1 is steeper as shown in Fig. 9. Logically, if $\theta_o$ is large, reference point R will encounter a large axial force change but within a large displacement, which means the force-displacement curve would be curvier ($F_o$ would increase slowly in this case) [6] [17] [20].

If we first take a look at the graphical results in Fig. 9 and Fig. 10, NEB-GM obviously is the most inaccurate one since it only considers bending. However, TBA-FEM and NTBA-GM take into account shear and axial strains, and their curves share the same trend with limited errors except for the turning points of the forces at left top corner of Fig. 9, which implies the following two points:

1. Model difference: in pure-compliant bistable mechanisms, the involved slender beams experience large axial loading, resulting in highly nonlinear buckling and post-buckling behaviors. Therefore, in this special case, the beam theory used NTBA-GM may not be able to capture the whole characteristics of the deformation during the operation of the bi-stable mechanisms. This means the linearly elastic assumption in considering bending, shear and axial strains may not be accurate enough to handle this scenario.

2. Numerical traps: as stated in [20], the mathematical essence of solving nonlinear buckling and post-buckling problems lies in multi solutions of ODEs. Supposing the model difference is small and can be neglected, the errors may possibly come from different numerical solutions that all satisfy the same governing equation. Practically speaking, for example, a 2nd-order post buckling

case and a 1st post buckling case may fall in a same numerical solution given a certain convergence criteria. It is very much possible since in Fig. 10, the corresponding results of potential energy from TBA-FEM and SD-FEM within the range of turning points of forces has very minor differences. This means, in a more general framework (for example, the principle of minimum potential energy in Lagrange mechanics), there may be several local minima just near the correct answer.

Then, we move to the comparison between the results of TBA-FEM and the proposed NTBA-GM. Mathematically speaking, TBA-FEM and NTBA-GM both follow the linearly elastic assumption in bending, shear and axial strains, which means the governing equations behind these two are basically the same. The only difference is that in TBA-FEM, linear Timoshenko beam theory is used via FEM to model the deflected beam and to locally approximate the solution whereas NTBA-GM directly uses a highly nonlinear BVP (BVP (13)) to model the deflected beam then to globally approximate the solution. Similarly, the gap between the results of TBA-FEM and NTBA-GM may possibly come from numerical traps stated above as well. Besides, as shown in Fig. 10, the 4 energy curves are closer to each other (in particular, the NTBA-GM curve almost covers that of SD-FEM), which also validates our opinion that the errors may possibly come from the multi-solution trap. However, this is a big topic out of the scope of this paper.

## 4.2 Compliant parallelograms

### 4.2.1 Geometric definition

In this section, we aim to model an ICB-based compliant parallelogram. To prove the feasibility of the proposed modeling strategy for ICBs of random shapes, we here define 2 irregular curves for the beam axes of 2 flexible beams, serving as the left compliant leg (Beam #1) and right compliant leg (Beam #2) of the studied compliant parallelogram (see Fig. 11). The geometries of the left beam and right beam are respectively defined in their corresponding local frames ($X_1-O_1-Y_1$ and $X_2-O_2-Y_2$) below: $\theta_{g1}(s) = \frac{\pi}{2} - 301.5929s^3$ and $\theta_{g2}(s) = \frac{\pi}{2} + 301.5929s^3$. Following Eq. (14), we can easily characterize the fabricated beam shapes. Logically, we can derive the initial curvature of the beam axes via: $\frac{d\theta_{g1}}{ds}(s) = -904.7787s^2 = \frac{1}{R_1(s)}$ and $\frac{d\theta_{g2}}{ds}(s) = 904.7787s^2 = \frac{1}{R_2(s)}$. In

<figure><img src="image_12.png" /><figcaption>Fig. 11: Schematic diagram of the studied ICB-based compliant parallelogram</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12: Detailed force diagram of the studied ICB-based compliant parallelogram</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13: Deflection range of the studied ICB-based compliant parallelogram</figcaption></figure>

Fig. 11, the ICB-based compliant parallelogram is composed of the mentioned 2 flexible beams and a rigid motion stage. Here, we are interested in the relationship between $F_o$ and the displacement $\Delta Y$ of the reference point R. Since the structure of the compliant parallelogram is symmetrical, we can simplify the modeling by working on the half of it (see Fig. 12).

## 4.2.2 Static modeling

Here, we proceed to model the simplified and equivalent mechanism presented in Fig. 12. Logically, according to BVP (13), Beam #1 can be modeled via the following BVP: Then, we can characterize the deflected shapes of Beam #1 and Beam #2 in their local coordinate systems respectively according to Eq. (14):

$$ \text{D.E.} \quad \frac{d^2\theta_{b1}}{ds^2} = -\frac{F_y}{EI} \frac{dx_1}{ds}(s) - \frac{F_x}{EI} \frac{dy_1}{ds}(s) - \frac{\frac{dR_1(s)}{ds}}{R_1(s)^2} $$

$$ \frac{1}{R_1(s)} = \frac{d\theta_{o1}}{ds}(s) = -904.7787s^2 $$

$$ \theta_{s1}(s) = \frac{k_t}{GA} (F_y \cos \theta_1(s) + F_x \sin \theta_1(s)) $$

$$ \theta_1(s) = \theta_{b1}(s) + \theta_{s1}(s) $$

$$ \text{B.C.} \quad \theta_1(0) = \theta_{g1}(0) = \frac{\pi}{2} $$

$$ \frac{d\theta_{b1}}{ds}(L_1) = \frac{M_o}{EI} + \frac{1}{R_1(L_1)}; \theta_1(L_1) = \theta_{g1}(L_1) = -\pi \quad (33) $$

Then the following geometric constraints should be satisfied:

<figure><img src="image_15.png" /><figcaption>Fig. 14: Force-displacement relationship of the studied ICB-based parallelogram mechanism</figcaption></figure>

$$ \text{fied: } \int_{0}^{L_1} \frac{dx_1}{ds}(s)ds = \int_{0}^{L_1} \cos \theta_{o1}(s)ds \text{ and } \int_{0}^{L_1} \cos \theta_{o1}(s)ds - \int_{0}^{L_1} \frac{dy_1}{ds}(s)ds = \Delta Y. \quad \text{The force equilibrium in the rigid motion stage is formulated: } F_y = \frac{E_o}{2}. $$

## 4.2.3 Numerical results

To verify the modeling methodology, we use solid-mechanics-based FEM as the reference. First, the geometric and material parameters chosen for the studied ICB is shown in the following: $fE = 200 \times 10^9$ Pa, $v = 0.30$, $w = 0.015$ m, $h = 0.005$ m, $L_1 = 0.25$ m, $k_t = \frac{12+11v}{10+10v}$, $G = \frac{E}{2(1+v)}$. To analyze the performance of the studied parallelogram mechanism, the relationship between the force $F_o$ and its corresponding displacement $\Delta Y$ plays an important role, which is plotted in Fig. 14. The deflection range of the studied parallelogram mechanism is graphically presented in Fig. 13.

## 4.2.4 Results analysis

Here we used a varying-initial-curvature case, so the initial radius of the beam is also varying within its beam length. Equivalently, there are infinitely different values of the initial radius to be tested here. As shown in Fig. 13 and Fig. 14, the proposed NTBA-GM has been both verified and proved feasible in predicting large deflection of ICBs as well as conducting mechanism synthesis. Clearly as presented in Fig. 14, the results of NTBA-GM and TBA-FEM almost cover each other since the models used both consider bending, shear and axial strains. The gap between the mentioned two and SD-FEM may mainly lie in the assumption made for all beam theories: plane stress assumption and the different convergence thresholds set up for the used numerical methods to solve the corresponding governing equations. In terms of the characteristics of the ICB-based parallelogram mechanism, the major advantage of this type is that the mechanism will not be over-constraint during its operation due to the involved ICBs. Then, if we reformulate the mechanism design problem into an optimization problem where the geometry of the ICB as the optimized variables, it will not be hard to achieve constant stiffness along its linear workspace or other desired performance objectives.
## 5 Conclusions

As a promising option for the existing mechanical applications, compliant mechanisms display many advantages due to their inherent characteristics, such as compliance and the increased accuracy. Obviously, the built-in flexible beams in CMs play a very important role in realizing the transfer of motion, force and energy. In this paper, we propose to adopt the modified Timoshenko beam theory including the considerations of bending, shear and axial strains to model large deflection of flexible beams and use weighted residual methods to solve the governing BVPs. To validate the modeling strategy, the numerical testing on a single beam and mechanism synthesis is presented along with detailed analysis. The results have proved the proposed modeling strategy feasible and straightforward in modeling CMs. In our future work, we will focus on deriving simple and straightforward 3-D beam models for designing and optimizing CMs.

## Acknowledgment

This work is partially supported by the National Natural Science Foundation of China (Grant No. 62073081), by the Project of Department of Education of Guangdong Province (Grant No. 2019KZDXM037), by the fund of Guangdong-Hong Kong-Macao Joint Laboratory for Intelligent Micro-Nano Optoelectronic Technology (No. 2020B1212030010), by the project ROBOCOP [ANR-19-CE19], by the project COSSROOTS [ANR-20-CE33], by the project Inventor (I-SITE ULNE, le programme d'Investissements d'Avenir, Métropole Européenne de Lille), France.

## References

[1] Larry L Howell. Compliant mechanisms. In *21st century kinematics*. Springer, London, 2013.

[2] D Farhadi Machekposhti, N Tolou, and JL Herder. A review on compliant joints and rigid-body constant velocity universal joints toward the design of compliant homokinetic couplings. *Journal of Mechanical Design*, 137(3), 2015.

[3] Nicolae Lobontiu. *Compliant mechanisms: design of flexure hinges*. CRC press, Boca Raton, 2020.

[4] Shorya Awtar. *Synthesis and analysis of parallel kinematic XY flexure mechanisms*. PhD thesis, Massachusetts Institute of Technology, 2003.

[5] Çağıl Merve Tanık, Volkan Parlaktaş, Engin Tanık, and Suat Kadioğlu. Steel compliant cardan universal joint. *Mechanism and Machine Theory*, 92:171–183, 2015.

[6] Jin Qiu, Jeffrey H Lang, and Alexander H Slocum. A curved-beam bistable mechanism. *Journal of microelectromechanical systems*, 13(2):137–146, 2004.

[7] Brian Trease and Sridhar Kota. Design of adaptive and controllable compliant systems with embedded actuators and sensors. *Journal of Mechanical Design*, 131(11), 2009.

[8] Ke Wu and Gang Zheng. Insight into numerical solutions of static large deflection of general planar beams

for compliant mechanisms. *Mechanism and Machine Theory*, 169:104598, 2022.

[9] Benliang Zhu, Xianmin Zhang, Hongchuan Zhang, Junwen Liang, Haoyan Zang, Hai Li, and Rixin Wang. Design of compliant mechanisms using continuum topology optimization: A review. *Mechanism and Machine Theory*, 143:103622, 2020.

[10] Olgierd Cecil Zienkiewicz, Robert Leroy Taylor, Robert Leroy Taylor, and Robert Lee Taylor. *The finite element method: solid mechanics, volume 2*. Butterworth-heinemann, 2000.

[11] Stephen P Timoshenko and James N Goodier. Theory of elasticity. 1951.

[12] Aimei Zhang and Guimin Chen. A comprehensive el LIptic integral solution to the large deflection problems of thin beams in compliant mechanisms. *Journal of Mechanisms and Robotics*, 5(2), 2013.

[13] Shorya Awtar, Alexander H Slocum, and Edip Sevinçer. Characteristics of beam-based flexure modules. *Journal of Mechanical Design*, 129(6), 2007.

[14] Guimin Chen and Ruiyu Bai. Modeling large spatial deflections of slender bisymmetric beams in compliant mechanisms using chained spatial-beam constraint model. *Journal of Mechanisms and Robotics*, 8(4), 2016.

[15] Hai-Jun Su. A pseudorigid-body 3r model for determining large deflection of cantilever beams subject to tip loads. *Journal of Mechanisms and Robotics*, 1(2), 2009.

[16] Ke Wu and Gang Zheng. A comprehensive static modeling methodology via beam theory for compliant mechanisms. *Mechanism and Machine Theory*, 169:104598, 2022.

[17] Guimin Chen and Fulei Ma. Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model. *Journal of Mechanical Design*, 137(2), 2015.

[18] GR Cowper. The shear coefficient in timoshenko's beam theory. 1966.

[19] Ke Wu and Gang Zheng. Solutions to large beam-deflection problems by taylor series and padé approximant. *Mechanism and Machine Theory*, 2022.

[20] Ke Wu and Gang Zheng. Theoretical analysis on nonlinear buckling, post-buckling of slender beams and bi-stable mechanisms. *Journal of Mechanisms and Robotics*, 14(3), 2022.

[21] Oliver A Bauchau and James I Craig. Euler-bernoulli beam theory. In *Structural analysis*, pages 173–221. Springer, 2009.

[22] AB Comsol. Structural mechanics module user's guide. COMSOL Multiphysics (TM) v, 5, 2012.

[23] Stephen Timoshenko. *Theory of elastic stability*. 2e. Tata McGraw-Hill Education, 1970.