Ke Wu
email: ke.wu@inria.fr

Gang Zheng
email: gang.zheng@inria.fr

Guimin Chen
email: guimin.chen@gmail.com

Extending Timoshenko Beam Theory for Large Deflections in Compliant Mechanisms

Compliant Mechanisms (CMs) have presented its inherently advantageous properties due to the fact that CMs utilize elastic deformation of the elementary flexible members to transfer motion, force and energy. Previously, the classic 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 accuracy, 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 numerically approximate the solution. We first focus on the singlebeam 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 contributions of the influences of shear and axial strains in beam deflections. Then, we prove the feasibility of the proposed modeling strategy via mechanism synthesis for a bi-stable mechanism 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 mechanical performances of the studied CMs.

Introduction 1.Compliant mechanisms

Compliant mechanisms (CMs) have been serving as a promising alternative for current mechanical applications due to their inherent advantages [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]. Conventional rigidbody 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 deformation of the built-in elementary flexible and compliant members [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]. This powerful concept has presented many desired characteristics over rigid ones: compliance due to the flexibility of the used elastic material [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]; monolithic design that results in increased motion precision with no concerns about errors from assembly [START_REF] Farhadi Machekposhti | A review on compliant joints and rigid-body constant velocity universal joints toward the design of compliant homokinetic couplings[END_REF], simplified manufacturing process [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] [2] [START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF], reduced weight [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] [3], cost reduction [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] [2] [START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF] and less maintenance [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] [START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF]. These mentioned desired characteristics have made CMs a gradually mature component in many mechanical applications, such as high-accuracy positioning motion stages [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF], compliant kinematic joints [START_REF] ¸agıl | Steel compliant cardan universal joint[END_REF], bi/tri-stable mechanisms [START_REF] Qiu | A curved-beam bistable mechanism[END_REF] and other more complex applications that integrate the former three [START_REF] Trease | Design of adaptive and controllable compliant systems with embedded actuators and sensors[END_REF].

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

As mentioned in the previous section, the built-in flexible members, which are normally straight beams and initially curved beams (ICBs) [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF], play an important role in realizing the functions of CMs. Therefore, accurately and efficiently modeling beams will be a common focus either for performance predictions or model-based optimization [START_REF] Zhu | Design of compliant mechanisms using continuum topology optimization: A review[END_REF]. There are two major aspects: first, the fundamental modeling theory should be properly chosen; second, and the corresponding 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 [START_REF] Cecil Zienkiewicz | The finite element method: solid mechanics[END_REF]. This is 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 highdimensional differential equation to govern the overall deformation of the studied object. 2. Euler-Bernoulli beam theory: in the current literature, Euler-Bernoulli beam theory [START_REF] Stephen | Theory of elasticity[END_REF] 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 [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF], and then Chen proposed a comprehensive elliptical solution for large beam-deflection problems for CMs [START_REF] Zhang | A comprehensive elliptic integral solution to the large deflection problems of thin beams in compliant mechanisms[END_REF]; 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 [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]; 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) [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF]; chained BCM [START_REF] Chen | Modeling large spatial deflections of slender bisymmetric beams in compliant mechanisms using chained spatial-beam constraint model[END_REF] and chained PRBM [START_REF] Su | A pseudorigid-body 3r model for determining large deflection of cantilever beams subject to tip loads[END_REF] 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 [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF] [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 [START_REF] Stephen | Theory of elasticity[END_REF]. 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 [START_REF] Chen | Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model[END_REF], 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 highdimensional 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 [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF].

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 [START_REF] Wu | A comprehensive static modeling methodology via beam theory for compliant mechanisms[END_REF]. It's straightforward and simple since it only considers the 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 [START_REF] Chen | Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model[END_REF]. For example, compliant parallelograms work under some stretching payload [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF], and the forward critical force needs to be precisely characterized in bi-stable mechanisms [START_REF] Chen | Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model[END_REF].

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 [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF]. Via this modeling strategy, we will be able to synthesize CMs as well as conduct model-based topology optimization for CMs more efficiently.

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

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:

d 3 y dx 3 = d 2 θ dx 2 = V (x) EI = dM(x) dx EI d 2 y dx 2 = dθ dx = M(x) EI dy dx = θ(x) = x 0 M(x) EI dx - k t V (x) GA = x 0 M(x) EI dx - k t dM(x) dx GA (1) with M(x) = M o + F x b + F y (a -x) (2)
where E denotes the Young's modulus; I denotes the second moment of inertia of the cross section (I = 1 12 h 3 w); A denotes the area of the cross section (A = wh); G denotes the shear modulus (G = E 2(1+ν) for isotropic materials); k t denotes Timoshenko shear coefficient [START_REF] Cowper | The shear coefficient in timoshenko's beam theory[END_REF] (k t = 12+11ν 10+10ν ) and ν is the Poisson's ratio of the material; d 2 y 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 dy dx is small, without considering geometric nonlinearity. Therefore, Eq. ( 2) can be simplified:

M(x) = M o + F y (L -x) (3)
where a ≈ L and b ≈ 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 curvature 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 deflections. Here, we proceed with M o being 0, which results in a very common beam-deflection case, beam-end concentrated loading: M(x) = F y (Lx). Rearranging this equation and Eq. ( 1), we will be able to arrive at the following:

dy dx = θ(x) = k t F y GA + F y L EI x - F y 2EI x 2 y(x) = k t F y GA x + F y L 2EI x 2 - F y 6EI x 3 (4)
Remark 1. In Eq. (1), 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.

k t V (x)
GA denotes the extra rotation of the Euler-Bernoulli beam under beam-end loading 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:

d 2 y dx 2 = dθ dx = M(x) EI = F y (L -x) EI dy dx = θ(x) = x 0 M(x) EI dx = θ(x) = F y L EI x - F y 2EI x 2 y(x) = x 0 F y L EI x - F y 2EI x 2 dx = F y L 2EI x 2 - F y 6EI x 3 (5)
As mentioned before, both linear Timoshenko beam theory and linear Euler Bernoulli beam theory are valid only if analyzing small deflections. Precisely, since the curva-

ture κ(x) = d 2 y dx 2 [1 + ( dy dx ) 2 ] -3/2
, then by assuming the Yaxis deflection of the beam is small, i.e. dy dx ≈ 0, we have:

κ(x) = d 2 y dx 2
Similarly, x is chosen as the independent variable due to the fact that the deformed beam axis still approximately 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 mentioned previous work above (linear Timoshenko beam theory and linear Euler Bernoulli beam theory) is all based on the small-deflection assumption where the linearized formulation of curvature is only valid for small deflections.

Geometrically nonlinear Timoshenko beam theory

with axial linearly elastic assumption Here, we would like to start with geometrically nonlinear 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 (θ b (s) and θ s (s)) of each cross section caused by bending and shear strains respectively. For the curvature of beams due to bending, we have the same linearly elastic assumption as Euler Bernoulli beam equation:

dθ b ds = M(s) EI(s) + 1 R(s) (6)
where

M(s) = M o + F y (a -x(s)) + F x (b -y(s)) (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:

θ s (s) = k t GA(s) (F y cos θ(s) + F x sin θ(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:

θ(s) = θ b (s) + θ s (s) (8)
If we take one step further to consider the linearly elastic elongation of each cross section, we have

ε a (s) = F y sin θ(s) -F x cos θ(s) EA(s) (9)
For each cross section, its infinitesimal length can be reformulated as (1 + ε a (s))ds. Therefore, we can arrive at the following by considering axial strain in deflected beams:

dx ds (s) = (1 + ε a (s)) cos θ(s); dy ds (s) = (1 + ε a (s)) sin θ(s) (10)
Differentiating Eq. ( 6) considering the relationships of Eqs. [START_REF] Trease | Design of adaptive and controllable compliant systems with embedded actuators and sensors[END_REF] [START_REF] Cecil Zienkiewicz | The finite element method: solid mechanics[END_REF], we can arrive at:

d 2 θ b ds 2 = dM(s) ds EI(s) -M(s)E dI(s) ds E 2 I(s) 2 + -dR(s) ds R(s) 2 = dM(s) ds EI - M(s) EI(s) 2 dI(s) ds + -dR(s) ds R(s) 2 (11) d 2 θ b ds 2 = dM(s) ds EI - M(s) EI(s) 2 dI(s) ds + -dR(s) ds R(s) 2 = - F y EI dx ds (s) - F x EI dy ds (s) + -dR(s) ds R(s) 2 (12)
Besides, there are 2 boundary conditions defined from the physical conditions:

θ b (0) = 0 and dθ b ds (L) = M o EI + 1 R(L)
. Rearranging Eqs. (8)(10)(12), we will be able to obtain the final formulation of the desired boundary value problem:

D.E. d 2 θ b ds 2 = - F y EI dx ds (s) - F x EI dy ds (s) + -dR(s) ds R(s) 2 θ s (s) = k t GA (F y cos θ(s) + F x sin θ(s)) θ(s) = θ b (s) + θ s (s) B.C. θ(0) = 0, and dθ b ds (L) = M o EI + 1 R(L) (13)
BVP [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF] 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 [START_REF] Wu | A comprehensive static modeling methodology via beam theory for compliant mechanisms[END_REF]. In BVP [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF], 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-Free end Pinned end Roller Fixed end Sliding Fig. 3: Different types of beam ends acterization of the deflected beam shape, we derive the following via Eqs. ( 10):

x(s) = s 0 dx ds (ξ)dξ; y(s) = s 0 dy ds (ξ)dξ (14)
As shown in BVP [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF] 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 (15)
with

E pt = L 0 1 2 EI( dθ ds ) 2 ds and E pa = L 0 1 2 EA( F y sin θ(s)-F x cos θ(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 dθ ds according to BVP [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF].

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.

Prerequisites

First, we proceed to nondimensionalize the Timoshenko beam governing equation [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF]. Therefore, the dependent variable θ and the independent variable s are nondimensionalized via:

θ = θ 1 = θ b + θ s 1 ; θ = 1 θ; θ b = 1 θb ; θ s = 1 θs ŝ = s L ; s = L ŝ ( 16
)
Note that θ, θ b and θ s are already non-dimensional but we still mark θ, θb and θs as the nondimensionalized ones for universal notation in this paper. Similarly, we have the same case for ε a (s) in Eq. ( 9) and Eq. ( 10):

d x d ŝ ( ŝ) = (1 + εa ( ŝ)) cos θ( ŝ); d ŷ d ŝ ( ŝ) = (1 + εa ( ŝ)) sin θ( ŝ) (17) where εa ( ŝ) = f y sin θ( ŝ)-f x cos θ( ŝ) 12
η 2 , with η being the slenderness ratio defined as η = h L . Then, we can nondimensionalize (13) using ( 16) accordingly:

D.E. d 2 θb d ŝ2 = -f y d x d ŝ ( ŝ) -f x d ŷ d ŝ ( ŝ) + -dr( ŝ) d ŝ r(s) 2 θs ( ŝ) = (12 + 11ν) 60 η 2 ( f y cos θ( ŝ) + f x sin θ( ŝ)) θ( ŝ) = θs ( ŝ) + θb ( ŝ) B.C. θ(0) = 0, and d θb d ŝ (1) = m o + 1 r(1) (18)
where

f x = F x L 2 EI ; f y = F y L 2 EI ; m o = M o L EI ; r( ŝ) = R(s) L ; η = h L ; ŝ ∈ [0, 1]; s ∈ [0, L].
In terms of characterizing the nondimensionalized beam shape, we can use the following nondimensionalized formulation according to Eq. ( 17): x( ŝ) = Remark 2. As long as the isotropic material used is settled, BVP (18) gradually converges to the following if η → 0:

D.E. d 2 θ d ŝ2 = -( f y cos θ( ŝ) + f x sin θ( ŝ)) + -dr( ŝ) d ŝ r(s) 2 B.C. θ(0) = 0 d θ d ŝ (1) = m o + 1 r(1) (19)
with the nondimensionalized formulation of characterizing the beam shape:

x( ŝ) = ŝ 0 cos θ( ξ)dξ; ŷ( ŝ) = ŝ 0 sin θ( ξ)dξ ( 20
)
The above BVP [START_REF] Wu | Solutions to large beamdeflection problems by taylor series and padé approximant[END_REF] with Eq. ( 20) is exactly the nondimensionalized governing equation of Euler Bernoulli beam theory [START_REF] Wu | Solutions to large beamdeflection problems by taylor series and padé approximant[END_REF]. We can easily draw the conclusion that higher 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 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:

D.E. d 2 θb d ŝ2 = -f y d x d ŝ ( ŝ) -f x d ŷ d ŝ ( ŝ) + -dr( ŝ) d ŝ r(s) 2 B.C. θ(0) = 0, and d θ d ŝ (1) = m o + 1 r(1) (21)
BVP (21) along with Eq. ( 17) is termed as modified geometrically nonlinear Euler Bernoulli beam theory incorporating axial elongation.

Numerical modeling through Galerkin method

As stated in BVP [START_REF] Cowper | The shear coefficient in timoshenko's beam theory[END_REF] and BVP [START_REF] Wu | Solutions to large beamdeflection problems by taylor series and padé approximant[END_REF], the governing equations of beam deflection are essentially a highly nonlinear BVP. Compared to the simple BVP [START_REF] Wu | Solutions to large beamdeflection problems by taylor series and padé approximant[END_REF] (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 [START_REF] Cowper | The shear coefficient in timoshenko's beam theory[END_REF] 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 [START_REF] Cowper | The shear coefficient in timoshenko's beam theory[END_REF].

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 contributor. Therefore, we can simplify the governing equation in BVP [START_REF] Cowper | The shear coefficient in timoshenko's beam theory[END_REF] by assuming θ b ≈ θ, equivalently θb ≈ θ. On the contrary, we closely look at BVP [START_REF] Cowper | The shear coefficient in timoshenko's beam theory[END_REF] again in the small-deflection cases. Obviously, cos θ b ≈ cos θ ≈ 1 and sin θ b ≈ sin θ ≈ 0 when the magnitude of θ is small. Therefore, we can still assume θ b ≈ θ in small-deflection cases. Then, we can rewrite Eqs. [START_REF] Chen | Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model[END_REF]:

d x d ŝ ( ŝ) ≈ d x d ŝ * ( ŝ) = (1 + ε * a ( ŝ)) cos θb ( ŝ) d ŷ d ŝ ( ŝ) ≈ d ŷ d ŝ * ( ŝ) = (1 + ε * a ( ŝ)) sin θb ( ŝ) (22)
where ε *

a ( ŝ) = f y sin θb ( ŝ)-f x cos θb ( ŝ) 12
η 2 . Now, we can logically rewrite the governing BVP into the following:

D.E. d 2 θb d ŝ2 = -f y d x d ŝ * ( ŝ) -f x d ŷ d ŝ * ( ŝ) + -dr( ŝ) d ŝ r(s) 2 θs ( ŝ) = (12 + 11ν) 60 η 2 ( f y cos θb ( ŝ) + f x sin θb ( ŝ)) B.C. θb (0) = 0, and d θb d ŝ (1) = m o + 1 r(1) (23)
Essentially speaking, we bypass the need of implicitly approximating θ s in BVP [START_REF] Cowper | The shear coefficient in timoshenko's beam theory[END_REF] to explicitly approximate θ s through BVP [START_REF] Timoshenko | Theory of elastic stability 2e[END_REF], which hugely decreases the complexity of approximation process. As proposed in [START_REF] Wu | Insight into numerical solutions of static large deflection of general planar beams for compliant mechanisms[END_REF], we have up to nearly 10 methods to treat this BVP. In this paper, we use Galerkin method to handle BVP [START_REF] Timoshenko | Theory of elastic stability 2e[END_REF], and a brief recap is presented here.

According to Weierstrass's first theorem, a given differentiable function for a random closed interval can be always approximated by a polynomial of a proper order. Therefore, we then define a function θb ( ŝ) which is an nth-order polynomial, to approximate the solution of the studied BVP [START_REF] Timoshenko | Theory of elastic stability 2e[END_REF]: θb ( ŝ) = ∑ n i=0 c i ŝi for i = 0, 1, ...n where there are n + 1 unknown coefficients to solve. Next, we can define an approximating residual function of BVP [START_REF] Timoshenko | Theory of elastic stability 2e[END_REF] utilizing the polynomial solution θb ( ŝ) as follows:

R e ( ŝ) = d 2 θb d ŝ2 + f y d x d ŝ * ( ŝ) + f x d ŷ d ŝ * ( ŝ) + dr( ŝ) d ŝ r(s) 2 (24)
The essence of weighted residual method is to force the residual Eq. (24) to zero in some average sense over the domain: 1 0 W i R e ( ŝ)d ŝ = 0. Here, we choose W i taken from the derivative of the approximating function θb ( ŝ) according to Galerkin method: W i ( ŝ) = ∂ θb ( ŝ) ∂c i = ŝi for i = 1, 2, 3...n + 1. Besides, θb ( ŝ) also needs to satisfy the boundary conditions defined in BVP [START_REF] Timoshenko | Theory of elastic stability 2e[END_REF]: θb (0) = 0 and d θb d ŝ (1) = m o + 1 r(1) . So far, we have altogether n + 3 equations to solve n + 1 polynomial coefficients c i (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 θb ( ŝ). 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 [START_REF] Wu | A comprehensive static modeling methodology via beam theory for compliant mechanisms[END_REF]. 2. Indirect way: we choose a globally initial guess for θb ( ŝ), such as trigonometric functions, which yields the corresponding c i for the θb ( ŝ). 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 [START_REF] Wu | Theoretical analysis on nonlinear buckling, post-buckling of slender beams and bi-stable mechanisms[END_REF].

As long as we have the approximate solution θb ( ŝ) available, we can explicitly obtain the approximation of θs (s) through the closed-form formulation below according to BVP [START_REF] Timoshenko | Theory of elastic stability 2e[END_REF]:

θs ( ŝ) = (12 + 11ν) 60 η 2 ( f y cos θb ( ŝ) + f x sin θb ( ŝ)) (25)
Finally, we can reach the modified approximate solution via:

θ( ŝ) = θb ( ŝ) + θs ( ŝ) (26)
Finally, replacing θ( ŝ) with θ( ŝ), we can just follow Eqs.

(17)(26) to characterize the deflected beam shape.

3.3 Deflection analysis on shear strain and axial elongation 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 models, Timoshenko beam models and solid mechanics. All the models are used to compare the results of the proposed modeling strategy (NTBA-GM), BVP [START_REF] Timoshenko | Theory of elastic stability 2e[END_REF] where solid mechanics 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 condition is schematically demonstrated in Fig. 4.

Besides, under this nondimensionalization framework, the analytical solutions of both linear Timoshenko and Euler Bernoulli beam theories Eq. ( 4) and Eq. ( 5) need to be nondimensionalized for more general comparisons as well. Eq. (29) Eq. ( 30)

Bending strain within small-deflection range [START_REF] Oliver | Euler-bernoulli beam theory[END_REF] NEB-GM Geometrically nonlinear Euler Bernoulli beam theory (Galerkin method) BVP [START_REF] Wu | Solutions to large beamdeflection problems by taylor series and padé approximant[END_REF] Bending strain within large-deflection range [START_REF] Wu | A comprehensive static modeling methodology via beam theory for compliant mechanisms[END_REF] EBA-FEM Linear Euler Bernoulli beam theory incorporating axial elongation (Finite element method)

See reference

Bending and axial strains within large-deflection range [START_REF] Ab Comsol | Structural mechanics module user's guide[END_REF] NEBA-GM Geometrically nonlinear Euler Bernoulli beam theory incorporating axial elongation (Galerkin method) BVP [START_REF] Oliver | Euler-bernoulli beam theory[END_REF] Bending and axial strains within large-deflection range N/A TB-AS Linear Timoshenko beam theory (Analytical solution)

Eq. ( 27) Eq. ( 28)

Bending and shear strains within small-deflection range [START_REF] Timoshenko | Theory of elastic stability 2e[END_REF] TBA-FEM Linear Timoshenko beam theory incorporating axial elongation (Finite element method)

See reference

Bending, shear and axial strains within large-deflection range [START_REF] Ab Comsol | Structural mechanics module user's guide[END_REF] NTBA-GM Geometrically nonlinear Timoshenko beam theory incorporating axial elongation (Galerkin method) BVP ( 23) Bending, shear and axial strains within large-deflection range N/A SD-FEM Solid mechanics (FEM)

See reference

All sorts of deformation within large-deflection range [START_REF] Ab Comsol | Structural mechanics module user's guide[END_REF] 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 nondimensionalization, Eq. ( 4) can be rewritten as the following:

d ŷ d x ( x) = θ( x) = 12 + 11ν 60 η 2 f y + f y x - 1 2 f y x2 ŷ( x) = 12 + 11ν 60 η 2 f y x + 1 2 f y x2 - 1 6 f y x3 (27)
To evaluate the beam-end rotation and vertical coordinate, we set x = 1:

d ŷ d x (1) = θ(1) = f y ( 12
+ 11ν 60 η 2 + 1 2 ) ŷ(1) = f y ( 12
+ 11ν 60 η 2 + 1 3 ) (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: where ŷ(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 η 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.

d ŷ d x = θ( x) = f y x - 1 2 f y x2 , and ŷ( x) = 1 2 f y x2 - 1 6 f y x3 (

Deflection analysis of a circularly curved beam

Similarly, we assume ν = 0.33,

f x = 0, f y = 0.2, m o = 0, r = 2 π , w L = 0.1 and η = h L = 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 ŷ(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 η increases, ŷ(1) increases as well, which numerically validates the theoretical conclusion drawn from BVP (18): larger η 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 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 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: bistable compliant mechanisms and compliant parallelograms under the proposed modeling framework BVP [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF] 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-beambased 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.

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 builtin flexible members are modeled BVP [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF]. 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 [START_REF] Timoshenko | Theory of elastic stability 2e[END_REF] is used here for dimensional modeling.

As shown in Fig. 7, the bi-stable mechanism is com- Here, we are interested in the relationship between F o and the displacement ∆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.

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:

D.E. d 2 θ b1 ds 2 = - F y1 EI dx 1 ds (s) - F x1 EI dy 1 ds (s) θ s1 (s) = k t GA (F y1 cos θ 1 (s) + F x1 sin θ 1 (s)) θ 1 (s) = θ b1 (s) + θ s1 (s) B.C. θ 1 (0) = θ o and dθ b1 ds (L 1 ) = M o1 EI (31)
For Beam #2, we have a similar BVP as follows:

0 0.
θ 2 (s) = θ b2 (s) + θ s2 (s) B.C. θ 2 (0) = θ 1 (L 1 ); θ 2 (L 2 ) = θ o and dθ b2 ds (L 2 ) = M o EI (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) = s 0 dx 1 dξ (ξ)dξ, y 1 (s) = s 0 dy 1 dξ (ξ)dξ, x 2 (s) = s 0 dx 2
dξ (ξ)dξ and y 2 (s) = s 0 dy 2 dξ (ξ)dξ. Besides, we also have the following forceequilibrium relationships:

F x = F x1 , F y = F y1 , F o = 2F y and M o1 = M o + F y (L r cos θ 1 (L 1 ) + x 2 (L 2 )) + F x (L r sin θ 1 (L 1 ) + y 2 (L 2 )
) and the geometric relationships:

(L 1 + L r + L 2 ) cos θ o = x 1 (L 1 ) + L r cos θ 1 (L 1 ) + x 2 (L 2 ) and ∆Y = (L 1 + L r + L 2 ) sin θ o -(y 1 (L 1 ) + L r sin θ 1 (L 1 ) + y 2 (L 2 ))
where ∆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).

Numerical results

To verify the modeling methodology, we use solidmechanics-based FEM as the reference. First, the geometric and material parameters chosen for the studied straight beam is shown in the following: E = 200 × 10 9 Pa; ν = 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; θ o = 0.1266 rad; k t = 12+11ν 10+10ν ; G = E 2(1+ν) ; 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 -∆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 ∆Y as well (see Fig. 10).

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 θ 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 θ 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) [START_REF] Qiu | A curved-beam bistable mechanism[END_REF] [17] [START_REF] Wu | Theoretical analysis on nonlinear buckling, post-buckling of slender beams and bi-stable mechanisms[END_REF].

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 postbuckling 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 [START_REF] Wu | Theoretical analysis on nonlinear buckling, post-buckling of slender beams and bi-stable mechanisms[END_REF], 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.

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: θ g1 (s) = π 2 -301.5929s 3 and θ g2 (s) = π 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: s) . In

dθ g1 ds (s) = -904.7787s 2 = 1 R 1 (s) and dθ g2 ds (s) = 904.7787s 2 = 1 R 2 (

Static modeling

Here, we proceed to model the simplified and equivalent mechanism presented in Fig. 12. Logically, according to BVP [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF], 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

Numerical results

To verify the modeling methodology, we use solidmechanics-based FEM as the reference. First, the geometric and material parameters chosen for the studied ICB is shown in the following:f E = 200 × 10 9 Pa, ν = 0.30, w = 0.015m, h = 0.005m, L 1 = 0.25m, k t = 12+11ν 10+10ν , G = E 2(1+ν) . To analyze the performance of the studied parallelogram mechanism, the relationship between the force F o and its corresponding displacement ∆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.

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

Fig. 1 :

 1 Fig. 1: Small (linear) deflection of Timoshenko beam compared to Euler-Bernoulli beam subjected to beam-end loading

Fig. 2 :

 2 Fig. 2: Large deflection of Timoshenko beam compared to Euler-Bernoulli beam under beam-end loading

  where x( ŝ) and ŷ( ŝ) 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.

Fig. 4 :loading conditions 1 . 1 rFig. 5 :Fig. 6 :

 41156 Fig. 4: Studied loading conditions

2 RFig. 7 :Fig. 8 : 7 .

 2787 Fig. 7: Schematic diagram of bi-stable compliant mechanisms

Fig. 11 :

 11 Fig. 11: Schematic diagram of the studied ICB-based compliant parallelogram

Fig. 12 :Fig. 13 :

 1213 Fig. 12: Detailed force diagram of the studied ICB-based compliant parallelogram

Fig. 14 :

 14 Fig. 14: Force-displacement relationship of the studied ICB-based parallelogram mechanism

Table 1 :

 1 Acronyms of different beam models

		Description	Governing equations	Model considerations	Refs
	EB-AS	Linear Euler Bernoulli beam theory (Analytical solution)

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 COSSEROOTS [ANR-20-CE33], by the project Inventor (I-SITE ULNE, le programme d'Investissements d'Avenir, Métropole Européenne de Lille), France.