# A Body-frame Beam Constraint Model

Ke Wu, Gang Zheng, Guimin Chen, Shorya Awtar

► To cite this version:

Ke Wu, Gang Zheng, Guimin Chen, Shorya Awtar. A Body-frame Beam Constraint Model. Mechanics and Machine Theory, 2024, 192, pp.105517. 10.1016/j.mechmachtheory.2023.105517. hal-04325275

HAL Id: hal-04325275

https://inria.hal.science/hal-04325275v1

Submitted on 5 Dec 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 CC BY 4.0 - Attribution - International License
# Highlights

## A Body-frame Beam Constraint Model

Ke Wu,Gang Zheng,Guimin Chen,Shorya Awtar

• 1 We propose a body-frame beam constraint model for intermediate-range deflections.

• 2 The proposed model is analytical, closed-form and concise.

• 3 The proposed model considers buckling initiation and intermediate-range buckling.

• 4 The proposed model considers different types of loading besides beam-end loading.

• 5 The proposed model considers continuous deflected beam shapes and the energy stored.
# A Body-frame Beam Constraint Model*

Ke Wu$^a$, Gang Zheng$^{a,*}$, Guimin Chen$^b$ and Shorya Awtar$^c$

$^a$Université de Lille, Inria, CNRS, Centrale Lille, UMR 9189 CRISTAL, Lille, F-59000, France

$^b$Xi'an Jiaotong University, Xi'an, Shaanxi, 710049, China

$^c$University of Michigan, Michigan, United States

## ARTICLE INFO

**Keywords:**

Geometrically Nonlinear Euler Bernoulli beam theory

Beam Constraint Model

Body frame (Frenet frame)

Linearization

## Abstract

Compliant Mechanisms (CMs) have presented their inherently advantageous properties since the nature of CMs is to utilize elastic deformation of the built-in elementary flexible members to transfer motion, force and energy. To predict the performances of CMs and optimize their designs, modeling CMs turns out to be a major concern in this field, particularly, such as high-precision flexure-based CMs. Reported in the literature, the famous Beam Constraint Model (BCM) expressed in Cartesian coordinate system serves as an effective tool to model the deflection of flexures. Inspired by BCM, in this paper, we take one step further to propose a Body-frame Beam Constraint Model (BBCM) to handle buckling and intermediate-range deflections under typical beam-end loading conditions. Numerical and analytical analysis also proves that BBCM can also provide reliable accuracy in complex loading conditions and complex beam geometries within small-to-intermediate-range deflections and may accurately characterize large deflections under some specific loading conditions as well.

# 1. Introduction

## 1.1. Compliant mechanisms

Compliant mechanisms (CMs), as a novel type of mechanical systems, utilize the elastic deformation of the built-in elementary flexible members to transfer motion, force and energy [1][2]. This novel and powerful concept has provided many desired properties for CM-based mechanical systems: compliance due to the flexibility of the used elastic material [1]; monolithic design that results in increased motion precision with no concerns about errors from assembly [3], simplified manufacturing process [1][3][4], reduced weight [1][4], cost reduction [1][3][4] and less maintenance [1][4]. These mentioned desired characteristics have made CMs a gradually mature component in many mechanical applications, such as compliant kinematic joints [3, 5, 6], high-accuracy positioning motion stages [7–9], bi/tri-stable mechanisms [10, 11] and other more complex applications that integrate the former three [12–14].

## 1.2. Commonly used theories and methods for modeling CMs

Modeling CMs plays a very important role in predicting the performances and optimizing the designs of CMs [15]. Here, the key is modeling large deflections of the flexible slender structures built in CMs. In the academic community of beam and rod theories [16], or even solid mechanics [17], modeling large 3D deflections considering all possible strains (like bending, shear and stretch) [16][17], as well as the corresponding methods [18][19][20][21] to solve the models (in the form of differential equations), have all been intensively studied from a theoretical perspective.

However, these theories and the corresponding implemented methods turn out to be too complex, unfriendly and inefficient for users from practical engineering, especially to DIY-design and model a specific CM [22][23], not to mention conducting model-based optimization. Therefore, to meet the needs of engineering, in the current literature, existing contributions in CM community provide many "more practical and engineering" solutions to this beam-deflection problem [23].

The most intuitively straightforward type of methods is called Pseudo-Rigid-Body Models (PRBMs) [1] which introduces the concept of rigid mechanisms and classic robotics, stating that a slender structure is assumed to be composed of a series of rigid links and rigid kinematic joints to behave like a flexible structure. In the literature, some

* This work is partially supported by the National Natural Science Foundation of China (Grant No. 62073081), by the project ROBOCOP [ANR-19-CE19], by the project COSSEROOTS [ANR-20-CE33], France.

* Corresponding author at: Université de Lille, Inria, CNRS, Centrale Lille, UMR 9189 CRISTAL, Lille, F-59000, France.
E-mail address: gang.zheng@inria.fr (G. Zheng).

ORCID(s):
PRBMs [24][25][26][27][28] are particularly designed for modeling flexure-based mechanisms, presenting remarkable contributions in this field. These methods, under the concept of rigid mechanisms, provide handy and simplified strategies to model these flexures via some feasible engineering assumptions. However, the mentioned methods also suffer from their limitations due to the assumptions made. For example, the proposed method in [24][26][27][28] are all based on constant-curvature assumption and logically only valid for the case of pure flexure-end moment and constant cross-sections; the PRBM proposed in [25] is still restricted to linear framework where a modified constant stiffness or compliance matrix is used to describe the relationship between forces and displacements. We have also noticed some PRBMs that are used to characterize large deflections of slender beams as well, such as Chained Pseudo-Rigid-Body Models (CPRBMs) under energy framework [29][30].

Besides, pioneering researchers dig into the essence those PRBMs aim to explain, such as geometrically nonlinear Euler Bernoulli beam theory (in the form of differential equations [15] written in body frame), to provide more theoretical results on beam-deflection problems bypassing the need of making engineering assumptions like PRBMs do [31][32][33]. The authors of [31] proposed an elegant analytical solution to large-deflection problems where several assumptions are made for the sake of the elegance, such as constant cross-sections and constant initial curvature. However, the proposed solution will still has to identify the sign of some functions if multi-solutions possibly exist, which might not be perfectly suitable for engineering. The contributors of [32] proposed a numerical solution in a very compact and clean form, to large-deflection problems of semicircular beams (with constant initial curvature), and based on the obtained results to facilitate the development of PRBMs for modeling planar compliant mechanisms. However, it will not be valid anymore if the sign of current curvature changes (Eq. (40) will not be valid if the sign $\frac{d\theta}{ds}$ changes). In other words, it will not be valid if inflection points exist, which is exactly the same problem reported in Section 4.3 of [22]. One more nice work should be noted that [33] starts from Lagrangian perspective to formulate the potential energy, and use Hamilton's principle to derive its Newtonian formulation for both dynamics and statics. This work is quite theoretical and still sticks to Cartesian coordinate system to describe the geometrical nonlinearity, resulting in the fact that the whole formulation is a bit more complex compared to the classic body-frame beam theory [15].

Essentially speaking, what we actually face in terms of modeling large deflection of slender structures via geometrically nonlinear Euler Bernoulli beam theory is a typical boundary value problem (BVP) of an ordinary differential equation (ODE) [22]. Therefore, our previous work proposes to use typical methods for solving BVPs to solve the mentioned problem. For example, we have used finite difference method [22], shooting method [22], weighted residual methods [22], Taylor series method and Padé approximant [34] to handle the original governing BVP. Besides, a comprehensive study under Newtonian framework is presented in [15] to consider different loading conditions in CMs, such as beam-end loading, distributed loading, distributed pressure, initially varying curvature of beams, varying crosssections of beams and so on. Recently, the authors of [35] extend Timoshenko beam theory for CMs to consider the shear along the beam when it is largely deflected. However, although the aforementioned BVP can be solved via numerical techniques due to its strong nonlinearity [22], there are still some researchers looking for closed-form solutions for more simplicity [7][36].

In Chapter 3 of Awtar's PhD thesis [7], he derived a simple, clear and also analytical formulation by assuming that the curvature of the deflected beam is small enough to transform the body-frame formulation into global Cartesian coordinate system, which is mathematically simpler and user-friendly. This methodology is termed as Beam Bonstraint Model (BCM), serving as a good tool of modeling intermediate-range deflections in high-precision flexure-based motion stages. It has been successfully and widely used in this field since its clear and precise formulation offers a sharp insight into flexure deflection and handy implementation. Guimin and Fulei extended Awtar's work by further taking care of the shear in the crosssections of beams [37]. Later on, some further contributions are made after the success of BCM: a more generalized beam constrained model for initially curved beams and beams with varying crosssections [38], the corresponding energy framework [38] and a chained Beam Constrained Model (CBCM) used for modeling large deflections [39][40].

Different from BCM, instead of transforming the body-frame formulation into global Cartesian coordinate system, in this paper, we propose a new Body-frame Beam Constraint Model (BBCM) where a closed-form formulation via linearization along body frame (Frenet frame) is derived for intermediate-range (post-)buckling and intermediate-range (deflection). It should be noted that many great contributions in mechanics also treat intermediate-range (post-)buckling and intermediate-range (deflection) in Cartesian coordinate system [41][42][43][44][45]. However, our proposed method explores another way to handle the mentioned problems. Under this strategy, first, we bypass the need of dealing with complex coordinate transformation for linearization. Particularly, this coordinate transformation
process will change the interval of the governing ordinary differential equations, which is bit tough to handle. Second, we can directly utilize the existing conclusions and results achieved in the mentioned classic geometrically exact beam and rod theories, such as energy formulation and beam-shape characterization. Finally, this manipulation provides more compact, clear and clean formulation without the need of choosing truncated series for a given accuracy compared to BCM.

## 1.3. Structure of the paper

In Section 2, we recall the classic geometrically nonlinear Euler Bernoulli beam theory as the prerequisites for deriving BBCM, and BCM for comparisons to the proposed BBCM. In Section 3, the detailed derivation and general formulation of BBCM is provided where the scope and differences to BCM are both theoretically discussed. Then, in Section 4, we apply BBCM to compliant mechanisms where intermediate-range buckling and intermediate-range deflection are studied where several typical cases are studied along with numerical results verified by solid-mechanics-based FEM (SD-FEM) and BCM as well. In particular, in the end of this section, the extended formulation of BBCM is also provided for more variety kinds of loading conditions. In the last section, we summarize the work presented in this paper as well as picture the perspective of our future work.

# 2. Recall of theory and prerequisites

## 2.1. Geometrically nonlinear Euler Bernoulli beam theory

The key of modeling compliant mechanisms is to model the built-in flexible structures. Therefore, we would like to recall the famous geometrically nonlinear Euler Bernoulli beam theory first [46][47][48][49][50]. Here, we need to notice the assumptions and conditions in this theory:

* The deflection is only caused by bending so axial stretch and shear are neglected.

* The material is homogeneous.

* The constitutive equation (1) is based on linear elasticity.

Taking a look at the constitutive equation(1), it implies that the curvature change $\frac{d\theta}{ds}(s)$ of the deflected slender beam at some point $s$, or equivalently point $(x(s), y(s))$ is proportional to the external moment $M(s)$ exerted at the same point $s$ or equivalently point $(x(s), y(s))$:

$$EI \frac{d\theta}{ds} = M(s), \quad s \in [0, L] \qquad (1)$$

where $\theta(s)$ refers to the rotation angle of the point $s$ with respect to X axis; $E$ is the Young's modulus of the material; $I$ denotes the second moment of inertia of the cross-section area; $L$ denotes the beam length; $F_x$, $F_y$ and $M_o$ are the constant forces and moment exerted at the beam-end point $s = L$. The loading conditions are shown in Fig. 1a. At $s = 0$ and $s = L$, the slender beam could be subjected to different boundary conditions as graphically demonstrated in Fig. 1b. Finally, we can end up with the following boundary value problem (BVP) in a general form:

$$\text{D.E.} \quad \frac{d^2\theta}{ds^2} = -\frac{F_y}{EI}(\cos\theta(s) + \frac{F_x}{F_y}\sin\theta(s)) \qquad (2)$$

$$\text{B.C.} \quad g(\theta(0), \theta'(0), \theta(L), \theta'(L)) = 0$$

where $g(\theta(0), \theta'(0), \theta(L), \theta'(L)) \in \mathbb{R}^2$ represents different independent boundary conditions defined at $s = 0$ and $s = L$. Note that (2) is the general formulation of a straight slender beam subjected by beam-end loading under different boundary conditions (Fig. 1b). Then, to characterize the deflected beam shapes, we have:

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

In terms of the elastic energy stored, we have:

$$E_p = \int_{0}^{L} \frac{1}{2} EI \left(\frac{d\theta}{ds}\right)^2 ds \qquad (4)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: General loading condition of a slender beam under different beam-end boundary conditions</figcaption></figure>

To simplify the studied BVP and also reveal the mathematical essence of the BVP, the governing equation (2) needs to be nondimensionalized first, as previously done in [7][51][36]:

$$ \text{D.E.} \quad \frac{d^2 \hat{\theta}}{ds^2} = -f_y \cos \hat{\theta}(\hat{s}) - f_x \sin \hat{\theta}(\hat{s}) \quad (5) $$

$$ \text{B.C.} \quad g(\hat{\theta}(0), \hat{\theta}'(0), \hat{\theta}(1), \hat{\theta}'(1)) = 0 $$

where

$$ f_x = \frac{F_x L^2}{EI}; \quad f_y = \frac{F_y L^2}{EI}; \quad m_o = \frac{M_o L}{EI}; \quad \hat{s} \in [0, 1] \quad (6) $$

## 2.2. Beam Constraint Model (BCM)

As stated in [52], BCM works in Cartesian coordinate system, which offers a clear and straightforward formulation to characterize the beam-end displacements. Then, let's recall the commonly-used formulations of BCM for straight beams with rectangular crosssections [7] written in non-dimensional form:

$$ \begin{gathered} \begin{bmatrix} f_y \\ m_o \end{bmatrix} = A \begin{bmatrix} \hat{y}(1) \\ \hat{\theta}(1) \end{bmatrix} - f_x B \begin{bmatrix} \hat{y}(1) \\ \hat{\theta}(1) \end{bmatrix} \\ \hat{x}(1) - 1 = -\frac{h^2}{12} f_x + \begin{bmatrix} \hat{y}(1) \end{bmatrix}^\top C \begin{bmatrix} \hat{y}(1) \\ \hat{\theta}(1) \end{bmatrix} - f_x \begin{bmatrix} \hat{y}(1) \\ \hat{\theta}(1) \end{bmatrix}^\top D \begin{bmatrix} \hat{y}(1) \\ \hat{\theta}(1) \end{bmatrix} \end{gathered} \quad (7) $$

where

$$ A = \begin{bmatrix} 12 & -6 \\ -6 & 4 \end{bmatrix}; \quad B = \begin{bmatrix} 1.2 & -0.1 \\ -0.1 & 2/15 \end{bmatrix}; \quad C = \begin{bmatrix} -0.6 & 1/20 \\ 1/20 & -1/15 \end{bmatrix}; \quad D = \begin{bmatrix} 1/700 & -1/1400 \\ -1/1400 & 11/6300 \end{bmatrix}; \quad (8) $$

Obviously, given any set of $f_x$, $f_y$ and $m_o$, to obtain beam-end coordinates $\hat{x}(1)$, $\hat{y}(1)$ and $\hat{\theta}(1)$, the formulations (7) along with (8) of BCM are closed-form:

$$ \begin{gathered} \begin{bmatrix} \hat{y}(1) \\ \hat{\theta}(1) \end{bmatrix} = (A - f_x B)^{-1} \begin{bmatrix} f_y \\ m_o \end{bmatrix} \\ \hat{x}(1) = -\frac{h^2}{12} f_x + \begin{bmatrix} \hat{y}(1) \end{bmatrix}^\top C \begin{bmatrix} \hat{y}(1) \\ \hat{\theta}(1) \end{bmatrix} - f_x \begin{bmatrix} \hat{y}(1) \end{bmatrix}^\top D \begin{bmatrix} \hat{y}(1) \\ \hat{\theta}(1) \end{bmatrix} + 1 \end{gathered} \quad (9) $$
Similarly, BCM also has a closed-form formulation for the stored energy [38]:

$$ \hat{E}_p = \frac{1}{2} \left( \frac{h^2}{12} \right) f_x^2 + \frac{1}{2} \left[ \hat{y}(1) \right]^T A \left[ \hat{y}(1) \right] + \frac{1}{2} f_x^2 \left[ \hat{y}(1) \right]^T D \left[ \hat{y}(1) \right] \quad (10) $$

Note that (9) only gives the formulation of the beam-end coordinates whereas according to [39], we can use the following to characterize its deflected beam shape:

$$ \begin{gathered} \hat{x}(\hat{s}) = \hat{s} - \frac{h^2 \hat{s}}{12} f_x - \left[ \frac{f_y}{m_o} \right]^T C^* \left[ \frac{f_y}{m_o} \right] \\ \left[ \begin{matrix} \hat{y}(\hat{s}) \\ \hat{\theta}(\hat{s}) \end{matrix} \right] = K^* \left[ \begin{matrix} f_y \\ m_o \end{matrix} \right] \end{gathered} \quad (11) $$

where matrix $C^*$ and matrix $K^*$ can be found in [39]. To dimensionalize the BCM results above, we have

$$ x(s) = \hat{x}(\hat{s})L; \quad y(s) = \hat{y}(\hat{s})L; \quad \theta(s) = \hat{\theta}(\hat{s}); \quad E_p = \frac{EI}{L} \hat{E}_p; \quad s \in [0, L], \quad \hat{s} \in [0, 1] \quad (12) $$

**Remark 1.** Note that the analytical solutions to BCM are in the form of transcendental functions that can be expressed as infinite series. In the final forms (7), (10) and (11), these transcendental functions are expanded and truncated.

### 3. Derivation of Body-frame Beam Constraint Model

As recalled in the last section, BCM has been serving as an effective tool to characterize intermediate-range deflections. BCM is formulated and implemented under Cartesian coordinate system. However, in this paper, we formulate the proposed beam model in body frame, namely **Body-frame Beam Constrain Model (BBCM)**. To start with, we can first approximate the right hand side of the ODE in (5) via the Maclaurin series of $\sin \hat{\theta}$ and $\cos \hat{\theta}$:

$$ \sin \hat{\theta} = \hat{\theta} - \frac{\hat{\theta}^3}{3!} + \frac{\hat{\theta}^5}{5!} - \frac{\hat{\theta}^7}{7!} + \frac{\hat{\theta}^9}{9!} - \dots; \quad \cos \hat{\theta} = 1 - \frac{\hat{\theta}^2}{2!} + \frac{\hat{\theta}^4}{4!} - \frac{\hat{\theta}^6}{6!} + \frac{\hat{\theta}^8}{8!} - \dots \quad (13) $$

If we just use the first terms of these series:

$$ \sin \hat{\theta} \approx \hat{\theta}; \quad \cos \hat{\theta} \approx 1 \quad (14) $$

we will end up with:

$$ \text{D.E.} \quad \frac{d^2 \hat{\theta}}{d \hat{s}^2} + f_x \hat{\theta}(\hat{s}) = -f_y \quad (15) $$

$$ \text{B.C.} \quad g(\hat{\theta}(0), \hat{\theta}'(0), \hat{\theta}(1), \hat{\theta}'(1)) = 0 $$

Obviously, the ODE of (15) is a typical non-homogeneous second-order linear ODE with constant coefficients. Therefore, its analytical solution always exists given the valid boundary conditions.

**Remark 2.** Note (15) is linearized from the original BVP (5) so it is able to locally work in the interval near the point of expansion where the linearized function (the truncated Maclaurin series) can relatively accurately approximate the original function. In other words, the value of $\theta$ should better be within or at least not be too far away from its radius of convergence around $\theta = 0$, which implies that this ODE only governs small-to-intermediate-range deflection of slender beams from the physical perspective.

To characterize the deflected beam shapes, we logically have the following:

$$ \hat{x}(\hat{s}) = \int_{0}^{\hat{s}} \cos(\hat{\theta}(\xi)) - f_{x}\beta d\xi; \quad \hat{y}(\hat{s}) = \int_{0}^{\hat{s}} \sin(\hat{\theta}(\xi))d\xi \quad (16) $$

where $\hat{x}(\hat{s})$ and $\hat{y}(\hat{s})$ are the nondimensionalized coordinates along the beam axis; $f_x\beta$ refers to the axial elastic stretch with $\beta = \frac{h^2}{12L^2}$ for rectangular cross sections and $\beta = \frac{r^2}{4L^2}$ for round cross sections.
**Remark 3.** In (16), we add a correction term (regarding the elastic stretch of the slender beam) along the global X axis by logically assuming that in small-to-intermediate-range deflections, the axial stretch has an insignificant but essential contribution to the deformation in global X axis. Note BCM also incorporates elastic stretching (see Eqs. 3.11 and 3.12 in Chapter 3 of Awtar's PhD thesis [7]).

In the above stated, we term (15) along with (16) as **Body-frame Beam Constraint Model (BBCM)**. In the following, we aim to handle two commonly-encountered types of beam deflections (buckling and general deflection both cause by beam-end loading) in CMs via manipulating the basic formulation of BBCM (15).

**Remark 4.** Here, we would like to clarify the different strategies to derive Beam BCM and BBCM.

* **BCM:** In Awtar's PhD thesis [7], Awtar transforms the formulation of Euler Bernoulli beam theory written in body frame to **Cartesian coordinate system**, and then linearizes the corresponding ordinary differential equation to find the closed-form solution. Therefore, in all previous work regarding BCM, Awtar and other contributors all stay in Cartesian coordinate system to promote this great work [37][38][51][53][54][55].

* **BBCM:** However, different from all previous work regarding BCM, we propose to still stay in **body frame** (Frenet frame), and directly linearize the governing ordinary differential equation of Euler Bernoulli beam theory in body frame to find the closed-form solution. To the best of our knowledge, in the current literature of compliant mechanisms, there is no existing work using our proposed method.

# 4. Application of Body-frame Beam Constraint Model to Compliant Mechanisms

In compliant mechanisms, there are two typical types of elastic deflection of flexible slender structures utilized to transfer motion, force and energy [1][23][37]:

* **(Post-)buckling:** constant-force (or torque) mechanisms [56], bi-stable mechanisms [37] and so on.

* **General deflection:** positioning stages [7][9], compliant grippers [57][58][59], compliant kinematic joints (hinges) [1] and so on.

In this section, we propose to handle the above two physical behaviors of a slender structure via utilizing the basic formulation of BBCM (15). As discussed in Remark 2, BBCM can only works within a small-range deflection due to the fact that (15) is the result of linearization. Therefore, BBCM can be utilized to deal with buckling initiation and intermediate-range post-buckling, and intermediate-range general deflection.

**Remark 5.** From the perspective of solving ordinary differential equations, (post-)buckling refers to multi-solutions to a boundary value problem (BVP) of an ordinary differential equation (ODE) since it may yield different buckling modes. On the contrary, general deflection only has one definite deformation, meaning that there is only one solution to its governing ordinary differential equation. Essentially, what we aim to deal with in the following is one or multi solutions to a BVP of an ODE.

## 4.1. Buckling initiation and intermediate-range post-buckling

The buckling initiation is studied here under different boundary conditions where a slender straight beam is only subjected by an axial load. Essentially, they all share the same governing equation derived from (15) with $f_y = 0$:

$$ \frac{d^2\hat{\theta}}{ds^2} + f_x\hat{\theta}(s) = 0 \qquad (17) $$

which is essentially the corresponding homogeneous ODE of the one in (15). Then, we can directly write down the general solution $\hat{\theta}_h(s)$ to (17):

$$ \hat{\theta}_h(s) = C_1 \cos(\sqrt{f_x} s) + C_2 \sin(\sqrt{f_x} s) \qquad (18) $$

Logically, with boundary conditions available, we will be able to find the exact analytical solution to (17). Mathematically speaking, the essence of buckling is all about bifurcation of differential equations no matter derived from
<figure><img src="image_3.png" /><figcaption>(a) Schematic diagram of a slender straight beam fixed at one end and exerted by an axial load</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Schematic diagram of a slender straight beam constrained by sliding passive kinematic joints at both ends</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Schematic diagram of a slender straight beam subjected by sliding clamped constraints at the both ends</figcaption></figure>

<caption>**Figure 2:** Buckling initiation and post-buckling modes of Case 1</caption>

<table><caption>**Table 1**
Buckling limits of Case 1</caption><thead><tr><th></th><th></th><th>1st</th><th>2nd</th><th>3rd</th><th>4th</th><th>5th</th><th>6th</th></tr></thead><tbody><tr><td rowspan="3">f<sub>x</sub><sup>*</sup></td><td>BBCM</td><td>(π/2)<sup>2</sup></td><td>(3π/2)<sup>2</sup></td><td>(5π/2)<sup>2</sup></td><td>(7π/2)<sup>2</sup></td><td>(9π/2)<sup>2</sup></td><td>(11π/2)<sup>2</sup></td></tr><tr><td>SD-FEM</td><td>2.4816</td><td>22.3370</td><td>62.0620</td><td>121.68</td><td>201.19</td><td>300.59</td></tr><tr><td>ER</td><td>0.57%</td><td>0.58%</td><td>0.61%</td><td>0.63%</td><td>0.66%</td><td>0.68%</td></tr></tbody></table>

solid mechanics or beam theories. Here, three types of buckling subjected to different boundary conditions are studied (see Fig. 2). Note the linearized ODE (17) is used here to characterize the buckling initiation (bifurcation of ODE (17)) since it is straightforward to analyze the buckling limits of different modes using (16) and (18). As stated in Remark 2, BBCM can only deal with intermediate-range deflections due to the linearization manipulation used so BBCM is only valid around the initiation areas. For exact buckling behaviors of slender beams, [60] is recommended for interested readers.

### 4.1.1. Case 1: a slender straight beam fixed at one end and exerted by an axial load at another end

As shown in Fig. 2a, the boundary conditions are

$$ \text{B.C. } \hat{\theta}_h(0) = 0; \quad \frac{d\hat{\theta}_h(1)}{d\hat{s}} = 0 \qquad (19) $$

which yields the analytical solution:

$$ \hat{\theta}_h(s) = C_2 \sin(\sqrt{f_x \hat{s}}); \quad f_x^* = \left[ \frac{\pi}{2} + (n-1)\pi \right]^2; \quad n = 1, 2, 3 \dots \qquad (20) $$

The detailed derivation and manipulation are demonstrated in Appendix.
<figure><img src="image_6.png" /><figcaption>**Figure 3:** Diagram of different buckling modes for Case 1</figcaption></figure>

<table><caption>**Table 2**
Buckling limits of Case 2</caption><thead><tr><th></th><th></th><th>1st</th><th>2nd</th><th>3rd</th><th>4th</th><th>5th</th><th>6th</th></tr></thead><tbody><tr><td rowspan="3">f<sub>x</sub><sup>*</sup></td><td>BBCM</td><td>&pi;<sup>2</sup></td><td>4&pi;<sup>2</sup></td><td>9&pi;<sup>2</sup></td><td>16&pi;<sup>2</sup></td><td>25&pi;<sup>2</sup></td><td>36&pi;<sup>2</sup></td></tr><tr><td>SD-FEM</td><td>9.8205</td><td>39.2950</td><td>88.4610</td><td>157.36</td><td>246.02</td><td>354.43</td></tr><tr><td>ER</td><td>0.50%</td><td>0.47%</td><td>0.41%</td><td>0.35%</td><td>0.29%</td><td>0.24%</td></tr></tbody></table>

**4.1.2. Case 2: a slender straight beam constrained by moving passive kinematic joints at both ends**
The boundary conditions demonstrated in Fig. 2b are

$$ \text{B.C.} \quad \frac{d\hat{\theta}_h(0)}{d\hat{s}} = 0; \quad \frac{d\hat{\theta}_h(1)}{d\hat{s}} = 0 \qquad (21) $$

which yields the analytical solution:

$$ \hat{\theta}_h(s) = C_1 \cos(\sqrt{f_x \hat{s}}); \quad f_x^* = (n\pi)^2; \quad n = 1, 2, 3... \qquad (22) $$

The detailed derivation and manipulation are demonstrated in Appendix.

**4.1.3. Case 3: a slender straight beam subjected by moving clamped constraints at the both ends**
In this Case, the boundary conditions as shown in Fig. 2c are

$$ \text{B.C.} \quad \hat{\theta}_h(0) = 0; \quad \hat{\theta}_h(\frac{1}{2}) = 0; \quad \hat{\theta}_h(1) = 0 \qquad (23) $$

which yields the analytical solution:

$$ \hat{\theta}_h(s) = C_1 \sin(\sqrt{f_x \hat{s}}); \quad f_x^* = (2n\pi)^2; \quad n = 1, 2, 3... \qquad (24) $$

Here, since we can not get enough information from the boundary conditions, it could be complex if still being stuck using the analytical formulation (24). We recommend directly using the original nonlinear strategy stated in [60] to characterize the exact buckling behaviors.
<figure><img src="image_7.png" /><figcaption>**Figure 4:** Bifurcation diagrams of Case 1</figcaption></figure>

### 4.1.4. Results analysis and discussions

In (20), (22) and (24), $f_x^*$ refers to the buckling limit and $n$ stands for the $n$th-order buckling mode. In (20), (22) and (24), the amplitudes $C_1$ and $C_2$ of the trigonometric functions are not defined so any value would satisfy the conditions, and they could be positive and negative, yielding Scenario 2 and Scenario 3 in Fig. 2a, Fig. 2b and Fig. 2c. Here, $C_1$ and $C_2$ are chosen as 1 for simplicity. However, we can still use (20), (22) and (24) along with (16) to characterize the buckling mode shapes of different orders, as shown in Fig. 3, Fig. 5 and Fig. 7. Here, BBCM and SD-FEM denote the results of the proposed Body-frame Beam Constraint Model and solid-mechanics-based finite element method (SD-FEM). SD-FEM is a very common module in most commercial FEM-based software frameworks, such as ANSYS, COMSOL and Abaqus. It means that to model deformable bodies, the software solves the corresponding governing differential equations arising from solid mechanics via finite element method. In the above cases, 20 hexahedral mesh elements are used for the linear buckling analysis (eigenvalue analysis), and logically planar bending, shear and stretch are all considered [17]. Besides, geometric nonlinearity is activated since intermediate-range deflections of beams are modeled, and the assumption of linear elasticity is made as well. Note that geometrical nonlinearity of the proposed BBCM comes from the consideration of bending strain happening in every single node defined along its beam axis (written in body frame), which is a typical characteristics of geometrically nonlinear beam and rod theories [61].

Similarly, the corresponding buckling limits of different orders (modes) are also available in (20), (22) and (24), which are verified by SD-FEM (See Tables. 1, 2 and 3). Note that BBCM, SD-FEM and ER denote the results of the proposed body-frame Beam Constraint Model, SD-FEM, and the errors with respect to SD-FEM.
<figure><img src="image_8.png" /><figcaption>**Figure 5:** Diagram of different buckling modes for Case 2</figcaption></figure>

<table><caption>**Table 3**
Buckling limits Case 3</caption><thead><tr><th></th><th></th><th>1st</th><th>2nd</th><th>3rd</th><th>4th</th><th>5th</th><th>6th</th></tr></thead><tbody><tr><td rowspan="3">f<sub>x</sub><sup>*</sup></td><td>BBCM</td><td>4π<sup>2</sup></td><td>16π<sup>2</sup></td><td>36π<sup>2</sup></td><td>64π<sup>2</sup></td><td>100π<sup>2</sup></td><td>144π<sup>2</sup></td></tr><tr><td>SD-FEM</td><td>39.94</td><td>159.82</td><td>359.61</td><td>638.74</td><td>995.7</td><td>1428</td></tr><tr><td>ER</td><td>1.15%</td><td>1.19%</td><td>1.20%</td><td>1.11%</td><td>0.88%</td><td>0.47%</td></tr></tbody></table>

Note that (69) and (77) should be regarded as semi-closed-form formulations since they also need a numerical process to get the value of $\alpha$ beforehand. However, the analytic formulations of (69) and (77) offer an intuitive understanding of deflected beam shape at least at small-to-intermediate-range deflections (see Figs. 4 and 6 where the numerical results are provided and compared to that of SD-FEM and Euler Bernoulli beam theory).

Taking a look at the two dash-line boxes in Fig. 4a and Fig. 4b, it can be concluded that from small to intermediate deflection range ($\hat{y}(1)$ is up to 0.4), BBCM can still capture relatively accurate results compared to that of solid mechanics. This is logical since BBCM is the linearized outcome of (5) with limited accuracy outside the radius of convergence. We can notice the same story in Case 2, specifically in the dash-line boxes in Fig. 6a and Fig. 6b. However, it is still promising since the closed-form formulations have provided a reliable prediction within the intermediate-range of deflection, and they can serve as the initial guesses for the original nonlinear governing equations stated in [60].

Overall speaking, BBCM can be used together with the numerical methodology proposed in [60]: first, we can analytically obtain the buckling limits and mode shapes as well as exact characterization of intermediate-range buckling via BBCM. Then, we can start from the pre-results of BBCM, such as buckling limits and mode shapes, to numerically handle the nonlinear buckling and post-buckling behaviors, which significantly speeds up the modeling process without both numerically and empirically guessing the buckling limits and mode shapes [60].

## 4.2. Intermediate-range deflection

This is the most common case in CMs, especially in compliant positioning mechanisms that require high-accuracy motions [1]. In the previous work, the BCM has been proved remarkably effective in terms of modeling up to
<figure><img src="image_9.png" /><figcaption>Figure 6: Bifurcation diagrams of Case 2</figcaption></figure>

intermediate-range deflections of slender straight beams (see Fig. 8). Here, in this section, we deal with the intermediate-range deflection via BBCM. The derivation, formulation comparison with BCM, feasibility analysis of BBCM, and mechanism synthesis are all provided accordingly.

4.2.1. Derivation

Starting from (15), the boundary conditions in this case yields:

$$
\text{B.C. } \hat{\theta}(0) = 0; \hat{\theta}'(1) = m_o \tag{25}
$$

Taking a look at (15), a typical non-homogeneous second-order linear ODE with constant coefficients, we normally
seek the analytical solution through the following procedure. We assume the solution $\hat{\theta}(\hat{s}) = \hat{\theta}_h(\hat{s}) + \hat{\theta}_p(\hat{s})$ where $\hat{\theta}_h(\hat{s})$
is the analytical solution to the corresponding homogeneous second-order linear ODE and $\hat{\theta}_p(\hat{s})$ is the extra linear
independent solution to the non-homogeneous part $-f_y$.

First, $\theta_p(s)$ can be directly solved:

$$
\hat{\theta}_p(\hat{s}) = -\frac{f_y}{f_x} \tag{26}
$$
<figure><img src="image_10.png" /><figcaption>Figure 7: Diagram of different buckling modes for Case 3</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: A slender straight beam subjected to general beam-end loading conditions</figcaption></figure>

Then, in terms of the corresponding homogeneous second-order linear ODE, we have

$$ \frac{d^2\hat{\theta}}{d\hat{s}^2} + f_x \hat{\theta}(\hat{s}) = 0 \qquad (27) $$

For this homogeneous second-order linear ODE (27), there are three scenarios depending on the sign of $f_x$. If $f_x > 0$, obviously, we can directly write down the general solution $\hat{\theta}_h(\hat{s})$ to (27):

$$ \hat{\theta}_h(\hat{s}) = C_1 \cos(\sqrt{f_x}\hat{s}) + C_2 \sin(\sqrt{f_x}\hat{s}) \qquad (28) $$

Therefore, the general solution can be then formulated as:

$$ \hat{\theta}(\hat{s}) = C_1 \cos(\sqrt{f_x}\hat{s}) + C_2 \sin(\sqrt{f_x}\hat{s}) - \frac{f_y}{f_x} \qquad (29) $$

Logically, we can determine the two constant coefficients and find the final exact solution (30) by satisfying (25).

$$ \hat{\theta}(\hat{s}) = \frac{\sin(\sqrt{f_x}\hat{s})}{\cos(\sqrt{f_x})} \left( \frac{\sin(\sqrt{f_x})f_y}{f_x} + \frac{m_o}{\sqrt{f_x}} \right) + \frac{\cos(\sqrt{f_x}\hat{s})f_y}{f_x} - \frac{f_y}{f_x}, \quad (f_x > 0) \qquad (30) $$
Then, if $f_x = 0$, rearranging (15), we will end up with the following:

$$ \frac{d^2\hat{\theta}}{d\hat{s}^2} = -f_y \qquad (31) $$

Satisfying (25), we will be able to find the analytical solution:

$$ \hat{\theta}(\hat{s}) = -\frac{1}{2} f_y \hat{s}^2 + (m_o + f_y) \hat{s}, \quad (f_x = 0) \qquad (32) $$

Finally, if $f_x < 0$, the general solution $\hat{\theta}_h(s)$ to (27) is:

$$ \hat{\theta}_h(\hat{s}) = C_1 e^{\sqrt{-f_x} s} + C_2 e^{-\sqrt{-f_x} s} \qquad (33) $$

Logically, the general solution can be then formulated as:

$$ \hat{\theta}(\hat{s}) = C_1 e^{\sqrt{-f_x} s} + C_2 e^{-\sqrt{-f_x} s} - \frac{f_y}{f_x} \qquad (34) $$

Then, we can determine the two constant coefficients by satisfying (25), and find the final exact solution (35). Now,

$$ \hat{\theta}(\hat{s}) = \frac{(f_y \sqrt{-f_x} e^{-\sqrt{-f_x}} + f_x m_o) e^{\sqrt{-f_x} \hat{s}} + (f_y \sqrt{-f_x} e^{-\sqrt{-f_x}} - f_x m_o) e^{-\sqrt{-f_x} \hat{s}}}{\sqrt{-f_x} (e^{\sqrt{-f_x}} + e^{-\sqrt{-f_x}}) f_x} - \frac{f_y}{f_x}, \quad (f_x < 0) \qquad (35) $$

we have the final analytical solution to (15), as well as the solution to BBCM for **intermediate-range deflections** of
slender straight beams subjected to beam-end loading conditions:

$$ \hat{\Gamma}(\hat{s}, f_x, f_y, m_o) = \begin{cases} \hat{\theta}(\hat{s}) = \frac{\sin(\sqrt{f_x}\hat{s})}{\cos(\sqrt{f_x})} \left( \frac{\sin(\sqrt{f_x})f_y}{f_x} + \frac{m_o}{\sqrt{f_x}} \right) + \frac{\cos(\sqrt{f_x}\hat{s})f_y}{f_x} - \frac{f_y}{f_x}, & (f_x > 0) \quad (30) \\ \hat{\theta}(\hat{s}) = -\frac{1}{2} f_y \hat{s}^2 + (m_o + f_y)\hat{s}, & (f_x = 0) \quad (32) \\ \hat{\theta}(\hat{s}) = \frac{(f_y \sqrt{-f_x} e^{-\sqrt{-f_x}} + f_x m_o) e^{\sqrt{-f_x} \hat{s}} + (f_y \sqrt{-f_x} e^{-\sqrt{-f_x}} - f_x m_o) e^{-\sqrt{-f_x} \hat{s}}}{\sqrt{-f_x} (e^{\sqrt{-f_x}} + e^{-\sqrt{-f_x}}) f_x} - \frac{f_y}{f_x}, & (f_x < 0) \quad (35) \end{cases} $$

with

$$ \hat{x}(\hat{s}) = \int_{0}^{\hat{s}} \cos(\hat{\Gamma}(\xi)) - f_x \beta d\xi; \quad \hat{y}(\hat{s}) = \int_{0}^{\hat{s}} \sin(\hat{\Gamma}(\xi)) d\xi \qquad (37) $$

where (36) is formulated in a piecewise form depending on the sign of $f_x$. In this BBCM formulation, (36) serves as the solution to the constitutive law whereas (37) characterizes the deformed beam shapes under beam-end loading conditions. Clearly, BBCM offers a purely analytic and closed-form formulation (Eq. (36) and Eq. (37)) so its modeling efficiency has been improved while its computational complexity has been decreased compared to numerical methods [15][22][60].

**Remark 6.** However, we are still interested in the continuity of the function $\hat{\Gamma}(\hat{s}, f_x, f_y, m_o)$ at $f_x = 0$. Therefore, we need to find out the limit of (36) while $f_x \to 0$ via L'Hôpital's rule. For $f_x \to 0$ and according to (30)(32)(35), we have

$$ \lim_{f_x \to 0^+} \hat{\Gamma}(\hat{s}, f_x, f_y, m_o) = \lim_{f_x \to 0^-} \hat{\Gamma}(\hat{s}, f_x, f_y, m_o) = -\frac{1}{2} f_y \hat{s}^2 + (m_o + f_y) \hat{s} \qquad (38) $$

Clearly, (38) proves the continuity of $\hat{\Gamma}(\hat{s}, f_x, f_y, m_o)$ at $f_x = 0$.
For nondimensional energy calculation, we have:

$$ \hat{E}_p = \int_0^1 \frac{1}{2} \left( \frac{\partial \hat{\Gamma}}{\partial \hat{s}} \right)^2 d\hat{s} + \frac{1}{2} \left( \frac{h^2}{12} \right) f_x^2 \quad (39) $$

Note that BCM also has a similar formulation for calculating the energy [38], and we have demonstrated the formulation of BCM in the next section.

### 4.2.2. Feasibility analysis of BBCM

In this section, we aim to analyze the scope of feasibility of the proposed BBCM via a straightforward enumeration strategy. Under the same boundary conditions (25), we use the original governing equation of geometrically nonlinear Euler Bernoulli beam theory (5) as the standard to compare the results from BBCM, equivalently (15) and also (36). In flexure-based CMs, the beam-end motions of the built-in flexures are the main contributors to the motions of these mechanisms. Therefore, we here use the following formulation to form the expression of the evaluated errors:

$$ e_r = \sqrt{e_{r_x}^2 + e_{r_y}^2 + e_{r_\theta}^2} \quad (40) $$

with

$$ e_{r_x}^2 = \left( \int_0^1 \cos(\hat{\Gamma}(\hat{s})) d\hat{s} - \int_0^1 \cos(\hat{\theta}(\hat{s})) d\hat{s} \right)^2; \quad e_{r_y}^2 = \left( \int_0^1 \sin(\hat{\Gamma}(\hat{s})) d\hat{s} - \int_0^1 \sin(\hat{\theta}(\hat{s})) d\hat{s} \right)^2; \quad e_{r_\theta}^2 = (\hat{\Gamma}(1) - \hat{\theta}(1))^2 \quad (41) $$

where $\hat{\Gamma}(\hat{s})$ and $\hat{\theta}(\hat{s})$ denote the solutions to (15) and (5) subjected to the boundary conditions (25).

Here, we consider $f_x$, $f_y$ and $m_o$ as the inputs of the model ((15) and (5)) where $f_x \in U_{f_x} = [\underline{f_x}, \overline{f_x}]$, $f_y \in U_{f_y} = [\underline{f_y}, \overline{f_y}]$ and $m_o \in U_{m_o} = [\underline{m_o}, \overline{m_o}]$, and $e_r$ as the corresponding output. We can then formulate the whole admissible input space: $U = U_{f_x} \times U_{f_y} \times U_{m_o}$ with $U_{f_x}, U_{f_y}, U_{m_o} \in \mathbb{R}$. Then, we sample $U_{f_x}, U_{f_y}$ and $U_{m_o}$ via a set of $N$ points, noted as $U_{f_x}^s, U_{f_y}^s$ and $U_{m_o}^s$ with fixed sampling lengths: $\frac{\overline{f_x}-\underline{f_x}}{N-1}, \frac{\overline{f_y}-\underline{f_y}}{N-1}$ and $\frac{\overline{m_o}-\underline{m_o}}{N-1}$, i.e., $U_{f_x} \approx U_{f_x}^s, U_{f_y} \approx U_{f_y}^s$ and $U_{m_o} \approx U_{m_o}^s$:

$$ U_{f_x}^s = \{ \underline{f_x}, \overline{f_x} + \frac{\overline{f_x} - \underline{f_x}}{N-1}, \dots, \underline{f_x} - \frac{\overline{f_x} - \underline{f_x}}{N-1}, \overline{f_x} \}; \quad U_{f_y}^s = \{ \underline{f_y}, \overline{f_y} + \frac{\overline{f_y} - \underline{f_y}}{N-1}, \dots, \underline{f_y} - \frac{\overline{f_y} - \underline{f_y}}{N-1}, \overline{f_y} \}; \quad (42) $$

$$ U_{f_x}^s = \{ \underline{m_o}, \overline{m_o} + \frac{\overline{m_o} - \underline{m_o}}{N-1}, \dots, \underline{m_o} - \frac{\overline{m_o} - \underline{m_o}}{N-1}, \overline{m_o} \} $$

and the input space can therefore be approximated by

$$ U \approx U^s = U_{f_x}^s \times U_{f_y}^s \times U_{m_o}^s \quad (43) $$

From an engineering point of view, the following values are assigned to the corresponding parameters in (42) to evaluate the error (40) distribution through the sampling strategy (or the enumeration framework):

$$ \underline{f_x} = -1; \quad \overline{f_x} = 1; \quad \underline{f_y} = -1; \quad \overline{f_y} = 1; \quad \underline{m_o} = -1; \quad \overline{m_o} = 1; \quad N = 0.025 \quad (44) $$

Then, we can obtain the following numerical results below in Table 4 where the results from BCM are compared. Note that the results from BCM follow the same evaluation framework as stated in (40) to (44). For clarification, let's take the first row of BBCM section in Table 4 for example: if $\hat{y}(1) < 0$, the errors of 95.66% of BBCM results are less than $5.5 \times 10^{-3}$ and the maximum error is $1.53 \times 10^{-2}$.

Obviously, in Table 4, under the same evaluation policy, BBCM maintains the same accuracy level as BCM in term of deflections within small-to-intermediate range.
<table><caption>**Table 4**
Accuracy evaluation results by BBCM and BCM</caption><thead><tr><th>Model</th><th colspan="3">Evaluation results</th></tr></thead><tbody><tr><td rowspan="3">BBCM</td><td>ŷ(1) &lt; 0.1;</td><td>95.66% e<sub>r</sub> &lt; 5.5 × 10<sup>-3</sup>;</td><td>e<sub>rmax</sub> = 1.53 × 10<sup>-2</sup></td></tr><tr><td>ŷ(1) &lt; 0.2;</td><td>95.48% e<sub>r</sub> &lt; 1.55 × 10<sup>-2</sup>;</td><td>e<sub>rmax</sub> = 4.20 × 10<sup>-2</sup></td></tr><tr><td>ŷ(1) &lt; 0.3;</td><td>95.33% e<sub>r</sub> &lt; 3.0 × 10<sup>-2</sup>;</td><td>e<sub>rmax</sub> = 8.47 × 10<sup>-2</sup></td></tr><tr><td rowspan="3">BCM</td><td>ŷ(1) &lt; 0.1;</td><td>96.05% e<sub>r</sub> &lt; 5.5 × 10<sup>-3</sup>;</td><td>e<sub>rmax</sub> = 1.49 × 10<sup>-2</sup></td></tr><tr><td>ŷ(1) &lt; 0.2;</td><td>95.80% e<sub>r</sub> &lt; 1.55 × 10<sup>-2</sup>;</td><td>e<sub>rmax</sub> = 4.10 × 10<sup>-2</sup></td></tr><tr><td>ŷ(1) &lt; 0.3;</td><td>94.66% e<sub>r</sub> &lt; 3.0 × 10<sup>-2</sup>;</td><td>e<sub>rmax</sub> = 8.52 × 10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 9:** Schematic diagram of a straight-beam-based compliant parallelogram</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10:** Numerical results of the studied compliant parallelogram by BBCM</figcaption></figure>

### 4.2.3. Mechanism analysis

In this section, we aim to validate the proposed BBCM to solve the constitutive equations of modeling CMs where the results of BCM are provided for comparisons. A typical CM, a **straight-beam-based compliant parallelogram** is studied for the demonstration of applying BBCM.

As shown in Fig. 9, the mechanism is composed of a rigid motion stage and two compliant beams (Beam #1 and Beam #2). K is the reference point. Other geometric information and loading conditions are also provided in Fig. 9. To statically model this mechanism, we need to follow the standard modeling procedure stated in solid mechanics [62].

#### 1. Constitutive equations
<figure><img src="image_14.png" /><figcaption>Figure 11: Numerical results of the studied compliant parallelogram by BCM</figcaption></figure>

For Beam #1, in coordinate system X₁-O₁-Y₁, we have

$$x_1(L) = L\hat{x}_1(1) = L \int_0^1 \cos(\hat{\Gamma}(\hat{s}, p_A, f_A, m_A)) - p_A\beta d\hat{s}; \quad y_1(L) = L\hat{y}_1(1) = L \int_0^1 \sin(\hat{\Gamma}(\hat{s}, p_A, f_A, m_A))d\hat{s} \quad (45)$$

with

$$p_A = \frac{P_A L^2}{EI}; \quad f_A = \frac{F_A L^2}{EI}; \quad m_A = \frac{M_A L}{EI} \quad (46)$$

For Beam #2, in coordinate system X₂-O₂-Y₂, we have

$$x_2(L) = L\hat{x}_2(1) = L \int_0^1 \cos(\hat{\Gamma}(\hat{s}, p_B, f_B, m_B)) - p_B\beta d\hat{s}; \quad y_2(L) = L\hat{y}_2(1) = L \int_0^1 \sin(\hat{\Gamma}(\hat{s}, p_B, f_B, m_B))d\hat{s} \quad (47)$$

with

$$p_B = \frac{P_B L^2}{EI}; \quad f_B = \frac{F_B L^2}{EI}; \quad m_B = \frac{M_B L}{EI} \quad (48)$$

## 2. Force equilibrium equations

$$F_K - F_A - F_B = 0; \quad P_K - P_A - P_B = 0 \\ M_K - M_A - M_B + F_K T/2 \sin \Phi - P_K T/2 \cos \Phi - F_B T \sin(\Phi) + P_B T \cos(\Phi) = 0 \quad (49)$$

## 3. Geometric compatibility equations

In the global coordinate system X-O-Y, we have

$$\begin{gathered} \begin{bmatrix} x_A \\ y_A \end{bmatrix} = \begin{bmatrix} L \\ T/2 \end{bmatrix}; \quad \begin{bmatrix} x_{A^*} \\ y_{A^*} \end{bmatrix} = \begin{bmatrix} x_1(L) \\ y_1(L) \end{bmatrix} + \begin{bmatrix} 0 \\ T/2 \end{bmatrix}; \quad \begin{bmatrix} x_B \\ y_B \end{bmatrix} = \begin{bmatrix} L \\ -T/2 \end{bmatrix}; \quad \begin{bmatrix} x_{B^*} \\ y_{B^*} \end{bmatrix} = \begin{bmatrix} x_2(L) \\ y_2(L) \end{bmatrix} - \begin{bmatrix} 0 \\ T/2 \end{bmatrix}; \\ AB = [(x_A - x_B)^2 + (y_A - y_B)^2]^{0.5} = T; \quad A^*B^* = [(x_{A^*} - x_{B^*})^2 + (y_{A^*} - y_{B^*})^2]^{0.5} = T; \\ |x_{A^*} - x_{B^*}| = T \sin \Phi; \quad |y_{A^*} - y_{B^*}| = T \cos \Phi; \quad \hat{\Gamma}(1, p_A, f_A, m_A) = \hat{\Gamma}(1, p_B, f_B, m_B) = \Phi; \\ \begin{bmatrix} x_K \\ y_K \end{bmatrix} = \frac{1}{2} \begin{bmatrix} x_A \\ y_A \end{bmatrix} + \frac{1}{2} \begin{bmatrix} x_B \\ y_B \end{bmatrix}; \quad \begin{bmatrix} x_{K^*} \\ y_{K^*} \end{bmatrix} = \frac{1}{2} \begin{bmatrix} x_{A^*} \\ y_{A^*} \end{bmatrix} + \frac{1}{2} \begin{bmatrix} x_{B^*} \\ y_{B^*} \end{bmatrix}; \quad \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix} = \begin{bmatrix} x_{K^*} \\ y_{K^*} \end{bmatrix} - \begin{bmatrix} x_K \\ y_K \end{bmatrix} \end{gathered} \quad (50)$$

where $[x_A \ y_A]^T$, $[x_{A^*} \ y_{A^*}]^T$, $[x_B \ y_B]^T$, $[x_{B^*} \ y_{B^*}]^T$, $[x_K \ y_K]^T$ and $[x_{K^*} \ y_{K^*}]^T$ are the coordinates of A, A*, B, B*, K and K*. Then, we can numerically solve Eq. (45) to (50) using Newton-Raphson method. Besides, $[\Delta x \ \Delta y \ \Phi]^T$ refers to the translational and rotational displacements of K, the reference point of the modeled mechanism.
<figure><img src="image_15.png" /><figcaption>Figure 12: The maximum-deflected configuration of the studied CM obtained by BBCM and SD-FEM</figcaption></figure>

We are interested in the relationships between applied loads $[F_K \ P_K \ M_K]^T$ and the corresponding displacements $[\Delta x \ \Delta y \ \Phi]^T$ both at the reference point K. The material properties and geometry information are provided in the following:

$$E = 200 \times 10^9 \text{ Pa}; w = 0.015 \text{ m}; h = 0.0045 \text{ m}; L = 0.25 \text{ m}; T = 0.2 \text{ m}; A = hw; I = \frac{wh^3}{12}; \beta = \frac{h^2}{12L^2}$$

The applied loads are set up as:

$$\begin{bmatrix} F_K \\ P_K \\ M_K \end{bmatrix} = iU = i \begin{bmatrix} 80 \text{ N} \\ 80 \text{ N} \\ 4 \text{ N.m} \end{bmatrix} \quad (i = 1, 2, 3, 4...10)$$

The numerical results are provided in Fig. 10 and Fig. 11. In both modeling process, the next-step initial guess is set up as the last-step answer for faster convergence [17]. It can be easily noticed that results from BBCM achieve the same accuracy level as BCM does where the maximum transverse deflection of the built-in flexible beams goes up to 10% of the beam length L.

Then, we can characterize the deflected beam shapes via (46) and (48) (see the maximum deflected configuration of the studied CM under 10U in Fig. 12), which serves as a handy tool to analyze the potential mechanical interference during the motion of the CM.

Besides, we can also calculate the stored elastic energy in the studied CM via BBCM:

$$E_p = \frac{EI}{L}(\hat{E}_{p1} + \hat{E}_{p2}) = \frac{EI}{L} \left[ \int_0^1 \frac{1}{2} \left( \frac{\partial \hat{\Gamma}(\hat{s}, p_A, f_A, m_A)}{\partial \hat{s}} \right)^2 d\hat{s} + \frac{1}{2} \left( \frac{h^2}{12} \right) p_A^2 + \int_0^1 \frac{1}{2} \left( \frac{\partial \hat{\Gamma}(\hat{s}, p_B, f_B, m_B)}{\partial \hat{s}} \right)^2 d\hat{s} + \frac{1}{2} \left( \frac{h^2}{12} \right) p_B^2 \right] \quad (51)$$

The results via BBCM and BCM are plotted in Fig. 13 where the errors of both are all below 4% along the motion of the studied CM.

### 4.2.4. General Body-frame Beam Constraint Model

As we have stated earlier, one of the key differences is the different coordinate systems used by BBCM and BCM. Being expressed in the body frame makes BBCM consistent with the fruitful contributions in rod [61] and beam [15][60] theories. This proves why BBCM offers standard formulations (or solutions) for characterization of
<figure><img src="image_16.png" /><figcaption>Figure 13: The stored elastic energy stored in the studied CM</figcaption></figure>

the deflected beam shape, energy calculation and buckling problems. Besides, we can easily explore more based on the previous contributions in rod and beam theories. For example, we can consider more different loading scenarios on the studied beam (see details in [15]). As shown in Fig. 14a, besides the beam-end loading, to take distributed moment $M^*$, distributed loading $q_y$ and initially varying curvature of the beam into account, based on BVP (2), we can reformulate the problem as a new BVP (52). Following the similar procedure such as (5) and (6) to nondimensionalize and linearize BVP (52), we have the final governing BVP (53) along (54), and the corresponding nondimensional graph is shown in Fig. 14b.

$$ \text{D.E.} \quad EI \frac{d^2\theta}{ds^2} = -(F_y \cos\theta(s) + F_x \sin\theta(s)) + \frac{dM^*}{ds} + q_y \cos\theta(s)(L-s) + EI \frac{d\kappa(s)}{ds} \quad (52) $$

$$ \text{B.C.} \quad \theta(0) = 0; \quad \theta'(L) = \frac{M_o + M^*(L)}{EI} + \kappa(L) $$

$$ \text{D.E.} \quad \frac{d^2\hat{\theta}}{ds^2} + f_x \hat{\theta} = -f_y + \frac{dm^*(\hat{s})}{d\hat{s}} + \hat{q}_y(\hat{s})(1-\hat{s}) + \frac{d\hat{\kappa}(\hat{s})}{d\hat{s}} \quad (53) $$

$$ \text{B.C.} \quad \hat{\theta}(0) = 0; \quad \hat{\theta}'(1) = m_o^* = m_o + m^*(1) + \hat{\kappa}(1) $$

$$ f_y = \frac{F_y L^2}{EI}; \quad f_x = \frac{F_x L^2}{EI}; \quad m_o(\hat{s}) = \frac{M_o L}{EI}; \quad m^*(\hat{s}) = \frac{M^*(s)L}{EI}; \quad \hat{q}_y(\hat{s}) = \frac{q_y(s)L^3}{EI}; \quad \hat{\kappa}(\hat{s}) = \kappa(s)L; \quad s \in [0, L]; \quad \hat{s} \in [0, 1] \quad (54) $$

It is noticed that (53) is a typical non-homogeneous second-order ODE, and $m^*(\hat{s})$, $\hat{q}_y(\hat{s})$ and $\hat{\kappa}(\hat{s})$ are all functions of $s$ in a general sense. Therefore, we assume the mentioned three variables are all polynomials:

$$ m^*(\hat{s}) = a_0 + a_1 \hat{s} + a_2 \hat{s}^2 + a_3 \hat{s}^3; \quad \hat{q}_y(\hat{s}) = b_0 + b_1 \hat{s}; \quad \hat{\kappa}(\hat{s}) = c_0 + c_1 \hat{s} + c_2 \hat{s}^2 + c_3 \hat{s}^3 \quad (55) $$

which suggests that $m^*(\hat{s})$, $\hat{q}_y(\hat{s})$ and $\hat{\kappa}(\hat{s})$ vary along the beam axis. Note that in theory, $q_y(s)$ denotes the general transverse distributed loading along the beam whereas in practice, it normally refers the gravity where $\hat{q}_y(\hat{s})$ is constant.

Similarly, the solution to (53) can be formulated as $\hat{\theta}(\hat{s}) = \hat{\theta}_h(\hat{s}) + \hat{\theta}_p(\hat{s})$ as well. Logically, we can derive the following:

$$ \hat{\theta}_p(\hat{s}) = A + Bs + Cs^2 \quad (56) $$
<figure><img src="image_17.png" /><figcaption>(a) Dimensional demonstration</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) Non-dimensional demonstration</figcaption></figure>

<caption>**Figure 14:** A slender straight beam subjected to complex loading conditions</caption>

where

$$A = \frac{(-f_y + a_1 + b_0 + c_1)f_x - 6a_3 + 2b_1 - 6c_3}{f_x^2}; \quad B = \frac{2a_2 - b_0 + b_1 + 2c_2}{f_x}; \quad C = \frac{3a_3 - b_1 + 3c_3}{f_x} \quad (57)$$
Therefore, similar to (36), we can arrive at the exact analytical solution (58).

$$
\hat{\Gamma}(\hat{s}) = \begin{cases}
C_1 \cos(\sqrt{f_x} s) + C_2 \sin(\sqrt{f_x} s) + A + Bs + Cs^2, & (f_x > 0) \\
C_1 = \frac{f_x a_1 + f_x b_0 + f_x c_1 - f_x f_y - 6a_3 + 2b_1 - 6c_3}{f_x^2} \\
C_2 = \frac{f_x^{\frac{3}{2}} \sin(\sqrt{f_x}) C_1 + m_o f_x - 2a_2 - 6a_3 + b_0 + b_1 - 2c_2 - 6c_3}{f_x^{\frac{3}{2}} \cos(\sqrt{f_x})} \\
\left(\frac{c_3}{4} + \frac{a_3}{4} - \frac{b_1}{12}\right)s^4 + \left(\frac{a_2}{3} - \frac{b_0}{6} + \frac{b_1}{6} + \frac{c_2}{3}\right)s^3 + \left(-\frac{f_y}{2} + \frac{a_1}{2} + \frac{b_0}{2} + \frac{c_1}{2}\right)s^2 \\
+ (m - a_3 - \frac{b_1}{6} - c_3 - a_2 - \frac{b_0}{2} - c_2 + f_y - a_1 - c_1)s, & (f_x = 0) \\
D_1 e^{\sqrt{-f_x} s} + D_2 e^{-\sqrt{-f_x} s} + A + Bs + Cs^2, & (f_x < 0) \\
D_2 = \frac{e^{\sqrt{-f_x} s} (e^{\sqrt{-f_x} s} (-f_y + a_1 + b_0 + c_1) (-f_x)^{\frac{3}{2}} + 6\sqrt{-f_x} (c_3 + a_3 - \frac{b_1}{3}) e^{\sqrt{-f_x} s} - f_x (f_x m_o - 2a_2 - 6a_3 + b_0 + b_1 - 2c_2 - 6c_3))}{(-f_x)^{\frac{5}{2}} (e^2 \sqrt{-f_x} + 1)} \\
D_1 = - \frac{D_2 f_x^2 + f_x a_1 + f_x b_0 + f_x c_1 - f_x f_y - 6a_3 + 2b_1 - 6c_3}{f_x^2}
\end{cases} \quad (58)
$$

**Remark 7.** Here, we can also consider shear strain in the studied flexure especially if the flexure is stubby. According to [35], we can just add shear influence into the model via

$$
\hat{\Gamma}_s(\hat{s}) = \hat{\Gamma}(\hat{s}) + \frac{12 + 11\nu}{5} \beta (f_y \cos \hat{\Gamma}(\hat{s}) + f_x \sin \hat{\Gamma}(\hat{s}) - \int_s^1 \hat{q}_y(\hat{s}) d\hat{s} \cos \hat{\Gamma}(\hat{s})) \quad (59)
$$

where $\nu$ is the Poisson's ratio of the material and $\beta = \frac{h^2}{12L^2}$. Besides the shear influence, we can also include the axial stretch in this general BBCM following [35]

$$
\begin{align}
\hat{x}(\hat{s}) &= \int_0^{\hat{s}} \left(1 + (f_y \sin \hat{\Gamma}_s(\hat{s}) - f_x \cos \hat{\Gamma}_s(\hat{s})) - \int_s^1 \hat{q}_y(\hat{s}) d\hat{s} \sin \hat{\Gamma}(\hat{s})\right) \beta \cos \hat{\Gamma}_s(\hat{s}) d\hat{s} \tag{60} \\
\hat{y}(\hat{s}) &= \int_0^{\hat{s}} \left(1 + (f_y \sin \hat{\Gamma}_s(\hat{s}) - f_x \cos \hat{\Gamma}_s(\hat{s})) - \int_s^1 \hat{q}_y(\hat{s}) d\hat{s} \sin \hat{\Gamma}(\hat{s})\right) \beta \sin \hat{\Gamma}_s(\hat{s}) d\hat{s}
\end{align}
$$

In the following, Scenario 1 and Scenario 2 are presented where the proposed BBCM (58) is used to further model distributed moment, distributed loading and varying initial curvature within small-to-intermediate-range deflections or even large-range deflections in some specific loading conditions. Finite element method (FEM) is used to verify the results where the loading conditions and all results have been nondimensionalized as provided in Table 5 and Table 6. In particular, the beam-end coordinates are presented where BBCM is compared with FEM. The error (ER) is defined in (61).

As demonstrated in Fig. 15 and Table. 5, BBCM can be used to consider distributed moment, distributed loading and varying initial curvature. The results also prove the reliable accuracy level of BBCM within small-to-intermediate-range deflections, as stated in Section 4.2.2.

Noticing that in BVP (52) and BVP (53), the nonlinear terms in the corresponding ODEs are the terms that have $F_x(f_x)$, $F_y(f_y)$ and $q_y(\hat{q}_y)$. Logically, if we exclude the mentioned three to only consider beams of initially varying curvature $\kappa(\hat{\kappa})$ subjected to beam-end moment $M_o(m_o)$ and distributed moment $M^*(m^*)$, the original studied ODE is already linear which has the exact analytic solutions. That means in this case, the exact analytic solutions to BBCM can be used to consider large-range deflections. The details of nondimensionalized loading conditions in Scenario 2 and the corresponding results of beam-end coordinates are presented in Table 6. In Fig. 16, the BBCM results perfectly capture FEM results as well. Therefore, it is logical to conclude that BBCM can be used to model large deflections in these specific scenarios with reliable accuracy.
$$ ER = \sqrt{\frac{(\hat{x}(1)_{BBCM} - \hat{x}(1)_{FEM})^2 + (\hat{y}(1)_{BBCM} - \hat{y}(1)_{FEM})^2 + (\hat{\theta}(1)_{BBCM} - \hat{\theta}(1)_{FEM})^2}{(\hat{x}(1)_{FEM})^2 + (\hat{y}(1)_{FEM})^2 + (\hat{\theta}(1)_{FEM})^2}} \quad (61) $$

<table><caption>**Table 5**
Details of nondimensionalized loading conditions in Scenario 1</caption><thead><tr><th>Loading</th><th>Undeformed</th><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>f<sub>x</sub></td><td>0</td><td>0.25</td><td>0.1</td><td>-0.1</td></tr><tr><td>f<sub>y</sub></td><td>0</td><td>0.25</td><td>0.1</td><td>-0.1</td></tr><tr><td>m<sub>o</sub></td><td>0</td><td>0.25</td><td>0.1</td><td>-0.1</td></tr><tr><td rowspan="4">m*(ŝ)</td><td>a<sub>0</sub></td><td>0</td><td>0</td><td>0.1</td><td>-0.1</td></tr><tr><td>a<sub>1</sub></td><td>0</td><td>0</td><td>0.1</td><td>-0.1</td></tr><tr><td>a<sub>2</sub></td><td>0</td><td>0</td><td>0.1</td><td>-0.1</td></tr><tr><td>a<sub>3</sub></td><td>0</td><td>0</td><td>0.1</td><td>-0.1</td></tr><tr><td rowspan="2">q̂<sub>y</sub>(ŝ)</td><td>b<sub>0</sub></td><td>0</td><td>0</td><td>1</td><td>0.3</td></tr><tr><td>b<sub>1</sub></td><td>0</td><td>0</td><td>2</td><td>0.3</td></tr><tr><td rowspan="4">k̂(ŝ)</td><td>c<sub>0</sub></td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>c<sub>1</sub></td><td>-4.5</td><td>-4.5</td><td>-4.5</td><td>-4.5</td></tr><tr><td>c<sub>2</sub></td><td>0.3</td><td>0.3</td><td>0.3</td><td>0.3</td></tr><tr><td>c<sub>3</sub></td><td>4.5</td><td>4.5</td><td>4.5</td><td>4.5</td></tr><tr><td rowspan="3">BBCM</td><td>x̂(1)</td><td>0.9961</td><td>0.9717</td><td>0.9930</td><td>0.9596</td></tr><tr><td>ŷ(1)</td><td>0.0000</td><td>0.2268</td><td>-0.0817</td><td>-0.2030</td></tr><tr><td>θ̂(1)</td><td>-0.0250</td><td>0.3817</td><td>-0.0060</td><td>-0.4424</td></tr><tr><td rowspan="3">FEM</td><td>x̂(1)</td><td>0.9961</td><td>0.9723</td><td>0.9932</td><td>0.9607</td></tr><tr><td>ŷ(1)</td><td>0.0007</td><td>0.2245</td><td>-0.0794</td><td>-0.1993</td></tr><tr><td>θ̂(1)</td><td>-0.0250</td><td>0.3769</td><td>-0.0034</td><td>-0.4369</td></tr><tr><td>ER</td><td></td><td>0.07%</td><td>0.50%</td><td>0.35%</td><td>0.62%</td></tr></tbody></table>

### 4.2.5. Discussions

In this section, we have demonstrated the detailed derivation of the solution to BBCM for deflections within an intermediate range. The accuracy level has been proved as the same as BCM via an enumeration strategy. Besides, its feasibility of modeling CMs has also been presented via modeling a typical compliant parallelogram. As stated previously, different from BCM expressed in the global coordinate system, BBCM is written in the body frame, which makes BBCM more practical and consistent with those previously extensive studies on beam and rod theories [15]. For example, in Section 4.2.4, we can directly extend BBCM to consider distributed moment, distributed loading and initially varying curvature of beams within intermediate-range deflections. Under these complex loading conditions, BBCM can also present reliable accuracy within intermediate-range deflections. In particular, it should be noticed that BBCM can accurately characterize large deflections of slender beams under some specific loading conditions.

## 5. Conclusions

In compliant mechanisms, modeling gradually becomes more important for performance prediction and model-based optimization (design). In this paper, we focus on modeling beams of intermediate-range deflections that commonly exist in CMs, specifically in high-precision CMs, such as positioning stages. Two representative cases are studied: intermediate-range (post-)buckling and intermediate-range deflections subjected to beam-end loading and geometrical conditions where semi-analytic and analytic solutions are provided. The accuracy has been compared and also verified by the famous BCM, and the feasibility in modeling CMs has been proved via a compliant parallelogram as well. Since BBCM still stays in the body frame which is consistent with classic rod and beam theories, we demonstrate extending BBCM to consider other complex loading conditions within intermediate-range deflections. Our
<figure><img src="image_19.png" /><figcaption>**Figure 15:** Graphical results of Scenario 1</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 16:** Graphical results of Scenario 2</figcaption></figure>

future work will include deriving the analytical closed-form solution to large-deflection problems of slender beams in
more complex loading conditions.

**Declaration of Competing Interest**

The authors declare that they have no known competing financial interests or personal relationships that could
have appeared to influence the work reported in this paper.
<table><caption>**Table 6**
Details of nondimensionalized loading conditions in Scenario 2</caption><thead><tr><th>Loading</th><th>Undeformed</th><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>f<sub>x</sub></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>f<sub>y</sub></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>m<sub>o</sub></td><td>0</td><td>1</td><td>-1</td><td>-1.5</td></tr><tr><td rowspan="4">m*(ŝ)</td><td>a<sub>0</sub></td><td>0</td><td>1</td><td>-0.5</td><td>-1</td></tr><tr><td>a<sub>1</sub></td><td>0</td><td>1</td><td>-0.5</td><td>-1</td></tr><tr><td>a<sub>2</sub></td><td>0</td><td>1</td><td>-0.5</td><td>-1</td></tr><tr><td>a<sub>3</sub></td><td>0</td><td>1</td><td>-0.5</td><td>-1</td></tr><tr><td rowspan="2">q̂<sub>y</sub>(ŝ)</td><td>b<sub>0</sub></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>b<sub>1</sub></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td rowspan="4">k̂(ŝ)</td><td>c<sub>0</sub></td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>c<sub>1</sub></td><td>2</td><td>2</td><td>2</td><td>2</td></tr><tr><td>c<sub>2</sub></td><td>3</td><td>3</td><td>3</td><td>3</td></tr><tr><td>c<sub>3</sub></td><td>4</td><td>4</td><td>4</td><td>4</td></tr><tr><td>BBCM</td><td>x̂(1)</td><td>0.2722</td><td>0.14770</td><td>0.8072</td><td>0.9413</td></tr><tr><td>BBCM</td><td>ŷ(1)</td><td>0.4562</td><td>0.2338</td><td>0.2902</td><td>-0.2590</td></tr><tr><td>BBCM</td><td>θ̂(1)</td><td>4.0000</td><td>7.0833</td><td>1.9583</td><td>0.4167</td></tr><tr><td rowspan="3">FEM</td><td>x̂(1)</td><td>0.2722</td><td>0.1496</td><td>0.8066</td><td>0.9412</td></tr><tr><td>ŷ(1)</td><td>0.4555</td><td>0.2357</td><td>0.2905</td><td>-0.2585</td></tr><tr><td>θ̂(1)</td><td>4.0000</td><td>7.0833</td><td>1.9583</td><td>0.4167</td></tr><tr><td>ER</td><td></td><td>0.02%</td><td>0.04%</td><td>0.03%</td><td>0.06%</td></tr></tbody></table>

# Appendix

## Detailed derivation of Section 4.1.1

To characterize the intermediate-range buckling shape, we need to follow the procedure below. Integrating (5), we will have the following:

$$ \frac{1}{2} \left( \frac{d\hat{\theta}}{ds} \right)^2 - f_x \cos \hat{\theta}(s) = C \qquad (62) $$

where C is an unknown constant which can be determined by the boundary conditions. Precisely,

$$ C = \frac{1}{2} \left( \frac{d\hat{\theta}}{ds}(1) \right)^2 - f_x \cos \hat{\theta}(1) = -f_x \cos \alpha \qquad (63) $$

where $\frac{d\hat{\theta}}{ds}(1) = 0$ and $\hat{\theta}(1) = \alpha$, and $\alpha$ is an unknown parameter to be determined. Then, we can reach

$$ \frac{1}{2} \left( \frac{d\hat{\theta}}{ds} \right)^2 - f_x \cos \hat{\theta}(s) = -f_x \cos \alpha \qquad (64) $$

Obviously, it is not hard to find that $\hat{\theta} = 0$ ($\alpha = \hat{\theta}(1) = 0$) is the solution to (64) (Scenario 1 in Fig. 2a). Besides, manipulating (64) drives us to obtain:

$$ \frac{d\hat{\theta}}{ds} = +\sqrt{2f_x(\cos \hat{\theta} - \cos \alpha)} \qquad (65) $$

if $\alpha > 0$ (Scenario 2 in Fig. 2a) and

$$ \frac{d\hat{\theta}}{ds} = -\sqrt{2f_x(\cos \hat{\theta} - \cos \alpha)} \qquad (66) $$
if $\alpha < 0$ (Scenario 3 in Fig. 2a). Taking Scenario 2 for example, we will have

$$f(\alpha) = \int_{0}^{\alpha^{-}} \frac{1}{\sqrt{2f_{x}(\cos \hat{\theta} - \cos \alpha)}} d\hat{\theta} - \int_{0}^{\alpha^{+}} d\hat{s} \qquad (67)$$

Obviously, we can therefore obtain the value of $\alpha$ via Newton-Raphson method. As long as $\alpha$ is available, the boundary conditions yields:

$$\text{B.C. } \hat{\theta}(0) = 0; \hat{\theta}(1) = \alpha; \hat{\theta}'(1) = 0 \qquad (68)$$

Then, we will be able to arrive at

$$\hat{\theta}(\hat{s}) = \frac{\alpha}{\sin \sqrt{f_x \hat{s}}} \sin(\sqrt{f_x \hat{s}}) \qquad (69)$$

where $f_x > \frac{\pi^2}{4}$ and this is only the formulation for the first order buckling of a slender straight beam. Then, using (16) yields the exact characterization of the deflected beam shapes.

### Detailed derivation of Section 4.1.2

Similarly, following the procedure like (62) to (67) which are originally stated in [60], we will be able to obtain the value of $\alpha$. Integrating (5), we will have the following:

$$\frac{1}{2} (\frac{d\hat{\theta}}{d\hat{s}})^2 - f_x \cos \hat{\theta}(\hat{s}) = C \qquad (70)$$

where $C$ is an unknown constant which can be determined by the boundary conditions. Precisely,

$$C = \frac{1}{2} (\frac{d\hat{\theta}}{d\hat{s}}(0))^2 - f_x \cos \hat{\theta}(0) = -f_x \cos \alpha \qquad (71)$$

where $\frac{d\hat{\theta}}{d\hat{s}}(0) = 0$ and $\hat{\theta}(0) = \alpha$, and $\alpha$ is an unknown parameter to be determined. Then, we can reach

$$\frac{1}{2} (\frac{d\hat{\theta}}{d\hat{s}})^2 - f_x \cos \hat{\theta}(\hat{s}) = -f_x \cos \alpha \qquad (72)$$

Obviously, it is not hard to find that $\hat{\theta} = 0$ ($\alpha = \hat{\theta}(1) = 0$) is the solution to (64) (Scenario 1 in Fig. 2b). Besides, manipulating (64) drives us to obtain:

$$\frac{d\hat{\theta}}{d\hat{s}} = -\sqrt{2f_x(\cos \hat{\theta} - \cos \alpha)} \qquad (73)$$

if $\alpha > 0$ (Scenario 2 in Fig. 2b) and

$$\frac{d\hat{\theta}}{d\hat{s}} = +\sqrt{2f_x(\cos \hat{\theta} - \cos \alpha)} \qquad (74)$$

if $\alpha < 0$ (Scenario 3 in Fig. 2b). Taking Scenario 2 for example, we will have

$$f(\alpha) = \int_{\alpha^1}^{\alpha^+} \frac{-1}{\sqrt{2f_x(\cos \hat{\theta} - \cos \alpha)}} d\hat{\theta} - \int_{\alpha^+}^{\alpha^+} d\hat{s} \qquad (75)$$

Logically, we can therefore obtain the value of $\alpha$ via Newton-Raphson method. As long as $\alpha$ is available, the boundary conditions yields:

$$\text{B.C. } \hat{\theta}(0) = \alpha; \hat{\theta}(1) = -\alpha; \hat{\theta}'(0) = 0; \hat{\theta}'(1) = 0 \qquad (76)$$

we will have:

$$\hat{\theta}(\hat{s}) = \alpha \cos(\sqrt{f_x \hat{s}}) \qquad (77)$$

where $f_x > \frac{\pi^2}{4}$ and this is only the formulation for the first order buckling of a slender straight beam. Similarly, using (16) yields the exact characterization of the deflected beam shapes.
References

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

[2] Larry L Howell. Introduction to compliant mechanisms. *Handbook of Compliant Mechanisms*, pages 1–13, 2013.

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

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

[5] Xin Dong, Mark Raffles, Salvador Cobos Guzman, Dragos Axinte, and James Kell. Design and analysis of a family of snake arm robots connected by compliant joints. *Mechanism and Machine Theory*, 77:73–91, 2014.

[6] Ç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.

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

[8] Shixun Fan, Hua Liu, and Dapeng Fan. Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity. *Mechanical Sciences*, 9(1):161, 2018.

[9] Ke Wu and Guangbo Hao. Design and nonlinear modeling of a novel planar compliant parallelogram mechanism with general tensural-compresural beams. *Mechanism and Machine Theory*, 152:1–23, 2020.

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

[11] Huy-Tuan Pham and Dung-An Wang. A constant-force bistable mechanism for force regulation and overload protection. *Mechanism and Machine Theory*, 46:899–909, 2011.

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

[13] Myeong-Gyu Song, Hyun-Woo Baek, No-Cheol Park, Kyoung-Su Park, Taeyong Yoon, Young-Pil Park, and Soo-Cheol Lim. Development of small sized actuator with compliant mechanism for optical image stabilization. *IEEE Transactions on Magnetics*, 46(6):2369–2372, 2010.

[14] Peng Qi, Chen Qiu, Hongbin Liu, Jian S Dai, Lakmal D Seneviratne, and Kaspar Althoefer. A novel continuum manipulator design using serially connected double-layer planar springs. *IEEE/ASME Transactions on Mechatronics*, 21(3):1281–1292, 2015.

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

[16] Stuart Antman. *Nonlinear Problems of Elasticity; 2nd ed*. Springer, Dordrecht, 1996.

[17] Singiresu S Rao. The finite element method in engineering. Butterworth-heinemann, 2017.

[18] Martin Braun and Martin Golubitsky. *Differential equations and their applications*, volume 2. Springer, 1983.

[19] Michael Renardy and Robert C Rogers. *An introduction to partial differential equations*, volume 13. Springer Science & Business Media, 2006.

[20] Harry Hochstadt. *Differential equations*. Courier Corporation, 2014.

[21] Edward L Ince. *Ordinary differential equations*. Courier Corporation, 1956.

[22] Ke Wu and Gang Zheng. Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms. *Mechanism and Machine Theory*, 172:104757, 2022.

[23] Mingxiang Ling, Larry L Howell, Junyi Cao, and Guimin Chen. Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey. *Applied Mechanics Reviews*, 72(3), 2020.

[24] Pier Paolo Valentini and Ettore Pennestrì. Second-order approximation pseudo-rigid model of leaf flexure hinge. *Mechanism and Machine Theory*, 116:352–359, 2017.

[25] Slaviša Šalinić and Aleksandar Nikolić. A new pseudo-rigid-body model approach for modeling the quasi-static response of planar flexure-hinge mechanisms. *Mechanism and Machine Theory*, 124:150–161, 2018.

[26] Mattia Cera, Marco Cirelli, Luca Colaiacovo, and Pier Paolo Valentini. Second-order approximation pseudo-rigid model of circular arc flexure hinge. *Mechanism and Machine Theory*, 175:104963, 2022.

[27] Pier Paolo Valentini, Marco Cirelli, and Ettore Pennestrì. Second-order approximation pseudo-rigid model of flexure hinge with parabolic variable thickness. *Mechanism and Machine Theory*, 136:178–189, 2019.

[28] Matteo Verotti. Analysis of the center of rotation in primitive flexures: Uniform cantilever beams with constant curvature. *Mechanism and Machine Theory*, 97:29–50, 2016.

[29] Mohui Jin, Zhou Yang, Collin Ynchausti, Benliang Zhu, Xianmin Zhang, and Larry L Howell. Large-deflection analysis of general beams in contact-aided compliant mechanisms using chained pseudo-rigid-body model. *Journal of Mechanisms and Robotics*, 12(3), 2020.

[30] Mohui Jin, Benliang Zhu, Jiasi Mo, Zhou Yang, Xianmin Zhang, and Larry L Howell. A cprbm-based method for large-deflection analysis of contact-aided compliant mechanisms considering beam-to-beam contacts. *Mechanism and Machine Theory*, 145:103700, 2020.

[31] Christian Iandiorio and Pietro Salvini. Elasto-kinematics and instantaneous invariants of compliant mechanisms based on flexure hinges. *Micromachines*, 14(4):783, 2023.

[32] Horacio Ahuett-Garza, Oscar Chaides, Pedro N Garcia, and Pedro Urbina. Studies about the use of semicircular beams as hinges in large deflection planar compliant mechanisms. *Precision Engineering*, 38(4):711–727, 2014.

[33] Sinwoo Jeong and Hong Hee Yoo. Flexibility modeling of a beam undergoing large deflection using the assumed mode method. *International Journal of Mechanical Sciences*, 133:611–618, 2017.

[34] Ke Wu and Gang Zheng. Solutions to large beam-deflection problems by taylor series and padé approximant for compliant mechanisms. *Mechanism and Machine Theory*, 177:105033, 2022.

[35] Ke Wu, Gang Zheng, and Guimin Chen. Extending timoshenko beam theory for large deflections in compliant mechanisms. *Journal of Mechanisms and Robotics*, 15(6):061012, 2023.

[36] Ke Wu. *Design, Modeling, Optimization and Control of Compliant Mechanisms*. PhD thesis, Centrale lille, 2023.
[37] Guimin Chen and Fulei Ma. Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model. *Journal of Mechanical Design*, 137(2), 2015.

[38] Shorya Awtar and Shiladitya Sen. A generalized constraint model for two-dimensional beam flexures: Nonlinear strain energy formulation. *Journal of Mechanical Design*, 135(3):031003, 2013.

[39] Fulei Ma and Guimin Chen. Modeling large planar deflections of flexible beams in compliant mechanisms using chained beam-constraint-model. *Journal of Mechanisms and Robotics*, 8(2), 2016.

[40] Shiladitya Sen and Shorya Awtar. A closed-form nonlinear model for the constraint characteristics of symmetric spatial beams. *Journal of Mechanical Design*, 135(3):031003, 2013.

[41] Samir Emam and Walter Lacarbonara. Buckling and postbuckling of extensible, shear-deformable beams: Some exact solutions and new insights. *International Journal of Non-Linear Mechanics*, 129:103667, 2021.

[42] Samir Emam and Walter Lacarbonara. A review on buckling and postbuckling of thin elastic beams. *European Journal of Mechanics-A/Solids*, 92:104449, 2022.

[43] Ali H Nayfeh and P Frank Pai. *Linear and nonlinear structural mechanics*. John Wiley & Sons, 2008.

[44] Teodor Atanackovic. *Stability theory of elastic rods*, volume 1. World Scientific, 1997.

[45] James F Doyle. *Static and dynamic analysis of structures: with an emphasis on mechanics and computer matrix methods*, volume 6. Springer Science & Business Media, 1991.

[46] Hans Irschik and Johannes Gerstmayr. A continuum mechanics based derivation of reissner's large-displacement finite-strain beam theory: the case of plane deformations of originally straight bernoulli-euler beams. *Acta Mechanica*, 206(1-2):1-21, 2009.

[47] Walter Lacarbonara. *Nonlinear structural mechanics: theory, dynamical phenomena and modeling*. Springer Science & Business Media, 2013.

[48] Milan Batista. Analytical treatment of equilibrium configurations of cantilever under terminal loads using jacobi elliptical functions. *International Journal of Solids and Structures*, 51(13):2308-2326, 2014.

[49] Christian Iandiorio and Pietro Salvini. Large displacements of slender beams in plane: Analytical solution by means of a new hypergeometric function. *International Journal of Solids and Structures*, 185:467-484, 2020.

[50] Stuart S Antman. Problems in nonlinear elasticity. *Nonlinear Problems of Elasticity*, pages 513-584, 2005.

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

[52] Shorya Awtar and Alexander H Slocum. Constraint-based design of parallel kinematic xy flexure mechanisms. 2007.

[53] Shorya Awtar and Alexander H Slocum. Closed-form nonlinear analysis of beam-based flexure modules. In *ASME 2005 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference*, pages 101-110. American Society of Mechanical Engineers Digital Collection, 2005.

[54] Shorya Awtar, Kevin Shimotsu, and Shiladitya Sen. Elastic averaging in flexure mechanisms: A three-beam parallelogram flexure case study. *Journal of Mechanisms and Robotics*, 2(4), 2010.

[55] Guimin Chen, Fulei Ma, Guangbo Hao, and Weidong Zhu. Modeling large deflections of initially curved beams in compliant mechanisms using chained beam constraint model. *Journal of Mechanisms and Robotics*, 11(1), 2019.

[56] Thanh Vu Phan, Huy-Tuan Pham, et al. Design and optimization of a large-stroke compliant constant-torque mechanism. *Journal of Technical Education Science*, (68):93-100, 2022.

[57] Dalibor Petković, Nenad D. Pavlović, Shahaboddin Shamshirband, and Nor Badrul Anuar. Development of a new type of passively adaptive compliant gripper. *Industrial Robot: An International Journal*, 40(6):610-623, 2013.

[58] Annem Narayana Reddy, Nandan Maheshwari, Deepak Kumar Sahu, and GK Ananthasuresh. Miniature compliant grippers with vision-based force sensing. *IEEE Transactions on Robotics*, 26(5):867-877, 2010.

[59] Dong Sun and James K Mills. Manipulating rigid payloads with multiple robots using compliant grippers. *IEEE/ASME transactions on mechatronics*, 7(1):23-34, 2002.

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

[61] Mordecai B Rubin and AH Cardon. Cosserat theories: shells, rods and points. solid mechanics and its applications, vol 79. *Appl. Mech. Rev.*, 55(6):B109-B110, 2002.

[62] Clive L Dym, Irving Herman Shames, et al. *Solid mechanics*. Springer, 1973.