A harmonic-based method for computing the stability of periodic oscillations of non-linear structural systems

In this paper, we present a validation on a practical example of a harmonic-based numerical method to determine the local stability of periodic solutions of dynamical systems. Based on Floquet theory and Fourier series expansion (Hill method), we propose a simple strategy to sort the relevant physical eigenvalues among the expanded numerical spectrum of the linear periodic system governing the perturbed solution. By mixing the Harmonic Balance Method and Asymptotic Numerical Method continuation technique with the developed Hill method, we obtain a purely-frequency based continuation tool able to compute the stability of the continued periodic solutions in a reduced computation time. This procedure is validated by considering an externally forced string and computing the complete bifurcation diagram with the stability of the periodic solutions. The particular coupled regimes are exhibited and found in excellent agreement with results of the literature, allowing a method validation.


Introduction
Determining the local stability of a dynamical system periodic solution is of primary interest in an engineering context since only stable solutions are experimentally encountered.Moreover, a change in the stability can lead to significant qualitative, and possibly dramatic, changes in the system response.A * Address all correspondence to this author.commonly-used method consists in computing the monodromy matrix of the linear periodic system governing the evolution of the perturbed solution and compare its eigenvalues (Floquet multipliers) moduli to one.Two well known numerical strategies are classically used.The first technique consists in a time-integration of the Jacobi matrix of the initial dynamical system over one period [1].In the second technique, the monodromy matrix is simply a by-product of a shooting continuation method [2].
The efficiency of those methods has been proved for years, but, inherently, the eigenvalues accuracy depends on the chosen time-step size, leading to possible long time computations.Furthermore, those time-domain approaches are not necessarily best suited for nonlinear periodic solution computed, for instance, with a harmonic balance method.In the latter, eigenvalues of the linear operator can naturally be computed with Hill's method also based on Fourier series expansion.This harmonic-based method deals with linear periodic systems and has been introduced by Hill one century ago for the determination of the lunar perigee [3].Since then, it has been applied to a wide range of physical problems like the computation of energy eigenvalues of the Schrödinger equation [4,5], the stability of limit cycles in electrical circuits [6] or the vibratory behavior of flexible rotating machines [7,8].However, although Hill's method often provides satisfactory results, the meaning of the computed eigenvalues is generally misunderstood and may lead to wrong results.Actually, this method is not of a trivial use since it requires to approx-imate the spectra of infinite-dimensional operators, i.e. to sort the most converged eigenvalues among all the numerical ones.
A recent work of the authors proposes a simple numerical strategy to sort these eigenvalues and therefore properly determinate the stability of periodic solutions of dynamical systems [9].A simple criterion, based on the energy distribution of the computed eigenvectors, is proposed and enables to extract the most converged eigenvalues (Floquet exponents) of the truncated problem, and therefore determinate the local stability of the solution.By mixing the standard harmonic balance method with the so-called asymptotic numerical method continuation technique based on the quadratic recast of the dynamical system [10], it becomes a straightforward and efficient procedure to continue the periodic solutions as well as their associated Floquet exponents.Whereas the case of a forced Duffing oscillator is proposed in [9] as a validation example, this article focuses on another one.
The main features of the developed harmonic-based continuation tool are first exposed.Then, the case of an externally excited string is considered.This problem is interesting since it extends the results obtained on the Duffing oscillator in two ways.Firstly, more than one mechanical degree of freedom is included.Secondly, since both transverse polarizations are considered in the model, one-to-one internal resonances are present, allowing a particularly complicated bifurcation diagram and showing several stable and unstable solutions as well as additional pitchfork bifurcations.Finally, this problem is qualitatively similar to the one of a perfect circular plate, where one to one resonances are observed between companion modes due to the symmetry of revolution.In the case of the string, out of plane motion are observed, which are equivalent to the travelling waves observed for the circular plate [11,12].

STABILITY OF PERIODIC SOLUTIONS General formulation
We consider continuous-time dynamical systems governed by the general equation where x is a N -dimensional state vector and f is a nonlinear Ndimensional vector field that depends on a control parameter λ.
In the following, f may explicitly depends on t (the system is non-autonomous) or not (the system is autonomous).We consider a periodic solution x 0 (t) of (1), with minimal period T , at the particular control parameter value λ = λ 0 .The stability of this periodic solution is studied by superimposing a small disturbance y (t): Substituting ( 2) into (1), assuming that f is at least twice continuously differentiable, expanding the result in a Taylor series about x 0 , and retaining only linear terms in the disturbance, we obtain: with The stability study of periodic solution x 0 consists in finding if the disturbance y(t) solution of (3) fades away or is amplified as t is increased.Since x 0 (t) is T -periodic in time, J(t) is also T -periodic in time considering its definition (4).Consequently, the system (3) is a linear system with periodic coefficients.The Floquet theory [13,14], used in the following, specifically deals with this kind of dynamical systems.

Floquet theory
The following developments are classical and can be found e.g. in the textbooks [1,15].The N -dimensional linear system (3) has N linearly independent solutions y n (t), so that any solution y(t) of (3) can be written: where c n are N constants that depend upon the initial conditions.First of all, we gather these fundamental solutions into a N × N matrix: By recalling that J(t) is T -periodic, one obtains so that y n (t + T ) verifies the system (6).Consequently, since the y n (t) are N independent fundamental solutions of (3), the y n (t + T ) can be expressed as linear combinations of the y n (t).
It thus exists a N × N constant matrix Φ, called the monodromy matrix, so that: This matrix can be used to study the stability of the periodic solution x 0 (t) since it maps a particular set of fundamental solution Y (t) at time t into their values at time t + T , thus defining a Poincaré map.The stability study is conducted by considering the eigenvalues ρ n of Φ, the so-called Floquet multipliers.
Considering again the fundamental solutions introduced above, every y n (t), for all n = 1, . . ., N , can be expressed in the Floquet form where p n (t + T ) = p n (t) is a T -periodic N -dimensional complex vector and α n is a complex number.Then, it follows from (9) and from the T -periodicity of p n that The above equation, compared to Eq. ( 8), proves that the α n are linked to the Floquet multipliers by, for all n = 1, . . ., N : The α n represent the so-called Floquet exponents.Whereas the Floquet multipliers are uniquely defined, the above equation shows that the α n are unique to within an additive integer multiple of 2iπ/T [14].This last result can also be viewed by replacing α n by α n + 2ikπ/T in (10).
Considering Eqs. ( 8), (9) or (10), the values of either the Floquet multipliers ρ n or the Floquet exponents α n can be used to determine the stability of periodic solutions.
If ℜ(α n ) < 0 (or |ρ n | < 1) for all n, all fundamental solutions y n converge toward zero as t is increased, so as disturbance y(t).The periodic solution is said to be asymptotically stable; If it exists a subscript n such that ℜ(α n ) > 0 (or |ρ n | > 1), the corresponding fundamental solutions increase exponentially, so as for some disturbances y(t).The periodic solution is unstable in this case.
The above statement stands for non-autonomous initial systems (1).If the initial system (1) is autonomous, one of the Floquet multipliers is always ρ n = 1.In this circumstances, only the set of remaining eigenvalues can provide information on the periodic solution stability.

Hill's method
Hill's method is a harmonic-based numerical approach used to determine the solutions of linear periodic systems like (3), which is used here to calculate the Floquet exponents α n .This method has been widely used together with perturbation methods to analytically obtain the stability of periodic solutions, for instance in the case of the Mathieu equation (see e.g.[16]).
The unknown periodic functions p n (t) introduced in ( 9) are expressed by the general Fourier series where the fundamental frequency ω reads ω = 2π/T and p k n are N dimensional complex vectors.By replacing (12) in the Floquet form ( 9), any fundamental solution y n (t) is an infinite sum of harmonic contributions Being T -periodic as well, the Jacobian matrix J (t) is also expanded in the infinite Fourier series where the J h are N × N matrices.By replacing the solution y n (t) and the Jacobian J (t) by their Fourier series in the linear periodic system (3), one obtains the following vector equation, for all n = 1, . . ., N : Since the sums have an infinite number of term, replacing superscript k by k − h in the left hand side term of the above equation does not changes it.It leads to, for all n = 1, . . ., N : (16) By applying the harmonic balance method to the above equation, i.e. by separately equating to zero each harmonic of (16) (for each value of k), the above equation can be rewritten in the following infinite-dimensional eigenproblem where H is the infinite-dimensional Hill matrix s is a complex number, q is an infinite dimensional vector and I is the identity matrix of appropriate size.By comparing ( 17) and ( 16), the eigenvalues and eigenvectors s and q are related to α n and p n by the following relations: valid for all l = 0, ±1, . . .± ∞ and n = 1, . . ., N .In the above equation, • T means the vector transpose of •.The above relations show that each Floquet exponent α n , for a given n, is associated to an infinite set of eigenvalues s l n , for all l = 0, ±1, . . .± ∞.Consequently, knowing the infinite set of s l n would enable to calculate the α n and then to assess the periodic solution stability.
However, for evident practical purposes and numerical computations, the infinite-dimensional problem ( 16)-( 18) is truncated to a finite dimension.Namely, matrix H is truncated to a N (2H + 1) × N (2H + 1) dimension, so that l = 0, ±1, . . ., ±H in (19).Let ŝl n and ql n denote the N (2H + 1) approximate eigensolutions of Eq. (17).If the α n are to be computed, two questions need to be solved: (1) do the (ŝ l n , ql n ) tend to (s l n , q l n ) as H tends to infinity and (2) for a given n, among all the computed ŝl n , which ones, as l is varied, can be used to compute α n .
Firstly investigated by Poincaré [17] more than one century ago, those questions were still arguing twenty years ago, in particular for the computation of Schrödinger equation [4].Recently, papers based on Fredholm theory provide rigorous proof of the convergence [18,19].
However, a delicate issue remains the choice of the eigensolutions of the truncated problem that can be used in a numerical procedure to evaluate the α n , n = 1, . . ., N .A simple idea consists in analysing the truncated problem eigenvectors distribution T and considering the one associated to l = 0: Since all the above developments are based on the assumption that the Fourier series ( 12) and ( 14) are convergent (p k n tends to zero as k tends to infinity), among all the computed eigenvectors, the components of the N ones associated with l = 0 (the q0 n ) have almost symmetric shapes since they involve only the p k n with |k| < H.For this reason, we assume that the q0 n converge faster to q 0 n than the others, for l = 0.With the above assumption, the numerical procedure to compute the N Floquet exponents α n is the following.The Hill matrix H, truncated to order H is assembled and its N (2H + 1) eigensolutions (ŝ l n , ql n ) computed.Then, the N ones associated to l = 0 are numerically selected by considering that the N eigenvectors q0 n are the N ones with the most symmetric shapes.The N associated eigenvalues are then the N Floquet exponents: n , n = 1, . . ., N .Their real part is then compared to 0 to assess the stability of the periodic solution.

THE HBM-ANM-HILL METHOD
One of the crucial point in the numerical method mentioned above is the expansion of the Jacobian J(t) in Fourier series (Eq.( 14)).Indeed, in the most general situation, the nonlinear evolution equation (1) shows nonlinearities of any kinds and the same is true for the Jacobian dependence on the periodic solution.Therefore, computation of its Fourier coefficients J h can be very cumbersome.
However, in the special case where the harmonic-balance method (HBM) combined with the asymptotic numerical method (ANM) is used to compute the Fourier coefficients of the periodic solution x 0 (as explained in [10]), the computation of the J h is rendered very efficient.The leading idea, also used for the ANM, is to systematically recast the dynamical system (1) in a quadratic polynomial form.From (1), the new quadratic system can be written as follows where c, l (•), q (•, •) are constant, linear and quadratic operators in z, respectively, and m (•) is a linear operator.The unkown vector z contains the original components of the state vector x and some new variables added to get the quadratic form.This procedure is fully explained in [10], which also shows that this quadratic recast can be applied to a large class of smooth systems with a few algebraic manipulations and a few additions of auxiliary variables.
A periodic solution of (21), of minimal period T = 2π/ω, is expanded into the Fourier series The HBM combined to the ANM continuation method enables to compute the z p 0 , p = −H, . . ., H and ω for various contiguous values of λ, detect possible bifurcations and compute different branches of solutions.For externally driven systems, the control parameter is often the forcing frequency λ = ω.Assuming the quadratic recast (21), for any λ = λ 0 , the Jacobian J (t) defined by ( 4) may be expressed in the form ) where J C , J L (•) and J Q (•, •) are constant, linear and quadratic N × N matrices with respect to z 0 .From then, it becomes straightforward to find the Fourier coefficients J h of J(t) as functions of those of z 0 (t).For each parameter λ 0 , the Hill matrix H is expressed in terms of the already computed z p 0 and ω and the stability of the periodic solution x 0 (t) is simply obtained following the previous Hill's method.
The selection of the N eigenvectors q 0 n among the N (2H + 1) eigenvectors of ( 17) is done by computing for each q l n , l = −H, . . ., H, n = 1, . . ., N , the median value of its components moduli.If q i is the i-th component of vector q of dimension Q = N (2H + 1), the median is i i|q i |/ i |q i |.The selected q 0 n are the N ones with the median value closest to Q/2, which are supposed to have the most symmetrical shape.Knowing those eigenvectors, the N corresponding eigenvalues are the α n , used for the stability study.

Theoretical model
The case of a forced hardening Duffing oscillator is presented in details in [9].Here, this example is extended to the case of an externally excited string.Many works have focused on non-linear vibrations of strings in the past and the possible destabilization of the in-plane motion giving birth to out of plane, possibly circular, motions, due to one to one internal resonances between the two transverse polarizations.Among others, the interested reader can refer to [20] for an early qualitative explanation of the out of plane motion, to [21,22] for a two degrees of freedom model analytically solved to obtain the solution branches as well as their stability, to [23] for experimental results.This problem shows qualitative similarities with nonlinear vibrations of circular plates, where the one to one internal resonance occurs between companion modes, giving birth to travelling waves [11,12].
Following [16], the nonlinear equations of motions of a vi-brating string can be written: In the above equations, x is the axial position, (v, w) are the two transverse displacements in the y and z directions, N 0 is the initial axial tension in the string, N is the axial tension added by the non-linear coupling with the transverse motion, f y and f z are the two components of the external force per unit length.L is the string length; ρ and E denote respectively the mass density and the Young's modulus of the string material; S and I are respectively the string cross section area and second moment.Finally, • and • ,x denote the first partial derivative of • with respect to time t and x, respectively.These equations are based on the classical von Kármán approximation that consists in keeping in the axial strain the first non-linear terms.The flexural rigidity effect has also been included.
In order to obtain a model with a minimum number of free parameters, Eqs.(24a-c) are transformed by introducing the following dimensionless parameters: with i = y, z. φ is the string cross section diameter.T 0 and f 0 are a time constant and a force constant, whose values (given below) are fixed so as to minimize the number of free parameters.One then obtains: where all the overbars have been dropped since all quantities are now dimensionless.The following parameters have been introduced: ε 1 is the ratio between flexural and axial stiffnesses and ε 2 is a measure of the non-linear terms amplitude.Those two parameters are often very small as compared to 1.
A solution to Eqs. ( 25a-c) is obtain by the following truncated modal expansion: where Φ k (x) is the k-th deformed shape of the string, solution of the linear part of Eqs.(25a-c) (i.e. with N = 0), identical for both y and z polarizations.It can be written: where Φ k (x) has been normalized so that and ω k is the k-th dimensionless angular frequency.By introducing Eqs. ( 26) into (25a-c), multiplying both sides of equations with Φ i (x) and using the orthogonality properties of the deformed shapes, one finally obtains, for all k = 1, . . ., K: In the above equations, a viscous damping term has been added, with modal damping factor ξ k for the k-th mode.It has also been chosen to introduce a small detuning between the frequencies of the two polarizations, by defining ωk = (1 + ε 3 )ω k , with ε 3 a given small parameter.The forcing terms in (29) are defined by: It must be noted that the same coefficient n 2 π 2 appears in (29a,b) and in (29c) because Φ k (0) = Φ k (1) = 0 and the special normalization of (28).

HBM-ANM-Hill formulation
Thanks to the above described mixed formulation (the unknowns are the two displacements v, w and the axial force N ), the problem to solve has only quadratic nonlinear terms and is naturally adapted to the HBM-ANM-Hill method [9,10].We focus on a harmonic forcing with frequency ω, on the y polarization only, concentrated at the middle of the string.It writes: The first step is to recast Eqs.(29a-c) into the following first order system: If we define the vector of unknowns z(t) = [p 1 . . .p K q 1 . . .q K r 1 . . .r K s 1 . . .s K N ] T , the mass, constant, linear and quadratic terms m( ż), c(ω, t), l(z, ω, t) and q(z, z, ω, t) are easily defined by comparing (31ae) to (21) (Since the system is non-autonomous, the bifurcation parameter is λ = ω).
The Jacobian used for the stability computation is the one associated to the dynamical system (31a-d), without equation (31e) which is static.It can be written in the form of (23): In the above definitions, all matrices are of size 2K × 2K, except J which size is 4K × 4K.0 is the 2K × 2K matrix filled with zeros and I is 2K × 2K identity matrix.

Results
The simulations presented below have been performed for a string of circular cross section, so that S = πφ 2 /4 and I =  Figure 1 shows the results of the HBM-ANM-Hill method applied to the string model of Eqs.(31).Only the first harmonic of the first modal coordinates p 1 (t) and q 1 (t) of both polarizations are shown as a function of the excitation frequency ω.The shapes of the branches are qualitatively conform to the ones obtained by perturbation methods for strings (see e.g.[16,21]) and for circular plates [11,12].By increasing the excitation frequency, a pitchfork (PF) bifurcation is first observed, giving birth to one unstable uncoupled branch (with no z-polarization, q 1 (t) = 0) and two stable coupled branches with both non-zero y and z polarizations.The stability has been determined by the Hill method.The saddle node (SN) bifurcations are precisely located at the turning points (the points with a vertical tangent), thus validating our method [24].The major interest of this numerical method, as compared to analytical perturbation methods, is that all unstable branches are automatically obtained and are found to be connected to one another.All simulations have been obtained by truncating the Fourier series to H = 5 harmonics.Another advantage of this numerical method is that the convergence, in term of harmonics number and modal truncation, can be precisely checked.Several simulations have been performed, by retaining up to H = 5 harmonics and K = 3 modes in the truncation.In this case, it has been found that H = 3 harmonics and K = 1 mode give a very accurate model, with almost zero truncation errors for this excitation amplitude range.This can be observed on Fig. 2, that shows that only the first harmonics of p 1 (t) and q 1 (t) have a significant amplitude.A conclusion is that no noticeable non-linear coupling is observed between the first mode (k = 1) and higher order modes (k > 1).
The stable coupled regime of the branches between PF and SN 1 bifurcations of Fig. 1 can be observed on Fig. 3. Two different branches (b 1 and b 2 of Fig. 1) share the same amplitude with different phases: the z-polarization of branch b 1 (qb 1 ) has a +π/2 phase difference as compared to the y-polarization (pb 1 ) whereas it is a −π/2 phase difference for branch b 2 (compare qb 2 and pb 2 ).This explains the almost circular orbits of Fig. 3 (they would be perfectly circular if both polarization had the same amplitude) and also that the string motion is clockwise for b 1 and counterclockwise for b 2 .
The performance of the Hill method for computing the stability is now addressed.Figure 4 enables to compare our stability results to those obtained with a classical time domain method.The latter consists in integrating (6) in time over one period with initial conditions Y (0) = I, where I is the N × N identity matrix.Then, the monodromy matrix is simply Φ = Y (T ) and the stability is obtained by comparing the Floquet multipliers ρ n moduli (its eigenvalues) to one.The time-integration is here performed by a fourth-order Runge-Kutta algorithm.As compared to the Hill method, one has to set an additional parameter: the number N s of sample points over the computed period for the numerical integration.Figure 4 shows |ρ n |, n = 1, . . .4, as a function of ω, computed with the Hill method as well as with the above described time-domain method, with various values of N s .It shows that the Hill method appears as a limit to the time domain method as N s is increased.Even N s = 6000 time integration points over the computed period are not sufficient to recover the right stability results and especially to place the saddle-node bifurcation points at the turning points of the solution branches.This large amount of integration points for the time-domain method is also related to an increased computational cost, as compared to the Hill method.Table 1 shows the needed running time to compute the continuation diagrams of Figure 4.The frequency method is by far the more efficient approach to obtain the continued stability map of the dynamical system.
Figure 5 highlights the relevance of our strategy to sort the most converged eigenvalues among all the computed ŝl n .Since for a given n, the real parts of the s l n are equal (ℜ(s l n ) = ℜ(α n ) for all l = 0, ±1, . . ., ±∞), it would be theoretically possible to use any of them to determine the periodic solution stability.However, among the truncated spectrum ŝl n , it may exists a set of eigenvalues with a real part that leads to erroneous stability predictions, for values of |l| close to H, which is clearly shown on Figure 5. Thus, a stability criterion based on the consideration of the real part of all eigenvalues would have led to erroneous results.Sorting the most symmetric eigenvectors ql n seems definitely a simple and correct method to follow the N more convergent physical eigenvalues among all the numerical ones.

CONCLUSION
In this article, a purely frequency-based numerical method to determine the local stability of periodic solutions of dynamical systems, based on computing the eigenvalues of the problem Hill matrix, has been presented.Its efficiency, in term of accuracy and computation time, has been shown when associated to a harmonic balance based continuation method.While the monodromy matrix is simply a by-product of a shooting continuation method in the time-domain approach, here, the Hill matrix is simply a by-product of a harmonic-balance continuation method in the frequency domain.The main interest of this method, which explains its efficiency, is that there is no need to switch from one domain to another for computing the stability.

FIGURE 2 . 1 FIGURE 3 .
FIGURE 2. AMPLITUDE OF THE H = 5 SIMULATED HAR-MONICS, FOR BOTH POLARIZATIONS WITH K = 3 MODES RE-TAINED IN THE TRUNCATION, AT THE SN1 POINT OF FIG. 1.

FIGURE 4 .
FIGURE 4. THE 4 FLOQUET MULTIPLIERS MODULI |ρn| AS FUNCTIONS OF THE EXCITATION FREQUENCY ω, CORRE-SPONDING TO THE SIMULATION OF FIG. 1 WITH K = 1 AND H = 5, WITHOUT THE MAIN UNSTABLE BRANCH.'•': ρn COMPUTED WITH THE HILL METHOD; '•': ρn COMPUTED WITH THE TIME DOMAIN DETERMINATION OF THE MON-ODROMY MATRIX, FOR VARIOUS NUMBERS Ns OF TIME SAM-PLES PER PERIOD (UP: Ns = 2000; BOTTOM: Ns = 6000).GRAY REGION CORRESPOND TO UNSTABLE PERIODIC SOLUTION.