<figure><img src="image_1.png" /></figure>

HAL
open science

# Vibration of mechanical systems with geometric nonlinearities: Solving Harmonic Balance Equations with Groebner basis and continuations methods

Aurélien Grolet, Fabrice Thouverez

► To cite this version:

Aurélien Grolet, Fabrice Thouverez. Vibration of mechanical systems with geometric nonlinearities: Solving Harmonic Balance Equations with Groebner basis and continuations methods. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722052

HAL Id: hal-01722052

https://hal.science/hal-01722052

Submitted on 2 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Vibration of mechanical systems with geometric nonlinearities: Solving Harmonic Balance Equations with Groebner basis and continuations methods

Aurélien GROLET<fnref n="1" />, Fabrice THOUVEREZ<fnref n="1" />

**Résumé** — This paper is devoted to the study of vibration of mechanical systems with geometric nonlinearities. After applying the harmonic balance method, one has to solve a system of multivariate polynomial equations whose solutions give the frequency component of the possible steady states. Computing solutions of HBM (harmonic balance method) equations for a particular frequency is usually done iteratively with Newton-Raphson methods resulting in only one solution depending on the first iterate. Here we intend to compute all solutions of HBM equations by using methods based on Groebner basis computation. This approach allows to reduce the complete system to an unique polynomial equation in one variable driving all solution of the problème. This way the procedure avoid the computation of multiple paths, as in homotopie techniques, which can be quickly to much time consuming. In addition continuation methods are used to extend the solution for multiple values of the frequency parameter. We apply those methods to a simple forced nonlinear dynamic system and give a representation of the multiple states possible versus frequency.

**Mots clés** — nonlinear vibration, harmonic balance method, Groebner basis, numeric continuation.

# 1 Introduction

This paper presents a method for computing multiple steady state solutions of dynamical systems with polynomial nonlinearities. The harmonic balance method (HBM) [1, 2, 3, 4] is used to transform the nonlinear differential equation into a set of nonlinear algebraic equation which solutions give the frequency component of the possible steady states. As we consider polynomial nonlinearities, the HBM equations consists in a system of polynomial equations.

Solving HBM equations is usually done for a fixed frequency by mean of Newton-like algorithms giving only one solution. The solution can then be extended using continuation algorithms [5] which follow the solution with relation to a parameter (frequency, curvilinear abscissa, damping, forces, ...). Bifurcations can be monitored, and several new branches of solutions can be computed from the bifurcation points depending on the bifurcation types [5, 6]. However, if there exists branches of solutions disconnected (in the sense that they are not arising from bifurcation points of the previous branches), or if the bifurcation points are out of the studied frequency range, continuation algorithms fail to detect them and new tool are needed to solve the polynomial system of equation induced by the HBM.

Finding all solutions of a multivariate polynomial system is a subject of particular interest since, systems of polynomial equations arise in many research fields. Most of the time, the system is solved by mean of Newton-Raphson algorithms [5], which gives only one solution depending on the initial iterate. If one wants to compute all solutions by mean of Newton Raphson, one has to discretized the initial condition space and try all initial iterates. Sometimes the algorithm do not converge, or the attraction pool for a particular solution is too narrow compared to the discretization used, leading to the fact that some solutions can be missed. Some authors has proposed [7] to modify the Newton Raphson iteration in a way such that all previously computed solution will not be computed again. Indeed, dividing each equation by a factor of the type $||x - x^*||$, where $x^*$ is a previously computed solution, introduces a singularity at each previously computed solution so that those solutions will not be computed again (supposing the singularity is strong enough). However, the more solution computed the more singularity

<footnote n="1">¹ LTDS, aurelien.grolet@ec-lyon.fr</footnote>

<footnote n="2">² LTDS, fabrice.thouverez@ec-lyon.fr</footnote>
in the system, which can lead to poor convergence rate or no convergence at all so that the modified algorithm outputs only a subset of all solutions.

Homotopy methods are an alternative which can actually compute all solutions of polynomial equation systems. The linear homotopy (see e.g [6]) relies on liking the roots of a polynomial system P to the known roots of a simpler system Q. The roots of Q are used as starting points for following the paths given by $H(x,t) = tP(x) + (1-t)Q(x)$ ($t \in [0, 1]$) by continuation methods on parameter t. System Q generally have more roots than system P so that some path will be divergent and their continuation will results only in a loss of time. Designing the polynomial Q may be more or less complex depending on the method used.

Finally algebraic methods such as Groebner basis [8], Rational univariate representation [9, 10] or matrix methods [11, 12, 13] exploit the relation between variables induced by the set of polynomial equation. Those method are mostly based on Groebner basis and computation in quotient algebra.

In this paper we propose to use Groebner basis to solve the set of HBM equations. We note that some attempts to use Groebner basis in the domain of structural mechanics has already been proposed in [14, 15] for the determination of multiple static equilibriums. Here, Groebner basis theory is used as a mean to derive solutions of polynomial equations through matrix method and eigenvalue problems. In order to reduce the computation time for the search of multiple solution, we propose to first find the multiple solutions of the undamped system and second recover the solutions of the damped system by applying a continuation on the damping parameter.

The paper is organized a follows : Section 2 described the application of the harmonic balance method, the methods for solving multivariate polynomial equations, and the continuation on damping. In section 3 we apply our method to a simple cyclic system with cubic nonlinearities.

# 2 Multiple steady states of systems with polynomial nonlinearities

## 2.1 Harmonic Balance Method

This section presents the application of the harmonic balance method [1, 2, 3, 4]. Let's consider the n dof nonlinear dynamic system given by :

$$M\ddot{u} + C\dot{u} + Ku + H(x, \dot{x}) = F(t) \quad (1)$$

where $u(t)$ is the vector of unknown of size $n$, M, C and K are respectively the mass, damping and stiffness matrices, $F(t)$ is the vector of excitation force which is assumed to be periodic with period $T = \frac{\omega}{2\pi}$, and finally $H(x, \dot{x})$ is the vector of nonlinear forces which are assumed to be polynomial. The HBM consist in searching the solution under the form of a truncated Fourier series up to H harmonics as follows :

$$u(t) = a_0 + \sum_{k=1}^{H} a_k \cos(k\omega t) + b_k \sin(k\omega t) \quad (2)$$

Then, Substituting Eq.(2) in Eq.(1) and projecting the resulting equation over the truncated Fourier basis leads to the following set of nonlinear algebraic equations :

$$Z(\omega)x + \tilde{H}(x) - \tilde{F} = P(x, \omega) = 0 \quad (3)$$

where $x = [a_0^T, a_1^T, b_1^T, ..., a_H^T, b_H^T]^T$ is the vector of unknown of size $n_h = n(2H + 1)$, Z is the matrix of dynamic stiffness computed as $Z = \text{diag}(K, (Z_k)_{1 \le k \le H})$ where $Z_k$ is given by

$$Z_k = \begin{bmatrix} K - (k\omega)^2 M & k\omega C \\ -k\omega C & K - (k\omega)^2 M \end{bmatrix} \quad (4)$$

Finally $\tilde{F}$ and $\tilde{H}(x)$ corresponds respectively to the excitation force and the nonlinear force in the frequency domain. Since we assumed the nonlinear force in the time domain H to be polynomial, the resulting nonlinear force in the frequency domain $\tilde{H}$ are also under polynomial form so that Eq.(3) is a system of $n_h$ polynomial equations. Fixing the frequency and solving the algebraic system in Eq.(3) then give the possible steady states of the system for a particular frequency $\omega$.
When the system is undamped, the HBM approximation only needs to take into account cosine terms, so that Eq.(3) simplifies to

$$Z_u(\omega)a + \tilde{H}(a) - \tilde{F} = P_u(a, \omega) = 0 \quad (5)$$

where $a = [a_0^T, a_1^T, \dots, a_H^T]^T$ is the vector of unknown of size $n_h = n(H + 1)$, and $Z_u = \text{diag}(K, (K - (k\omega)^2 M)_{1\le k\le H})$.

## 2.2 Continuation methods on the damping parameter

This section described the continuation method used to transform a solution of the undamped system in Eq.(5) into a solution of the damped system in Eq.(3). Let's consider a solution $a^{(u)}$ of system in Eq.(5) for a particular frequency $\omega_u$, this solution will be used as a starting point for the continuation on damping parameter. Substituting $C$ by $\varepsilon C$ in Eq.(4) where $\varepsilon$ is a new parameter going from 0 to 1, results in the following form for Eq.(3) :

$$Z(\omega, \varepsilon)x + \tilde{H}(x) - \tilde{F} = P(x, \omega, \varepsilon) = 0 \quad (6)$$

For $\varepsilon = 0$ the system is undamped and the vector $x_0 = [a^{(m)}, 0_H]$ (of size $(n(2H + 1))$) is solution of Eq.(6), and for $\varepsilon = 1$ the system is fully damped. At first sight, a simple way to apply continuation would be to set $\omega$ to the fixed value $\omega_u$ and apply a sequential continuation on parameter $\varepsilon$. However, as the damping increase with $\varepsilon$, the frequency of the solution may be shifted so that the continuation stop prematurely. In order to overcome this drawback, $\omega$ is also considered as a variable and a continuation process with two parameters $\omega$ and $\varepsilon$ based on a predictor/corrector is applied. The prediction/correction process is then applied until $\varepsilon = 1$, giving a solution $(x^{(d)}, \omega^{(d)})$ of the damped system.

During the continuation a step length management is applied by monitoring the number of iteration needed in the correction process, if this number is greater than a selected threshold, the length of the predicted vector is decrease and the correction process starts again. For some starting points, it may happen that the solution disappears as the damping increase leading to an unfinished continuation. This situation can be checked by monitoring for exemple the step length size, if this length is to small and $\varepsilon \ne 1$ then we assume that the solution has disappeared due to too high damping. As a by product of the continuation method, if the solution has disappeared, we have the maximum damping value for this solution to exist (i.e the value of $\varepsilon$ for the last converged correction). This information can be useful in design strategies for removing critical solutions.

## 2.3 Solving multivariate polynomial systems with groebner basis

This section presents the use of Groebner basis for solving multivariate polynomial system $P$ such as the HBM equations in Eqs.(3,5). First, note that polynomial system $P$ defines an ideal $I$ over the ring $\mathbb{C}[x]$ of complex valued polynomials in the variable $x$. The concept of Groebner basis for polynomial ideals has been introduced in the 60's by B.Buchberger during its PhD under the direction of professor W.Grobner [16]. A Grobner basis is a particular kind of generating subset $G$ for the ideal $I$ with the additional property that the division of any element of ring $\mathbb{K}[x]$ by $G$ gives a unique remainder [17]. This definition extends the concept of euclidean division for multivariate polynomials, and allows for calculation in the quotient space $\mathbb{Q} = \mathbb{K}[x]/I$. Buchberger has proposed an algorithm for computing Groebner basis in his thesis which has been improved several time, particularly by J.C. Faugère [18, 19, 20].

When computed with an elimination order, Groebner basis are under a special triangular which can be used directly for finding zeros of $P$. Indeed, the last equation will be an univariate polynomial which can be solved to give values for $x_n$. Then, the penultimate equation (a two-variate polynomials) is solved giving values for $x_{n-1}$, etc... The procedure is then repeated until values for $x_1$ has been found. However it is known that computing Groebner basis relatively to an elimination order can be very time consuming as the number of variables grows.

In the case of non-elimination order (such as total degree ordering), Groebner basis are computed faster, but they no longer are in triangular form. However, they can be used to compute a monomial basis $\mathcal{B}$ for the quotient space $\mathbb{Q}$, which in turns is used to build special matrices called multiplication matrices [11, 8]. It can be shown that the eigenvalues of those multiplication matrices are related to the
zeros of the polynomial system [11]. Solving polynomial system then corresponds to solve eigenvalue problems (or equivalently univariate polynomials, like with elimination order), and simply read the solution in selected eigenvectors. The detailed description of this method can be found in [11, 12, 13]. In this study, we use symbolic computation softwares that already include packages related to Groebner basis and multiplication matrices computation such as the 'Groebner' package of the Maple software, so that the computation of multiple solutions will corresponds to solving eigenvalue problems.

## 3 Application on a simple exemple

This section presents the application of methods described in section 2 on a simple dynamic system with cubic nonlinearities. Let's consider the following n dof cyclic dynamic system :

$$m\ddot{u}_i + c\dot{u}_i + (k + 2k_c)u_i - k_c u_{i-1} - k_c u_{i+1} + k_{nl}u_i^3 = F_i(t), \quad 1 \le i \le n \quad (7)$$

with convention $u_{n+1} = u_1$. Parameters m, c, k are respectively the mass, damping and stiffness coefficient, $k_c$ and $k_{nl}$ are the coupling stiffness and the nonlinear stiffness coefficient, and finally $F_i(t)$ is the excitation force which is assumed to be under the form $F_i(t) = f_i \cos(\omega t)$.

In the remaining of this paper we will consider a system with only n = 6 dofs. This number can seem pretty small, but it indeed sufficient for the system to develop a large number of solution. Moreover we will consider a special case of forcing where all dof are forced with the same amplitude, i.e. $f_i = f$ for $1 \le i \le n$. Finally, the numerical value used through this example are the following :

$$m = 1, c = 0.1, k = 1, k_c = 1, k_{nl} = 1, f = -1 \quad (8)$$

In accordance with the results of a linear analysis, we choose to study our system in the range of frequency [0, 1]. The particular frequency for the search of multiple solution will be set to 0.47 Hz.

### 3.1 Solution of the undamped system

In order to find solutions of the undamped system, the HBM is applied with only one harmonics resulting in a set of $n_h = 6$ algebraic equations with variable $a = [a_1, ..., a_6]^T$ given by :

$$(k + 2k_c - \omega^2 m)a_i - k_c a_{i-1} - k_c a_{i+1} + \frac{3}{4}k_{nl}a_i^3 - f_i = 0, \quad 1 \le i \le n_h \quad (9)$$

Solutions of Eq.(9) are then sought using methods described in section 2.3. In this example, it turns out that Groebner basis computation is not necessary since the system in Eq.(9) is already in a Groebner basis form for any total degree monomial ordering, and the dimension of the quotient space $\mathbb{Q}$ is maximum being equals to $3^{n_h}$ [21, 8]. Taking this facts into account, we only need to build the multiplication matrices $M_{a_i}$ of size $3^{n_h} = 729$, and solve the different eigenvalue problems. Multiplication matrices are build for the particular frequency $\omega_m = 3 \text{ rad.s}^{-1}$ with help of the Maple v15 software. Each multiplication matrices has a total of 217 real eigenvalues. We noted that the eigenvectors corresponding to the real eigenvalues are not perfectly symmetric, this fact may be due round off error when converting rational numbers into floating point numbers in the multiplication matrices. The resulting eigenvectors are thus refined using the HBM equation of the undamped system with one harmonic. Merging the different results leads to a total of 254 different solutions for the undamped system. In order to avoid unnecessary continuation du to cyclic symmetry, we choose to sort the solution by family, a family $F_a$ being defined by $b \in F_a$ either if $b = a$, or $b$ is a circular permutation of $a$. Applying this sorting method leads to only 46 families of solution for the undamped system, corresponding to motions going from weakly to strongly localized. Finally, the 46 solutions of the undamped system are refined using an HBM approximation with 5 harmonics.

### 3.2 From undamped to damped

Here the 46 solutions of the undamped system are used as starting points for the continuation algorithm described in section 2.2 in order to derive solutions for the damped system. The continuation is
carried out on the damped HBM system derived with $H = 5$ harmonics. From the 46 solutions, only 23 converged to a solution of the fully damped system. The method of damping continuation is illustrated in Fig.1 : first a solution of the undamped system is computed (green curve), then the continuation on the damping parameter is applied to derive a solution of the damped system (black curve) and finally the solution of the damped system is continued (magenta curve).

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Illustration of the continuation on the damping parameter</figcaption></figure>

## 3.3 Continuation of the damped solutions

Finally the 23 converged solutions for the damped system are continued using a classical arc-length continuation. In order to avoid unnecessary continuations, for each starting point we test if this starting point belong to a previously computed curve of solution by computing the distance of the starting point to the curve. The same strategy is used for terminating the continuation of closed curve solutions (i.e. we monitor the distance between the current point and the starting point). Using this starting point management, only 6 curves of solutions are continued over the 23 starting points, leading to 6 families of solutions for the damped system. In addition, stability of solution is evaluated using the so called monodromy matrix [5].

The resulting FrF diagram is depicted on Fig.2 in a classical frequency-amplitude plot, and on Fig.3 in an energy-frequency plot.

We see that there exists a small interval of stability for some different solutions in the range [0.44, 0.48] Hz. The stability of a particular stable solution of the 3rd family (red solution) is illustrated by numerical integration in Fig.4. After integrating the system over 5000 periods with starting point taken from the 3rd family at $\omega = 0.45Hz$ one can still see the accordance between the HBM solution and the integrated solution.

# 4 Conclusion

This paper propose a method to compute multiple steady states of nonlinear dynamic equation with polynomial nonlinearities. The harmonic balance method is used to transform the set of nonlinear differential equations into polynomial equations, which are solved by mean of computation in the induced quotient algebra. In order to decrease the computation time, we proposed to work on the undamped system first. This has several advantages : (i) decrease the number of variables, and (ii) put the polynomial
<figure><img src="image_3.png" /><figcaption>Fig. 2 – Frequency Amplitude plot of the 6 family of solution for the damped system (◦:stable, ·:unstable).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Energy Frequency plot of the 6 family of solution for the damped system (◦:stable, ·:unstable).</figcaption></figure>

system under Groebner basis form so that no Groebner basis computation are needed. Once the undamped system has been totally solved, solutions of the damped system are recovered by using a continuation method on the damping parameter leading to a full frequency response diagram. The proposed method has been applied on a simple example.

In this paper solutions of polynomial systems are computed by solving eigenvalue problems related to multiplication matrices. Clearly the major limitation of this method is the exponential growth of the multiplication matrix size with the number of variables, which considerably limits the number of harmo-
<figure><img src="image_5.png" /><figcaption>Fig. 4 – Comparison between temporal integration and HBM solutions for a stable point of the 3rd family at ω = 0.45Hz (· : temporal integration ○ : HBM)</figcaption></figure>

nics to be used in the harmonic balance method. To maintain the number of variables sufficiently small
one should consider model reduction.

Références

[1] G. Groll and D.J. Ewins. The harmonic balance method with arc-length continuation in rotor stator contact problems. *Journal of sound and vibration*, 241 (2) :223–233, 2001.

[2] T.M. Cameron, J.H. Griffin, and P. Grant. An alternating frequency time domain method for calculating the steady state response of nonlinear dynamic systems. *Journal of applied mechanics*, 56 (1) :149–154, 1989.

[3] Cochelin B. and Vergez C. A high order purely frquency based harmonic balance formulation for continuation of periodic solutions. *Journal of sound and vibration*, 324 :243–262, 2009.

[4] Grolet A. and Thouverez F. On a new harmonic selection technique for harmonic balance method. *Mechanical Systems and Signal Processing*, 30:43–60, 2012.

[5] A.H. Nayfey and B. Balanchandran. *Applied nonlinear dynamics*. Wiley-Interscience, 1995.

[6] Sarrouy E., Grolet A., and Thouverez F. Global and bifurcation analysis of a structure with cyclic symmetry. *International Journal of Nonlinear Mechanics*, 46:727–737, 2011.

[7] C.S Liu, W. Yeih, and Atluri Satya. An enhanced fictuous time integration method for nonlinear algebraic equation with multiple solutions. *CMESy*, 59(3):301–323, 2010.

[8] D. Cox, J. Little, and D. O’Shea. *Ideals, varieties and algorithms (2nd edition)*. Springer, 1996.

[9] F. Rouillier. Solving zero dimensinal polynomial system through the rational univariate representation. AAECC, 9 :433–461, 1999.

[10] K. Ouchi, J. Kayser, and J.M. Rojas. The exact rational univariate representation and its application. *preprint elsevier*, 2003.

[11] H.J. Stetter W.Auzinger. An elimination algorithm for the computation of all zeros of a system of multivariate polynomial equations.

[12] Moller H.M and Tenberg R. Multivariate polynomial system solving using intersections of eigenspaces. *Journal of symbolic computation*, 32 :513–531, 2001.

[13] Moller H.M and Stetter H.J. Multivariate polynomial equations with multiple zeros solved by matrix eigen-problems. *Numerische Mathematik*, 70 :311–329, 1995.

[14] VanDervort R.W. *Geometrically nonlinear analysis of rectangular orthotropic plates using Groebner basis*. PhD thesis, Tennessee Technological University, 2009.

[15] Shanmugasundaram A. *An application of the method of Groebner bases to a geometrically nonlinear free vibration analysis of composite plates*. PhD thesis, Tennessee Technological University, 2009.
[16] B. Buchberger. *An algorithm for finding the basis element of residue class ring of a zero dimensional polynomial ideal*. PhD thesis, J. Kepler University, 1965.

[17] B. buchberger. Groebner basis, a short introduction for systems theorists.

[18] J.C. Faugere, P. Gianni, and D. Mora. Efficient computation of zero-dimensional groebner basis by change of ordering. *Journal of symbolic computation*, 16 :329–344, 1993.

[19] J.C. Faugere. A new efficient algorithm for computing groebner basis (f4). 2002.

[20] J.C. Faugere. A new efficient algorithm for computing groebner basis without reduction to zero (f5). 2004.

[21] Hanzon B. and Jibetean D. Global minimization of a multivariate polynomial using matrix methods. *Journal of global optimisation*, 27 :1–23, 2003.