Path Analysis for Hybrid Rigid–Flexible Mechanisms

: Hybrid rigid–ﬂexible mechanisms are a type of compliant mechanism that combines rigid and ﬂexible elements, being that their mobility is due to rigid-body joints and the relative ﬂexibility of bendable rods. Two of the modeling methods of ﬂexible rods are the Cosserat rod model and its simpliﬁcation, the Kirchhoff rod model. Both of them present a system of differential equations that must be solved in conjunction with the boundary constraints of the rod, leading to a boundary value problem (BVP). In this work, two methods to solve this BVP are applied to analyze the inﬂuence of external loads in the movement of hybrid compliant mechanisms. First, a shooting method (SM) is used to integrate directly the shape of the ﬂexible rod and the forces that appear in it. Then, an integration with elliptic integrals (EI) is carried out to solve the workspace of the compliant element, considering its buckling mode. Applying both methods, an algorithm that obtains the locus of all possible trajectories of the mechanism’s coupler point, and detects the buckling mode change, is developed. This algorithm also allows calculating all possible circuits of the mechanism. Thus, the performance of this method within the path analysis of mechanisms is demonstrated.


Introduction
A compliant body is one whose motion depends on its geometry, its material properties and the location and magnitude of the applied forces. If a body of this kind belongs to a mechanism, it is known as a compliant element. When a mechanism is fully composed of compliant elements, e.g., slender rods, it is named a compliant mechanism, and it gains all or part of its mobility, thanks to the relative flexibility of those compliant elements. On another note, if the mechanism combines rigid and compliant elements, it takes the name of a hybrid compliant or hybrid rigid-flexible mechanism [1].
Regarding applications, such compliant-and, especially, hybrid-compliant-mechanisms are an alternative to robotic systems of rigid bodies [2], or even cable robots [3] used in tasks where there is a human-machine interface. Additionally, such hybrid-flexible mechanisms are a plausible alternative to sub-systems subjected to impacts, such as four-link legs in bio-inspired mobile multipod robots [4]. Finally, its elasticity, if appropriately tuned, can play a beneficial role in the dynamic characteristics of systems used in high-speed link motions, reducing the need for balancing [5].
For the complete kinematic characterization of a linkage with rigid bars, all the circuits of the mechanism should be obtained [6]. Each circuit is the set of all possible orientations of links that can be calculated without disconnecting any of the joints [7]. If the linkage needs to be disassembled to move from one position to another, these positions lie on different circuits. These mechanism's circuits do not depend on the input link chosen. In [8,9] is presented the circuit analysis for the Watt and Stephenson -like six-bar mechanisms. It should be noted that, in the case of compliant mechanisms, it is not needed to disassemble any link to obtain the different circuits but to deform a flexible rod through an external load.
To solve the movement of these hybrid rigid-flexible mechanisms, it is needed to obtain the change of shape of compliant elements and the loads that these bodies withstand.
The derivation of the mathematics structure of this section draw on Antman's work [13] (chs. 4 and 8) on nonlinear elasticity. Nevertheless, the nomenclature used here are those used in the robotics community, in the line of the work of Caleb [22].
The framework needed to describe the shape of a slender rod in space includes three main aspects that, coupled, generate the nonlinear system of differential equations that has to be solved in order to obtain the relationship between force and deformation. These are as follows: a kinematic definition of the deformed rod, material constitutive laws and static equilibrium equations.

Kinematics
A three-dimensional parametric Cartesian curve p(s) ∈ R 3 that links the centroids of each transverse section and the orthonormal rotation matrix R(s) ∈ SO(3) (SO(3) is the special orthogonal subgroup in three dimensions, SO(3) = {R ∈ R 3x3 | R T R = I, and det(R) = 1}) that orients a local frame attached to the section are used to define the deformed shape of a rod. The principal axes of each section are usually named the xand y-axes. The z-axis is perpendicular to the cross-section and tangent to the deformed shape; see Figure 1. A scalar reference arc-length parameter s is used to position and orient each cross-section. This parameter is within the finite interval s ∈ [0, L], being that L is the length of the rod in the initial state. Due to the assumption that the sections are not distorted, the whole deformation of the flexible bar can be described by mapping from s to a homogeneous rigid-body transformation, T(s) ∈ SE(3) (SE (3) is the special Euclidean subgroup in three dimensions, with p(s) ∈ R 3 and R(s) ∈ SO (3)).
where we will often refer to T(s) as a "frame" hereafter. For the unloaded pre-curved state of the rod, the magnitudes described above are designated with the subscript • .
A priori, R • (s) could be assigned arbitrarily. However, one can establish conventions governing the assignment of reference orientations such that the mapping form R • (s) to R(s) has an easily interpretable meaning in terms of material strains. As mentioned above, we have chosen to assign reference orientation such that the z axis of the local frame is parallel to the tangent curve to the reference curve (see Figure 1). Then, the following applies: where e 1 , e 2 and e 3 are used for the standard basis vectors of the local frame: [1 0 0] T , [0 1 0] T and [0 0 1] T , respectively. The superscript is used to indicate differentiation with respect to the arc length s. Now, let us consider the derivatives of p(s) and R(s) with respect to s. The position and orientation evolve along the arc length according to the rates of change, linear v(s) ∈ R 3 and angular u(s) ∈ R 3 . These rates of change are defined in the local frame, and obtained from the above mentioned derivatives, defined in the global frame, with the use of R(s) T as follows: where U(s) is a skew-symmetric matrix of the following form: Because U(s) is defined through three independent values, u x (s), u y (s) and u z (s), it can be expressed with a vector u(s) = u x (s) u y (s) u z (s) T such that u = U ∨ and U = u, where the ∨ operator denotes conversion of an element of so(3) (so (3) is the Lie algebra of SO(3)) to its corresponding element in R 3 . The inverse operation, denoted by , maps R 3 to so(3), so that u ∨ = u. In the following, u(s) is used whenever a vector form is needed to express U.
When v(s) and u(s) are those of the reference, i.e., the unloaded shape of the rod, they are designated as v • and u • . Because of the convention of the rotation matrix stated in (2), an initially straight rod has v • = [0 0 1] T and u • = 0.
The material strain can be defined by the variation of v(s), ∆v(s) = v(s) − v • (s), and u(s), ∆u(s) = u(s) − u • (s). The local elongation is represented by ∆v z ; values above 1 correspond to extension, values below 1 to compression and values equal to 1 mean that the length of the curve has not changed. The shear deformation along the local axes is depicted by ∆v x and ∆v y . The x and y components of ∆u(s), ∆u x and ∆u y , measure the bending along the local xand y-axes. Concerning ∆u z , it collects the value of the torsion.

Material Constitutive Relationships
The way the material behaves (it can be elastic or plastic, for instance) affects the constitutive relationships between strain and stress. However, the Cosserat rod model only is valid if the cross-sections do not deform. ∆v(s) and ∆u(s) can be used to define the stress field since it is obtained from the strain field. Additionally, ∆v(s) and ∆u(s) can also be used to express internal force n(s) and moment m(s), due to the fact that they are calculated from the stress field.
Hereafter, the linear deformation v is related to internal forces by a linear constitutive law, which is also used to relate the angular deformation u to the internal moment. Diagonal stiffness matrices, i.e., K SE (s) and K BT (s), allow expressing the aforementioned magnitudes in the local frame. It is worth noting that the subscript SE indicates shear and extension, and BT, bending and torsion. n loc (s) = K SE (s)∆v(s) (5) m loc (s) = K BT (s)∆u(s) (6) with where A is the area of cross section, E is Young's modulus, G is the shear modulus, I x and I y are the second moment of area of the rod cross section about the principal axes x and y, and J z is the polar moment about the local axis z. The linear constitutive law applied here is commonly used in mechanics for the small strain range. This happens for most metals and for polymers in the range of small strain. The rotation matrix is used to change n(s) and m(s) from the local axes to the fixed global axes.

Static Equilibrium Equations
The forces and moments on an infinitesimal piece of the rod, [s, s + ds], are considered to obtain the static equilibrium equations. The internal force and moment exerted on the piece by the material of [0, s) are named −n(s) and −m(s). On the other side of the portion, the material of (s + ds, L] exerts a force n(s + ds) and a moment m(s + ds). In addition, distributed forces f(s) and moments l(s) can also be applied on the interval ds.
On the one hand, the equilibrium of forces yields the following: −n(s) + n(s + ds) + f(s)ds = 0 Expanding the magnitudes in s + ds with the Taylor series until the element of first order, simplifying terms, and dividing by ds gives the following: On the other hand, a balance of moments, expanding variables in s + ds with the Taylor series until the element of first order, simplifying terms, and taking into account (12), yields the following: So, the nonlinear ordinary differential equations for the equilibrium of a Cosserat rod describing the evolution of the internal force n(s) and moment m(s) along arc length are as follows: m (s) + p (s) × n(s) + l(s) = 0

Static Model
If the kinematic, constitutive, and equilibrium laws are combined, a static model of the flexible rod can be achieved.
Thus, the Cosserat static model for a flexible rod can be obtained using m and n as state variables, and combining (3), (14) and (15) to obtain a set of differential equations. If it is considered that the rod is deformed only by bending and torsion, the Cosserat rod model is simplified, leading to the Kirchhoff model. In the reference rod's configuration, shear and extension deformations do not usually appear, i.e., v = v • = e 3 , where e 3 is a unit vector tangent to the deformed shape. In this way, a new system of differential equations can be stated:

Kirchhoff Static Model Equations for the Planar Case
The three-dimensional deformation of a rod in space is defined by (16). However, if the rod only deforms within a plane, the model can be further simplified and a shorter formulation is obtained. This has an impact on the computational cost when the resolution method is implemented within an algorithm.
The rod that appears in Figure 2 deforms within the plane OXY. The local frame attached to each cross-section is oriented accordingly, so note that expressions in the previous chapter vary slightly now. To orientate each cross-section the rotation matrix is used and it depends on the angle θ = θ(s): The derivative of R with respect to arc length s is as follows: and allows finding vector u in the following: Thus, only the component u z of vector u is not null: This makes sense, as there is no bending moment with respect to the local y-axis and no torsion (expressed through u x ).
Taking (20) into (16) considering this time v = v • = e 1 , and taking n and m as state variables, the expression of the system of nonlinear ordinary differential equations in (16) yields the following: If no distributed force and moment (f and l) are applied along the rod, we obtain the following: Hence, in this particular case, the internal force n is constant. Such a rod with no distributed loads is quite common in the applications.
The system of differential equations in (22) can be solved by using a direct integration method, such as the 4th order R-K, or a method based on elliptic integrals.

Direct Integration Using Fourth-Order Runge-Kutta Method
Typical problems in robotics are the forward kinematic problem (FKP) and the inverse kinematic problem (IKP). In the first one, the mechanism inputs are known, and the position and orientation of the end-effector are unknown. In contrast, in the IKP, the end-effector is wanted to describe a given trajectory, so its position is known, but the needed inputs to accomplish the task must be calculated.
When a flexible rod is part of a mechanism, usually, it is under external loads coming from that mechanism. Here, the aim should be to calculate how these loads deform the rod. So, it could be considered as a kind of forward kinematic (FK) position problem because the position and orientation of the extreme in which force and moment are applied are not known.
So focusing on the FKP of a rod that is clamped at the proximal end s = 0 and subjected to a known force n ext and moment m ext at the distal end s = L (see Figure 2), the rod undergoes a deformation that depends on these boundary conditions and the geometric and mechanical characterization of the rod.
If the effect of gravity is not taken into account, there is not distributed force along the length. Therefore, the shape and the load status of the rod are defined by (22).
This takes the form of a BVP, where the independent variable is s in [0, L], and the vector of dependent variables, y, and the non-linear function, f(s, y), are the following: The boundary conditions are at s = 0 of the kinematic kind, and at s = L of the load kind. Regarding the kinematic variables, the position and orientation at s = 0 are data, and unknown at s = L. In relation to internal forces, and because no distributed force is considered, its value does not change along the length of the rod. At s = 0, and because of the static equilibrium, the internal force takes the same value that of the external load, n ext , at s = L, but in opposite direction (although the internal forces could be taken out of the integration, in this work, they are included because they are needed to solve the FKP of the whole mechanism). With respect to the internal moment, the value is a datum at s = L (in fact, m(L) = m ext ), and unknown at s = 0.
Then, y at both ends are the following: being Dirichlet conditions.
Since there are three unknowns at the initial of the range of integration s = 0, three guess values (one for each variable) are needed: guess = n x guess n y guess m z guess T . Once these values are chosen, the initial value problem is solved using the 4th R-K method.
where y is a vector that contains the unknown variables at each s, function f is given, s is the independent variable, and subscript 0 indicates values at initio.

Integration Using Elliptic Integrals
The system of differential equations in (22) can be further developed looking for an analytical integration. If only end-point load is applied (expressed in terms of magnitude R and orientation ψ as in Figure 2), we can express the internal force as the following: For constant E and I, with a stress-free reference straight rod, considering the wellknown Bernoulli-Euler law, and upon substitution of (26) and (27) into (22), we obtain the following: Therefore, the system of differential Equation (22) can be simplified to the following: Elliptic integrals is the classical mathematical tool to solve (29) because of its rapid computation. The following equations are used in this work to solve for the coordinates of an arbitrary point along the beam, and they are found in [21] (ch. 4).
where the functions F(k, φ) and E(k, φ) are the incomplete elliptic integrals of the first and second kind, respectively. Moreover, the following holds: Thus, the value of the end-tip force R in terms of k and φ can be obtained from (32) as: In addition, the bending moments along the rod are given by the following: The functions F(k, φ) and E(k, φ) depend on two parameters: • The non-dimensional parameter k It is known as the modulus of the function and it can vary between −1 and 1. In this application, k corresponds roughly but non-linearly to the magnitude of the force R.

•
The variable φ It is called the amplitude of the elliptic integral. It is measured in radians, and it varies continuously along the beam from φ 1 on the left edge to φ 2 on the right. These amplitudes, φ 1 and φ 2 , are defined by the boundary conditions of the rod extremes. The variable φ is related to the beam angle θ by the following relation: Hence, (29) is integrated to yield (30), (31), and (33), i.e., a parametric system that uniquely defines a deformed rod in terms of parameters k and ψ for known boundary conditions and a given buckling mode.

Boundary Conditions for a Clamped-Pinned Rod
The studied mechanism in this work has a clamped-pinned rod (as in Figure 2) that substitutes one of the rigid bars of the four-bar linkage.
The first extreme clamped implies a given slope θ 1 for the rod at s = 0; hence, using (35), we can state the following: A pinned end-tip implies that no moment is acting at that point. Upon application on (34), we obtain the following: This means that the end-tip is a point of null curvature, then an inflection point of the deflected curve. Additionally, the second part of the equation implies the following: being that q is the order of the Buckling Mode shape. Note that feasible values of k may depend on ψ as, for example, is the case of the limiting value of k min = ||cos ψ−θ 1 2 || that is necessary to produce real solutions in (36). This condition is indicated with the gridded area in Figure 3a. Therefore k is in the range k = (k min , 1) or k = (−1, −k min ). Provided that range for k, φ 1 will have values in [−π/2, π/2]. On the other hand, φ 2 will have discrete values π/2, 3π/2, and so on, for buckling modes 1, 2 · · · , respectively. The variable φ will vary continuously from φ 1 to φ 2 . As shown in Figure 3, if we set a duple of values to k and ψ, inside the feasible range (Figure 3a), we will obtain the deformed shape of the rod for every buckling mode chosen (Figure 3b).
From a practical point of view, the data will come from physical quantities that should be used to find those parameters k and ψ. For example, if end-tip force is given, R and ψ, along with boundary conditions, (33) should be used to find k, and then the deformed shape with (30) and (31). If the end-tip position, X P and Y P , is given for some boundary conditions, again the objective will be to find the corresponding parameters k and ψ from (30) and (31). In most occasions, this process is numerical.
The different problems that can be stated regarding the shape of a rod accomplishing given conditions may produce a multiplicity of solutions for each buckling mode considered. For clamped-pinned bars those multiple solutions do not have to be symmetrical, but refer to completely different configurations. For example, given a same end-point position, we can obtain different solutions for given buckling modes as shown in  The curves plotted in Figure 5 were calculated for the first buckling mode of the rod. To obtain the workspace for other buckling modes, it is necessary to execute the calculation again, changing φ 2 in (38), for the desired mode. For instance, Figure 6a shows the workspace for positive k and buckling mode 2 of the rod, and Figure 6b represents that for negative k and the same mode. As for the orange lines in Figures 5 and 6, they come from the calculation of the rod's end-tip position for the higher limit of the range of k and varying ψ. The green lines are calculated with the lower limit of the interval of k and varying ψ. To close the workspace area, the purple lines are plotted. They are related to the calculation with the higher limit of the interval of ψ and varying k. In contrast, the cyan lines derive from the computation with the lower limit of the interval of ψ and varying k.

Closed-Loop Hybrid Mechanisms
Once the fundamentals of elastic rods are understood, this type of element can be used in conjunction with rigid bars to create new hybrid mechanisms. These mechanisms are characterized by combining rigid and flexible bars to perform a task, such as describing a trajectory or performing a given force. In this paper, one of the mechanisms proposed in [21] (ch. 12), the four-bar linkage with one flexible clamped-pinned rod (see Figure 7), is the subject of study. The dimensional data and the description of each bar are given in Table 1. Each angle showed in Figure 7 is explained in Table 2.  Angle Description α Input angle, which positions the crank θ 1 Orientation of the rigid coupler BC θ 2 Relative angle between BC and EF segments. It has a constant value of 90º through the cycle.
The flexible segment CD is a slender rod, made of nitinol, with a circular cross-section; its properties are shown in Table 3. Table 3. Properties of the flexible rod.

General Equations for the 4-Bar Linkage
The problem to solve is the FKP since the inputs to the mechanism, i.e., the α angle and the external forces, F ext , and moment M ext , are known, but the position and orientation of the coupler BC and the shape of the flexible rod CD are unknown.
A shooting method that combines the direct integration with a 4th order R-K method (explained in Section 2.3) and a N-R scheme is used.
The resolution begins from a known position of the mechanism for an α angle of 0 degrees.
In the next step, a given α is introduced in the calculation. Thus, the crank AB rotates and the coordinates of B point can be calculated. Then, the coordinates of point C must be obtained to position the coupler BC. However, this calculus requires a deeper analysis because the rod CD deforms, so one needs to integrate (22) to obtain the position of C point.
Here, the SM is used to solve the position of the C point and, thus, that of the mechanism, for each step. This method is focused on the minimization, with the aid of the N-R scheme, of a residue function that represents all the constraints of the mechanism, and is given by (39).
where m * z,C is the bending moment at C, ∑ M z B the set of all moments about B produced by the external loads, and the superscript * indicates that the variables are not the exact solutions; instead, they are the approximate solutions obtained numerically in each step of the SM.
The first two components of (39) are restrictions of the loads kind, and state that the moments about points B and C must be null because in those points there are hinges. The other restriction is of the geometric kind and imposes that the distance between points B and C is the length of the rigid coupler BC.
To find the solution of the mechanism's position, the SM needs two inputs: the residue function (39) and a vector of guess values, guess = n x guess n y guess m z guess T , to start the integration of the rod CD. These guess values are needed to solve the BVP shown in Section 2.3. The 4th order R-K method explained in that section can be particularized for the case of the rod CD if m ext is equal to 0 since there is a hinge in the point C.
For the first α step, these guess variables can take arbitrary values. Despite this, previous calculations were made in this work to obtain guess values that improve the convergence of the SM. For the next α steps, the vector of guess values is composed by the n x n y m z T values of the solution of the previous α step.
With these two inputs to the SM, the process to obtain the position of point C and the shaped of rod CD is as follows: 1.
For each guess, apply 4th order R-K to obtain y * = x * C y * C θ * C n * x,C n * y,C m * z,C T ; see Section 2.3.

2.
Evaluate the residue (39). To do this, x * C and y * C should be used to calculate θ 1 , and θ 1 , n * x,C and n * y,C to solve ∑ M z B . m * z,C from the direct integration of the rod CD is directly used to compare its value with 0.

3.
If all terms are below a tolerance, exit: the final solution is obtained. 4.
If not, calculate variation of residue with guess: J = d residue d guess .

6.
Start the loop again.
The results are obtained when the three components of the residue are below a tolerance. With this, the position of the mechanism is obtained for the current α step, and the results are the guess values for the next α step.

Load-Dependent Paths
Applying the aforementioned SM, the trajectory of point F can be calculated, but it depends on external loads because these affect the deformed shape of the flexible rod CD. Depending on the load applied, point C takes different positions and, in consequence, the position and the orientation of the coupler BC varies. This influence of the loads in the path described by point F can be seen in Figure 8 where a set of trajectories for different load cases is plotted.
If these trajectories are plotted as a function of the moment that the external forces and moments generate in point B, it results in Figure 9. In Figure 9a, the trajectories plotted for constant external moments applied to the coupler with no external forces result in level curves. However, if external forces are applied, the relative position between points B and F varies during the movement, so the moment of external loading about B takes different values. For this reason, the trajectories produced are not flat when forces are applied. Despite this, they belong to the same three-dimensional surface that contains the previous level set in Figure 9a, as can be seen in Figure 9b.
Hence, by applying this method to solve the FKP of mechanisms of this kind, it is possible to obtain the trajectories of the end-effector points. Moreover, this method predicts the three-dimensional surface in which the paths under different load cases are.

Buckling Mode Change
Next, if the SM is combined with the integration using EI, explained in Section 2.4, a powerful algorithm that solves the FKP and shows the workspace of the flexible rod CD is obtained. In this way, it is possible to know in which buckling mode rod CD is working at each instant of the work cycle. The process should follow these steps:

1.
Integrate the deformed shape using the SM explained in Section 3.1. By doing this, the distribution of internal moments along the bar is obtained.

2.
Search for sign changes of the internal moments along the bar. Each sign change means a point of curvature change. The number of points of changes of the curvature corresponds to the buckling mode of the flexible rod. For instance, if there are two curvature changes, the rod is in the second buckling mode. It must be borne in mind that, in the case of this mechanism, there is always a point of null moment, i.e., null curvature, in the distal end of the slender rod CD because of the hinge.

3.
Once the buckling mode is known, solve the workspace of the rod CD, applying the concepts of Section 2.4.2, and plot it.
An example of the application of this algorithm is shown in Figures 10-15 (an animation is attached in the Supplementary Materials: Video S1), in which a complete turn of the crank AB is simulated while an external horizontal force of 5 N at point F is applied. Taking into account that the points of null moment of the rod CD are plotted in red, the flexible rod CD is working in its first buckling mode when the cycle starts (see Figure 10). Moreover, k parameter of the EI that allows calculating the workspace is positive. The rod CD is working in this mode (see Figure 11) until it reaches an α of 159 degrees, where its buckling mode changes from 1 to 2 (see Figure 12). In addition, k also changes from positive values to negatives ones. The rod CD remains in its second buckling mode (see Figure 13) until the input α acquires a value of 285 degrees (see Figure 14). From this pose on, the rod CD works in its first buckling mode, with negatives k, until the crank AB completes the turn (see Figure 15).     If the trajectory of point F is plotted against the moment of external loading about B for this load case, Figure 16 is obtained. It can be seen that the path is contained in the locus exposed on the previous Section 3.2.
It is worth noting that a cut appears on the surface for moments less than −2 Nm. The potential of this method to obtain the locus in which trajectories can lie for different load cases is highlighted here. Looking at the shadowed surface of Figure 16, it can be detected that load cases that generate moments about B less than −2 Nm cannot produce complete closed trajectories of the point F.
Thus, it is demonstrated that the SM proposed in Section 3.1 can detect and solve buckling mode changes as long as the trajectory of the end-effector point is contained in its corresponding locus.

Multiplicity of Mechanism Circuits
Besides buckling mode changes, the method here developed also can find all mechanism circuits. The process to obtain all circuits is the following:

1.
Integrate the deformed shape with a 4th order R-K method in points along a circumference with the center at point B and a radius equal to the length of bar BC. In this way, all potential solutions of rod CD that satisfy geometric conditions of the mechanism are calculated. Once the forces in the distal end of the rod CD are known, continue with the next step.

2.
Find the points of the circumference where there is a sign change in the summation of moments about B. These points are the potential circuits of the mechanism. 3.
Apply a SM to refine the potential initial solutions. The function to minimize is (39).

4.
Run the SM explained in Section 3.1 to solve the cycle for each calculated circuit.
An example of the application of this process is shown in Figures 17-22. For the case of the mechanism analyzed in this paper, three circuits appear when external forces of F ext,x = −4 N and F ext,y = −4 N, and a external moment of M ext,z = 4 Nm are applied on point F. In the initial position of the first circuit, the rod CD is straight (Figure 17), and completes a turn without changing the buckling mode, but with a sign change of EI's parameter k. As for the second circuit, the rod CD starts and carries out the whole cycle in its second buckling mode ( Figure 19). Finally, if the mechanism moves on its third circuit, the rod CD works in its first buckling mode ( Figure 21). The complete trajectories for the first, second and third circuits are plotted in Figures 18, 20 and 22, respectively. In these figures, it can be seen that each circuit generates a specific trajectory.
For a better understanding, Videos S2-S4 are available as Supplementary Materials.

Algorithm's Flowchart
To summarize all steps of the entire algorithm, its flowchart is given in Figure 23.  Figure 23. Flowchart of the entire algorithm.

Discussion
In this paper, the performance of the SM that combines the 4th order R-K method with a numerical minimization method, i.e., N-R scheme, to solve the position of hybrid compliant mechanisms is shown. This method not only allows solving the FKP, but also allows predicting the locus of all possible paths of the mechanism's coupler under different load cases. In this manner, this method can be used in the synthesis process of mechanisms to know whether the concerned mechanism can describe the desired trajectory supporting the external loading.
Furthermore, the SM is combined with the use of EI to obtain the workspace of the flexible rod, allowing detecting changes of buckling modes.
Future work could look for an algorithm able to detect abrupt changes or loops on the deformed shape of the compliant element. If this information is known, a highly detailed knowledge of the mechanism motion can be acquired. All this would allow designing mechanisms able to foresee the possible problems that could appear in the prototyping phase.
Further research can be conducted in the multiplicity of circuits of this kind of mechanism, applying the method here proposed. An ambitious goal could be an algorithm to obtain the loci of all trajectories for different load cases and all the mechanism circuits since, as can be seen in Figures 17-22, diverse circuits lead to different trajectories' loci. Detecting the influence of external loading in these circuits also should be studied.

Conclusions
The contribution of this paper is an algorithm that combines two numerical methods (a shooting method and elliptic integrals) to solve the movement of a hybrid rigid-flexible mechanism, taking into account the influence of the external load in that movement. The external load can affect the coupler trajectory by changing the deformed shape of the slender rod since, depending on the load, the flexible rod can work in different buckling modes, or even the whole mechanism can operate in different circuits. For this reason, the proposed algorithm is also developed to detect the circuit and buckling mode in which the mechanism works.
With this work, the potential of this method to study the behavior of compliant mechanisms is demonstrated. Its potential is such that it can be used, for instance, in the modeling process of humanoid robots, such as the one that appears in [23].

Acknowledgments:
The authors wish to acknowledge the financial support received from the Spanish government through the Ministerio de Economía y Competitividad (Project DPI2015-64450-R and PID2020-116176GB-100 (MINECO/FEDER, UE)), and the support for the research group through Project Ref. IT949-16, provided by the Departamento de Educación, Política Lingüística y Cultura from the regional Basque Government.

Conflicts of Interest:
The authors declare no conflict of interest.