Locomotion Dynamics for Bio-inspired Robots with Soft Appendages: Application to Flapping Flight and Passive Swimming

Abstract

In animal locomotion, either in fish or flying insects, the use of flexible terminal organs or appendages greatly improves the performance of locomotion (thrust and lift). In this article, we propose a general unified framework for modeling and simulating the (bio-inspired) locomotion of robots using soft organs. The proposed approach is based on the model of Mobile Multibody Systems (MMS). The distributed flexibilities are modeled according to two major approaches: the Floating Frame Approach (FFA) and the Geometrically Exact Approach (GEA). Encompassing these two approaches in the Newton–Euler modeling formalism of robotics, this article proposes a unique modeling framework suited to the fast numerical integration of the dynamics of a MMS in both the FFA and the GEA. This general framework is applied on two illustrative examples drawn from bio-inspired locomotion: the passive swimming in von Karman Vortex Street, and the hovering flight with flexible flapping wings.

Appendices

Appendix 1: Model of External Forces for the Two Examples

In this appendix, we provide the models of external forces used in the numerical simulation of flapping-wing flight and in the swimming simulation of the dead fish. The first model is an adaptation of the model of Dickinson et al. (1999), which was proposed in Porez et al. (2014), while the second is an extension of the Lighthill swimming model proposed in Candelier et al. (2011), Porez et al. (2014).

1.1 1.1 Model of the Aerodynamical Forces Exerted on a Flapping Wing

Due to the complex topology of the flow surrounding the flapping wings of an insect, the detailed modeling of aerodynamic forces involved in flapping-wing flight is a challenging issue, which extends beyond the scope of this article (Liu 2009). In Porez et al. (2014), we adopted the following simplified analytical model of the density of aerodynamic forces per unit of span length l,

$$\begin{aligned} \forall X\in [0,l]:\; F_\mathrm {aero}(X) = F_\mathrm {steady}(X) - \left( \mathcal {M}_f \frac{\partial \eta }{\partial t}\right) (X) + \left( \mathrm {ad}_{\eta }^\mathrm {T}(\mathcal {M}_f \eta )\right) (X),\nonumber \\ \end{aligned}$$

(40)

where, in accordance with our Cosserat beam modeling, the four terms of (40) describe densities of aerodynamic forces per unit of wing span. This model is directly inspired from the model of Dickinson et al. (1999), which was initially established using a dynamically scaled rigid model of a fruit fly’s wing. Once multiplied by dX (X being the span-wise abscissa), (40) represents a model of the aerodynamic force component applied on each chordwise strip (or blade, as used in aerodynamics). Such a model was applied in Whitney and Wood (2010) to represent the dynamics of flexible flapping wing. On the right hand side of (40), we have 1$^\circ$) the quasi-steady forces which model the leading edge vortex effect, together with the rotational circulation preventing the lift from decreasing when the wing quickly turns on itself at the end of a stroke (Dickinson et al. 1999),⁸ 2$^{o}$) the effect of added mass, in which $\mathcal {M}_f(X)$ represents the $6\times 6$ tensor of added mass of the X-cross section. In addition, note that the quasi-steady forces can be detailed as $F_\mathrm{steady}= \mathrm{Ad}_h^\mathrm{T}.^\mathrm{a}F_\mathrm{steady}$, where $^\mathrm{a}F_\mathrm{steady}$ is the wrench of aerodynamic forces expressed in a frame aligned with the incident flow centered on the aerodynamical center.⁹ This frame is hereafter referred to as the aerodynamical frame (indexed $^\mathrm{a}$), which is related to the cross-sectional frame (attached to the leading edge) by the rigid (homogeneous) transformation h(X). In this aerodynamical frame of reference, the steady forces can be represented as follows, $^\mathrm{a}F_\mathrm{steady}=(L_\mathrm{stat}^\mathrm{T}+D_\mathrm{stat}^\mathrm{T},0^\mathrm{T})^\mathrm{T}$ with $L_\mathrm{stat}(X)$ and $D_\mathrm{stat}(X)$ the lift and drag forces exerted on the X-cross section.¹⁰ This forces are of the form

$$\begin{aligned} \left\{ \begin{array}{rcl} L_\mathrm{stat} &{}=&{} (1/2)\;\rho _\mathrm{air}\;c(X)C_{L,\mathrm{stat}}\; |V_\mathrm{a}(X)|^{2},\\ D_\mathrm{stat} &{}=&{} (1/2)\;\rho _\mathrm{air}\;c(X)C_{D,\mathrm{stat}}\; |V_\mathrm{a}(X)|^{2}, \end{array}\right. \end{aligned}$$

(41)

where $\rho _\mathrm{air}$ is the density of the air, c(X) is the X-cross-sectional cord length, and $|V_\mathrm{a}(X)|$ is the norm of the velocity of the aerodynamical center of the blade of spanwise abscissa X. In addition, $C_{L,\mathrm{stat}}$ and $C_{D,\mathrm{stat}}$ are empirical dimensionless coefficients which, in the case of a fruit fly (Dickinson et al. 1999), are of the form

$$\begin{aligned} \left\{ \begin{array}{rcl} C_{L,\mathrm{stat}} &{}=&{} 1.8 \sin (2 \beta ),\\ C_{D,\mathrm{stat}} &{}=&{} 1.92-1.55 \cos (2 \beta ), \end{array}\right. \end{aligned}$$

(42)

where $\beta (X)$ is the angle of attack of the X-cross section.¹¹ Finally, the effect of added mass is defined as the reaction due to the acceleration of the mass of the fluid surrounding each blade of the wings. In the case of 40, the term $\mathrm{ad}_{\eta }^\mathrm{T}(\mathcal {M}_f\eta )$ is the added mass term due to Coriolis-centrifugal accelerations.

1.2 1. 2 Model of hydrodynamical forces exerted on a swimming body

In the second example, we consider the case of a neutrally buoyant fish (-like robot) with elliptic cross sections. The external forces are hydrodynamic forces locally exerted on the cross sections of the fish (-like robot). Following the works of Lighthill (1960, 1971), several refinements of the Large Amplitude Elongated Body Theory (LAEBT) have been proposed the recent years (Boyer et al. 2010; Candelier et al. 2011, 2013; Porez et al. 2014). In particular, the LAEBT has been extended to model the external forces acting on a carangiform fish swimming in an ambient flow defined by its field of velocity $v_{f}$ in any point of the ambient space. Based on the integration of the pressure field (given by the unsteady Bernoulli equation) in the curvilinear coordinates of the Cosserat beam prolonged to the nearby flow, the model can be detailed as follows (Candelier et al. 2013):

$$\begin{aligned}&\forall X\in ]0,l[: F_{ext}(X) = - \mathcal {M}_f \frac{\partial \eta _{r}}{\partial t} +\mathrm {ad}^\mathrm {T}_{\eta _{r}} (\mathcal {M}_f \eta _{r})-F_\mathrm {flow}, \end{aligned}$$

(43)

$$\begin{aligned}&F_\mathrm {ext}(l) =-(\mathcal {M}_f \eta _{r})\left( e_{1}^{T}\eta _{r}\right) - \frac{1}{2}\left( \eta _{r}^T \mathcal {M}_f \eta _{r}\right) e_1, F_\mathrm {ext}(0)=0, \end{aligned}$$

(44)

with $e_1= (1 \ 0 \ 0 \ 0 \ 0 \ 0)^T$ and $\mathcal {M}_{f}(X)$ the $6\times 6$ added mass tensor of the X-cross section¹² in is body frame F(X). Since F(X) is centered on the cross-sectional mass center, $\mathcal {M}_{f}(X)$ takes the form:

$$\begin{aligned} \mathcal {M}_f=\left( \begin{array}{cc} m_{f} &{} 0 \\ 0 &{} I_{f} \\ \end{array} \right) , \end{aligned}$$

(45)

where we introduced the matrices of linear and angular added inertia expressed in the cross-sectional frames centered on the elliptic cross-sectional centroid:

$$\begin{aligned} m_{f}=\left( \begin{array}{ccc} 0&{} 0 &{} 0 \\ 0&{} m_{f2} &{} 0\\ 0&{} 0 &{} m_{f3} \\ \end{array} \right) , \text { }I_{f}=\left( \begin{array}{ccc} I_{f1}&{} 0 &{} 0 \\ 0&{} 0 &{} 0\\ 0&{} 0 &{} 0 \\ \end{array} \right) \end{aligned}$$

(46)

In (43)–(44), $\eta _{r}=(V_{r}^{T},0^{T})^{T}$ is the relative velocity of the cross section’s frame with respect to the fluid¹³, where $V_{r}(X)=V(X)-V_{f}(X)$ with V(X) and $V_{f}(X)=R^{T}(X)v_{f}(O(X))$ being the linear velocity of the beam and the fluid (without the beam) both expressed at the center O(X) of the X-cross section in the basis of F(X). With this conventional notation, the two first terms in Eq. (43) (as in the case of the flapping wing) represent the forces exerted by the fluid mass laterally accelerated by each cross section. In contrast to the flapping wing, this hydrodynamic model accounts for the vorticity shed in the wake (Candelier et al. 2013) through the second term of (44). Finally, in the case of the planar swimming, the last term of (43) takes the form $F_\mathrm{flow}=(f_\mathrm{flow}^{T},m_\mathrm{flow}^{T})^{T}$, where:

$$\begin{aligned}&f_\mathrm {flow} =-m_{b} R^{T}\gamma _{f}+R^{T} \left( \left( \nabla _{X}v_{f}\right) \left( m_{f}V_{r}\right) \right) , \end{aligned}$$

(47)

$$\begin{aligned}&m_\mathrm {flow} =\left( 0,0,\widetilde{m}_{f}\partial V_{r,2}/\partial t\right) ^{T}, \end{aligned}$$

(48)

with $m_{b}(X)$ the mass of the X-cross section, and $\gamma _{f}$ and $\nabla _{X}v_{f}$ the acceleration and the gradient of the velocity field¹⁴ of the ambient flow (without the body), both evaluated at the center O(X) of the cross section and expressed in the axes of the fixed frame ${F}_{e}$. The field of hydrodynamic couples exerted along the beam given by (48) only depends on the cross section variations since $V_{r,2}=(0,1,0)V_{r}$ and $\widetilde{m}_{f}=(m_{f}/4)(a\partial b/\partial X+3b \partial a/\partial X$ with a(X) and b(X) the length of the two axes of the elliptic X-cross section (see Fig. 3b).

Appendix 2: Newton–Euler Model of the MAV of Fig. 1a

The purpose of this appendix is to provide the detail of the notations used in the generalized Newton–Euler model in (17)–(20). First, note that all quantities referring to purely rigid movements were introduced prior to (17)–(20). In particular, $g_j$, $\eta _j$, $\dot{\eta }_j$ denote the associated transformation, velocity and acceleration, respectively, of the $\mathcal {B}_j$-body frame, now defined as the $j\mathrm{th}$ floating frame. In addition, $\mathcal {M}_j$, $\mathcal {F}_{j}$ are the $6\times 6$ inertia matrix and the $6\times 1$ vectors of Coriolis-centrifugal and external forces, applied to the flexible body when it is locked in its current configuration $r_{\mathrm{e}j}$. The additional variables used in (17)–(20) are related to the elastic degrees of freedom. In (17), $m_j$ is the matrix of modal inertia, $M_j^\mathrm{T}$ is the matrix of inertia coupling between rigid reference accelerations and modal accelerations, and $Q_{j}$ is the $N_{ej}\times 1$ vectors of modal inertia (Coriolis-centrifugal), external and internal (restoring elastic) forces, while $\varPhi _j(l_{j})=(\varPhi _{dj}(l_{j})^\mathrm{T},\varPhi _{rj}(l_{j})^\mathrm{T})^\mathrm{T}$ is the matrix of linear and angular shape functions expressed at $X=l_j$, (i.e., at the last cross section of $\mathcal {B}_j$). In (18)-(20), $R_{g_{j,i}}$ is the $6\times 6$ matrix transforming a twist from the basis of ${F}_i$ to the basis of ${F}_j$, leaving the origin unchanged, and $A_j$ is the unit twist supporting the joint axis. Consequently, $A_i\dot{r}_i = g_{ai}^{-1}\dot{g}_{ai}$ accounts for the joint velocity, and $\zeta _i$ contains the corresponding accelerations ($\ddot{r}_i$), as well as the residual velocity-dependent accelerations resulting from the time differentiation of 19. Similarly, $\varPhi _i(l_{i})\dot{r}_{\mathrm{e}i} = \dot{g}_{\mathrm{e}i} g_{\mathrm{e}i}^{-1}$ represents the $6\times 1$ (elastic) velocity vector introduced by the body $\mathcal {B}_i$ which precedes $\mathcal {B}_j$ along the structure’s branch. The interested reader is referred to Boyer and Coiffet (1996) for the detailed expression of these matrices in the general case. In the special case of the MAV described in Sect. 3, because the two flexible wings $\mathcal {B}_j$, $j=1$, 2, are terminal bodies, $\varPhi _j(l_{j})$ plays no role in the model. For a pure torsional deformation (here described with one torsional mode $\varphi _{t1}$ of modal coordinate $\theta _1$ per wing), the inertia matrix and the vector of Coriolis-centrifugal (and external) forces in the wings’ generalized NE model (17) can be computed with the exact kinematics (i.e., not linearized) of the wings in the floating frame. The resulting expressions are the following

$$\begin{aligned} \mathcal {M}_j= & {} \left( \begin{array}{cc} {\mathrm{m}}_j 1_{3\times 3} &{} m\widehat{s}_j^\mathrm{T}\\ m\widehat{s}_j &{} J_j \\ \end{array}\right) \text { ,}\nonumber \\ \mathcal {F}_{j}= & {} \left( \begin{array}{c} {\mathrm{m}}_j(V_j \times \varOmega _j) - \varOmega _j \times (\varOmega _j \times ms_j)+f_{\mathrm{ext,}} j\\ ms_j \times (V_{j} \times \varOmega _j)-\varOmega _j\times (J_j \varOmega _j) + b_{j} + m_\mathrm{ext,j} \\ \end{array}\right) , \end{aligned}$$

(49)

where $\mathrm{m}_j =\mu _0$ is the total mass of $\mathcal {B}_j$, $b_{j}=\dot{\theta }_1\left( 0,-\mu _{2\tau ,1}\varOmega _{jz},\mu _{2\tau ,1}\varOmega _{jy}\right) ^\mathrm{T}$, and the matrices of first and second momenta of inertia are obtained from

$$\begin{aligned} ms_j = \left( \begin{array}{c} \mu _{x0} \\ 0 \\ 0 \\ \end{array}\right) , \quad J_j = \left( \begin{array}{ccc} \mu _{2} &{} 0 &{} 0 \\ 0 &{} \mu _{xx0}+0.5\mu _{2} &{} 0 \\ 0 &{} 0 &{} \mu _{xx0}+0.5\mu _{2} \end{array}\right) .\nonumber \end{aligned}$$

In addition, we have

$$\begin{aligned} m_j = \mu _{2\tau \tau ,11}, \quad M_j =\left( \; 0,\; 0,\; 0,\; \mu _{2\tau ,1},\; 0,\; 0 \;\right) , \quad Q_j= - k_{\tau ,11}\theta _{1} + Q_{\mathrm{ext}, j}.\nonumber \end{aligned}$$

Computing the above expressions requires the value of the following inertial parameters ($\rho$ denotes the mass per unit of beam volume, X denotes the span-wise beam coordinate),

$$\begin{aligned} \begin{array}{lll} \mu _0 = \int _\mathrm{span}\rho A \mathrm{d}X, &{} \mu _{x0} = \int _\mathrm{span}X\rho A\mathrm{d}X, &{}\mu _{xx0} = \int _\mathrm{span}X^2\rho A\mathrm{d}X, \\ \mu _2 = \int _\mathrm{span}\rho I_{\rho }\mathrm{d}X, &{} \mu _{2\tau ,1} = \int _\mathrm{span}\rho I_{\rho }\varphi _{t1}\mathrm{d}X , &{}\mu _{2\tau \tau ,11} = \int _\mathrm{span}\rho I_{\rho } \varphi _{t1} \varphi _{t1} \mathrm{d}X, \end{array}\nonumber \end{aligned}$$

as well as that of the modal stiffness,

$$\begin{aligned} k_{\tau ,11} = \int _\mathrm{span}GI_{\rho }\varphi '_{t1}\varphi '_{t1}\mathrm{d}X\nonumber , \end{aligned}$$

where G is the twisting modulus of the beam, A and $I_{\rho }$ its cross-sectional area and polar inertia, respectively. These parameters being time independent, they can be computed once initially, then used for all simulation times. Finally, we need to compute the external forces $f_{\mathrm{ext,}j}$, $m_{\mathrm{ext,}j}$, and $Q_{\mathrm{ext,}j}$, which are exerted by gravity and the flow surrounding the wing. The gravity force field can be directly included within the algorithm as a net overall acceleration added to $\dot{\eta }_0$ given by (22).¹⁵ The model of external forces required by the algorithm thus reduces to that of aerodynamical forces, which are primarily given by

$$\begin{aligned} \begin{array}{rcl} f_{\mathrm{ext}, j} &{}=&{} \int _\mathrm{span}f_\mathrm{aero}(X)\; dX, \\ m_{\mathrm{ext}, j} &{}=&{} \int _\mathrm{span}(\left( X,0,0\right) ^\mathrm{T}\times f_\mathrm{aero}(X)+m_\mathrm{aero}(X))\; dX,\\ Q_{\mathrm{ext}, j} &{}=&{} \int _\mathrm{span}(m_{\mathrm{aero}} \varphi _{t1})(X) dX, \end{array} \end{aligned}$$

(50)

where $(f_\mathrm{aero}^\mathrm{T}, m_\mathrm{aero}^\mathrm{T})(X)^\mathrm{T}= F_\mathrm{aero}(X)$ is defined by (40). Additional analytical details could be provided regarding these forces. However, in the simulation discussed in Sect. 7, only the steady component of (40) is used. Meanwhile, the value of $f_{\mathrm{ext,}j}$, $m_{\mathrm{ext,}j}$, and $Q_{\mathrm{ext,}j}$ is obtained from a numerical integration based on the Gauss quadrature method. This space integration is performed at each iteration of the algorithm, using the current values of $V_\mathrm{a}$ and $\beta$ required by the model of steady aerodynamic forces (41)–(42), along with the exact expressions (i.e., not linearized) of the position and velocity of the wing material points in the floating frame.

Appendix 3: Mixed Dynamics Algorithm for Soft Locomotion Systems

The algorithm is structured in three recursions. It uses a Boolean variable $b_o(j)=1$ or $b_o(j)=0$ depending on whether the $j\mathrm{th}$ joint’s movement or torque is imposed. Note that if the joint is the result of the discretization of the beam’s strain field, we have $b_o(j)=0$. The term $\zeta _{j}$ of (20) is of the form $\zeta _j= A_j \ddot{r}_j + \nu _j$, where $\nu _j$ is a velocity-dependent acceleration term deduced from the time differentiation of (18). It is known for any value of $b_o(j)$, whereas $\zeta _j$ is known only when $b_o(j)=1$. For the sake of concision, we use the simplified notations $g_{ei}(l_{i})=g_{ei}$ and $\varPhi _{i}(l_{i})=\varPhi _{i}$.

• First Forward Recursion

  $$\begin{aligned}&\mathbf{For } \,\,j=1 \,\mathbf{ to } j=n \mathbf{, compute }\nonumber \\&\quad g_j = g_i g_{\mathrm{e}i} g_{\mathrm{a}j}, \end{aligned}$$

  (51)
  $$\begin{aligned}&\quad \eta _j = \mathrm{Ad}_{g_{j,i}} \eta _i + R_{g_{j,i}} \varPhi _i \dot{r}_{ei} + A_j \dot{r}_j, \end{aligned}$$

  (52)
  $$\begin{aligned}&\quad \mathrm{Ad}_{g_{i,j}},\, R_{g_{i,j}},\, \mathcal {M}_j,\, m_{j},\, M_j,\, \mathcal {F}_j,\, Q_j \end{aligned}$$

  (53)
  $$\begin{aligned}&\quad \text{ If } b_o(j)=1,\; \text{ compute } \zeta _j,\end{aligned}$$

  (54)
  $$\begin{aligned}&\quad \text{ If } b_o(j)=0,\; \text{ compute } \nu _j.\\&\mathbf{End \, for. }\nonumber \end{aligned}$$

  (55)

• Backward Recursion

  $$\begin{aligned}&\mathbf{For }\, j=n \,\mathbf{ to } \,j=1 \,\mathbf{, compute }\nonumber \\&\quad \text{ If } b_o(j)=1\nonumber \\&\quad \varUpsilon = \widetilde{\mathcal {M}}^+_j,\end{aligned}$$

  (56)
  $$\begin{aligned}&\quad \gamma = \zeta _j,\end{aligned}$$

  (57)
  $$\begin{aligned}&\quad \text{ Else }\nonumber \\&\quad k = A^\mathrm{T}_j \widetilde{\mathcal {M}}^+_j A_j, \end{aligned}$$

  (58)
  $$\begin{aligned}&\quad E = 1_{6\times 6} - k^{-1} \left( A_j A^\mathrm{T}_j\right) \widetilde{\mathcal {M}}^+_j, \end{aligned}$$

  (59)
  $$\begin{aligned}&\quad \varUpsilon = \widetilde{\mathcal {M}}^+_j E,\end{aligned}$$

  (60)
  $$\begin{aligned}&\quad \gamma = E \nu _j + k^{-1} A_j \left( \tau _j + A^\mathrm{T}_j \widetilde{\mathcal F}^+_j\right) ,\end{aligned}$$

  (61)
  $$\begin{aligned}&\quad \text{ End } \text{ if }\nonumber \\&\quad \mathcal {M}^+_i =\mathcal {M}^+_i + \mathrm{Ad}^\mathrm{T}_{g_{j,i}} \varUpsilon \mathrm{Ad}_{g_{j,i}}, \end{aligned}$$

  (62)
  $$\begin{aligned}&\quad {M^+_i}^\mathrm{T}= {M^+_i}^\mathrm{T}+ \mathrm{Ad}_{g_{j,i}}^\mathrm{T}\varUpsilon R_{g_{j,i}} \varPhi _i, \end{aligned}$$

  (63)
  $$\begin{aligned}&\quad m^+_i = m^+_i + \varPhi ^\mathrm{T}_i R_{g_{i,j}} \varUpsilon R_{g_{j,i}} \varPhi _i, \end{aligned}$$

  (64)
  $$\begin{aligned}&\quad \mathcal F^+_i = \mathcal F^+_i + \mathrm{Ad}^\mathrm{T}_{g_{j,i}} \left( \widetilde{\mathcal F}^+_j - \widetilde{\mathcal {M}}^+_j \gamma \right) , \end{aligned}$$

  (65)
  $$\begin{aligned}&\quad Q^+_i = Q^+_i + \varPhi ^\mathrm{T}_i R_{g_{i,j}} \left( \widetilde{\mathcal F}^+_j - \widetilde{\mathcal {M}}^+_j \gamma \right) , \end{aligned}$$

  (66)
  $$\begin{aligned}&\quad \widetilde{\mathcal {M}}^+_i = \mathcal {M}^+_i - {M^+_i}^\mathrm{T}{m^+_i}^{-1} M^+_i, \end{aligned}$$

  (67)
  $$\begin{aligned}&\quad \widetilde{\mathcal F}^+_i = \mathcal F^+_i - {M^+_i}^\mathrm{T}{m^+_i}^{-1} Q^+_i, \\&\mathbf{End \,for. }\nonumber \end{aligned}$$

  (68)

• Second Forward Recursion

  $$\begin{aligned}&\mathbf{For } \,j=1 \mathbf{to }\, j=n \,\mathbf{, compute }\nonumber \\&\quad \mu = \mathrm{Ad}_{g_{j,i}} \dot{\eta }_i + R_{g_{j,i}}\varPhi _i \ddot{r}_{ei} + \nu _j,\end{aligned}$$

  (69)
  $$\begin{aligned}&\quad \text{ If } b_o(j)=1\nonumber \\&\quad \tau _j = A^\mathrm{T}_j \left( \widetilde{\mathcal {M}}^+_j \dot{\eta }_j - \widetilde{\mathcal F}^+_j\right) ,\end{aligned}$$

  (70)
  $$\begin{aligned}&\quad \text{ Else } \nonumber \\&\quad \ddot{r}_j = k^{-1} \left( \tau _j - A^\mathrm{T}_j \left( \widetilde{\mathcal {M}}^+_j \mu - \widetilde{\mathcal F}^+_j\right) \right) ,\end{aligned}$$

  (71)
  $$\begin{aligned}&\quad \text{ End } \text{ if }\nonumber \\&\quad \dot{\eta }_j = \mu +A_j\ddot{r}_{j},\end{aligned}$$

  (72)
  $$\begin{aligned}&\quad \ddot{r}_{ej} = {m^+_j}^{-1} \left( Q^+_j - M^+_j \dot{\eta }_j\right) ,\\&\mathbf{End \,for. }\nonumber \end{aligned}$$

  (73)

These three recursions are initialized as follows. The first recursion is initialized by the current state of the reference body $(g_0,\eta _0)$, the backward recursion is initialized by $\widetilde{\mathcal {M}}_{j}^{+}=\widetilde{\mathcal {M}}_{j}=\mathcal {M}_j -M_j^\mathrm{T}m_j^{-1}M_j$ and $\widetilde{\mathcal F}_j^+=\widetilde{\mathcal F}_j=\mathcal F_j-M_j^\mathrm{T}m_j^{-1}Q_j$, for all j, while the last recursion is initialized by the current reference state $(g_0,\eta _0)(t)$ and the current reference accelerations $\dot{\eta }_0(t)$, computed after the second recursion from $\dot{\eta }_0 = \widetilde{\mathcal {M}}_0^{+-1}\widetilde{\mathcal F}_0^+$. Finally, in the case that $\mathcal {B}_j$ is a rigid body, all modal matrices of (17)–(20) can be removed and $\widetilde{\mathcal {M}}_j^+=\mathcal {M}_j^+$, $\widetilde{\mathcal F}_j^+=\mathcal F_j^+$. In particular, when applying the algorithm to the GEA, all bodies are rigid and the algorithm is reduced to a rigid version given in Featherstone (2008). Finally, this algorithm only requires derivation of the generalized Newton–Euler model in a preliminary stage, all remaining calculations can be automatically performed (either numerically or symbolically).