Gang Zheng • Ke Wu

Gang Zheng

Insight into numerical solutions of static large deflection of general planar beams for Compliant Mechanisms

Keywords:

We formulated the large beam-deflection problem as a general nonlinear BVP.

• We reviewed the commonly used methods for modeling beam deflection in CMs.

• We introduced efficient numerical methods to solve beam deflection problems.

Introduction 1.Compliant mechanisms

As a unique and novel branch of robotic systems, Compliant Mechanisms (CMs) have been presenting their inherent advantages in many existing engineering applications over the last few decades [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]. Unlike conventional rigid-body mechanisms, CMs are most likely to be designed monolithic, serving as devices to transfer motion, force, and energy only through the deformation of the built-in elementary flexible components. According to the state of the art, CMs have demonstrated many meaningful advantages over conventional rigid-body ones, such as increased motion precision [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][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] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF], reduced weight [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF][START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF], cost reduction [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[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][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]. Due to these mentioned merits, CMs have been integrated into many engineering applications, serving as a type of promising alternative to conventional rigid mechanisms [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[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]. Besides, slender straight beams and ICBs are mostly utilized as the elementary flexible components in CMs: straight beams present relatively ideal degrees of freedom (DoF) and degrees of constraint (DoC), serving as a suitable choice for 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], due to their large deflection range however with a relatively small strain range. Therefore, the efficient methods for modeling nonlinear (intermediate or large) deflection of the mentioned two types of slender beams are highly needed for conducting static synthesis of a whole CM. Fortunately, several impressive and great contributions have been completed by the pioneering researchers, which are reviewed in the next section.

Literature review about static modeling of slender beams in the field of CMs

To properly design and then optimize flexure-based CMs, accurate and efficient modeling approaches are of great importance. Logically, the mentioned static modeling process highly depends on that of the involved flexible components as the thick parts are normally treated/assumed as rigid bodies [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF][START_REF] Ling | Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey[END_REF]. In this short review section, we only take into account the commonly-used modeling methods in the research community of CMs.

Indeed, the key essence behind large-deflection scenarios of general planar beams lies in geometrically nonlinear Euler-Bernoulli beam theory, which is essentially a BVP of an ODE. Howell introduced elliptical integrals to handle this problem [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF], which has been the most accurate solution to this beam theory since then. For straightforward implementation, Howell and his co-authors then introduced a more designer-friendly approximation method called pseudo-rigid-body model (PRBM), serving as an efficient modeling method in the early stages of the design process, which has proved its effectiveness in many engineering cases [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]. Much work has been completed regarding its derivatives afterwards [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][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][START_REF] Venkatasubramanian | Pseudo-rigid-body models of initially-curved and straight beams for designing compliant mechanisms[END_REF] for its widespread use. Besides, according to the current literature of CMs, the most commonly-used numerical beam model termed as BCM [START_REF] Awtar | Characteristics of beam-based flexure modules[END_REF], which was proposed in [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF], has been serving as an effective modeling tool in beam-deflection problems since then. Essentially, BCM is also derived from geometrically nonlinear Euler-Bernoulli beam theory by making some approximations on the curvature of the modeled beams [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF]. Its derivatives [START_REF] Chen | Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model[END_REF][START_REF] Ma | Modeling large planar deflections of flexible beams in compliant mechanisms using chained beam-constraintmodel[END_REF] have also been quite popular primarily due to the success of the original BCM. Furthermore, Finite Element Method (FEM) is a widely used method for numerically solving differential equations arising in engineering [START_REF] Singiresu | The finite element method in engineering[END_REF] and in this case, beam-deflection problems can be theoretically described using the governing equations behind solid mechanics and beam theory respectively and logically solved by FEM. It is reliable so many researchers from different research areas tend to use it as a means of verifying their own mathematical models. (Note that here FEM specifically refers to the reliable and accurate post-optimized algorithm built in commercial software or researchers' own coding). In 2020, a comprehensive CM review has been completed in [START_REF] Ling | Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey[END_REF] where the authors 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 method 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 algorithms actually introduce the concept of FEM where the summation of relatively smaller deformation of all discretized elements results in a feasible analysis of larger deformation (geometric nonlinearity). Therefore, we can logically consider all chained BCM [START_REF] Ma | Modeling large planar deflections of flexible beams in compliant mechanisms using chained beam-constraintmodel[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][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][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][START_REF] Venkatasubramanian | Pseudo-rigid-body models of initially-curved and straight beams for designing compliant mechanisms[END_REF] as FEM-related or chained algorithms. They all have been proved to be effective options for model modeling large deflection of slender straight beams and ICBs.

Although the above mentioned methods have already been successfully applied to different applications in practice, these modeling methods may still have their limitations. We would like to evaluate them from the following aspects: physical theories used; analytical or numerical; accuracy within its scope; deformation range; computational expense; beam shape characterization; practicability for multiple existing engineering scenarios; challenges in statically synthesizing large-deflection compliant mechanisms; feasibility for mechanism optimization.

1 Physical theories used: FEM: solid mechanics, Euler Bernoulli beam theory, Timoshenko beam theory; Elliptical integral method: Euler Bernoulli beam theory; (chained-)BCM: Euler Bernoulli beam theory; (chained-)TBCM: Timoshenko beam theory; (chained-)PRBM (mostly Euler Bernoulli beam theory). a Euler Bernoulli beam theory: it assumes that the deformation of the beam only comes from bending deflection, and it also does not consider axial elastic elongation or compression. b Timoshenko beam theory: it assumes that the deformation of the beam comes from both bending deflection and the shear strain of every cross-section, and it also does not consider axial elastic elongation or compression. c Solid mechanics: it considers more complex deformation of any deformable body so in the case studied using solid mechanics to analyze beam-deflection problems is supposed to be the most accurate one. This is why solid-mechanics-based FEM has been widely used as a verification method in the academic field of mechanical engineering.

Remark 1. (Chained-)BCM and (chained-)TBCM both consider the axial elastic stretching by adding an engineering approximation in their governing equations (see Eq. 3.11 and Eq 3.12 in Chapter 3 of Awtar's PhD thesis [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF]) although they are based on Euler Bernoulli beam theory and Timoshenko beam theory.

2 Analytical or numerical: FEM, (chained) BCM, and chained PRBM are numerical methods whereas elliptical integral method is analytical. One of PRBM's main advantages compared to elliptical integral method or other numerical methods is that it simplifies the nonlinear beam-deflection problem into a very straightforward one that can be solved by closed-form formulations [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF][START_REF] Chen | Finding the optimal characteristic parameters for 3r pseudo-rigid-body model using an improved particle swarm optimizer[END_REF], which helps users to design CMs in early stages (essentially, PRBM tends to serve as an engineering or empirical approximation rather than an exact analytical solution to the beam model). Chained PRBM [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][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] Venkatasubramanian | Pseudo-rigid-body models of initially-curved and straight beams for designing compliant mechanisms[END_REF]] logically provides better accuracy however at the cost of transforming the closed-form approximation into a numerical one. 3 Accuracy within its scope: Commercial FEM-based software is reliable to trust as most of them have passed the standard modeling tests proposed by International Association Engineering Modeling (NSFEMS), which is why researchers tend to use it as a means of verification [START_REF] Singiresu | The finite element method in engineering[END_REF]. In CM modeling, a huge portion of designs employ slender beams as the elementary flexible components so any method based on geometrically nonlinear Euler-Bernoulli beam theory (it is chosen here since we can ignore shear deformation on the cross section of slender beams) can do the trick. Obviously, elliptical integral method provides the most accurate solution in this case. PRBM loses some accuracy due to its own assumptions made (it is limited in terms of estimating larger end slope of slender beams [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[END_REF]). Besides, BCM can accurately capture intermediate deformation whereas it cannot handle large-deflection scenarios due to its assumptions made. 4 Deformation range: FEM, elliptical integral method, chained BCM, and (chained) PRBM are valid analyzing large deformation whereas BCM is valid for intermediate deformation.

Remark 2. In particular, nonlinear buckling is a type of large deflection for slender beams. Logically, BCM and TBCM cannot handle nonlinear buckling [START_REF] Awtar | Synthesis and analysis of parallel kinematic XY flexure mechanisms[END_REF]. Chained-BCM and chained-TBCM are used to analyze bistable mechanisms (with post-buckling involved [START_REF] Chen | Kinetostatic modeling of fully compliant bistable mechanisms using timoshenko beam constraint model[END_REF][28]) so is chained-PRBM. There's also work just using PRBM to handle nonlinear buckling but its accuracy is compensated by classic beam theory [START_REF] Sönmez | A compliant bistable mechanism design incorporating elastica buckling beam theory and pseudo-rigid-body model[END_REF].

Main contributions

Clearly, the existing methods have already proved to be effective in many used cases under some reasonable assumptions related to practices. However, those methods have as well their own application limitations, for example, elliptical integrals cannot handle non-tip-loading cases and BCM requires small vertical deflection, etc. Compared to those existing results, the main contributions of this paper are as follows.

• Formulation: we formulate the large beam-deflection problem as a nonlinear BVP of an ODE by applying Euler-Bernoulli Beam theory. Note that such a formulation can be used as well to model more general and complex scenarios by considering distributed forces, moments, and so on. However, for the sake of simplicity, this paper focuses only on two commonly-used flexible members in CMs [START_REF] Ling | Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey[END_REF]: straight beams and initially curved beams (ICBs).

• Numerical methods: we apply several useful numerical methods (including finite difference method, shooting method, weighted residual methods, and so on) to solve the formulated BVPs that describe large beamdeflection problems. We have presented detailed derivations of each method regarding solving the BVPs and provided in-depth discussions about the advantages and differences of these methods. Numerical results are presented accordingly.

• Proving feasibility in modeling CMs: we have proved the feasibility of the mentioned numerical methods in modeling 2 representative CMs: straight-beam-based cross-axis flexural pivots [START_REF] Larry | Handbook of compliant mechanisms[END_REF] and ICB-based hinges [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]. All numerical methods have been used and compared in modeling the 2 CMs with numerical results provided accordingly.

Modeling of large deflection via geometrically nonlinear Euler-Bernoulli Beam Theory

In this paper, all of our proposed methods aim to model large beam-deflection problems via geometrically nonlinear Euler-Bernoulli beam theory. Assuming that, for any point along the beam axis, its corresponding curvature d d is proportional to the exerted moment ( ) at (named as elastic assumption), i.e.,

d d = ( ) (1)
where is the elastic modulus and is the second moment of area of the beam's cross-section. If we consider the cross section varies along the beam axis (i.e., the second moment of inertia of the cross section ( ) is a function of ), we will have:

= ( ) ( )
Note that ( ) is formulated based on the centroidal axis, and it 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 and the beam length is defined along it as well. Besides, the centroidal axis lies in a plane so only planar beams are studied in this paper. If the slender beam has an initial curvature 1 ( ) where we consider the varying curvature by defining 1 ( ) as a function of , we will logically arrive at:

= ( ) ( ) + 1 ( ) (2)
Here, ( ) is defined as the exerted moment at any arbitrary point on the beam, noted as ( , ( ), ( )) due to all external loads on the beam, such as beam-tip loads, distributed loads, distributed moment and pressure along the beam body. In this paper, we only focus on two representative examples: beam-tip loads on constant-cross-section straight beams and constant-cross-section ICBs of constant curvature (i.e., ( ) and ( ) are both constant noted as and in the following) respectively since these two loading scenarios frequently happen in CMs. Therefore, we can derive the corresponding ( ) according to the general beam-tip loading shown in Fig. 1:

( ) = ( -( )) + ( -( )) + (3)
where an arbitrary force (which can be decomposed into a vertical force and a horizontal force ) and an arbitrary moment at its beam end ( , ) or equivalently along the beam axis. Obviously, we can also consider many other loading scenarios, such as flexible beams with varying cross-section, varying curvature, pressure and so on, and it is evident that all those common loading cases can be easily formulated like (3) by recalculating its right-hand terms with correct applied moments at .

Rearranging Eq. ( 2) and Eq. ( 3), we will end up with:

d d = ( -( )) + ( -( )) + + (4)
Then, conducting differentiation on Eq. ( 4) with respect to results in the following ODE:

d 2 d 2 = - (cos ( ) + sin ( )) (5)
As stated above, we will focus on the following two cases since they commonly-used in the design process of CMs as standards parts:

Case 1: Straight beams:

1 = 0;
Case 2: Circular beams (constant-curvature beams): 1 = ( is an arbitrary constant real number).

For these two cases, different numerical methods will be detailed in Section 5 to solve large beam deflection problems, based on which two CMs (straight-beam-involved one and circular-beam-involved one) will be analyzed in Section 6. Remark 3. Note that there is a special type of beam-deflection scenarios called nonlinear buckling where the studied beam is only subjected to an axial force shown in Fig. 2. Its governing equation can be rearranged from BVP (6) as follows:

D.E. d 2 d 2 = - (cos ( ) + sin ( )) = - sin ( ) B.C. (0) = 0 d d ( ) = + 1 = 0 (7)
where = 0, = 0 and 1 = 0. Physically speaking, if reaches the buckling limit, the straight slender beam will experience nonlinear buckling, possibly resulting in 3 different (deformed) beam shapes shown in Fig. 2.

Mathematically speaking, it is normal that BVPs may have multiple solutions given some specific conditions. Normally, in applied mathematics, the concept of bifurcation theory [START_REF] Chow | Methods of bifurcation theory[END_REF] needs to be introduced to solve these BVPs where some techniques are used, such as adjusting initial guesses, determining the boundary conditions by making most of the information of the governing equation, perturbation method and so on. This topic is out of the scope of this paper. However, a comprehensive solution to handle nonlinear buckling problems of slender beams numerically has been recently proposed in [START_REF] Wu | Theoretical analysis on nonlinear buckling, post-buckling of slender beams and bi-stable mechanisms[END_REF].

Recall of some useful numerical methods

In this paper, we will adopt different numerical methods which are necessary to solve beam-deflection problems via geometrically nonlinear Euler-Bernoulli beam theory, the majority of which are based on Weierstrass's first theorem and Newton-Raphson method. Runge-Kutta method and Euler method are also involved in this paper. Therefore, we first proceed with a quick theory recap on them hereafter for a better demonstration of the rest of the paper.

Numerical method to approximate continuous function

Recall that this paper aims to find a continuous solution of the geometrically nonlinear Euler-Bernoulli beam model [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF], but its exact analytical solution is very hard to find. Therefore, we bypass the difficult task to approximate the exact analytical solution using polynomials as basis functions, which is based on the employment of Weierstrass's first theorem. This theorem implies that, instead of finding the exact solution to BVP (6), we can find some specific polynomials that can quite accurately approximate the exact solution.

Numerical method to solve algebraic equation

There are several methods to handle nonlinear equations or systems of equations reported in the literature. In this paper, we would like to proceed with a commonly used method called Newton-Raphson method for these problems. Here, a brief recap is presented to solve a general nonlinear equation ( ) = 0 ∈ ℝ ×1 with ∈ ℝ ×1 ( ≥ ), which can be realized by the following iterative equation:

= -1 - ( ) -1 | = -1 ( -1 )
where represents the th iterative solution with 0 being a given initial guess value. After the convergence (in practice the convergence criteria is set as || +1 -|| ≤ where is a pre-defined threshold), we can finally obtain the solution to ( ) = 0. The dimension of the Jacobian matrix ( ) is ℝ × . Obviously, the major computation of this algorithm lies in calculating the pseudo inverse ( > ) or the inverse ( = ) of ( ) [42][43]. Therefore, the higher values of and of the studied matrix, which leads to higher dimension of the matrix, would result in more computational time for computing the pseudo inverse or the inverse of this matrix.

Numerical methods to approximate ODE with known initial conditions

Given a general nonlinear ODE

d d = ( , )
for ∈ [ 0 , ] with known initial condition ( 0 ) = 0 . It is difficult or impossible to find out its analytical solution if ( , ) is highly nonlinear. Numerical methods are then needed to be used to approximate its solution. Euler method and Runge-Kutta method are two commonly used ones, which will be recalled hereafter. Consider = ( ) is the solution to the above ODE with ( 0 ) = 0 where 0 being the known initial condition.

Euler's method: Denote ℎ as the fixed sampling step of the space [ 0 , ] for , i.e., +1 = + ℎ , and then we proceed with +1 = + ( , )ℎ for 0 ≤ ≤ -1. Finally, we will end up with ( +1 ) ≈ +1 .

Runge-Kutta method: Similarly, Runge-Kutta method is also a popular choice for numerically solving ODE with known initial conditions. The formula of the fourth Runge-Kutta method is presented below. Similarly, ℎ is defined as the step size. Then, we proceed with the following formula:

1 = ℎ ( , ), 2 = ℎ ( + ℎ 2 , + 1 2 ), 3 = ℎ ( + ℎ 2 , + 2 2 ), 4 = ℎ ( + ℎ , + 3 ) +1 = + 1 6 ( 1 + 2 2 + 2 3 + 4 )
which computes the approximate solution as ( +1 ) ≈ +1 .

Applying different numerical methods to solving beam deflection problems

In this section, numerical methods are explained in terms of solving geometrically nonlinear Euler-Bernoulli beam model [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF], and the detailed derivation is presented accordingly, where numerical methods presented in Section 3 will be applied. Note that apart from Non-iterative IVP method (NIM) stated in Section 4.3, the other methods can directly deal with inflection points. After having the numerical solution of BVP (6), i.e. Θ( ), we can then characterize the deformed beam shape via:

( ) = ∫ 0 cos(Θ( ))d ; ( ) = ∫ 0 sin(Θ( ))d ; (8)
for ∈ [0, ] with being the length of the deformed beam.

Finite difference method

Finite difference method defines sampling points (0

= 0 < Θ 1 < Θ 2 < Θ 3 < Θ 4 ... < Θ -1
) as unknowns.

Then by using those points to approximate

d 2 d 2 ≈ Θ +1 -2Θ +Θ -1
[ ∕( -1)] 2 , the BVP ( 6) can be written into the following discrete manner:

Θ 2 -2Θ 1 + 0 [ ∕( -1)] 2 + (cos Θ 1 + sin Θ 1 ) = 0 Θ +1 -2Θ + Θ -1 [ ∕( -1)] 2 + (cos Θ + sin Θ ) = 0, for 2 ≤ ≤ -2 Θ -1 -Θ -2 ∕( -1) - -1 ∕( -1) + (cos Θ -1 + sin Θ -1 ) = 0 (9)
where the first and the last equation take into account the two boundary conditions presented in [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF]. Clearly, there are altogether equations available and unknowns to be solved. Then, we can numerically handle ( 9), which contains -1 equations, using Newton-Raphson method presented in Section 3 to obtain the values of Θ for 0 ≤ ≤ -1, equivalently the approximation Θ( ) of the solution ( ) to BVP [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF]. Besides, we take one step further to add an iterative strategy for automatically choosing the number of sampling points : if |Θ ( ) -Θ -1 ( )| ≤ where is the th option ( = 1, 2, 3...) for and Θ ( ) is the approximation under sampling points, then is settled: = ; if not, = +1 = + 1. This strategy enables the method adaptive to solve the studied system (9) without manually setting up .

Shooting method

As the term implies, shooting method is essentially a trial method via mathematical methods to direct its guessing to reach the corrective answer. Here, we first need to transform BVP ( 6) into an IVP (initial value problem) which can be directly solved using Euler's method or Runge-Kutta method presented in Section 3:

D.E.
( ) = ( ) -( ( ) + 1 ( ) ) = 0 (11)
We can solve this type of nonlinear equations via utilizing the presented Newton-Raphson method in Section 3. After the convergence, we can then obtain the desired initial condition for IVP [START_REF] Sam Han | Design, simulation, and fabrication of a quadstable monolithic mechanism with x-and y-directional bistable curved beams[END_REF]. Therefore, ( ) can be approximated via Euler's method or Runge-Kutta method presented in Section 3 for IVP [START_REF] Sam Han | Design, simulation, and fabrication of a quadstable monolithic mechanism with x-and y-directional bistable curved beams[END_REF].

Non-iterative IVP method

Inspired by [START_REF] Bs Shvartsman | Large deflections of a cantilever beam subjected to a follower force[END_REF], we propose an IVP-based non-iterative method for BVP (6) by using all the boundary and ODE information. After a straightforward manipulation of (6), we have

d 2 d 2 = d d ( 2 2 ) = - (cos + sin ) (12)
where denotes the curvature d d . Then conducting integration of ( 12) constrained by the B.C. of BVP [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF], we obtain

= d d = ± 2 (cos -cos ) + 2 (sin -sin ) + ( + 1 ) 2 0.5 (13)
Again, integrating (13), we will finally arrive at:

= ∫ 0 ± 2 (cos -cos ) + 2 (sin -sin ) + ( + 1 ) 2 -0.5 d (14)
Obviously, we can numerically compute the value of easily via using Newton-Raphson method. The sign of d d should not change when 0 ≤ ≤ in Eq .( 13), which means this method does not consider the scenarios where inflection points exist along the deformed beam. Besides, the sign of d d needs to be determined before solving [START_REF] Ling | Kinetostatic and dynamic modeling of flexure-based compliant mechanisms: a survey[END_REF], either + or -. In practice, if the numerical solution of ( 14) can not be found using Newton-Raphson method given a proper convergence criteria, rechecking the sign of d d is always essential. Then, BVP ( 6) can be transformed into:

D.E. d 2 d 2 = - (cos ( ) + sin ( )) B.C. ( ) = d d ( ) = + 1 (15)
where the boundary conditions are all given at the end of beam = . Therefore, if we consider the above ODE with an initial condition defined at = , i.e., a special IVP, then we can directly handle it to approximate ( ) via either Euler's method or Runge-Kutta method presented in Section 3 without touching the iterative process of shooting method.

Weighted residual methods

As reported in the literature, weighted residual method [START_REF] Hatami | Weighted residual methods: principles, modifications and applications[END_REF] is an approximation technique for solving differential equations. According to Theorem 1, we can then define a function Θ( ), which is an th-order polynomial, to approximate the solution of BVP (6), i.e.,

Θ( ) = ∑ =0 (16)
where there are + 1 unknown coefficients ( = 0, 1, 2, 3... ). With the above chosen polynomial solution Θ( ), an approximation residual can be defined as follows:

( ) = d 2 Θ d 2 + (cos Θ( ) + sin Θ( )) ≠ 0 (17)
In terms of satisfying the boundary conditions defined in (6), we then have the corresponding constraints expressed via Θ( ) as follows:

B.C. Θ(0) = 0 dΘ d ( ) = + 1 (18)
Therefore, we logically arrive at:

Θ(0) = ∑ =0 0 = 0 dΘ d ( ) = ∑ =0 -1 = + 1 (19)
The essence of weighted residual method is to force the residual [START_REF] Chen | Finding the optimal characteristic parameters for 3r pseudo-rigid-body model using an improved particle swarm optimizer[END_REF] to zero in some average sense over the domain:

∫ ( )d = 0 ( = 1, 2, 3... -1) (20)
where is the weight function; -1 refers to the least number of weight functions, equivalently the least number of the equations derived from (20) (for instance, we could have + 1 equations derived from [START_REF] Yu | 5r pseudo-rigid-body model for inflection beams in compliant mechanisms[END_REF] in least square method in Section 4.4.3); is the corresponding integral domain. Logically, we can arrive at deducing the + 1 unknown coefficients defined in ( 16) via at least -1 + 2 = + 1 equations from [START_REF] Yu | A novel model of large deflection beams with combined end loads in compliant mechanisms[END_REF] and [START_REF] Yu | 5r pseudo-rigid-body model for inflection beams in compliant mechanisms[END_REF]. Such a problem can be solved by applying Newton-Raphson method presented in Section 3. Note that in [START_REF] Yu | A novel model of large deflection beams with combined end loads in compliant mechanisms[END_REF] and [START_REF] Yu | 5r pseudo-rigid-body model for inflection beams in compliant mechanisms[END_REF], there are at least 3 equations available so we have + 1 ≥ 3, equivalently ≥ 2, which means 2th or higher-order polynomials should be used to approximate the solution of BVP [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF].

In order to apply weighted residual methods, it is necessary to initialize a trial guess for ( 16) to start its iterative computation which might play an important role in the convergence and local minima. In this context, we can have two different ways to set the initial values of :

1. Direct way: We randomly initialize the values of , which yields a random shape of Θ( ). 2. Indirect way: We choose a globally initial guess for Θ( ), such as trigonometric functions, which yields the corresponding for the Θ( ).

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.

Similarly, we can also implement the iterative strategy for automatically choosing the order of pre-set polynomial : if |Θ ( ) -Θ -1 ( )| ≤ where is the th option ( = 1, 2, 3...) for and Θ ( ) is the approximation using a polynomial of th order, then is settled: = ; if not, = +1 = + 1. This strategy enables the method adaptive to solve the studied system of equations ( 19) and ( 20) without manually setting up .

Another important issue related to weighted residual methods is how to choose the weight function. Intuitively, different choices of those weight functions lead to different implementations of the weighted residual methods with different performances. The following explains some popularly used techniques to choose weight functions.

Collocation method

As a commonly-used type of weighted residual method, 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]. In this method, the weight functions are taken from the family of Dirac function in the domain, i.e.

( ) = ( -). The Dirac function has the following property:

( ) = ( -) = 1 = 0 otherwise ( 21
)
where is termed as collocation points. Substituting ( 16), ( 17) and ( 21) into [START_REF] Yu | 5r pseudo-rigid-body model for inflection beams in compliant mechanisms[END_REF] results in the forcing of the residual to zero at specific points, i.e.,

∫ 0 ( ) ( )d = ( ) = d 2 Θ d 2 ( ) + (cos Θ( ) + sin Θ( )) = d 2 ( ∑ =0 ) d 2 ( ) + (cos( ∑ =0 ) + sin( ∑ =0 )) ( = 1, 2, 3... -1). (22)
where = and ( 22) contains -1 equations. Considering 2 extra boundary conditions, there are altogether + 1 equations to solve + 1 polynomial coefficients ( = 0, 1, 2, 3... ).

Subdomain method

Subdomain method can be regarded as a modification of the presented collocation method. The idea of such a method is to force the weighted residual to zero not only at fixed collocation points in the domain but also over various subsections of the domain. To implement this method, the weight functions are set to the same integers, and the integral over the entire domain is broken into ( -1) subdomains that are sufficient to evaluate all unknown coefficients defined in [START_REF] Yu | A pseudo-rigid-body 2r model of flexural beam in compliant mechanisms[END_REF]. Substituting ( 16) and ( 17) into (20), we will have:

∫ -1 ( -1) -1 ( )d = ∫ -1 ( -1) -1 [ d 2 Θ d 2 ( ) + (cos Θ( ) + sin Θ( ))]d = ∫ -1 ( -1) -1 [ d 2 ( ∑ =0 ) d 2 + (cos( ∑ =0 ) + sin( ∑ =0 ))]d ( = 1, 2, 3... -1). (23)
where ( 23) contains -1 equations. Considering 2 extra boundary conditions, there are altogether + 1 equations to solve + 1 polynomial coefficients ( = 0, 1, 2, 3... ).

Least square method

In this method, the continuous summation of all the squared residuals is minimized:

min 0 , 1 , 2 ... ( ), with ( ) = ∫ 0 2 ( )d s.t. Θ( ) = dΘ d ( ) = + 1 (24)
where 0 , 1 , 2 ... are the polynomial unknown coefficients defined in [START_REF] Yu | A pseudo-rigid-body 2r model of flexural beam in compliant mechanisms[END_REF]. In order to achieve a minimal of this scalar function, the derivatives of ( ) with respect to all the unknown parameters must be zero:

( ) = 2 ∫ 0 ( ) ( ) d = 0 ( = 0, 1, 2, 3... ) (25)
Comparing with [START_REF] Yu | 5r pseudo-rigid-body model for inflection beams in compliant mechanisms[END_REF], the weight functions in this method are chosen as:

( ) = ( ) ( = + 1 = 1, 2, 3... + 1) (26)
We will then end up with:

∫ 0 ( ) ( ) = ∫ 0 [ d 2 Θ d 2 ( ) + (cos Θ( ) + sin Θ( ))] [ d 2 Θ d 2 ( ) + (cos Θ( ) + sin Θ( ))]d = ∫ 0 [ d 2 ( ∑ =0 ) d 2 + (cos( ∑ =0 ) + sin( ∑ =0 ))] [ d 2 ( ∑ =0 ) d 2 + (cos( ∑ =0 ) + sin( ∑ =0 ))]d ( = 1, 2, 3... + 1) (27)
where ( 27) contains + 1 equations. Considering 2 extra boundary conditions, there are altogether + 3 equations to solve + 1 polynomial coefficients ( = 0, 1, 2, 3... ).

Galerkin method

In Galerkin method, instead of using the derivative of the residual with respect to the unknown defined in ( 16), the derivative of the approximating function is used, and the weight functions are chosen as follows:

( ) = Θ( ) ( = + 1 = 1, 2, 3... + 1) (28)
Hence [START_REF] Yu | 5r pseudo-rigid-body model for inflection beams in compliant mechanisms[END_REF] will yields

∫ 0 Θ( ) ( )d = ∫ 0 [ ∑ =0 ] [ d 2 ( ∑ =0 ) d 2 + (cos( ∑ =0 ) + sin( ∑ =0 ))]d ( = 1, 2, 3... + 1) (29)
where ( 29) contains + 1 equations. With 2 extra boundary conditions, there are altogether + 3 equations to solve + 1 polynomial coefficients ( = 0, 1, 2, 3... ).

Method of moments

In this method, the weight functions are chosen from the family of polynomials:

( ) = , = 1, 2, 3... -1 (30)
Rearranging (20) will yields

∫ 0 ( ) ( )d = ∫ 0 [ d 2 Θ( ) d 2 + (cos Θ( ) + sin Θ( ))]d = ∫ 0 [ d 2 ( ∑ =0 ) d 2 + (cos( ∑ =0 ) + sin( ∑ =0 ))]d ( = 1, 2, 3... -1) (31)
where [START_REF] Chen | Modeling large spatial deflections of slender bisymmetric beams in compliant mechanisms using chained spatial-beam constraint model[END_REF] contains -1 equations. With 2 extra boundary conditions, there are altogether + 1 equations to solve + 1 polynomial coefficients ( = 0, 1, 2, 3... ).

Rayleigh Ritz method

In this method, the principle of integration by parts is used:

∫ 0 ( )d = ∫ [ d 2 Θ( ) d 2 + (cos Θ( ) + sin Θ( ))]d = [ dΘ( ) d ] 0 -∫ 0 d d dΘ( ) d d + ∫ 0 (cos Θ( ) + sin Θ( ))d = 0 = [ d( ∑ =0 ) d ] 0 -∫ 0 d d d( ∑ =0 ) d d + ∫ 0 (cos( ∑ =0 ) + sin( ∑ =0 ))d = 0 (32)
for 1 ≤ ≤ + 1, where the weight functions are chosen in the following manner:

= Θ( ) ( = + 1 = 1, 2, 3... + 1) (33)
Note that (32) contains + 1 equations. With 2 extra boundary conditions, there are altogether + 3 equations to solve + 1 polynomial coefficients ( = 0, 1, 2, 3... ).

Weighted residual methods from global to local

From Section 4.4.1 to Section 4.4.6, all the mentioned weighted residual methods are handled in a global manner. However, these methods can be all implemented in a local manner for higher accuracy and approximating complex solutions of more BVPs, such as BVPs' solutions of highly frequent vibration (where in the case of beam-deflection problems, it refers to high-order buckling phenomena). Unlike the procedure explained in Section 4.4.1 to Section 4.4.6 from a global perspective, the following explains how BVP ( 6) is handled using one of the mentioned weighted residual methods, i.e. collocation method in a local manner [START_REF] Lawrence F Shampine | Solving boundary value problems for ordinary differential equations in matlab with bvp4c[END_REF].

Generally speaking, by dividing the interval [0, ] into + 1 points 0 = 0 < 1 < ⋯ < = , collocation method approximates the solution of the above BVP via linear combination of well-chosen basis functions by satisfying D.E. and B.C. at all collocation points. Precisely, denote by Θ( ) the approximate solution of BVP [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF], and suppose that a cubic polynomial ( ) is chosen to approximate the solution of BVP [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF]

on each sub-interval [ -1 , ], i.e., Θ( ) = ∑ =1 ( ) where
( ) = ∑ 3 =0 , ∀ ∈ -1 , 0, otherwise (34)
with for 0 ≤ ≤ 3 being the 4 unknown polynomial coefficients on the th interval. In total, we have 4 unknown parameters to be identified. Clearly, Θ( ) should satisfy the D.E. in BVP ( 6) for all collocations points, and this yields the following 3 algebraic equations to be fulfilled:

d 2 Θ d 2 ( -1 ) + (cos Θ( -1 ) + sin( -1 )) = d 2 ( ∑ 3 =0 ) d 2 ( -1 ) + (cos( 3 ∑ =0 -1 ) + sin( 3 ∑ =0 -1 )) = 0 d 2 Θ d 2 ( ̄ -1 ) + (cos Θ( ̄ -1 ) + sin( ̄ -1 )) = d 2 ( ∑ 3 =0 ) d 2 ( ̄ -1 ) + (cos( 3 ∑ =0 ( ̄ -1 ) ) + sin( 3 ∑ =0 ( ̄ -1 ) )) = 0 d 2 Θ d 2 ( ) + (cos Θ( ) + sin( )) = d 2 ( ∑ 3 =0 ) d 2 ( ) + (cos( 3 ∑ =0 ) + sin( 3 ∑ =0 )) = 0 (35)
where ̄ -1 = -1 + 2 . Note that the approximated solution Θ( ) 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.,

-1 ( -1 ) = ( -1 ), ∀ ∈ [2, ] (36)
Hence we have the above -1 algebraic equations due to the continuity conditions. Together with the 2-dimensional B.C. in [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF], we finally obtain a set of algebraic equations with dimension 3 + -1 + 2 = 4 + 1, which will be used to compute those 4 unknowns parameter defined in [START_REF] Wu | Efficient spatial compliance analysis of general initially curved beams for mechanism synthesis and optimization[END_REF]. So far, we have demonstrated the 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]. In common weighted residual methods (see Section 4.4.1 to Section 4.4.6), they are all implemented in a global manner, which means a single global polynomial Θ( ) ∈ 1 is used to approximate the solution of (6). However, the solution of the locallyperformed collocation method only ensures the continuity of Θ( ), equivalently Θ( ) ∈ 0 , which ends up with non-differentiable piecewise polylines for the solution Θ( ). Apparently, this non-differentiability of Θ( ) contradicts the nature of beam deflection since ( ) is always continuous and differentiable, equivalently ( ) ∈ 1 [START_REF] Peter Boresi | Advanced mechanics of materials[END_REF]. To treat such a problem, we take one step further to modify the collocation method via ensuring the differentiability of Θ( ), i.e., Θ( ) ∈ 0 . By imposing the following derivative continuity conditions at all collocation points:

d -1 d ( -1 ) = d d ( -1 ), ∀ ∈ [2, ] (37)
This extra added constraints of derivative continuity leads to a differentiable global approximation of the solution of the studied BVP [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF], bypassing the need of conducting curve fitting process afterwards. It is clear that (37) yields another -1 equations, ending up with altogether 4 + 1 + -1 = 5 algebraic equations available. In summary, the collocation method enables us to solve the 4 unknowns based on the 4 + 1 equations from ( 6), ( 35) and ( 36), or 5 equations from ( 6), [START_REF] Jin | Design of compliant mechanisms using a pseudo-rigid-body model based topology optimization method[END_REF] [START_REF] Sun | Fem-based mechanics modeling of bio-inspired compliant mechanisms for medical applications[END_REF] and [START_REF] Sun | Automatic design in matlab using pde toolbox for shape and topology optimization[END_REF] which can be numerically solved via Newton-Raphson method presented in Section 3.

Numerical results and detailed analysis of single-beam modeling

As we have stated in Section 2, BVP [START_REF] Fan | Design and development of a novel monolithic compliant xy stage with centimeter travel range and high payload capacity[END_REF] is the key we need to handle, which yields two categories of beamdeflection problems: large deflection of planar slender straight beams and circular beams (see Fig. 3 and Fig. 14). These two scenarios are commonly utilized in CMs to provide desired motions according to the current literature [START_REF] Larry | Compliant mechanisms. In 21st century kinematics[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][START_REF] Lobontiu | Compliant mechanisms: design of flexure hinges[END_REF], and they are easy to promote in the standard manufacturing process.

In the following, we study large deflection problems of a planar slender straight beam and a circular beam under a series of loading cases. All numerical methods in Section 4 are used respectively to solve the mentioned beamdeflection problems, followed by FEM verification. Comparisons are made among all the used numerical methods and FEM with detailed analysis provided afterwards.

Large deflection problems of planar slender straight beams

In this section, we solve one single large deflection problem of a planar slender straight beam using all numerical methods mentioned in Section 4 respectively. To start with, the geometric and material parameters chosen for the studied straight beam (see Fig. 3) is shown in the following: To both verify the accuracy and efficiency of these numerical methods, we use beam elements in FEM where the information of mesh elements is presented here:

Edge element: 19; Mesh vertices: 20

Since the shear strain is not considered in Euler-Bernoulli beam model, is logically not involved in our numerical algorithms.

Large deflection problems of slender circular beams

In this section, we solve one single large deflection problem of a planar circular beam using all numerical methods mentioned in Section 4 respectively. Similarly, the geometric and material parameters chosen for the studied circular beam (see Fig. 14) is shown in the following:

Numerical results and discussions

Generally speaking, all proposed methods that are used to solve beam-deflection problems in Section 4 present higher accuracy and higher efficiency compared to FEM. In Figs. 4 to 13 and Fig. 15 to 24, the graphical results of deformed beam shapes are provided using the methods discussed in Section 4 respectively, followed by FEM verification. Likewise, in Tables 1 to 10 and 11 to 20, the beam-end coordinates are also presented respectively and Note that 6th-order polynomials are used to approximate the solutions.

Table 7: Beam-end coordinates by GM

Finite difference method

The graphical results of deformed beam shapes and beam-end coordinates are presented in Figs. equations need to be solved using Newton-Raphson method. 2. Deflection of the slender circular beam: as noted in Table 11, the maximum error reaches 1.29%, and the maximum computational expense is around 3 s. Besides, 1280 sampling points are used which means that 1279 equations need to be solved using Newton-Raphson method.

Compared with FEM, FDM presents reliable accuracy and better efficiency in general. In particular, solving deflection problems of the studied slender circular beam is more time-consuming than that of the studied slender circular beam. This is because more sampling points defined lead to more unknowns or equations to be solved using Newton-Raphson method (see Section 3.2).

Shooting method

The graphical results and beam-end coordinates are presented in Figs. Note that Runge-Kutta method is used, and the predefined threshold for convergence is set as 1 -20 .

Table 12:

Beam-end coordinates by SH 1. Deflection of the slender straight beam: as noted in Table 2, the maximum error reaches 0.09%, and the maximum computational expense is around 0.91 s. Besides, Runge-Kutta method is used in each iteration, and the predefined threshold for convergence is set as 1 -20 . the predefined threshold for convergence is set as 1 -20 .

Deflection of the slender circular beam: as noted in

According to the results shown in Tables 2 &12, SH has been proved feasible and more efficient in terms of solving beam-deflection problems following the framework of Euler-Bernoulli beam theory. Note that 6th-order polynomials are used to approximate the solutions.

Table 16: Beam-end coordinates by LSM

Non-iterative IVP method

The graphical results and beam-end coordinates are presented in Figs. 6 & 17 and Tables 3 & 13 where noniterative IVP method is denoted by NIM.

1. Deflection of the slender straight beam: as noted in Table 3, the maximum error reaches 0.09%, and the maximum computational expense is around 0.30 s. Obviously, there is no inflections points in Case Note that 4th-order polynomials are used to approximate the solutions.

Table 18: Beam-end coordinates by MOM Case 3 and Case 4. In Case 1 and Case 2, d d keeps positive along [0, ] so the positive sign + is chosen for Eq. ( 13). On the contrary, in Case 3 and Case 4, d d keeps negative along [0, ] so the negative signis chosen for Eq. ( 13). Finally, Runge-Kutta method is used to approximate the solution of the studied IVP [START_REF] Su | A pseudorigid-body 3r model for determining large deflection of cantilever beams subject to tip loads[END_REF]. 2. Deflection of the slender circular beam: as noted in Table 13, the maximum error reaches 0.16%, and the [START_REF] Su | A pseudorigid-body 3r model for determining large deflection of cantilever beams subject to tip loads[END_REF], which is not suitable for synthesizing/modeling complex compliant mechanisms where the deformed shapes of slender beams are difficult to predict, especially the existence of inflection points and the sign of d d of each deformed slender beams.

Collocation method

The graphical results and beam-end coordinates are presented in Figs. 7 & 18 and Tables 4 & 14 where collocation method is denoted by CMG.

1. Deflection of the slender straight beam: as noted in Table 4, the maximum error reaches 2.68%, and the maximum computational expense is around 0.22 s. Besides, 6th-order polynomials are used to approximate the solution. 2. Deflection of the slender circular beam: as noted in Table 14, the maximum error reaches 1.28%, and the maximum computational expense is around 0.23 s. Besides, 6th-order polynomials are used to approximate the solution.

As explained in Section 4.4.1, CMG satisfies the constraints at all defined collocation points (see Eq. ( 22)) and the 2 given boundary conditions [START_REF] Zhu | Pseudo-rigid-body model for the flexural beam with an inflection point in compliant mechanisms[END_REF]. Therefore, the nature of CMG does not take into account forcing the weighted residuals to be 0 over the whole domain [0, ], which could possibly lead to some relatively large errors (residuals) at some non-collocation points, such as the beam end ( = ). Obviously, higher-order polynomials result in more accurate approximation since more collocation points will be used to satisfy the constraints as well as obtain the polynomial coefficients.

Subdomain method

The graphical results and beam-end coordinates are presented in Figs. 8 & 19 and Tables 5 & 15 where subdomain method is denoted by SUM.

1. Deflection of the slender straight beam: as noted in Table 5, the maximum error reaches 0.09%, and the maximum computational expense is around 0.37 s. Besides, 6th-order polynomials are used to approximate the solution. 2. Deflection of the slender circular beam: as noted in Table 15, the maximum error reaches 0.16%, and the maximum computational expense is around 0.34 s. Besides, 6th-order polynomials are used to approximate the solution.

Different from CMG, SUM needs to satisfy the constraints where the residuals need to be forced to 0 over several subdomains. This method logically presents better accuracy than CMG, and the results shown in Table 5 and 15 have also proved the conclusions.

Least square method

The graphical results and beam-end coordinates are presented in Figs. 9 & 20 and Tables 6 &16 where least square method is denoted by LSM.

1. Deflection of the slender straight beam: as noted in Table 6, the maximum error reaches 0.09%, and the maximum computational expense is around 0.50 s. Besides, 6th-order polynomials are used to approximate the solution. 2. Deflection of the slender circular beam: as noted in Table 16, the maximum error reaches 0.16%, and the maximum computational expense is around 0.50 s. Besides, 6th-order polynomials are used to approximate the solution.

Essentially speaking, LSM reformulates the form of forcing residuals to 0 into a minimization problem where the corresponding weight functions are determined. However, LSM often suffers from falling into local minimal due to its nature, especially dealing with large deflection problems of circular beams.

Galerkin method

The graphical results and beam-end coordinates are presented in Figs. 10 & 21 and Tables 7 & 17 where Galerkin method is denoted by GM.

1. Deflection of the slender straight beam: as noted in Table 7, the maximum error reaches 0.20%, and the maximum computational expense is around 0.31 s. Besides, 6th-order polynomials are used to approximate the solution. 2. Deflection of the slender circular beam: as noted in Table 17, the maximum error reaches 0.22%, and the maximum computational expense is around 0.21 s. Besides, 6th-order polynomials are used to approximate the solution.

In GM, the weight function are chosen via following [START_REF] Wu | Theoretical analysis on nonlinear buckling, post-buckling of slender beams and bi-stable mechanisms[END_REF]. GM both presents reliable accuracy and good convergence compared to CMG and LSM.

Method of moments

The graphical results and beam-end coordinates are presented in Figs. 11 & 22 and Tables 8 & 18 where method of moment is denoted by MOM.

1. Deflection of the slender straight beam: as noted in Table 8, the maximum error reaches 0.34%, and the maximum computational expense is around 0.30 s. Besides, 4th-order polynomials are used to approximate the solution. 2. Deflection of the slender circular beam: as noted in Table 18, the maximum error reaches 0.19%, and the maximum computational expense is around 0.14 s. Besides, 4th-order polynomials are used to approximate the solution.

MOM only needs 4th-order polynomials to properly approximate the solution, which ends up a smaller number of unknowns and equations to be solved compared to CMG, SUM, LSM, and GM. This would be good news for saving computational cost when using Newton-Raphson method.

Rayleigh Ritz method

The graphical results and beam-end coordinates are presented in Figs. 12 & 23 and Tables 9 &19 where Rayleigh Ritz method is denoted by RRM.

1. Deflection of the slender straight beam: as noted in Table 9, the maximum error reaches 0.20%, and the maximum computational expense is around 0.41 s. Besides, 4th-order polynomials are used to approximate the solution. 2. Deflection of the slender circular beam: as noted in Table 19, the maximum error reaches 0.24%, and the maximum computational expense is around 0.14 s. Besides, 4th-order polynomials are used to approximate the solution.

Similarly, RRM only needs 4th-order polynomials to properly approximate the solution, which means fewer unknowns and equations need to be solved compared to CMG, SUM, LSM and GM. Therefore, this would be good news for saving computational cost when using Newton-Raphson method as well.

Collocation method implemented in a local manner

The graphical results and beam-end coordinates are presented in Figs. 13 & 24 and Tables 10 & 20 where collocation method implemented in a local manner is denoted by CML.

1. Deflection of the slender straight beam: as noted in Table 10, the maximum error reaches 0.10%, and the maximum computational expense is around 0.41 s. Besides, 3 collocation points and 2 sub-intervals are defined. 2. Deflection of the slender circular beam: as noted in Table 20, the maximum error reaches 0.64%, and the maximum computational expense is around 0.36 s. Besides, 3 collocation points and 2 sub-intervals are defined.

As stated in Section 4.5, CML is proposed for complex BVPs that have solutions of highly frequent vibration. The corresponding example in beam-deflection problems refers to high-order buckling phenomena (specifically, higherorder buckling needs more sub-intervals and collocation points defined). Obviously, CML can also do its duties well in terms of accurately and efficiently solving common large beam-deflection problems according to the results shown in Figs. [START_REF] Kumar Rai | Synthesis of path generating compliant mechanisms using initially curved frame elements[END_REF] & 24 and Tables 10 &20.

Summary for the presented methods

In this section, we aim to summarize modeling beam-deflection problems using the 10 presented methods in terms of efficiency, accuracy, and feasibility respectively.

1. Efficiency: Compared to FEM using 1-D beam elements, all the 10 methods present better efficiency in modeling beam-deflection problems. Among these 10 methods, FDM is the most time-consuming one due to its modeling nature. It essentially transforms the BVP into a very high-dimensional system of equations compared to weighted residual methods, which results in more computational expense. 2. Accuracy: All the presented 10 methods have been proved reliable in modeling beam-deflection problems by FEM. In particular, collocation methods will lose some accuracy if the number and the property of pre-set polynomials are not properly chosen, which is explicitly reflected in the results of CMG and CML (see Table 4 & Table 14 and Table 10 & Table 20). However, the above mentioned is all about numerical errors, and we always have means to reduce them. For example, we can always increase the order of the pre-set polynomials defined in all weighted residual methods and decrease the predefined threshold of Newton Raphson method to reduce numerical errors. In terms of the used constitutive model, we here take one step further to analyze the accuracy compared to actual beam deflection in reality. As we have stated, all the presented methods are based on Euler-Bernoulli beam theory (Eq. ( 1)) where some assumptions are made. This constitutive model ignores the shear deformation on cross sections of the beam and assumes that the studied beam is inextensible. To handle the above 2 mentioned problems, Timoshenko beam theory needs to be used to consider the shear deformation of beam cross sections, and Eq. ( 1) also needs to be modified to consider the axial extension or compression. 3. Feasibility: As stated before, NIM can not handle beam-deflection problems where inflection points are involved, and LSM often suffers from falling into local minimal solving relatively complex and large beamdeflection problems, such as circular beams. Apart from NIM and LSM, the rest proposed 8 methods present full feasibility in terms of solving large beam-deflection problems.

Mechanism synthesis via the 10 proposed methods

In this section, we aim to model 2 representative compliant mechanisms to prove the feasibility of the proposed 10 methods: straight-beam-based cross-axis compliant revolute joint [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] (see Fig. 25) and circular-beam-based compliant revolute joint (see Fig. 36). To synthesize the 2 mentioned compliant mechanisms, we have to follow the 3 modeling principles of mechanics [START_REF] Pieter | Advanced strength of materials[END_REF] The above 3 relationships determine the essence of modeling deformable bodies regardless of their linear or nonlinear properties. Essentially, we need to solve a nonlinear system of equations that are composed of the above three types of relationships. Therefore, we have two strategies to implement the programming process. The first one is packing the constitutive equations into a function, and calling the function in the iterative process of Newton-Raphson method. This strategy turns out to be straightforward in programming but strongly complex in calculation. The second strategy is directly coding the constitutive equations in the entire nonlinear system of equations and solving it using Newton-Raphson method. The latter has been proved a little more complex in coding but more efficient in computation. Finally, the modeling results are compared and verified by FEM results (where solid mesh elements are used) with corresponding errors presented accordingly.

Straight-beam-based cross-axis compliant revolute joint

In this section, a straight-beam-based cross-axis compliant revolute joint is modeled where the presented 10 methods are used to handle the constitutive beam model (Euler-Bernoulli beam theory) one by one. As shown in Fig. 25, the revolute joint is composed of two slender straight beams (beam #1 and beam #2) and a rigid link. Each of the are the coordinates of A, A * , B, B * , K and K * . Then, we can numerically solve Eq. ( 38) to (41) using Newton-Raphson method. Besides, Δ Δ Φ refers to the translational and rotational displacements of K, the reference point of the modeled revolute mechanism. Particularly, the 10 proposed methods in Section 4 are used respectively to handle the constitutive equations ( 38) and [START_REF] Larry | Handbook of compliant mechanisms[END_REF], therefore ending up with 10 final modeling results to be compared and analyzed in the next section.

Numerical results

Modeling the studied revolute mechanism is all about finding out the relationships between applied loads and the corresponding displacements Δ Δ Φ both at the reference point K. The material properties and geometry information are provided in the following: The results of and Δ Δ Φ are provided via using the 10 proposed methods in Section 4 respectively, followed by FEM verification (see Fig. 26 to Fig. 35). Obviously, every method proposed in Section 4 presents high accuracy compared to solid-mechanics-based FEM with errors all below 0.5%. In all these numerical methods, the next-step initial guess is set up as the last-step answer for faster convergence [START_REF] Singiresu | The finite element method in engineering[END_REF]. Particularly, in Fig. 28b, no inflection points exist in the 2 deformed beams of the compliant revolute joint so NIM can still work in this case. The graphical results by FEM are shown in Fig. 37a where the fabricated and the 10 -deformed mechanism are presented respectively.

Circular-beam-based compliant revolute joint

In this section, a circular-beam-based compliant revolute joint is studied where the presented 10 methods are used to handle the constitutive beam model (Euler-Bernoulli beam theory) accordingly. The revolute joint is composed of two slender circular beams (beam #1 and beam #2) and a rigid link shown in Fig. 36. Each of the beams is fixed at one end respectively (O 1 and O 2 ) and rigidly attached to a rigid link AB on the other end. Similarly, , , and ℎ refer to the length, radius, width, and thickness of the 2 involved elementary beams respectively. The global coordinate frame X-O-Y is established at the mid-point of O 1 and O 2 . A combined load composed of , and is applied at the reference point K (the mid-point of AB), resulting in a deformed shape of the studied revolute mechanism shown in Fig. 36.

, and are the forces and moment exerted at the tip A of beam #1 by the rigid link. , and are the forces and moment exerted at the tip B of beam #2 by the rigid link. Besides, A * , B * and K * refer to the points after deformation accordingly. Φ is the rotational angle of the revolute joint. 2 local coordinate systems are set up for each flexible beam (see X 1 -O 1 -Y 1 and X 2 -O 2 -Y 2 in Fig. 36). Note that we assume O 1 O 2 =AB = for this revolute joint.

, and are the coordinates of A, A * , B, B * , K and K * . Then, we can numerically solve Eq. ( 42) to (45) using Newton-Raphson method. Besides, Δ Δ Φ refers to the translational and rotational displacements of K, the reference point of the modeled revolute mechanism. In particular, 9 out of the 10 proposed methods (except NIM) in Section 4 are used respectively to handle the constitutive equations ( 42) and ( 43), therefore ending up with 9 final modeling results to be compared and analyzed in the next section as well. In Section 5.3.3, we have found out that it's easy to generate inflection points in deformation of circular beams. Therefore, we don't use NIM to solve the constitutive equations (42)(43) of the slender circular beams.

Numerical results

Modeling the revolute mechanism is all about studying the relationships between the applied loads and the resulting displacements. The material properties and geometry information are provided in the following:

Figure 1 :

 1 Figure 1: General beam-tip loading condition of a thin straight beam

with 1

 1 are utilized. Besides, for ODE[START_REF] Gräser | On the influence of the flexure hinge orientation in planar compliant mechanisms for ultra-precision applications[END_REF] we have two boundary conditions: the first one is (0) = 0 implying that the tangent angle at = 0 is zero since it is fixed, and the second condition is d being the constant initial curvature of the studied flexible beam (circular beams) which derives from Euler-Bernoulli beam assumption at = . Summarizing what we have stated above, the beam-deflection problem is essentially a boundary value problem (BVP) of ODE (

Figure 2 :

 2 Figure 2: Nonlinear buckling of a slender straight beam

Theorem 1 (

 1 Weierstrass's first theorem). If ( ) is a given continuous function for ≤ ≤ where and represent the lower bound and upper bound of , and if is an arbitrary positive quantity, then there exists an approximating polynomial ( ) = ∑ =0 with ∈ ℝ and ∈ ℤ such that | ( ) -( )| < for ≤ ≤ .

  in the th iteration. Then, we can define IVP[START_REF] Sam Han | Design, simulation, and fabrication of a quadstable monolithic mechanism with x-and y-directional bistable curved beams[END_REF] as a function ( ) where obviously is the input and d d ( ) is the output of the function. Therefore, what we actually aim at solving is:

= 200 × 10 9

 9 Pa; = 0.3; = 0.01 m; ℎ = 0.004 m; = 0.2 m;In terms of its beam-tip loading, it is shown as follows:

= 200 × 10 9 Figure 3 :Figure 4 :

 934 Figure 3: The geometry of the studied slender straight beam and its general beam-tip loading

Figure 7 :

 7 Figure 7: Graphical results by CMG

Figure 13 :

 13 Figure 13: Graphical results by CML

Figure 14 :

 14 Figure 14: The geometry of the studied slender circular beam and its general beam-tip loading

Figure 19 :

 19 Figure 19: Graphical results by SUM

Figure 23 :

 23 Figure 23: Graphical results by RRM

  [START_REF] Robert R Archer | An introduction to mechanics of solids[END_REF][51]: a) Constitutive relationships b) Force equilibrium or force compatibility relationships c) Geometric equilibrium or geometric compatibility relationships

Figure 25 :

 25 Figure 25: Constitutive, force equilibrium and compatibility relationships of a straight-beam-based cross-axis compliant revolute joint and its schematic geometry

=

  200 × 10 10 Pa; = 0.01 m; ℎ = 0.004 m; = 0.2

=

  200 × 10 10 Pa; = 0.01 m; ℎ = 0.004 m; = 0.1 m; = ; = ℎ 3 12 ; = 0.25 m;

Table 1 :

 1 Beam-end coordinates by FDM

	To both verify the accuracy and efficiency of these numerical methods, we use beam elements in FEM where the
	element information is presented here:

Edge element: 44; Mesh vertices: 45

Since the shear strain is not considered in Euler-Bernoulli beam model, is logically not involved in our numerical algorithms.

Table 2 :

 2 Beam-end coordinates by SH

				The deformed beam shape by NIM
				The deformed beam shape by FEM
		0.15						Case	Case 1 Case 2	Case 3	Case 4
								(N)	1000	500	-500	-1000
		Case 1		Case 2			Loading	(N)	1000	500	-500	-1000
		0.1						(N.m)	60	30	-30	-60
								( ) (m) 0.01491 0.09098	0.1739	0.1556
							NIM	( ) (m)	0.1598	0.1527 -0.08810 -0.1113
	Y (m)	0.05					FEM	( ) (rad) 2.3171 ( ) (m) 0.01492 0.09105 1.6512 ( ) (m) 0.1599 0.1527	-0.8061 -1.1085 0.1739 0.1556 -0.8817 -0.1114
				Undeformed shape		( ) (rad) 2.3185	1.6512	-0.8066 -1.1092
		0
								( )	0.07%	0.08%	0.00%	0.00%
							ER	( )	0.06%	0.00%	0.08%	0.09%
								( )	0.06%	0.00%	0.06%	0.06%
		-0.05			Case 3		CE (s)	NIM FEM	0.30 13.00	0.17 12.00	0.13 7.00	0.17 8.00
		-0.1		Case 4			Note that Runge-Kutta method is used.
		0	0.05	0.1	0.15	0.2	0.25
				X (m)
	Figure 6: Graphical results by NIM

Table 3 :

 3 Beam-end coordinates by NIM

Table 4 :

 4 Beam-end coordinates by CMG

				The deformed beam shape by SUM
				The deformed beam shape by FEM	Case	Case 1 Case 2	Case 3	Case 4
		0.15						(N)	1000	500	-500	-1000
							Loading	(N)	1000	500	-500	-1000
		Case 1		Case 2				(N.m)	60	30	-30	-60
		0.1						( ) (m) 0.01491 0.09098	0.1739	0.1556
							SUM	( ) (m)	0.1598	0.1527 -0.08810 -0.1113
								( ) (rad) 2.3175	1.6513	-0.8061 -1.1084
	Y (m)	0.05		Undeformed shape	FEM	( ) (m) 0.01492 0.09105 ( ) (m) 0.1599 0.1527 ( ) (rad) 2.3185 1.6512	0.1739 -0.8817 -0.1114 0.1556 -0.8066 -1.1092
		0						( )	0.07%	0.08%	0.00%	0.00%
							ER	( )	0.06%	0.00%	0.08%	0.09%
								( )	0.04%	0.00%	0.06%	0.07%
		-0.05			Case 3		CE (s)	SUM FEM	0.30 13.00	0.25 12.00	0.24 7.00	0.37 8.00
							Note that 6th-order polynomials are used to approximate the
		-0.1		Case 4			solution.
		0	0.05	0.1	0.15	0.2	0.25
				X (m)
	Figure 8: Graphical results by SUM

Table 5 :

 5 Beam-end coordinates by SUM

				The deformed beam shape by LSM
				The deformed beam shape by FEM	Case	Case 1 Case 2	Case 3	Case 4
		0.15						(N)	1000	500	-500	-1000
							Loading	(N)	1000	500	-500	-1000
		Case 1		Case 2				(N.m)	60	30	-30	-60
		0.1						( ) (m) 0.01491 0.09098	0.1739	0.1556
							LSM	( ) (m)	0.1598	0.1527 -0.08810 -0.1113
								( ) (rad) 2.3177	1.6512	-0.8061 -1.1085
	Y (m)	0.05		Undeformed shape	FEM	( ) (m) 0.01492 0.09105 ( ) (m) 0.1599 0.1527 ( ) (rad) 2.3185 1.6512	0.1739 -0.8817 -0.1114 0.1556 -0.8066 -1.1092
		0						( )	0.07%	0.08%	0.00%	0.00%
							ER	( )	0.06%	0.00%	0.08%	0.09%
								( )	0.03%	0.00%	0.06%	0.06%
		-0.05			Case 3		CE (s)	LSM FEM	0.46 13.00	0.40 12.00	0.39 7.00	0.50 8.00
							Note that 6th-order polynomials are used to approximate the
		-0.1		Case 4			solution.
		0	0.05	0.1	0.15	0.2	0.25
				X (m)
	Figure 9: Graphical results by LSM

verified by FEM. In particular, computational expenses (denoted by CE) and errors (denoted by ER) with respect to FEM are also presented accordingly in Tables 1 to 10 and 11 to 20.

Table 6 :

 6 Beam-end coordinates by LSM

				The deformed beam shape by GM
				The deformed beam shape by FEM	Case	Case 1 Case 2	Case 3	Case 4
		0.15						(N)	1000	500	-500	-1000
							Loading	(N)	1000	500	-500	-1000
		Case 1		Case 2				(N.m)	60	30	-30	-60
		0.1						( ) (m) 0.01495 0.09099	0.1739	0.1556
							GM	( ) (m)	0.1598	0.1527 -0.08810 -0.1113
								( ) (rad) 2.3148	1.6502	-0.8062 -1.1092
	Y (m)	0.05		Undeformed shape	FEM	( ) (m) 0.01492 0.09105 ( ) (m) 0.1599 0.1527 ( ) (rad) 2.3185 1.6512	0.1739 -0.8817 -0.1114 0.1556 -0.8066 -1.1092
		0						( )	0.20%	0.07%	0.00%	0.00%
							ER	( )	0.06%	0.00%	0.08%	0.09%
								( )	0.16%	0.06%	0.05%	0.00%
		-0.05			Case 3		CE (s)	GM FEM	0.22 13.00	0.19 12.00	0.19 7.00	0.31 8.00
		-0.1		Case 4
		0	0.05	0.1	0.15	0.2	0.25
				X (m)
	Figure 10: Graphical results by GM

Table 8 :

 8 Beam-end coordinates by MOM

				The deformed beam shape by RRM
				The deformed beam shape by FEM	Case	Case 1 Case 2	Case 3	Case 4
		0.15						(N)	1000	500	-500	-1000
							Loading	(N)	1000	500	-500	-1000
		Case 1		Case 2				(N.m)	60	30	-30	-60
		0.1						( ) (m) 0.01495 0.09099	0.1739	0.1556
							RRM	( ) (m)	0.1598	0.1527 -0.08810 -0.1113
								( ) (rad) 2.3148	1.6502	-0.8062 -1.1092
	Y (m)	0.05		Undeformed shape	FEM	( ) (m) 0.01492 0.09105 ( ) (m) 0.1599 0.1527 ( ) (rad) 2.3185 1.6512	0.1739 -0.8817 -0.1114 0.1556 -0.8066 -1.1092
		0						( )	0.20%	0.07%	0.00%	0.00%
							ER	( )	0.06%	0.00%	0.08%	0.09%
								( )	0.16%	0.06%	0.05%	0.00%
		-0.05			Case 3		CE (s)	RRM FEM	0.41 13.00	0.31 12.00	0.30 7.00	0.31 8.00
							Note that 4th-order polynomials are used to approximate the
		-0.1		Case 4			solutions.
		0	0.05	0.1	0.15	0.2	0.25
				X (m)
	Figure 12: Graphical results by RRM

Table 9 :

 9 Beam-end coordinates by RRM

Table 10 :

 10 Beam-end coordinates by CML

Table 11 :

 11 Beam-end coordinates by FDM

					The deformed beam shape by SH
		0.3			The deformed beam shape by FEM		Case	Case 1	Case 2 Case 3	Case 4
										(N)	600	-150	-300	-600
		0.25	Case 4		Case 3			Loading	(N)	600	-150	-300	-600
										(N.m)	20	-5	-10	-20
		0.2							SH	( ) (m) -0.07488 0.1280 ( ) (m) 0.1882 0.2091 0.08049 -0.05971 0.2863 0.3045
						Case 2			( ) (rad)	3.4016	2.3505	1.0960	0.2706
	X (m)	0.1 0.15				Undeformed shape	FEM	( ) (m) -0.07487 0.1278 ( ) (m) 0.1882 0.2092 0.08060 -0.05969 0.2862 0.3045 ( ) (rad) 3.4032 2.3518 1.0967 0.2710
										( )	0.01%	0.16%	0.03%	0.00%
		0.05			Case 1			ER	( ) ( )	0.00% 0.05%	0.05% 0.06%	0.14% 0.06%	0.03% 0.15%
		0							CE (s)	SH FEM	1.20 12.00	0.80 12.00	0.77 28.00	1.07 67.00
		-0.05
			-0.05	0	0.05	0.1	0.15	0.2
					Y (m)
	Figure 16: Graphical results by SH

Table 12

 12

		0.3	Case 4 Box 1		The deformed beam shape by NIM The deformed beam shape by FEM	Case	Case 1	Case 2 Case 3	Case 4
		0.25			Case 3			Loading	(N) (N)	600 600	-150 -150	-300 -300	-600 -600
									(N.m)	20	-5	-10	-20
		0.2							( ) (m) -0.07488 0.1280	NA	NA
	X (m)	0.15	Box 2		Case 2 Undeformed shape	NIM FEM	( ) (m) ( ) (rad) ( ) (m) -0.07487 0.1278 0.1883 0.2091 3.4028 2.3506 ( ) (m) 0.1882 0.2092 0.08060 -0.05969 NA NA NA NA 0.2862 0.3045
		0.1							( ) (rad)	3.4032	2.3518	1.0967	0.2710
									( )	0.01%	0.16%	NA	NA
		0.05							ER	( )	0.05%	0.05%	NA	NA
					Case 1			( )	0.01%	0.05%	NA	NA
		0							CE (s)	NIM FEM	0.22 12.00	0.35 12.00	NA 28.00	NA 67.00
		-0.05							Note that Runge-Kutta method is used.
			-0.05	0	0.05	0.1	0.15	0.2
					Y (m)
	Figure 17: Graphical results by NIM

, the maximum error reaches 0.16%, and the maximum computational expense is around 1.20 s. Besides, Runge-Kutta method is used in each iteration, and

Table 13 :

 13 Beam-end coordinates by NIM

					The deformed beam shape by CMG
		0.3			The deformed beam shape by FEM		Case	Case 1	Case 2 Case 3	Case 4
										(N)	600	-150	-300	-600
		0.25	Case 4		Case 3			Loading	(N) (N.m)	600 20	-150 -5	-300 -10	-600 -20
										( ) (m) -0.07391 0.1269	0.2865	0.3044
		0.2							CMG	( ) (m)	0.1902	0.2093 0.08022 -0.06017
						Case 2			( ) (rad)	3.3790	2.3584	1.0965	0.2665
	X (m)	0.15				Undeformed shape	FEM	( ) (m) -0.07487 0.1278 ( ) (m) 0.1882 0.2092 0.08060 -0.05969 0.2862 0.3045 ( ) (rad) 3.4032 2.3518 1.0967 0.2710
		0.1
										( )	1.28%	0.70%	0.10%	0.03%
									ER	( )	1.06%	0.05%	0.47%	0.80%
		0.05								( )	0.71%	0.28%	0.02%	1.66%
		0			Case 1			CE (s)	CMG FEM	0.11 12.00	0.08 12.00	0.12 28.00	0.23 67.00
									Note that 6th-order polynomials are used to approximate the
		-0.05							solutions.
			-0.05	0	0.05	0.1	0.15	0.2
					Y (m)
	Figure 18: Graphical results by CMG

Table 14 :

 14 Beam-end coordinates by CMG

Table 15 :

 15 Beam-end coordinates by SUM

					The deformed beam shape by LSM
		0.3			The deformed beam shape by FEM	Case	Case 1	Case 2 Case 3	Case 4
									(N)	600	-150	-300	-600
		0.25	Case 4		Case 3		Loading	(N) (N.m)	600 20	-150 -5	-300 -10	-600 -20
									( ) (m) -0.07486 0.1280	0.2863	0.3045
		0.2						LSM	( ) (m)	0.1883	0.2091 0.08050 -0.05971
						Case 2		( ) (rad)	3.4027	2.3506	1.0967	0.2706
	X (m)	0.15			Undeformed shape	FEM	( ) (m) -0.07487 0.1278 ( ) (m) 0.1882 0.2092 0.08060 -0.05969 0.2862 0.3045 ( ) (rad) 3.4032 2.3518 1.0967 0.2710
		0.1							( )	0.01%	0.16%	0.03%	0.00%
								ER	( )	0.05%	0.05%	0.12%	0.03%
									( )	0.01%	0.05%	0.00%	0.15%
		0.05			Case 1			CE (s)	LSM FEM	0.45 12.00	0.14 12.00	0.41 28.00	0.50 67.00
		0
		-0.05
			-0.05	0	0.05	0.1	0.15	0.2
					Y (m)
	Figure 20: Graphical results by LSM

Table 17 :

 17 Beam-end coordinates by GM

					The deformed beam shape by MOM
		0.3			The deformed beam shape by FEM		Case	Case 1	Case 2 Case 3	Case 4
										(N)	600	-150	-300	-600
		0.25	Case 4		Case 3			Loading	(N) (N.m)	600 20	-150 -5	-300 -10	-600 -20
										( ) (m) -0.07473 0.1280	0.2863	0.3045
		0.2							MOM	( ) (m)	0.1884	0.2091 0.08053 -0.05970
						Case 2			( ) (rad)	3.4025	2.3507	1.0961	0.2706
	X (m)	0.15				Undeformed shape	FEM	( ) (m) -0.07487 0.1278 ( ) (m) 0.1882 0.2092 0.08060 -0.05969 0.2862 0.3045 ( ) (rad) 3.4032 2.3518 1.0967 0.2710
		0.1
										( )	0.19%	0.16%	0.03%	0.00%
									ER	( )	0.11%	0.05%	0.09%	0.02%
		0.05								( )	0.02%	0.05%	0.05%	0.15%
		0			Case 1			CE (s)	MOM FEM	0.14 12.00	0.10 12.00	0.11 28.00	0.13 67.00
		-0.05
			-0.05	0	0.05	0.1	0.15	0.2
					Y (m)
	Figure 22: Graphical results by MOM

Table 19 :

 19 Beam-end coordinates by RRMNote that 3 collocation points and 2 sub-intervals are defined.

					The deformed beam shape by CML
		0.3			The deformed beam shape by FEM		Case	Case 1	Case 2 Case 3	Case 4
		0.25	Case 4		Case 3			Loading	(N) (N)	600 600	-150 -150	-300 -300	-600 -600
										(N.m)	20	-5	-10	-20
		0.2								( ) (m) -0.07488 0.1281	0.2863	0.3044
						Case 2		CML	( ) (m)	0.1881	0.2091 0.08053 -0.05931
	X (m)	0.15				Undeformed shape	FEM	( ) (rad) ( ) (m) -0.07487 0.1278 3.4028 2.3503 ( ) (m) 0.1882 0.2092 0.08060 -0.05969 1.0959 0.2706 0.2862 0.3045
		0.1								( ) (rad)	3.4032	2.3518	1.0967	0.2710
										( )	0.01%	0.23%	0.03%	0.03%
		0.05							ER	( )	0.05%	0.05%	0.09%	0.64%
					Case 1				( )	0.01%	0.06%	0.07%	0.15%
		0							CE (s)	CML FEM	0.36 12.00	0.11 12.00	0.10 28.00	0.24 67.00
		-0.05
			-0.05	0	0.05	0.1	0.15	0.2
					Y (m)
	Figure 24: Graphical results by CML

Table 20 :

 20 Beam-end coordinates by CML maximum computational expense is around 0.35 s. In Case 1 and Case 2, there is no inflection points in the exists. As shown in Fig.17, it is not hard to observe one inflection point in Box 1 and at least one inflection point in Box 2. In the end, Runge-Kutta method is used to approximate the solution of the studied IVP[START_REF] Su | A pseudorigid-body 3r model for determining large deflection of cantilever beams subject to tip loads[END_REF].As stated in Section 4.3, NIM can not handle beam-deflection problems where inflection points are involved. According to the graphical results shown in Figs.6 & 17, it is obvious that the largely-deformed circular beams tend to produce inflection points. Besides, we always need to determine the sign of d d to calculate the beam-end angle ( defined in Section 4.3) first before approximating the solution of the studied IVP

	deformed beam shapes, and that NIM can not find the solutions of Case 3 and Case 4 since the sign of d keeps positive along [0, ] so the positive sign + is chosen for Eq. (13). Note d d changes along [0, ] where at d least one inflection point ( d d = 0)

2) Force equilibrium 3) Geometric compatibility relationships The applied loads are set up as:

The results of and Δ Δ Φ are provided via using 9 proposed methods in Section 4 respectively (except NIM), followed by FEM verification (see Fig. 38 to Fig. 46). It is noticed that all methods used present high accuracy compared to solid-mechanics-based FEM with maximum errors around 2%. Likewise, in all these numerical methods, the next-step initial guess is set up as the last-step answer for faster convergence [START_REF] Singiresu | The finite element method in engineering[END_REF]. The graphical results by FEM are shown in Fig. 37b where the fabricated and the 10 -deformed mechanism are presented respectively.

Error analysis

This section aims to analyze the origin of errors when using the proposed 10 methods. According to Figs. 26 to 35, the errors are all below 0.5% using the proposed 10 methods to solve straight-beam-deflection problems. The maximum errors in terms of modeling circular-beam-based revolute joints are around 2%. It's obvious to conclude that all these used methods have proved to be effective modeling CMs. However, we still need to analyze where these errors are generated and how we can reduce them. To start with, the reason why we use solid-mechanics-based FEM as the reference method is that it considers all sorts of deformation in any deformable object.

Numerical errors

The first type of errors is called numerical errors, which comes from the essence of numerical methods. The aim of numerical methods is to approximate the solution to the studied mathematical models so they would inevitably lead to some numerical errors. However, we have many ways to decrease these errors. For instance, we can increase the order of the pre-set polynomials defined in all weighted residual methods; we can decrease the predefined threshold of Newton Raphson method's convergence criteria in all weighted residual methods, NIM, SH, and NIM.

Model errors

As long as the numerical errors have been decreased to their minimal, what is left is the model errors due to the non-perfect constitutive model. In our paper, Euler Bernoulli beam theory (Eq. ( 1)) is adopted which is based on some assumptions: ignoring the shear deformation on cross sections of the beam and assuming that the studied beam is inextensible. Logically, model errors always exist compared to solid-mechanics-based FEM since it takes care of all sorts of deformation in any deformable object. 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.

Conclusions

Compliant mechanisms have been serving as a promising option in designing mechanical systems. The constitutive equations of the elementary flexible members, such as slender straight beams and slender circular beams are studied in this paper. The planar beam-deflection problems are formulated as BVPs of corresponding ODEs using the most classic Euler-Bernoulli beam theory. Then, the 10 proposed numerical methods are used respectively to solve planar large deflection of slender straight beams and slender circular beams with promising results and detailed analysis provided. These methods have also been proved to be effective in terms of synthesizing CMs where 2 representative compliant revolute joints are modeled. Our future work will focus on modeling spatial deflection of slender beams for studying spatial CMs.

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.

Abbreviations

Nomenclature

Young's Modulus The actuation force along Y axis exerted at the reference point K of the studied mechanism

The payload along X axis exerted at the reference point K of the studied mechanism

The force exerted at beam #1 end A along Y axis of its local coordinate system X 1 -O 1 -Y 1

The force exerted at beam #1 end A along X axis of its local coordinate system X 1 -O 1 -Y 1

The moment exerted at beam #1 end of its local coordinate system X 1 -O 1 -Y 1

The force exerted at beam #2 end B along Y axis of its local coordinate system X 2 -O 2 -Y 2

The force exerted at beam #2 end B along X axis of its local coordinate system X 2 -O 2 -Y 2

The moment exerted at beam #2 end of its local coordinate system The loading exerted at the reference point K of the studied mechanism