# Modeling and frequency domain analysis of nonlinear compliant joints for a passive dynamic swimmer

###### Abstract

In this paper we present the study of the mathematical model of a real life joint used in an underwater robotic fish. Fluid-structure interaction is utterly simplified and the motion of the joint is approximated by Düffing’s equation. We compare the quality of analytical harmonic solutions previously reported, with the input-output relation obtained via truncated Volterra series expansion. Comparisons show a trade-off between accuracy and flexibility of the methods. The methods are discussed in detail in order to facilitate reproduction of our results. The approach presented herein can be used to verify results in nonlinear resonance applications and in the design of bio-inspired compliant robots that exploit passive properties of their dynamics. We focus on the potential use of this type of joint for energy extraction from environmental sources, in this case a Kármán vortex street shed by an obstacle in a flow. Open challenges and questions are mentioned throughout the document.

## 1 Introduction

How much of the diverse behavior we observe in animals is a direct expression of the dynamics of the individual’s body? In the last two decades, many characteristics of animal locomotion on land were successfully linked to the mechanical properties of the legs. The springy behavior observed in the trajectories of the center of mass during running, walking and jumping can be explained by the mechanical properties of the limbs and its tunning (Farley et al. (1993); Farley (1996); Ferris et al. (1998); Dickinson et al. (2000); Kerdok et al. (2002); Moritz and Farley (2003); Roberts and Azizi (2011)). The hypothesis that running, walking and jumping is tuned to the resonance of the underlying mechanical system is strongly supported by experimental evidence and by machines constructed based on this idea, the passive dynamic walkers (Alexander (1990); Ahlborn and Blake (2002); McGeer (1990); Thompson and Raiber (1989); Collins et al. (2005)). The consequence of such setting is energy efficient performance and alleviation of the controller.

Due to similarities between running and swimming (Bejan and Marden (2006); Kokshenev (2010)), it is not surprising (but not less exiting) that efficient locomotion in fluids was reported to relay on the dynamics of the body of the animal. Living trouts have been observed to exploit the energy in the flow they inhabit to reduce their swimming efforts (Liao et al. (2003)). Later, euthanized trouts performed passive self-propulsion when placed in the von Kármán vortex street shed by an obstacle in a flow (Liao (2007)). The experimental results are supported by mathematical models, analytical and numerical (Eldredge and Pisani (2008); Kanso and Newton (2009); Alben (2009)). These models do not fully agree with each other, however this is expected due to the mathematical complexity of the interaction between structures and fluids, of which the passive case is the worst scenario: unprescribed motion of the interface boundary. This alone represents an open challenge for the mathematical modeling community. Lest the challenge remains unsolved for too long, robotic researchers design their machines with less fluid-dynamic rigor, pursuing the first passive dynamic swimmer. The final objective is to build a swimming machine that can perform at least as well as fish, and at comparable power ratings (Harper et al. (1997); Lauder et al. (2007)).

For a robot to extract the energy in the surroundings, its mechanical properties have to be tuned to the environmental energy storage. Stated this way, the problem is one of energy harvesting, were nonlinear properties are believed to be beneficial (Cottone et al. (2009)). Therefore, we need to pin down the resonance characteristic of the actuators and joints to be used in the robot that, as their biological counterparts, are generally nonlinear. As if difficulty was lacking, the study of resonance of nonlinear systems has suffer a very slow developmental process that started early in the 1960’s and has yet not overcome its infancy. In a technical report from 1958 by Brilliant (1958) we read:

Sometimes nonlinearity is avoided, not because it would have an undesired effect in practice, but simply because its effect cannot be computed.

Nowadays the situation is not completely different. However, we have more powerful computers, new simulation methods and some novel uses of classical tools promise to open the path ahead Peng et al. (2007); Vakakis et al. (2009).

In this paper we present the study of the joint of a robot fish from a mathematical point of view. The amplitude of the oscillations of the joint in response to periodic forcing is studied. In section 2 we briefly introduce the real mechanical device and we move to its mathematical model in section 3. In the same section the model for fluid-structure interaction is briefly described. Approximated solutions methods are introduced in section 4. Results are presented in section 5. Finally, we close the paper in section 6 with a discussion on the implications of the results and the relevance of the approach to the design of robots and the test of controllers based on resonances.

## 2 A simple compliant joint

To extract energy from the environment, the robotic platform has to be optimally driven by interaction forces. In this case, by the interaction between the rigid body of the robot and the surrounding turbulent flow. The feat can not be performed if the angle trajectories of the joints are fully prescribed by the controller. Therefore, the joints of our robot are compliant, i.e. the motion of the joint is not only defined by the controller, but also by external actions. In Figure 1a we see the details of the joints of the robot fish used in Ziegler et al. (2011). Each joint behaves as a rotational spring. The restoring torque is generated when the relative angle between the two connected bodies is not zero. The force producing the torque is given by the extension of a linear spring fixed to the first body. The spring is connected via an inelastic thread to an appendage of the second body. When the deflection angle is zero, the extension of the spring is minimum as well as the force it exerts. We call tension to this minimum force value and it is referred with the letter . Measurements of the torque for , are given in Figure 1b. The values of the parameters used throughout this paper are given in Table 1. The two parameters are distances that can be seen in the figures. The elastic constant of the linear spring is and denotes the moment of inertia of the joint around the axis of rotation. The linear specific damping coefficient of the joint is denoted with . The parameters and correspond to the amplitude and frequency of the external forcing, respectively.

Name | Value |
---|---|

## 3 Mathematical model

A geometrical representation of the joint is given in Figure 2. The parameters , and are fixed at construction time and always (Table 1). The tension in the spring at its shortest length (), is the controlled parameter and a servomotor can change it dynamically.

The torque applied to the bodies connected to the joint is thus,

(1) |

Where we have defined the parameters and . The formula is obtained by calculating the deformation of the spring as a function of the deflection angle of the appendage (). The deformation is inside parenthesis in the numerator of equation (1) and when multiplied by the stiffness , it gives the force due to the deformation of the linear spring. The outer factors come from the product of the force and the moment arm. Note that the the torque is linear in the controlled input .

For the third order Taylor expansion gives

(2) |

Where

(3) | |||||

(4) |

And the equation of motion of the deflection angle is

(5) |

Where denotes the moment of inertia around the axis of rotation. The specific damping is given by , and we have defined , . is the specific net effect of all other external torques acting on the joint. The approximating equation of motion is the well studied Düffing’s equation.

To quantify the error introduced by the approximation, we calculated the angle at which the difference between the torque produced by equation (1) and equation (2) is equal to a reference error given by , where is a reasonable resolution for a force sensor working in a range. These angles are plotted in Figure 3 for different values of the tension.

### 3.1 Hardening, linear and softening spring

As can be seen from equation (2), modulates the intensity of the cubic non-linear term. This term vanishes when

(6) |

rendering a linear spring for the angles where the third order approximation is valid. For bigger values of the spring will be softening () and for smaller values it will be a hardening spring () (Fig. 4). However, the full expression of the torque (Eq. (1)) contains higher order terms. Hence, the linear behavior will be even more evident if the higher order terms cancel each other. In Figure 4 curves of torque versus angle for several values of are shown. In particular we show the curve for which up to seventh order nonlinearities give a minimum contribution (found via optimization ) together with the curve at . These illustrates the power of the actuation chose, since we can control the dynamical properties of a virtual rotational spring at the joint (more details in Ziegler et al. (2011)).

### 3.2 Forcing model

As explained before, the external torques acting on the joint are due to fluid-structure interaction. This kind of interaction still poses a great challenge for the mathematical modeling community. In the search for simplified models of the forces generated in the interaction between flexible bodies and turbulent flows, we found the work of Kanso and Newton (2009) and Alben (2009) instructive. Using equation (3.8) given by Alben (2009), we can calculate the pressure difference on the boundary of a slender body in a vortex street. If the width of the body is bigger than the separation among vortices and it is placed in the middle of the vortex street, the forces on the body can be approximated by a sine function , where is time and is an initial offset that sets the balance between the sine and cosine behavior of the forcing. The frequency is proportional to the speed of the flow plus the vorticity of the vortices (assumed to be equal for each vortex) scaled by a factor that depends on the geometrical properties of the wake. The amplitude is proportional to the density of the fluid and the square of the vorticity of each vortex. With a few additional assumptions, the torque acting on the joint can be made proportional to this force. Here we adopt this over-simplified forcing model to avoid diverting the attention of the reader from the core ideas of our work. In this manner we postpone a detailed study with a more elaborated forcing model.

## 4 Solution methods

In this section we present two independent methods to estimate the amplitude of oscillations of the joint under periodic forcing. The first method turns out to be excellent for estimating the amplitude of the first harmonic. The second method is based on a Volterra series expansion of the Duffing equation. It allows to study the response of the joint to more general forcing conditions, but the tension range for which it is valid is smaller.

### 4.1 Harmonic solutions of Düffing’s equation

Under periodic forcing , equation (5) has been extensively studied (see Holmes and Rand (1976); Luo and Han (1997) and references therein). Following these analyses, we show here how to maximize the amplitude of the periodic response of the joint by tunning the tension parameter.

The key point of the analysis in Luo and Han (1997), is that we search for the amplitude of solutions of (5) that are periodic (this rules out sub-harmonics, supra-harmonics and chaotic motion), . Under this assumption it can be shown that the amplitude of such solutions is given by the roots of the polynomial,

(7) |

The roots of this polynomial can be obtained analytically using a computer algebra system as Maxima (2009) and they establish the relation between the amplitude of the oscillations and the parameters of the equation, . In the case at hand, we have and , therefore . For a given value of the parameters, only one root corresponds to the observed amplitude (there are unstable amplitudes). This implies, that is not enough to look at the roots, but we must also check their stability. This adds some complexity to the evaluation of the results that adds to the limitation of pure harmonic inputs.

### 4.2 Volterra series expansion

When the amplitude of the forcing is sufficiently small, the behavior of Duffing’s oscillator in the neighborhood of the origin can be described by polynomial Volterra functionals(Theorem 3.1 of Rugh (1981)). This means that the angle of the oscillations can be described by an expansion of the form

(8) |

where is the order of the expansion and each term is the multi-dimensional integral

(9) |

where is the forcing signal and we wish to determine the kernel functions . However, for the objective at hand, it is more useful to calculate the Fourier transform of these functionals. We have not found this process described in the literature, hence we describe it succinctly. We used the inverse Fourier transform definition on the input signal and substitute it into the functional definition Eq. (9). By inverting the order of integration and splitting variables, we make the kernel transform appear. Then, we compare this expression to the formal definition of the inverse Fourier transform of and isolate the desired result. This process can be applied for the first three orders, but it gets cumbersome for higher ones. Therefore, before showing the results, we will briefly introduce some notation used here to simplify the presentation.

An ordered set of arguments will be denoted . In general we have,

For example, the fifth order kernel will be written

where . For the integration variable in multiple integrals we will write . For multiple summations over indexes we will write . Having defined the notation, we continue with our presentation.

For zero initial conditions, , and working in the frequency domain, each kernel can be obtained recursively based on lower order kernels (details of the calculation are given in Peyton-Jones and Billings (1989)). The first order kernel derived from the model in Eq. (5), is equivalent to the frequency response function (also known as transfer function) of a linear second order system,

(10) |

The following equations present the Volterra kernels up to seventh order using recurrent relations,

(11) |

Due to the symmetry of the eq. (5), all even order kernels are zero. Note that these kernels are valid for any second order cubic oscillator. The recursive formulas were obtained using the probing method of Peyton-Jones and Billings (1989), which proceeds as follows. We start from a polynomial nonlinear ordinary differential equation (as Eq. (5)) and assume the response can be represented by Volterra series, i.e. Eq. (8). We substitute this into the system’s differential equation, and equate similar terms. This last step produces equations relating different order kernels and the expressions are always recursive (we obtain Eqs. (11) in our case). The factorial in the recursive relations appear from the multinomial expansion of nonlinear polynomial terms.

We proceed to calculate the Fourier transform of the output given by Eq. (9) subject to harmonic forcing. We will use Eqs. (11) in order to obtain an explicit relation between the parameters of our model and the amplitude of the first harmonic of the output.

We recall that the Fourier transform of a general bandwith limited periodic function is a scaled Dirac Comb (Shah function, Impulse train, Dirac train, etc.) (see Schetzen (1980)),

(12) |

where is a positive integer, is the complex amplitude of the -th harmonic and is the impulse function (Dirac delta function). The Fourier transform of each of the output terms in (8) is

(13) |

where we have used Eq. (12) and and,

(14) |

The integral in (13), has an interesting geometric interpretation in terms of convolutions in hyperplanes, we refer the interested reader to Lang and Billings (1996). Therein, the output frequency range of nonlinear systems that are representable by Volterra series is analytically calculated. Additionally, in Eq.(16) of that paper the frequency spectrum of the output signal is represented as the superposition of contributions from the nonlinearities. In the case studied herein, the input has only one frequency, therefore , , with the imaginary unit. Replacing these values in all the equations and using the relations in Eq. (11) we obtain the desired result,

(15) |

Where is given by Eq. (10) with . The amplitude of the oscillation is obtanied by taking the double of the modulus of the complex number .

## 5 Results

Figure (a)a shows the amplitude of periodic oscillations in the plane, according to Eq. (7). The line of maxima is shown with dashes. This line describes the value of the tension that produces maximum amplitude for a given forcing frequency. In Fig. (b)b we plot the logarithm of the relative difference between the amplitudes given by (7) and (15). The two approximation differ for regions of low frequency and low tension where the system is most nonlinear (Luo and Han (1997)). To compare with the observed amplitude (obtained by simulating (5) without any approximation), we extracted the curves of amplitude against tension for (vertical dotted lines in Fig. 5). These curves are shown in Figure 6 together with the amplitude calculated using Eq. (15) corresponding to the Volterra kernel expansion. For the frequency, both models predict the simulated amplitude accurately. For the frequency, the amplitude predicted by Eq. (7) drifts away from the simulated value for lower tensions. The amplitude calculated using the Volterra expansion diverges for low tensions at this frequency.

In the same figure we show the amplitudes measured in simulations using the exact expression for the torque (1) (Octave (Eaton (2002)) function ode45, absolute and relative tolerances, ). The amplitude of the oscillations is obtained using the absolute value of the analytic signal of the angle (Octave function hilbert). The agreement between the approximated results and the simulated ones is noteworthy.

These results show that the approximation from Luo and Han (1997) can be used for a tension controller designed for a joint of this kind, to obtain periodic responses with the same frequency as the forcing. The Volterra approximation fails when the system becomes strongly nonlinear, however these expansion can be used for more general responses or when the inputs have a more complicated frequency spectrum (as in a realistic scenario).

## 6 Discussions and conclusion

Though the quote from Brilliant (1958) remains valid, we used the knowledge about Düffing’s equation to understand analytically a compliant joint. This gives a corner stone to study more complicated setups (e.g. chain of joints, as in the robotic fish) and define the robust engineering design of robots with the desired resonance properties. Additionally, the results showed here provide a reference solution to the problem of finding the right pretension for a given external forcing, related to the idea of adaptive controllers. For a forcing with a sufficiently slow varying frequency, the tension could be adjusted to maintain the amount of energy transferred into the system to its maximum possible value, i.e. keep the system close to the line of maxima in Fig. (a)a. With the results herein the adjustment could be done by direct calculations. However, other methods like the frequency oscillators(AFO) (Buchli et al. (2006)) were used to achieve a similar objective. Applications of the latter has been reported in Buchli and Ijspeert (2008), however, due to lack of theoretical results on the resonance frequency of the nonlinear platform studied therein, the theoretical solution to the problem was unknown and a grounded evaluation of the performance was not possible. Approximations as the one presented here could be used to verify the results generated by the adaptive oscillators.

In our results we exclude any information about oscillations with frequencies different from that of the forcing since that is the case where Eq. (7) is valid. The higher order harmonics in the system’s output could be determined and compared with the simulated results. Traditional harmonic balance method can not be used to find these results. Additionally, a more complete panorama of the frequency response of the system could be obtained with a nonlinear analysis method as the NOFRF (proposed in Lang and Billings (2005) and showcased in Peng et al. (2007)), which is based on Volterra series expansions as the one presented in this work. NOFRF give the response of the system to a given set of inputs, therefore the natural ensemble of forcing signals must be known to acquire useful information. This ensemble has to be compiled from (or modeled based on) fluid-structure interactions data. Comparing the present results with the information provided by NOFRF will be the objective of our next work. Proved that such analysis is helpful for the design of compliant robots, we will extend the analysis to biped and quadruped robots. In that situation the forcing comes from the ground reaction forces, that are generally not smooth (impacts) and contain a rich frequency spectrum. We consider those scenarios as highly complex problems to be addressed with a mature toolbox of methods.

Concluding, we have provided a method that allows the construction of a controller that would maximize the harvesting of environmental energy, under the circumstances defined herein. The method will be extended to cover more realistic situations than the simplified periodic forcing, including forcing with wider frequency spectrum. Additionally, the same methodology presented here can be used to verify results obtained using heuristic methods as AFO or neural networks. Our results can be easily validated and we invite our colleagues to do it.

## Acknowledgements

The authors would like to thank Cecilia Tapia Siles for sharing insights and information about passive swimmers.

#### Funding.

The research leading to these results has received funding from the European Community’s Seventh Framework Programme FP7/2007-2013-Challenge 2 –Cognitive Systems, Interaction, Robotics– under [grant agreement No 248311-AMARSi]. RB is founded by the CAPES Foundation, Ministry of Education of Brazil. MZ received funding for this work from the SNSF [project no. 122279] (From locomotion to cognition).

#### Author contributions.

JPC, RB and ZQL worked on the mathematics of the models and in the simulations. MZ designed and built the joint studied in this work. All authors collaborated for the production of the manuscript.

## References

- Ahlborn and Blake (2002) Ahlborn, B. K. and Blake, R. W. (2002). Walking and running at resonance. Zoology 105, 165–74.
- Alben (2009) Alben, S. (2009). On the swimming of a flexible body in a vortex street. J. Fluid Mech. 635, 27.
- Alexander (1990) Alexander, R. M. (1990). Three Uses for Springs in Legged Locomotion. Int. J. Robot. Res. 9, 53–61.
- Bejan and Marden (2006) Bejan, A. and Marden, J. H. (2006). Unifying constructal theory for scale effects in running, swimming and flying. J. Exp. Biol. 209, 238–48.
- Brilliant (1958) Brilliant, M. B. (1958). Theory of the analysis of nonlinear systems. Technical Report 345 MIT - Research laboratory of electronics.
- Buchli and Ijspeert (2008) Buchli, J. and Ijspeert, A. J. (2008). Self-organized adaptive legged locomotion in a compliant quadruped robot. Autonomous Robots 25, 331.
- Buchli et al. (2006) Buchli, J., Righetti, L. and Ijspeert, A. J. (2006). Engineering entrainment and adaptation in limit cycle systems : From biological inspiration to applications in robotics. Biol. Cybern. 95, 645–64.
- Collins et al. (2005) Collins, S., Ruina, A., Tedrake, R. and Wisse, M. (2005). Efficient bipedal robots based on passive-dynamic walkers. Science 307, 1082–5.
- Cottone et al. (2009) Cottone, F., Vocca, H. and Gammaitoni, L. (2009). Nonlinear Energy Harvesting. Phys. Rev. Lett. 102, 80601.
- Dickinson et al. (2000) Dickinson, M. H., Farley, C. T., Full, R. J., Koehl, M. A. R., Kram, R. and Lehman, S. (2000). How Animals Move: An Integrative View. Science 288, 100–106.
- Eaton (2002) Eaton, J. W. (2002). GNU Octave Manual. Network Theory Limited, http://www.octave.org.
- Eldredge and Pisani (2008) Eldredge, J. D. and Pisani, D. (2008). Passive locomotion of a simple articulated fish-like system in the wake of an obstacle. J. Fluid Mech. 607.
- Farley (1996) Farley, C. T. (1996). Leg stiffness and stride frequency in human running. J. Biomech. 29, 181–186.
- Farley et al. (1993) Farley, C. T., Glasheen, J. and McMahon, T. A. (1993). Running springs: speed and animal size. J. Exp. Biol. 185, 71–86.
- Ferris et al. (1998) Ferris, D. P., Louie, M. and Farley, C. T. (1998). Running in the real world: adjusting leg stiffness for different surfaces. Proc. R. Soc. Lond. [Biol] 265, 989–94.
- Harper et al. (1997) Harper, K., Berkemeier, M. and Grace, S. (1997). Decreasing the energy costs of swimming robots through passive elastic elements. In Proc. ICRA 1997 vol. 3, pp. 1839 –1844,.
- Holmes and Rand (1976) Holmes, P. and Rand, D. (1976). The bifurcations of duffing’s equation: An application of catastrophe theory. J. Sound Vib. 44, 237–253.
- Kanso and Newton (2009) Kanso, E. and Newton, P. K. (2009). Passive locomotion via normal-mode coupling in a submerged spring-mass system. J. Fluid Mech. 641, 205.
- Kerdok et al. (2002) Kerdok, A. E., Biewener, A. A., McMahon, T. A., Weyand, P. G. and Herr, H. M. (2002). Energetics and mechanics of human running on surfaces of different stiffnesses. J. Appl. Physiol. 92, 469–78.
- Kokshenev (2010) Kokshenev, V. B. (2010). Key principle of the efficient running, swimming, and flying. EPL (Europhysics Letters) 90, 48005.
- Lang and Billings (1996) Lang, Z. Q. and Billings, S. A. (1996). Output frequency characteristics of nonlinear systems. International Journal of Control 64, 1049–1067.
- Lang and Billings (2005) Lang, Z.-Q. and Billings, S. a. (2005). Energy transfer properties of non-linear systems in the frequency domain. International Journal of Control 78, 345–362.
- Lauder et al. (2007) Lauder, G. V., Anderson, E. J., Tangorra, J. and Madden, P. G. A. (2007). Fish biorobotics: kinematics and hydrodynamics of self-propulsion. J. Exp. Biol. 210, 2767–80.
- Liao (2007) Liao, J. C. (2007). A review of fish swimming mechanics and behaviour in altered flows. Philos. T. Roy. Soc. B 362, 1973–93.
- Liao et al. (2003) Liao, J. C., Beal, D. N., Lauder, G. V. and Triantafyllou, M. S. (2003). Fish exploiting vortices decrease muscle activity. Science 302, 1566–9.
- Luo and Han (1997) Luo, A. C. J. and Han, R. P. S. (1997). A quantitative stability and bifurcation analyses of the generalized duffing oscillator with strong nonlinearity. J. Franklin I. 334, 447–459.
- Maxima (2009) Maxima (2009). Computer Algebra System. Version 5.18.1. http://maxima.sourceforge.net/.
- McGeer (1990) McGeer, T. (1990). Passive Dynamic Walking. Int. J. Robot. Res. 9, 62–82.
- Moritz and Farley (2003) Moritz, C. T. and Farley, C. T. (2003). Human hopping on damped surfaces: strategies for adjusting leg mechanics. Proc R Soc Lond [Biol] 270, 1741–6.
- Peng et al. (2007) Peng, Z., Lang, Z.-Q. and Billings, S. (2007). Resonances and resonant frequencies for a class of nonlinear systems. J. Sound Vib. 300, 993–1014.
- Peyton-Jones and Billings (1989) Peyton-Jones, J. C. and Billings, S. A. (1989). Recursive algorithm for computing the frequency response of a class of non-linear difference equation models. International Journal of Control 50, 1925–1940.
- Roberts and Azizi (2011) Roberts, T. J. and Azizi, E. (2011). Flexible mechanisms: the diverse roles of biological springs in vertebrate movement. J. Exp. Biol. 214, 353–361.
- Rugh (1981) Rugh, W. J. (1981). Nonlinear system theory : the Volterra/Wiener approach. Baltimore, London, Johns Hopkins University Press.
- Schetzen (1980) Schetzen, M. (1980). The Volterra and Wiener theories of nonlinear systems. Wiley-Interscience.
- Thompson and Raiber (1989) Thompson, C. M. and Raiber, M. H. (1989). Passive Dynamic Running pp. 74–83. Berlin/Heidelberg: Springer-Verlag.
- Vakakis et al. (2009) Vakakis, A. F., Gendelman, O. V., Bergman, L. A., McFarland, D. M., Kerschen, G. and Lee, Y. S. (2009). Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems, vol. 156, of Solid mechanics and its applications. Springer Netherlands, Dordrecht.
- Ziegler et al. (2011) Ziegler, M., Hoffmann, M., Carbajal, J. P. and Pfeifer, R. (2011). Varying Body Stiffness for Aquatic Locomotion. In Proc. ICRA 2011, Shanghai, China.