Ke Wu

Gang Zheng

A Comprehensive Static Modeling Methodology via Beam Theory for Compliant Mechanisms ?

Keywords:

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.

Introduction 1.Compliant mechanisms

Compliant Mechanisms (CMs) have gradually become a promising alternative for designing mechanical applications, attracting more and more research attention over the last few decades [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]. Conventional mechanisms are able to transfer motion, force and energy via the cooperative operation of the involved rigid links and kinematic pairs. However, CMs can present the same function only depending on the deformation of elementary flexible members [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF][START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF]. Among the existing designs, CMs have demonstrated several more desirable features than conventional rigidbody mechanisms: simplified manufacturing process [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF][START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF][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]; lighter in weight [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF][START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF]; less maintenance needed [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF][START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF] and better motion precision [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF]. The above mentioned advantageous characteristics have paved the way for CMs to be integrated in a series of mechanical designs, such as high-accuracy positioning motion stages [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF][START_REF] Gräser | On the influence of the flexure hinge orientation in planar compliant mechanisms for ultra-precision applications[END_REF][START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF][START_REF] Wu | Design and nonlinear modeling of a novel planar compliant parallelogram mechanism with general tensuralcompresural beams[END_REF], compliant kinematic joints [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][START_REF] Dong | Design and analysis of a family of snake arm robots connected by compliant joints[END_REF][START_REF] Çagıl Merve Tanık | Steel compliant cardan universal joint[END_REF], bi/tri-stable mechanisms [START_REF] Qiu | A curved-beam bistable mechanism[END_REF][START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF][START_REF] Sam Han | Design, simulation, and fabrication of a quadstable monolithic mechanism with x-and y-directional bistable curved beams[END_REF][START_REF] Pham | A constant-force bistable mechanism for force regulation and overload protection[END_REF][START_REF] Brian T Edwars | A pseudo-rigid-body model for functionally binary pinnedpinned segments used in compliant mechanisms[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][START_REF] Song | Development of small sized actuator with compliant mechanism for optical image stabilization[END_REF][START_REF] Qi | A novel continuum manipulator design using serially connected double-layer planar springs[END_REF]. In the current literature, the majority of studies on CMs utilize straight beams and ICBs as the elementary flexible members: straight beams provide relatively ideal DoF and DoC, being a suitable choice of high-accuracy positioning CMs [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF][START_REF] Gräser | On the influence of the flexure hinge orientation in planar compliant mechanisms for ultra-precision applications[END_REF][START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF][START_REF] Wu | Design and nonlinear modeling of a novel planar compliant parallelogram mechanism with general tensuralcompresural beams[END_REF]; ICBs are commonly used in large-deflection CMs, such as bi/tri-stable CMs [START_REF] Qiu | A curved-beam bistable mechanism[END_REF][START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF][START_REF] Sam Han | Design, simulation, and fabrication of a quadstable monolithic mechanism with x-and y-directional bistable curved beams[END_REF][START_REF] Pham | A constant-force bistable mechanism for force regulation and overload protection[END_REF][START_REF] Brian T Edwars | A pseudo-rigid-body model for functionally binary pinnedpinned segments used in compliant mechanisms[END_REF] and path-generating CMs [START_REF] Kumar Rai | Synthesis of path generating compliant mechanisms using initially curved frame elements[END_REF], since ICBs have a large deflection range along with a relatively small strain range. To conduct static synthesis of a whole CM, the efficient methods for modeling nonlinear (intermediate or large) deflection of the mentioned two types of slender beams are therefore highly needed. Fortunately, several impressive and great contributions have been completed by the pioneering researchers, which are reviewed in the next section.

Nonlinear modeling for beam deflection problems 1.2.1. Commonly used modeling methods for CMs

In terms of modeling slender beams, many practically useful methods have been reported in the literature [START_REF] Ling | Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey[END_REF]. In 1940s, the authors of [START_REF] Hj Barten | On the deflection of a cantilever beam[END_REF] and [START_REF] Ke Bisshopp | Large deflection of cantilever beams[END_REF] explored using elliptical integrals to analytically solve large beam deflection problems where only vertical beam-end loading is considered. Howell introduced this elliptical-integral method to analytically solve the geometrically nonlinear Euler Bernoulli beam equation (essentially a BVP of an ODE) where a more general beam-end loading is considered, and this has been serving as the most accurate solution to this beam theory so far [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]. For easy-to-follow implementation, Howell and his co-workers then introduced a more designerfriendly method called pseudo-rigid-body model (PRBM) to analyze large deformation of slender straight beams and ICBs, which is essentially a simplified model via introducing the concept of rigid-body mechanism [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]. This original PRBM enables designers to handle efficient modeling in early stages of the design process, which has proved its effectiveness in many circumstances. Another remarkable contribution called Beam Constraint Model (BCM) for slender straight beams was first proposed in Dr. Shorya Awtar's PhD thesis [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF] along with his related contribution [START_REF] Awtar | A generalized constraint model for two-dimensional beam flexures: Nonlinear strain energy formulation[END_REF] for ICBs, serving as an effective tool for analyzing straight flexible members in CMs within intermediate deformation range. Dr. Chen Guimin then extended this work and proposed Timoshenko Beam Constraint Model (TBCM) for the purpose of characterizing stubby beams (which takes the influence of shear force into account) [START_REF] Chen | Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model[END_REF].

Considering large beam-deflection scenarios, Finite Element Analysis (FEA) [START_REF] Singiresu | The finite element method in engineering[END_REF] is a reliable choice since it has already presented its superiority in accuracy regarding modeling any type of deformation and been widely used as a means for verifying other modeling methods. In particular, FEA through total Lagrange approach serves as a feasible choice to analyze the nonlinear dynamic and static deformation of studied objects [START_REF] Thomas | Finite element model for dynamic analysis of timoshenko beam[END_REF][26] [START_REF] Bathe | Finite element formulations for large deformation dynamic analysis[END_REF], which can be a powerful tool to analyze the properties of CMs as well. Besides, the comprehensive review presented in [START_REF] Ling | Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey[END_REF] has mentioned another type of modeling method: chained algorithm. Here, we would like to recall its original definition in Howell's classic work [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]: "The chained algorithm requires discretization of the object being modeled into beam elements and each element is analyzed in succession, the chained algorithm is similar to Finite Element Analysis and lends itself well to CM analysis, the chained algorithm is so named because it requires discretization of the object being modeled into beam elements and analyzes each element in succession". This definition clearly implies that chained algorithm actually introduces the concept of FEA where the summation of relatively smaller deformation of all discretized elements results in feasible analysis of larger deflection (geometric nonlinearity). Therefore, we can logically consider Chained Beam Constraint Model (CBCM), nR-PRBM and some multiple-segment modeling method [START_REF] Anil | A general and efficient multiple segment method for kinetostatic analysis of planar compliant mechanisms[END_REF] as FEA-related methods or chained algorithms: 3R-PRBM are available [START_REF] Su | A pseudorigid-body 3r model for determining large deflection of cantilever beams subject to tip loads[END_REF][START_REF] Yu | A pseudo-rigid-body 2r model of flexural beam in compliant mechanisms[END_REF][START_REF] Chen | Finding the optimal characteristic parameters for 3r pseudo-rigid-body model using an improved particle swarm optimizer[END_REF] (three rotational pivots are involved), and more nR-PRBM (n > 3) serve as feasible options [START_REF] Zhu | Pseudo-rigid-body model for the flexural beam with an inflection point in compliant mechanisms[END_REF][START_REF] Yu | A novel model of large deflection beams with combined end loads in compliant mechanisms[END_REF][START_REF] Yu | 5r pseudo-rigid-body model for inflection beams in compliant mechanisms[END_REF][START_REF] Jin | A cprbm-based method for large-deflection analysis of contact-aided compliant mechanisms considering beam-to-beam contacts[END_REF][START_REF] Jr | A 3-d chain algorithm with pseudo-rigid-body model elements[END_REF]. Similarly, nR-PRBM are also proposed for ICBs as well [START_REF] Venkatasubramanian | Pseudo-rigid-body models of initially-curved and straight beams for designing compliant mechanisms[END_REF]. Chained beam constraint model (CBCM) has been proved to be an effective choice for modeling large deformation of straight beams [START_REF] Chen | Modeling large spatial deflections of slender bisymmetric beams in compliant mechanisms using chained spatial-beam constraint model[END_REF][START_REF] Ma | Modeling large planar deflections of flexible beams in compliant mechanisms using chained beam-constraintmodel[END_REF] as well as ICBs [START_REF] Chen | Modeling large deflections of initially curved beams in compliant mechanisms using chained beam constraint model[END_REF].

The above mentioned work has been a great success due to their effectiveness in modeling CMs, however some of them may suffer from their own inherent limitations. For example, although FEA and chained algorithms are accurate and reliable, they are based on the premise of increasing computational expense. In other words, they are high dimensional models from the viewpoint of mathematics. Obviously, BCM and TBCM are specifically used for modeling intermediate range of deflection, which are perfect choices for flexure-based CMs with a relatively small stroke such as highly accurate positioning stages, but logically not ideal options for large-deflection CMs.

Modeling via beam theory

To have a better understanding of modeling slender beams, we need to step back to focus on the classic Euler-Bernoulli beam equation:

EI d✓ ds = M(s) ( 1
)
where E is the elastic modulus and I is the second moment inertia of the cross-section area, which assumes the linear relationship between the beam curvature d✓ ds and the exerted moment M at (x, y) or equivalently s along the beam axis. For a general loading condition displayed in Fig. 1, a general slender beam is in its deformed configuration under an arbitrary force F o and an arbitrary moment M o at its beam end (a, b) or equivalently L along the beam axis.

There are some implicit assumptions made in classic Euler-Bernoulli beam equation: Eq. 1 assumes that beams are inextensible (L is constant) and ignores the influence of shear deformation of cross sections. Besides, it only takes care of planar deflection of slender beams. Therefore, our work will follow the assumptions accordingly.

For easier implementation, F o is decomposed into a horizontal force F x and a vertical force F y . Given the above configuration, the moment M at any point (x, y) along the beam can be formulated as follows:

EI d✓ ds = M(x, y) = M(s) = F y (a * x(s)) + F x (b * y(s)) + M o (2)
Then, conducting differentiation on Eq. 2 with respect to s results in the following equation: for ODE [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]. Summarizing what we have stated above, the beam-deflection problem can be described by:

d 2 ✓ ds 2 = * F y EI (cos ✓(s) + F x F y sin ✓(s)) (3)
D.E. d 2 ✓ ds 2 = * F y EI (cos ✓(s) + F x F y sin ✓(s)) B.C. ✓(0) = 0 d✓ ds (L) = M o EI (4)
which is essentially a nonlinear BVP of ODE (3).

Methods to solve BVP

We would like to emphasize that analytic and exact solution of the above BVP is highly recommended in practice, like the proposed elliptical integral method [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF] , which however is semi-analytical since closed-form analytic solution cannot be obtained when computing the involved integrals in such a method. If analytic exact solution is not possible or too complex, we might then seek its approximation, either via the approximation (linearization) of nonlinear terms in this BVP (4), or via the approximation of the solution of this BVP.

Linearization of nonlinear terms in BVP

It is clear that the BVP (4) is nonlinear, and thus some extra assumptions can be imposed for the purpose of linearly approximating those nonlinear terms of (4) in Cartesian coordinates. Precisely, since

d✓ ds = d 2 y dx 2 [1 + ( dy dx ) 2 ] 3_2 (5)
then by assuming the Y-axis deflection of the beam a f L 10 , i.e. dy dx ˘0, we have:

d✓ ds = d 2 y dx 2 (6)
Rearranging Eq. 2 via the simplified curvature definition Eq. 6, we can finally arrive at Eq. 7 where BCM is originally derived from (see details in Section 3 of [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF]): [START_REF] Wu | Design and nonlinear modeling of a novel planar compliant parallelogram mechanism with general tensuralcompresural beams[END_REF] which is much easier to handle as it is now linear. Consequently, BCM and TBCM are both limited to a f L 10 due to this assumption.

EI d 2 y dx 2 = F y (a * x) + F x (b * y) + M o

Numerical methods to approximate the solution of BVP

Another technique to handle nonlinear BVPs is to approximate its solution via numerical methods. For BVP (4) (general beam-end loading) and its similar forms, many academic contributions have been completed. Ref. [START_REF] Banerjee | Large deflection of cantilever beams with geometric non-linearity: Analytical and numerical approaches[END_REF][START_REF] Sung | A nonlinear shooting method for two-point boundary value problems[END_REF] recommended shooting method (where the BVP is transformed into an initial value problem (IVP) under the concept of shooting method, and the IVP is solved using Runge-Kutta method). Ref. [START_REF] Dado | A new technique for large deflection analysis of non-prismatic cantilever beams[END_REF] studied a non-prismatic cantilever beam under general end loading along with uniformly distributed vertical and horizontal loading (the authors solved the ODE via introducing optimization concept to the minimize the equation residual until it's infinitely near 0 where the optimization process is functioned by Newton' method). Ref. [START_REF] Gao | Optimization algorithm-based approach for modelling large deflection of cantilever beam subjected to tip load[END_REF] recently proposed Particle-swarm-optimization (PSO)-based algorithm to handle the BVP. It is worth noting that Ref. [START_REF] Bailey | Nonlinear flexure coupling elements for precision control of multibody systems[END_REF] used finite-difference method implementing a three-stage Lobatto IIIa formula to solve the BVP, which has demonstrated the effectiveness of the proposed method from an applied mathematics point of view.

In the literature, there are also many simplified loading scenarios of beam-deflection problems solved by numerical methods. Ref. [START_REF] Mutyalarao | Large deflections of a cantilever beam under an inclined end load[END_REF] investigated a slender straight beam under an inclined load/follower load at the beam end (given a specified tip-angle (0) in ( [START_REF] Mutyalarao | Large deflections of a cantilever beam under an inclined end load[END_REF]. Fig. 1), and then the authors approximates the integral ( [START_REF] Mutyalarao | Large deflections of a cantilever beam under an inclined end load[END_REF]. Eq. 18) using Gauss-Legendre quadrature method to find the needed loading); Ref. [START_REF] Beléndez | Large and small deflections of a cantilever beam[END_REF] and Ref. [START_REF] Ghuku | A theoretical and experimental study on geometric nonlinearity of initially curved cantilever beams[END_REF] focused on straight beams and ICBs subjected to a vertical load at the beam end respectively; Ref. [START_REF] Kumar | Techniques based on genetic algorithms for large deflection analysis of beams[END_REF] proposed a genetic-algorithm-based shooting method for the case of an in-plane load and a transversal load exerted at the beam tip; although the authors of [START_REF] Bs Shvartsman | Large deflections of a cantilever beam subjected to a follower force[END_REF] studied the same case of a follower load subjected at the beam end, they artistically transformed the shooting-methodbased BVP into a direct IVP, which largely saves computational expenses; Ref. [START_REF] Lee | Large deflections of cantilever beams of non-linear elastic material under a combined loading[END_REF] utilized shooting method to solve a combined-loading case where a slender beam is subjected to a vertical concentrated load at the beam tip and a uniformly distributed load along its body; noticeably, Ref. [START_REF] González | Stiffness of a curved beam subjected to axial load and large displacements[END_REF] investigated the case of axially extensible ICBs under a horizontal load at the beam tip since the majority of studies regarding modeling beams assume axial inextensibility.

Contributions of the paper

As discussed above in Section 1.2, the large beam-deflection problems of a slender beam can be basically formulated as nonlinear BVPs, such as (4), and analytic (or semi-analytic) solutions have been proposed only for the uniform slender beam with simple beam-end loading. For more general loading scenarios (such as large deflection with varying cross sections, distributed bending force/moment, etc.), no analytic solution has been reported in the literature, and conventional method is to use chained algorithm to numerically get approximated solutions. Compared to those mentioned results, the main contributions of this paper are as follows:

1. Formulation: we formulate various planar and static beam-deflection problems as a general nonlinear BVP including varying cross sections, distributed body force, beam-end loading and so on. 2. Numerical method: we propose a modified collocation method, serving as a quite universal method to approximate the solution of the deduced nonlinear BVPs of beam-deflection problems. Besides, the hot topic of dealing with inflection points [START_REF] Zhang | A integral to the large deflection problems of thin beams in compliant mechanisms[END_REF] is considered and solved in this method. 3. Solving constitutive models: we apply the proposed method to study 8 loading scenarios (considering distributed loading, pressure, constant-curvature ICBs, varying-curvature ICBs, beams with varying cross sections and so on) that may happen in CMs. To the best of our knowledge, no existing beam models could cover all these scenarios. Moreover, the proposed collocation method has proved to be more efficient than FEA (using 1-D beam elements). 4. Proving feasibility in modeling CMs: we have proved the feasibility of the proposed method in conducting static synthesis of CMs by modeling four representative ones (such as compliant legs and compliant parallelograms).

Comprehensive modeling methodology and numerical solution

In this section, we first derive the general ODE for the large-deflection problem of slender beams that represent as many loading scenarios as possible in planar CMs, and then explain the proposed modified collocation method in detail.

Modeling of large-deflection of slender beams with general configuration

In this paper, we aim to extend the previous work, from the designing-CM-oriented point of view, by taking into account more loading scenarios on flexible slender beams that exist in CMs. Precisely, for a slender beam with uniformly varying cross sections depicted in Fig. 2, it is assumed to be subjected to general beam-end loading where F x and F y are the horizontal force and the vertical force respectively decomposed from an arbitrary force F o ; M o is an arbitrary moment exerted at the beam end as well; q x (s) and q y (s) denote the horizontal and vertical distributed loads along the beam body (a typical example is gravity where the sum of q x (s) and q y (s) determines the direction of the gravity, and the gravity could also be q x (s) or q y (s) itself alone); q n (s) denotes the distributed load whose direction is always normal to the beam axis (a typical example is pressure exerted on the side surface of the beam); M < (s) is the distributed bending moment along the beam body (for example, conducting polymer actuators can be a potential targeted spot [START_REF] Farajollahi | Characterization and dynamic charge dependent modeling of conducting polymer trilayer bending[END_REF] for this loading scenario); L denotes the beam length. With all the above defined loads, we then proceed to derive the BVP for the general loading of a general slender beam shown in Fig. 2. From Eq. (1), we can have:

F x F y F o M o (x(s), y(s), s) q y q x M (s) * q n (x(L), y(L), L) X Y
d✓ ds = M(s) EI(s)
where we consider the cross section varies along the beam axis (the second moment of inertia of the cross section

I(s) is a function of s). If the slender beam has an initial curvature 1 R(s)
where we consider the varying curvature by

defining 1 R(s)
as a function of s, we will arrive at:

d✓ ds = M(s) EI(s) + 1 R(s) (8)
To properly define the geometric property of cross sections and the beam length, we would like to clarify the exact beam axis we refer to. As is stated in [START_REF] Dennis | Advanced engineering mathematics[END_REF], the geometric property of cross section is described by the second moment of inertia of the cross section I(s), and I(s) is formulated based on the centroidal axis. Besides, the centroidal axis is an axis that passes through the centroid of the cross section [START_REF] Dennis | Advanced engineering mathematics[END_REF]. Therefore, the cross section varies along the centroidal axis. In this paper, the centroidal axis lies in a plane, which means only planar beams are considered in this work. Similarly, its beam length L refers to the length of the centroidal axis as well.

According to the definition of M(s) in Eq. ( 1), M(s) is the exerted moment at any arbitrary point on the beam, noted as (s, x(s), y(s)) due to all external loads on the beam. Therefore, M(s) can be formulated as:

M(s) =F y (x(L) * x(s)) + F x (y(L) * y(s)) + M o + M < (s) * L s q y (s)[x( ) * x(s)]d + L s q x (s)[y( ) * y(s)]d + L s q n (s)(x( ) * x(s)) cos ✓( )d + L s q n (s)(y( ) * y(s)) sin ✓( )d (9)
where (L, x(L), y(L)) are the coordinates of the beam end. Rearranging Eq. ( 8) and Eq. ( 9), we will end up with:

EI(s) d✓ ds =F y (x(L) * x(s)) + F x (y(L) * y(s)) + M o + M < (s) * L s q y (s)[x( ) * x(s)]d + L s q x (s)[y( ) * y(s)]d + L s q n (s)(x( ) * x(s)) cos ✓( )d + L s q n (s)(y( ) * y(s)) sin ✓( )d + EI(s) R(s) ( 10
)
where E denotes the Young's modulus of the chosen isotropic material; d✓ ds denotes the curvature of the beam after deformation. We would like to emphasize that ( 10) is nonlinear whose right-hand term contains 2 unknown independent constants x(L) and y(L). So, in order to remove these constant dependencies, by applying dx ds = cos ✓ and dy ds = sin ✓, we calculate the derivative of ( 10) with respect to s, which yields:

E dI ds d✓ ds + EI(s) d 2 ✓ ds 2 = * (F y cos ✓(s) + F x sin ✓(s)) + dM < (s) ds + E( dI ds (s) R(s) * I(s) dR ds (s) R(s) 2 ) + q y (s) cos ✓(s)(L * s) + dq y ds x(s)(L * s) + dq y ds s L r 0 cos ✓(⇠)d⇠dr + q x (s) sin ✓(s)(s * L) + dq x ds y(s)(s * L) + dq x ds L s r 0 sin(✓(⇠))d⇠dr + q n (s) cos ✓(s) s L cos ✓(⇠)d⇠ + dq n ds x(s) s L cos ✓(⇠)d⇠ + dq n ds L s r 0 cos ✓(⇠)d⇠ cos ✓(r)dr + q n (s) sin ✓(s) s L sin ✓(⇠)d⇠ + dq n ds y(s) s L sin ✓(⇠)d⇠ + dq n ds L s r 0 sin ✓(⇠)d⇠ sin ✓(r)dr (11)
subjected to the following two boundary conditions:

B.C. ✓(0) = 0 d✓ ds (L) = M o EI(L) + 1 R(L) (12)
In summary, for a slender beam under general loading, we can describe its large deflection via [START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF] and [START_REF] Sam Han | Design, simulation, and fabrication of a quadstable monolithic mechanism with x-and y-directional bistable curved beams[END_REF], which is a highly nonlinear BVP. Obviously, it is complex or even impossible to deduce the associated analytic solution, and therefore we will present a modified collocation method in the next subsection to approximate solution of the deduced BVPs.

Solution approximation via modified collocation method

As reported in the literature of applied mathematics, collocation method has already been a widely accepted option of solving two-point BVPs, and it essentially conducts collocation with several continuous polynomial functions to approximate the target unknown function under given boundary conditions [START_REF] Russell | A collocation method for boundary value problems[END_REF] [START_REF] Lawrence F Shampine | Solving boundary value problems for ordinary differential equations in matlab with bvp4c[END_REF]. To recall it, we would like to briefly explain how BVP ( 13) is handled using collocation method in a local manner.

We first arrange the BVP [START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF] into a more succinct one for easier demonstration when manipulating it:

D.E. d 2 ✓ ds 2 = f (s, d✓ ds , ✓...) B.C. ✓(0) = 0 d✓ ds (L) = M o EI(L) + 1 R(L) (13)
where f represents the nonlinear function in [START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF].

Generally speaking, by dividing the interval [0, L] into i + 1 points 0 = s 0 < s 1 < 5 < s i = L (named as collocation points), collocation method approximates the solution of the above BVP via linear combination of wellchosen basis functions by satisfying D.E. and B.C. at all collocation points. Precisely, denote by ⇥(s) the approximate solution of ( 13), and suppose that a cubic polynomial p k (s) is chosen to approximate the solution ✓ of ( 13

) on each sub-interval [s k*1 , s k ], i.e., ⇥(s) = ≥ i k=1 p k (s)
where

p k (s) = < b k0 + b k1 s + b k2 s 2 + b k3 s 3 , ≈s À ⌅ s k*1 , s k ⇧ 0, otherwise (14)
with b kj for 0 f j f 3 being the 4 unknown polynomial coefficients on the kth interval. In total, we have 4i unknown parameters to be identified. Clearly, ⇥(s) should satisfy the D.E. in ( 13) for all collocations points, and this yields the following 3i algebraic equations to be fulfilled:

d 2 ⇥ ds 2 (s k*1 ) = f (s k*1 , d⇥ ds (s k*1 ), ⇥(s k*1 )) d 2 ⇥ ds 2 ( Ñ s k k*1 ) = f ( Ñ s k k*1 , d⇥ ds ( Ñ s k k*1 ), ⇥( Ñ s k k*1 )) d 2 ⇥ ds 2 (s k ) = f (s k , d⇥ ds (s k ), ⇥(s k )) (15)
where

Ñ s k k*1 = s k*1 + s k 2 .
Note that the approximated solution ⇥(s) is a combination of a series of local cubic polynomials, therefore additional equation should be imposed to guarantee the continuity at all collocation points, i.e.,

p k*1 (s k*1 ) = p k (s k*1 ), ≈k À [2, i] ( 16
)
Hence we have the above i * 1 algebraic equations due to the continuity conditions. Together with the 2-dimensional B.C. in [START_REF] Pham | A constant-force bistable mechanism for force regulation and overload protection[END_REF], we finally obtain a set of algebraic equations with dimension 3i + i * 1 + 2 = 4i + 1, which will be used to compute those 4i unknowns parameter b kj defined in [START_REF] Brian T Edwars | A pseudo-rigid-body model for functionally binary pinnedpinned segments used in compliant mechanisms[END_REF].

So far, we have demonstrated the common collocation method performed in a local manner [START_REF] Lawrence F Shampine | Solving boundary value problems for ordinary differential equations in matlab with bvp4c[END_REF]. It can be easily noticed that only continuity at all collocation points are satisfied but the derivatives at these collocation points might not exist, which ends up with non-differentiable piecewise polylines for the solution ⇥(s). Obviously, this result contradicts the nature of beam deflection as ✓(s) is always continuous and differentiable. To handle such a problem, we propose to take one step further to ensure the differentiability to modify the collocation method, i.e. ⇥(s) À C 1 , by imposing the following derivative continuity conditions at all collocation points:

dp k*1 ds (s k*1 ) = dp k ds (s k*1 ), ≈k À [2, i] (17)
This added constraints of derivative continuity results in a differentiable global approximation of the solution of the studied BVP, bypassing the need of conducting curve fitting process afterwards.

It is clear that [START_REF] Qi | A novel continuum manipulator design using serially connected double-layer planar springs[END_REF] yields another i * 1 equations, ending up with altogether 4i + 1 + i * 1 = 5i algebraic equations available. In summary, the collocation method enables us to solve the 4i unknowns b kj based on the 4i + 1 equations from ( 13), ( 15) and ( 16), or 5i equations from ( 13), [START_REF] Trease | Design of adaptive and controllable compliant systems with embedded actuators and sensors[END_REF] [START_REF] Song | Development of small sized actuator with compliant mechanism for optical image stabilization[END_REF] and [START_REF] Qi | A novel continuum manipulator design using serially connected double-layer planar springs[END_REF] which can be numerically solved via Newton-Raphson method.

It is worth noting that numerical tools would always need an initial guess to start its iterative computation which might play an important role in the convergence and local minima [START_REF] Hull | Comparing numerical methods for ordinary differential equations[END_REF] [START_REF] Tsai | An introduction of numerical methods[END_REF]. In our context, we can have two different ways to set the initial values of b kj :

1. Direct way: We randomly initialize the values of b kj , which yields a random shape of ⇥(s). 2. Indirect way: We choose a globally initial guess for ⇥(s), which yields the corresponding b kj for every cubic polynomial p k (s).

The direct way is suitable for the scenario where the problem tends to have only one single solution, i.e., it contains only global minima. 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, and it is suitable to treat local minima case. In this paper, the direct way of initializing guesses is used for all loading scenarios.

Once we obtain all suitable b kj for every cubic polynomial p k (s), then we can arrive at a proper approximation ⇥(s) to ✓(s). Then, the deformed beam shape (x(s) and y(s)) can be solved using:

x(s) = s 0 cos(⇥(⇠))d⇠; y(s) = s 0 sin(⇥(⇠))d⇠; for s À [0, L].
Here, we would like to analyze the proposed modified collocation method further in terms of following topics: choosing the number of collocation elements i, solving inflection-point-involved beam-deflection problems and the difference compared to FEA respectively.

1. we would like to briefly discuss the influence of i in terms of computational time and accuracy.

Computational time: as we have stated above, the total number of equations to be solved is 5i, and the number of unknowns is 4i. Let us formulate the system of equations as F (x) = 0 À R 5iù1 with x À R 4iù1 . Then, the numerical solution can be obtained via Newton-Raphson method:

x m = x m*1 * 4 )F (x) )x 5 *1 x=x m*1 F (x m*1 )
where x m represents the mth iterative solution with x 0 being a given initial guess condition. After the convergence (in practice the convergence criteria is set as x m+1 * x m  f ✏ where ✏ is a predefined threshold), we can finally obtain the solution x to F (x) = 0. The dimension of the Jacobian matrix )F (x) )x is R 5iù4i . Obviously, the major computation of this algorithm lies in calculating the pseudo inverse of )F (x) )x . Therefore, the higher value of i of the studied matrix, which leads to higher dimension of the matrix, would result in more computational time for computing the inverse of this matrix.

Accuracy: mathematically, the larger value of i is chosen, the more cubic polynomials are used to approximate the solution of each subpart, which logically can reach the solution of the model more accurately. However, this is not essential. Collocation method is just one of the numerical methods to solve BVPs so it always relies on a global governing ODE (which in this case is Eq. ( 1)). Therefore, it is possible to have a precise approximated solution of such an ODE by even using one relatively low order of polynomial. In this paper, according to the testing results (see Table 1 to Table 8), i = 2 is accurate enough and used for all loading scenarios. It has been shown that a larger i may not be necessary to get high accuracy, which could inevitably introduce more unnecessary computational cost.

2. In beam-deflection problems, dealing with inflection points are always a hot area in this field [START_REF] Zhang | A integral to the large deflection problems of thin beams in compliant mechanisms[END_REF] [START_REF] Kimball | Modeling of flexural beams subjected to arbitrary end loads[END_REF]. We would like to emphasize that this proposed modified collocation can efficiently handle the existence of inflection points along the unformed beam. According to [START_REF] Kimball | Modeling of flexural beams subjected to arbitrary end loads[END_REF] and [START_REF] Zhang | A integral to the large deflection problems of thin beams in compliant mechanisms[END_REF], the essence of inflection points in a beam is about the sign change of d✓ ds . If no inflection point exists, the slope is either monotonically increasing or decreasing, d✓ ds must always be positive or negative, respectively. If an inflection point exists, the sign of d✓ ds changes once along the range of independent variable s À [0, L]. Mathematically speaking, if d✓ ds (s < ) = 0 and

d 2 ✓ ds 2 (s < ) ë 0, then s = s < is an inflection point.
In [START_REF] Zhang | A integral to the large deflection problems of thin beams in compliant mechanisms[END_REF], inflection points have been intensively explained.

Given a mathematical BVP that describes large deflection of slender beams, this problem becomes essential when we use elliptical integrals to obtain analytical solution of this BVP. The main reason is that such an analytical solution depends on the sign of d✓ ds . However, our proposed method is numerical by using polynomials to numerically approximate the solution of the BVP. Therefore, the inflection points will be implicitly treated by the proposed numerical method to seek the solution of Eqs. (15)(16)(17).

3. The discretization idea of the proposed collocation method looks similar to that of classic FEA but they are essentially two different methods to solve beam-deflection problems. FEA does not have a global governing equation for beam deflection problems, thus it always needs to discretize the studied object into fine mesh elements, and apply either Euler-Lagrange equation or virtual work principal to these mesh elements to obtain the dynamical or static model of the studied beam. However, collocation method is a commonly-used numerical methods to solve BVPs, and therefore it always relies on a global governing ODE (Eq.( 1)). As a result, it is likely that the solution of the BVP can be approximated by a low-order polynomial with precision).

Application to model large-deflection of various slender beams in CMs

In this section, we aim to apply the proposed methodology to model 8 beam-deflection scenarios that exist in CMs, and then solve it using the proposed modified collocation method. We would like to first present the different BVP formulations rearranged from [START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF] and [START_REF] Sam Han | Design, simulation, and fabrication of a quadstable monolithic mechanism with x-and y-directional bistable curved beams[END_REF], which are used to model the mentioned 8 loading scenarios. We then provide the numerical results afterwards to prove the effectiveness of the proposed modeling methodology.

Formulations of 8 different scenarios in the proposed framework

From Section 3.1.1 to Section 3.1.7, beams with rectangular cross section are studied. As graphically displayed in Fig. 3a, its initial cross-section parameters are set as w(0) = 0.01 m and h(0) = 0.004 m respectively, and the length of the beam L = 0.2 m The cross-section parameters at any random point along the beam axis are defined as functions of s: w(s) and h(s). In Section 3.1.8, we focus on a round-cross-section beam whose geometric definition is shown in Fig. 3b. Similarly, the diameter of the initial cross section is fixed to d(0) = 0.005 m, and the diameter at any random point along the beam axis is defined as a function of s: d(s). It should be noted that L, w(s), h(s) and d(s) are all defined based on the beam axis and the beam axis could be straight or pre-curved due to the initial curvature R(s) *1 of beams as provided in Eq. 10. The material chosen is structural steel with its Young's modulus E = 200 ù 10 9 Pa and Poisson's ratio ⌫ = 0.3. As the shear strain is not considered in Euler-Bernoulli beam model, ⌫ is logically not involved in our numerical algorithms. In the following, we will start from the most basic one to more complicated loading scenarios that often happen in CMs as depicted in Fig. 4.

General loading of a slender straight beam

Modeling the pure beam-end loading of a slender straight beam as shown in Fig. 4a yields the following simplifications for the BVP [START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF] as: R(s) *1 = 0; M < (s) = 0; q x (s) = 0; q y (s) = 0; q n (s) = 0; h(s) = h(0); w(s) = w(0); I = h(0) 3 w(0) 12 ;

w (0) which is the most commonly studied scenario BVP (4) in CM research community, such as flexure-based highprecision CMs [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF]. The related numerical results, by applying the proposed modified collocation method, are presented in Fig. 5 and Table 1.

h (0) s L Fixed base w (s) h (s)

General loading of a slender straight beam with distributed moment along its beam axis

Modeling general loading of a slender straight beam (Fig. 4a) with distributed moment along its beam axis yields the following simplifications for the BVP [START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF] as: R(s) *1 = 0; M < (s) = 40 * 100s; q x (s) = 0; q y (s) = 0; q n (s) = 0; h(s) = h(0); w(s) = w(0); I = h(0) 3 w(0) 12 ;

Then, we will have the simplified version of the BVP (11): This scenario often happens in modeling compliant beam-like conducting polymer actuators [START_REF] Farajollahi | Characterization and dynamic charge dependent modeling of conducting polymer trilayer bending[END_REF] which is a novel multi-disciplinary study in the field of CMs. The related numerical results, by applying the proposed modified collocation method, are presented in Fig. 6 and Table 2.

Scenario 2: D.E. d 2 ✓ ds 2 = * F y EI (cos ✓(s) + F x F y sin ✓(s)) * dM < (s) ds EI B.C. ✓(0) = 0 d✓ ds (L) = M o + M < (L) EI (19)

General loading of a slender straight beam considering uniformly distributed loading

In some circumstances, designing CMs needs to consider the uniformly distributed loading on their bodies, among which one of the most representative cases is considering gravity of CMs [START_REF] Dennis | Advanced engineering mathematics[END_REF] [START_REF] Radaelli | Gravity balanced compliant shell mechanisms[END_REF]. Here, we aim to model the uniformly distributed loading on a slender straight beam (see the exact geometry in Fig. 4a) along with the general beam-end loading, which yields the following simplifications from the BVP (11): R(s) *1 = 0; M < (s) = 0; q x (s) = q x ; q y (s) = q y ; q n (s) = 0; h(s) = h(0); w(s) = w(0); I = h(0) 3 w(0) 12 ;

Therefore, we can come to:

Scenario 3: D.E. d 2 ✓ ds 2 = * F y EI (cos ✓(s) + F x F y sin ✓(s)) + q y cos ✓(s)(L * s) EI + q x sin ✓(s)(s * L) EI B.C. ✓(0) = 0 d✓ ds (L) = M o EI (20)
The related numerical results, via the proposed modified collocation method, are presented in Fig. 7 and Table 3.

3.1.4. General loading of a slender straight beam along with pressure on its body Fluid (gas or liquid)-driven CMs have not been intensively explored according to the literature but some representative studies on this topic have been completed [START_REF] Zentner | On the classification of compliant mechanisms[END_REF][64] [START_REF] Linß | Pneumatically driven compliant structures based on the multi-arc principle for the use in adaptive support devices[END_REF]. In this section, a slender straight beam subjected to pressure along its body and beam-end loading is modeled, which leads to the following conditions for the BVP (11): R(s) *1 ; M < (s) = 0; q x (s) = 0; q y (s) = 0; q n (s) = q n ; h(s) = h(0); w(s) = w(0); I = h(0) 3 w(0) 12 ;

Logically, we can arrive at:

Scenario 4: D.E. d 2 ✓ ds 2 = * F y EI (cos ✓(s) + F x F y sin ✓(s)) + q n cos ✓(s) î s L cos ✓(⇠)d⇠ EI + q n sin ✓(s) î s L sin(⇠)d⇠ EI B.C. ✓(0) = 0 d✓ ds (L) = M o EI ( 21
)
The related numerical results, via the proposed modified collocation method, are presented in Fig. 8 and Table 4.

General loading of a slender straight beam of uniformly varying cross sections

Flexible beams of varying cross sections present more complex but more diverse mechanical responses, which enriches the selection of elementary members and more creative design possibilities for CMs [START_REF] Lobontiu | Design of circular cross-section corner-filleted flexure hinges for three-dimensional compliant mechanisms[END_REF]. In this section, we aim to model a slender straight beam of uniformly varying cross sections (see the exact beam geometry in Fig. 4b) under general beam-end loading, which results in the following simplifications for the BVP (11): 4 12 ;

R(s) *1 = 0; M < (s) = 0; q x (s) = 0; q y (s) = 0; q n (s) = 0; h(s) = h(0)(1 * 4s); w(s) = w(0)(1 * 4s); I = h(0) 3 w(0)(1 * 4s)
We will have the following based on the mentioned simplifications: The related numerical results, via the proposed modified collocation method, are presented in Fig. 9 and Table 5.

General loading of a slender initially curved beam with a constant curvature

ICBs hugely expand the design scope of CMs they have a large deflection range along with a relatively small strain range, which is perfectly suitable for large-scale movement mechanisms [START_REF] Qiu | A curved-beam bistable mechanism[END_REF] [START_REF] Verotti | Effect of initial curvature in uniform flexures on position accuracy[END_REF]. Here, we model an ICB with a constant radius of curvature R = 0.1 m (see the exact geometry in Fig. 4c). The simplifications made here for the BVP [START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF] are as follows:

M < (s) = 0; q x (s) = 0; q y (s) = 0; q n (s) = 0; h(s) = h(0); w(s) = w(0); R(s) = R; I = h(0) 3 w(0) 12 ;
Therefore, the BVP we actually deal with for this scenario is as follows:

Scenario 6: D.E. d 2 ✓ ds 2 = * F y EI (cos ✓(s) + F x F y sin ✓(s)) B.C. ✓(0) = 0 d✓ ds (L) = M o EI + 1 R (23)
The related numerical results, via the proposed modified collocation method, are presented in Fig. 10 and Table 6.

3.1.7. General loading of a slender initially curved beam with a varying curvature

We have discussed ICBs with a constant curvature in the last section. Here in this section, we move onto a slightly different loading scenario: a slender ICB with a varying curvature (see Fig. 4d) under general beam-end loading [START_REF] Wu | Efficient spatial compliance analysis of general initially curved beams for mechanism synthesis and optimization[END_REF], which leads to the following assumptions for the BVP (11): M < (s) = 0; q x (s) = 0; q y (s) = 0; q n (s) = 0; h(s) = h(0); w(s) = w(0); R(s) = 1 0.1 + 1000s 2 ; I = h(0) 3 w(0) 12 ;

Then, we obtain the simplified version of the BVP (11):

Scenario 7: D.E. d 2 ✓ ds 2 = * F y EI (cos ✓(s) + F x F y sin ✓(s)) * dR ds R(s) 2 B.C. ✓(0) = 0 d✓ ds (L) = M o EI + 1 R(L) (24)
The related numerical results, via the proposed modified collocation method, are presented in Fig. 11 and Table 7.

Combined loading of a slender initially curved beam of round cross section and varying curvature

In this section, we aim to model a relative complex scenario that may happen to CM design process to verify both the accuracy of the BVP (11) that describes the complex loading scenario and the feasibility of the method that's used to solve it. Also, this is to show that slender beams with different cross sections can fit in the methodology. Supposing a slender ICB of round cross section and varying curvature (see Fig. 4e) is subjected to general beam-end loading and uniformly distributed loading along its body, here are the conditions assumed for the BVP [START_REF] Gerson | Large displacement low voltage multistable micro actuator[END_REF] regarding this scenario: M < (s) = 0; q x (s) = q x ; q y (s) = q y ; q n (s) = 0; d(s) = d(0); I = ⇡d(0) 4 64 ; R(s) = 1 0.1 + 1000s 2 ;

which leads to the following:

Scenario 8: D.E. d 2 ✓ ds 2 = * F y cos ✓(s) + F x sin ✓(s) EI * dR ds R(s) 2 + q y (s) cos ✓(s)(L * s) + q x (s) sin ✓(s)(s * L) EI B.C. ✓(0) = 0 d✓ ds (L) = M o EI(L) + 1 R(L) (25)
The related numerical results, via the proposed modified collocation method, are presented in Fig. 12 and Table 8.

Discussion on numerical results

In this section, we present and discuss the numerical results of the 8 loading scenarios (Section 3.1.1 to Section 3.1.8) via applying the modified collocation method. The exact deformed beam shapes and the beam-end coordinates were modeled both compared and verified by FEA results with geometric nonlinearity activated.

In FEA, we used beam elements (Euler-Bernoulli beam), and the mesh information of each loading scenario is presented below: 1 to 8, the graphical results of deformed beam shapes and the numerical results of beam-end coordinates are presented respectively. In each loading scenario, the results obtained are compared and verified by reliable FEA results in terms of beam shapes (Figs. 5 to 12) and beam-end coordinates (Tables 1 to 8) where the high accuracy and efficiency of the proposed method have been proved. In Tables 1 to 8, MCM denotes the modified collocation method; ER denotes the error of MCM compared to FEA (ER = FEA*MCM FEA ); CE denotes the computational time.

Application to model representative CMs

As noted in the current literature, Awtar first systematically studied compliant parallelograms [START_REF] Awtar | A generalized constraint model for two-dimensional beam flexures: Nonlinear strain energy formulation[END_REF][61] [START_REF] Awtar | Closed-form nonlinear analysis of beam-based flexure modules[END_REF]. In Awtar's pioneering work of modeling these compliant parallelograms, all flexible beams are limited to relatively small deflection range for highly accurate positioning CMs, which has proved to be one of most important contributions in recent years. With the increasing interest of large-deflection-based CMs in the research community [START_REF] Ling | Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey[END_REF], we here proceed to work on these CMs using the modeling methodology for elementary flexible members proposed in this paper to further prove its feasibility of synthesizing CMs. We would like first to start with the lumped compliant legs (straight-beam-based and ICB-based) to the typical compliant parallelograms (also straight-beam-based and ICB-based respectively). Here, we present modeling the mentioned CMs in terms of much larger deflection range following the modeling principles stated in Awtar's PhD work, which is a typical formulation in mechanics [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF]: a) Constitutive relationships b) Force equilibrium or force compatibility relationships c) Geometric equilibrium or geometric compatibility relationships

The above three relationships determine the essence of modeling deformable bodies regardless of its linear or nonlinear properties. Finally, the modeling results are compared and verified by FEA results with corresponding errors presented accordingly. Note that solid mesh element (under the framework of solid mechanics) are used in this FEA since we want to consider all types of deformation in beams, such as the axial extension of beams, to verify if our proposed modeling methodology is feasible in mechanism synthesis.

Modeling lumped compliant legs

In this section, by applying the proposed modeling methodology, we aim to model two types of compliant legs: straight-beam-based ones and ICB-based ones.

Straight-beam-based lumped compliant leg

To model the compliant leg as shown in Fig. 14a, three types of equations need to be developed which describe constitutive relationships, force equilibrium and geometric compatibility relationships respectively as demonstrated in Fig. 13a. As displayed Fig. 14a, the compliant leg can be characterized by the following geometric parameters: w, h and L i (i = 1, 2) denote depth, width and length of the beam, and W denotes the distance between the center points of the two beams in one leg.

From the perspective of the global coordinate system X-Y-Z, F x , F y and M z are the forces exerted at the reference point R of the leg, and x, y and ↵ are the corresponding translational and angular displacements. Similarly, f x1 , f y1 , m z1 , x 1 , y 1 , ✓ 1 are defined in its local coordinate system X 1 -Y 1 -Z 1 and f x2 , f y2 , m z2 , x 2 , y 2 , ✓ 2 are defined in its local coordinate system X 2 -Y 2 -Z 2 respectively as displayed in Fig. 13a. Logically, we can have the following relationships via synthesizing the above defined variables:

1) Constitutive relationships

If we consider Scenario 1 in Section 3.1.1 as an equation of f xi , f yi , m zi , x i , y i , ✓ i , E, w, h, L i (see Fig. 13a and Fig. 14a, i = 1, 2) , we then can formulate the constitutive relationships as follows:

f sb (f x1 , f y1 , m z1 , x 1 , y 1 , ✓ 1 , E, w, h, L 1 ) = 0; f sb (f x2 , f y2 , m z2 , x 2 , y 2 , ✓ 2 , , E, w, h, L 2 ) = 0; (26)
where f sb = 0 denotes the constitutive equation of straight beams.

2) Force equilibrium The deformed beam shape by MCM The deformed beam shape by FEA 3) Geometric compatibility relationships

f x2 cos ✓ 1 * f y2 sin ✓ 1 = F x ; f y2 cos ✓ 1 + f x2 sin ✓ 1 = F y ; m z1 = m z2 * f x2 y 2 + f y2 (W * 0.5L 1 + 0.5L 2 * x 2 ); M z = m z2 ; f y1 = F y ; f x1 = F x ; (27)
x =0.5L 1 + 0.5L 2 * (L 1 * x 1 + W cos ✓ 1 * 0.5L 1 cos ✓ 1 + 0.5L 2 cos ✓ 1 * x 2 cos ✓ 1 * y 2 sin ✓ 1 ) + W ; y = y 1 + (W * 0.5L 1 * 0.5L 2 ) sin ✓ 1 + L 2 sin ✓ 1 * x 2 sin ✓ 1 + y 2 cos ✓ 1 ; ↵ =✓ 1 + ✓ 2 ; (28)
Therefore, given the geometric and material parameters E, w, h, L 1 , L 2 , W and any three of F x , F y , M z , x, y and ↵, the rest three can be obtained by solving Eqs. 26 to 28, which leads to the following:

f sbl (F x , F y , M z , x, y, ↵, E, w, h, L 1 , L 2 , W ) = 0 (29)
Logically, Eq. 29 is a system of nonlinear equations (Eqs. 26 to 28) formulated in a matrix form. Here, f sbl = 0 explains the force-displacement relationship of the leg's reference point R (see Fig. 13a). Therefore, we can use Newton's method or any other similar solvers to handle it.

It is worth noting that we can also consider f sbl = 0 is the constitutive equation of the leg when synthesizing more complex CMs instead of painstakingly developing the model from the most basic beam constitutive equations like the one in [START_REF] Wu | Design and nonlinear modeling of a novel planar compliant parallelogram mechanism with general tensuralcompresural beams[END_REF].

Here are the geometric and material parameters chosen for the compliant leg: E = 69 ù 10 10 Pa; w = 0.01 m; h = 0.001 m;

L 1 = L 2 = 0.1 m; W = 0.2 m; (30)
The loading conditions are:

F x = *1 N; F y = 1, 2...10 N; M z = 0 N.m;
Here, the relationships of F y against x, F y against y and F y against ↵ are explored where the results (AM denotes the results from the proposed solution and FEA logically denotes the ones from finite element analysis) along with based leg is also going to be developed upon constitutive relationships, force equilibrium and geometric compatibility relationships as shown in Fig 13b . Note that semi-circular ICBs are chosen for idea demonstration (see Fig. 14b) where the compliant leg can be characterized by the following geometric parameters: w, h and D i (i = 1, 2) denote depth, width and diameter of the beam curve, and W denotes the distance between the center points of the two beams in one leg.

As demonstrated in Fig 13b ,F x , F y and M z are the forces exerted at the reference point R of the leg, and x, y and ↵ are the corresponding translational and angular displacements, which are defined in the global coordinate system X-Y-Z. Similarly, f x1 , f y1 , m z1 , x 1 , y 1 , ✓ 1 are defined in its local coordinate system X 1 -Y 14b, i = 1, 2) , we then can arrive at the following constitutive relationships:

f icb (f x1 , f y1 , m z1 , x 1 , y 1 , ✓ 1 , E, w, h, D 1 ) = 0; f icb (f x2 , f y2 , m z2 , x 2 , y 2 , ✓ 2 , , E, w, h, D 2 ) = 0; (31)
where f icb = 0 denotes the constitutive equation of ICBs where in this case semi-circular ICBs are referred (see Fig. 14b).

2) Force equilibrium

f x2 cos ✓ 1 * f y2 sin ✓ 1 = F x ; f y2 cos ✓ 1 + f x2 sin ✓ 1 = F y ; m z1 = m z2 * f x2 y 2 + f y2 (W * 0.5L 1 + 0.5L 2 * x 2 ); M z = m z2 ; f y1 = F y ; f x1 = F x ; (32)
3) Geometric compatibility relationships

x =0.5L 1 + 0.5L 2 * (L 1 * x 1 + W cos ✓ 1 * 0.5L 1 cos ✓ 1 + 0.5L 2 cos ✓ 1 * x 2 cos ✓ 1 * y 2 sin ✓ 1 ) + W ; y = y 1 + (W * 0.5L 1 * 0.5L 2 ) sin ✓ 1 + L 2 sin ✓ 1 * x 2 sin ✓ 1 + y 2 cos ✓ 1 ; ↵ =✓ 1 + ✓ 2 ; (33)
Logically, given the geometric and material parameters E, w, h, D 1 , D 2 , W and any three of F x , F y , M z , x, y and ↵, the rest three can be obtained by solving Eqs. 31 to 33, which leads to the following: Essentially, Eq. 34 is a system of nonlinear equations (Eqs. 31 to 33) formulated in a matrix form, and f icbl = 0 explains the force-displacement relationship of the leg's reference point R (see Fig. 13b). Obviously, many math tools can handle it in an elegant and efficient manner, such as Newton's method. Note that f icbl can be considered as a constitutive equation of the ICB-based lumped compliant leg when modeling a more complex CM.

f icbl (F x , F y , M z , x, y, ↵, E, w, h, D 1 , D 2 , W ) = 0 (34) F x F y R R * Δx Δy α f x1 θ 1 f y1 δx 1 δy 1 f x2 m z2 θ 2 f y2 δx 2 δy 2 θ 1 X 1 Y 1 Z 1 X 2 Y 2 Z 2 X Y Z m z1 M z (a) Straight-beam-based lumped compliant leg X Y Z R Δx Δy F y α F x R * f x1 f y1 θ 1 δx 1 δy 1 X 1 Y 1 Z 1 X 2 Y 2 Z 2 f x2 m z2 f y2 θ 2
Here are the geometric and material parameters chosen for the compliant leg: The loading conditions are:

F x = *1 N; F y = 1, 2...10 N; M z = 0.5 N.m;
The relationships of F y against x, F y against y and F y against ↵ are plotted in Fig. 16a to Fig. 16c where AM results are compared and also verified by FEA results along with the acceptable errors. Besides, the deformed shapes of the compliant leg evaluated in FEA are depicted in Fig. 16d to graphically present its deformation.

Modeling compliant parallelograms

In this section, we aim to model two types of compliant parallelograms that integrate straight beams and ICBs respectively based on the work of Section 4.1.

Straight-beam-based compliant parallelogram

The first studied parallelogram is based on two straight-beam-based lumped compliant legs as shown in Fig. 17a (where the parameter T is introduced to characterize the distance between the two compliant legs). To model it, its constitutive relationships, force equilibrium and geometric compatibility relationships need to be considered.

As depicted in Fig 17a,F X , F Y and M Z denote the forces exerted at the reference point R of the parallelogram, and X, Y and denote the X-axis translational displacement, Y-axis translational displacement and Z-axis angular displacement. Similarly, F xi , F yi and M zi are exerted at the reference points A and B of the two compliant legs respectively. x i , y i and ↵ i denote the X-axis translational displacement, Y-axis translational displacement and Z-axis angular displacement (i = 1 or 2). Then, the model can be developed through the following:

1) Constitutive relationships

Here, Eq. 29 is considered as the elementary equation of constitutive relationships (see Fig. 14a and Fig. 17a):

f sbl (F x1 , F y1 , M z1 , x1, y1, ↵ 1 , E, w, h, L 11 , L 21 , W 1 ) = 0; f sbl (F x2 , F y2 , M z2 , x2, y2, ↵ 2 , E, w, h, L 12 , L 22 , W 2 ) = 0; (36)
where L ji denotes the length of jth straight beam in the ith compliant leg, and W i stands for the distance between the built-in two straight beams in the ith compliant leg (i = 1, 2 and j = 1, 2). The relationships of F Y against X, F Y against Y and F Y against are plotted in Fig. 19a to Fig. 19c where AM results are compared and also verified by FEA results along with the corresponding errors (the maximum is around 7.5%). Fig. 19d demonstrates the deformed shapes of the compliant parallelogram evaluated in FEA.

B A R T X Y Z A B * ΔX F x1 F y1 M z1 A B * * F x2 F y2 R F X F Y A * B * ΔY M Z M z2 M z1 M z2 F x1 F y1 F x2 F y2 Δy 1 Δy 2 α 1 α 2 Δx 1 Δx 2 (a) Straight-beam-based lumped compliant parallelogram A X Y Z * F x1 F y1 M z1 A B * * F x2 F y2 R F X F Y A * ΔY M Z M z2 M z1 F x1 F y1 Δy 1 α 1 Δx 1 R B B * M z2 F x2 F y2 Δy 2 α 2 Δx 2 A B ΔX T (b) ICB-based lumped compliant parallelogram

Error analysis

In this section, we aim to analyze error sources that exist in the final synthesis of the 4 studied CMs (compliant legs and compliant parallelograms) in Section 4.1 and 4.2. As is shown in Fig. 15, Fig. 16, Fig. 18 and Fig. 19, the errors are around 2% to 8%. However, in Table 1 to 8, most of the errors are within 1% or even smaller. This is logical since the errors that exist in final synthesis are the summation of the error of each built-in flexible member in the studied CMs. Here, we focus on the 2 major error sources of modeling each flexible members: model errors and numerical errors.

Model errors

The first type of errors comes from the inaccurate constitutive model, which in this case is due to the assumptions made in Euler Bernoulli beam equation (Eq. 1). This constitutive model ignores the shear deformation on cross sections of the beam and assumes that the studied beam is inextensible. To solve the above 2 mentioned problems, our future work will focus on Timoshenko beam theory which includes the shear deformation of beam cross sections, and Eq. 1 also needs to be modified to consider the axial extension or compression.

Numerical errors

The second type of errors is called numerical errors, which is due to the essence of numerical methods. As we have stated, the modified collocation method is essentially a numerical method to approximate the solution to the studied BVPs (which in this paper refer to Eq. 18 to Eq. 25). Therefore, this proposed method would also result in some numerical errors. However, we can always increase the number of collocation points and the order of the pre-set polynomials defined in the modified collocation method and decrease the predefined threshold of Newton Raphson method's convergence criteria to reduce numerical errors.

Conclusions

Compliant Mechanisms display many desired properties for mechanical applications depending on elastic deformation of the involved compliant members. This paper introduces a comprehensive methodology to model general slender beams undergoing planar large deflection via directly solving the BVPs of an ODE that is derived from the most classic Euler-Bernoulli beam theory. 8 scenarios regarding a wide variety of loading conditions and beam geometry have been studied, and the synthesis of 4 compliant mechanisms has been conducted as well, which have been both verified by FEA. Our future work will include conducting CM synthesis using the proposed solutions and proposing a universal methodology for CM design from a system-optimization point of view.

Figure 1 :

 1 Figure 1: General loading condition of a general slender straight beam

Figure 2 :

 2 Figure 2: General loading condition of a general slender beam

Figure 3 :

 3 Schematic geometry definition of general slender beams Then, we will have the simplified version of the BVP (11): Scenario 1: D.E.

(a) Scenarios 1 to 4 : 4 :

 44 straight beam with constant rectangular cross sections (b) Scenario 5: straight beam with uniformly varying rectangular cross sections (c) Scenario 6: ICB with constant curvature and constant rectangular cross sections (d) Scenario 7: ICB with uniformly varying curvature and constant rectangular cross sections (e) Scenario 8: ICB with uniformly varying curvature and constant round cross sections Figure The exact beam geometries of the studied slender beams in 8 loading scenarios

Scenario 5 :

 5 F y cos ✓(s) + F x sin ✓(s))

Scenario 1 :

 1 15 edge elements; Scenario 2: 20 edge elements; Scenario 3: 15 edge elements; Scenario 4: 15 edge elements; Scenario 5: 459 mesh elements; Scenario 6: 21 edge elements; Scenario 7: 22 edge elements; Scenario 8: 27 edge elements; From Figs. 5 to 12 and Tables

Figure 5 :

 5 Figure 5: Scenario 1 results of beam shape

Figure 6 :

 6 Figure 6: Scenario 2 results of beam shape

Figure 7 :

 7 Figure 7: Scenario 3 results of beam shape

Figure 8 :

 8 Figure 8: Scenario 4 results of beam shape

Table 5 :Figure 10 :

 510 Figure 10: Scenario 6 results of beam shape

  ICB-based lumped compliant legFigure 13: Constitutive, force equilibrium and compatibility relationships of compliant legs

E = 69 ù 1 WFigure 14 :

 69114 Figure 14: Lumped compliant legs

Figure 17 :

 17 Constitutive, force equilibrium and compatibility relationships of compliant parallelograms In Fig 17b, F X , F Y and M Z denote the forces exerted at the reference point R of the parallelogram, and X, Y and denote the X-axis translational displacement, Y-axis translational displacement and Z-axis angular displacement respectively. Similarly, F xi , F yi and M zi are exerted at the reference points A and B of the two compliant legs respectively. x i , y i and ↵ i denote the X-axis translational displacement, Y-axis translational displacement and Z-axis angular displacement (i = 1 or 2). Then, the model can be developed through the following equations: 1) Constitutive relationships

Figure 18 : 41 ) 3 ) 19 :

 1841319 Figure 18: Results of the straight-beam-based compliant parallelogram

Table 1 :

 1 Scenario 1 results of beam-end coordinates

			Case	Case 1 Case 2	Case 3	Case 4
	Loading	F x (N) F y (N) M o (N.m)	900 300 40	225 75 10	*225 *75 *10	*900 *300 *40
	MCM i g 2		x(L) (m) y(L) (m) ✓(L) (rad)	0.03183 0.1908 0.1588 0.05357 *0.02754 *0.05968 0.1976 0.1889 2.1804 0.4636 *0.2504 *0.6067
	FEA		x(L) (m) y(L) (m) ✓(L) (rad)	0.03188 0.1908 0.1589 0.05353 *0.02755 *0.05973 0.1976 0.1880 2.1813 0.4632 *0.2505 *0.6070
	ER		x(L) y(L) ✓(L)	0 .16% 0 .06% 0 .04%	0.00% 0.07% 0.09%	0.00% 0.04% 0.04%	0.49% 0.08% 0.05%
	CE (s)	MCM (i = 2) FEA	0.31 14	0.26 4	0.27 3	0.60 6

Table 2 :

 2 Scenario 2 results of beam-end coordinates

			Case	Case 1 Case 2	Case 3	Case 4
	Loading i g 2	F x (N) F y (N) M o (N.m)	900 300 40	225 75 10	*225 *75 *10	*900 *300 *40
	MCM		x(L) (m) y(L) (m) ✓(L) (rad)	0.03183 0.1908 0.1588 0.05357 *0.02754 *0.05968 0.1976 0.1889 2.1804 0.4636 *0.2504 *0.6067
	FEA		x(L) (m) y(L) (m) ✓(L) (rad)	0.03188 0.1908 0.1589 0.05353 *0.02755 *0.05973 0.1976 0.1880 2.1813 0.4632 *0.2505 *0.6070
	ER		x(L) y(L) ✓(L)	0 .16% 0 .06% 0 .04%	0.00% 0.07% 0.09%	0.00% 0.04% 0.04%	0.48% 0.08% 0.05%
	CE (s)	MCM (i = 2) FEA	0.31 14	0.26 4	0.27 3	0.60 6

Table 3 :

 3 Scenario 3 results of beam-end coordinates

Table 4 :

 4 Scenario 4 results of beam-end coordinates

	0.3415 *0.5643

Table 6 :

 6 Scenario 6 results of beam-end coordinates

	Y (m)	0.04 0.06 0.08 0.12 0.14 0.1	Case 1		The deformed beam shape by MCM Case 2 The deformed beam shape by FEA Undeformed shape Case 3	Loading MCM i g 2 FEA	Case F x (N) F y (N) M o (N.m) x(L) (m) y(L) (m) ✓(L) (rad) x(L) (m) y(L) (m) ✓(L) (rad)	Case 1 Case 2 Case 3 Case 4 900 225 *225 *900 300 75 *75 *300 40 10 *10 *40 0.01751 0.08919 0.1541 0.1827 0.1038 0.1011 0.05503 *0.01285 4.0390 3.1592 2.2097 1.3040 0.01758 0.08923 0.1541 0.1826 0.1038 0.1012 0.05503 *0.01281 4.0384 3.1586 2.2105 1.3069
		0 0.02					Case 4	ER	x(L) y(L) ✓(L)	0 .40% 0 .00% 0 .01%	0.16% 0% 0.06%	0.40% 0.00% 0.04%	0.05% 0.31% 0.22%
		-0.02						CE (s)	MCM (i = 2) 0.87 FEA 9	0.23 4	0.45 4	0.53 16
		-0.04
			0	0.05	0.1	0.15	0.2
					X (m)
	Figure 11: Scenario 7 results of beam shape

Table 7 :

 7 Scenario 7 results of beam-end coordinates

  1 -Z 1 and f x2 , f y2 ,

	Y (m)	0.02 0.06 0.1 0.12 0.14 0.08 0.04	Case 1		The deformed beam shape by MCM Case 2 Undeformed shape The deformed beam shape by FEA Case 3	Loading MCM i g 2 FEA	Case F x (N) F y (N) M o (N.m) q x (N_m) q y (N_m) x(L) (m) y(L) (m) ✓(L) (rad) x(L) (m) y(L) (m) ✓(L) (rad)	Case 1 Case 2 Case 3 Case 4 675 225 *225 *675 225 75 *75 *225 30 10 *10 *30 2500 1250 *1250 *2500 5000 2500 *2500 *5000 *0.007232 0.08018 0.1519 0.1836 0.08681 0.09606 0.06986 0.03081 4.4327 3.3358 2.2057 1.4453 *0.007239 0.08060 0.1520 0.1835 0.8670 0.09579 0.06964 0.03069 4.4320 3.3313 2.2053 1.4473
		-0.02 0					Case 4	ER	x(L) y(L) ✓(L)	0 .10% 0 .13% 0 .02%	0.52% 0.28% 0.14%	0.07% 0.32% 0.02%	0.05% 0.39% 0.14%
		-0.04	0	0.05	0.1	0.15	0.2	CE (s)	MCM (i = 2) FEA	1.79 13	0.50 6	0.46 12	0.41 12
					X (m)
	Figure 12: Scenario 8 results of beam shape

Table 8 :

 8 Scenario 8 results of beam-end coordinates m z2 , x 2 , y 2 , ✓ 2 are defined in its local coordinate system X 2 -Y 2 -Z 2 respectively as displayed in Fig.13b. Then, we are ready to develop the mentioned three types of relationships for modeling the ICB-based compliant leg:1) Constitutive relationshipsHere, we consider Scenario 6 in Section 3.1.6 as an equation of f xi , f yi , m zi , x i , y i , ✓ i , E, w, h, D i (see Fig.13band Fig.

2) Force equilibrium

3) Geometric compatibility relationships

Logically, given the geometric and material parameters E, w, h, L ji , W i , T and any three of F X , F Y , M Z , X, Y and , the rest three can be obtained by solving Eqs. 36 to 38.

Here below are the geometric and material parameters chosen for the parallelogram: E = 69 ù 10 10 Pa; w = 0.01 m; h = 0.001 m; The relationships of F Y against X, F Y against Y and F Y against are plotted in Fig. 18a to Fig. 18c where AM results are compared and also verified by FEA results along with the corresponding errors (the maximum is around 6%). Fig. 18d demonstrates the deformed shapes of the compliant parallelogram evaluated in FEA to graphically present its deformation.

ICB-based compliant parallelogram

In this section, we aim to study ICB-based compliant parallelograms which have two ICB-based lumped compliant legs as shown in Fig. 17b (where the parameter T is introduced to characterize the distance between the two compliant legs). Likewise, modeling this type of parallelograms requires considering constitutive relationships, force equilibrium and geometric compatibility relationships as well.

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.