Comparison of ANM and Predictor-Corrector Method to Continue Solutions of Harmonic Balance Equations

In this work we apply and compare two numerical path continuation algorithms for solving algebraic equations arising when applying the Harmonic Balance Method to compute periodic regimes of nonlinear dynamical systems. The first algorithm relies on a predictor-corrector scheme and an Alternating Frequency-Time approach. This algorithm can be applied directly also to non-analytic nonlinearities. The second algorithm relies on a high-order Taylor series expansion of the solution path (the so-called Asymptotic Numerical Method) and can be formulated entirely in the frequency domain. The series expansion can be viewed as a high-order predictor equipped with inherent error estimation capabilities, which permits to avoid correction steps. The second algorithm is limited to analytic nonlinearities, and typically additional variables need to be introduced to cast the equation system into a form that permits the efficient computation of the required high-order derivatives. We apply the algorithms to selected vibration problems involving mechanical systems with polynomial stiffness, dry friction and unilateral contact nonlinearities. We assess the influence of the algorithmic parameters of both methods to draw a picture of their differences and similarities. We analyze the computational performance in detail, to identify bottlenecks of the two methods.


Introduction
Harmonic Balance (HB) permits the efficient approximation of periodic solutions of nonlinear ordinary differential equations. Often, we want to determine the solution as a function of a free parameter. To this end, numerical path continuation is commonly applied. Continuation provides higher robustness and efficiency as compared to simply computing the solution for a sequence of equidistant parameter values. Perhaps more importantly, continuation allows us to overcome turning points and therefore to capture multiple solutions for a single parameter value. The purpose of this work is to enlighten the strengths and weaknesses of two popular methods for continuing solutions of the HB equations. HB uses a truncated Fourier series as approximation ansatz. Substitution into the ordinary differential equation system gives a residual term, which is then made orthogonal to the Fourier basis function (Fourier-Galerkin projection). This corresponds to requiring that the Fourier coefficients of the residual are zero, for those harmonics retained in the ansatz. When we apply this method to the equation of motion of a mechanical system, with generalized coordinates q q q, subjected to periodic forcing, we obtain the algebraic equation system, Herein, Ω is the angular excitation frequency,q q q is the vector of Fourier coefficients of the approximation for q q q(t), andf f f nl andf f f ex are the Fourier coefficients of the nonlinear forces and the external forces, respectively. S S S is dynamical stiffness matrix representing linear internal forces proportional to q q q,q q q, andq q q, where overdot denotes derivative with respect to time. S S S is block diagonal (different harmonics are decoupled in the linear case).

Alternating-Frequency-Time scheme with Predictor-Corrector Continuation (AFT-PreCo)
The AFT scheme approximates the Fourier coefficientsf f f nl by sampling q q q andq q q at a certain number of equidistant time instants N along one period, evaluating f f f nl (q q q,q q q) in the time domain, and applying the discrete Fourier transform to the samples of the nonlinear force. Hence, the nonlinear forces are represented by N samples. N affects the accuracy of the procedure and has to be selected properly. For polynomial forces, the sampling procedure is exact beyond a certain N.
Otherwise,f f f nl is an approximation and contains aliasing errors. The great advantage of the AFT procedure is that it can be easily applied to smooth as well as non-smooth nonlinear forces. The downside of the AFT scheme is that it is not straight-forward to efficiently compute derivatives of order higher than one. The AFT scheme is commonly combined with PreCo continuation. Suppose we want to determineq q q for a range of the free parameter Ω . Starting from an initial solution X X X 0 with X X X := q q q T Ω T , we can do a prediction, where X X X 1 is the unit tangent to the solution path at the point X X X 0 and α is the step length of the prediction. X X X pre PreCo will usually not satisfy Eq. (1). Hence, correction steps are necessary, commonly in the form of Newton iterations, to improve the estimate until ||r r r(X X X)|| < ε for a given tolerance ε. Note that the equation system is under-determined. This is often resolved by considering an additional equation that ensures that the new solution point lies at a certain distance α from X X X 0 (arc-length continuation). In the PreCo, the step length is a crucial parameter: A too small step length leads to spurious computation effort, a too large step length may lead to many and costly correction iterations or even divergence. The step length is commonly adjusted automatically with the intent to achieve a desired number of correction iterations. However, it is advisable to choose reasonable upper and lower bounds, to avoid overlooking important features of the solution path and getting stuck near branching points, respectively. These bounds are usually chosen based on experience.

Classical Harmonic Balance with Asymptotic Numerical Method (cHB-ANM)
The ANM can be interpreted as high-order predictor (order p), which is a Taylor series expansion of the solution path, around X X X 0 , in the arc length α, X X X pre ANM = X X X 0 + αX X X 1 + . . . + α p X X X p .
The step size α is determined such that the estimated error, based on the range of utility of the Taylor series, remains below the tolerance ε, with the intent to completely avoid correction iterations. For the ANM, the algebraic equation system is commonly recast into quadratic form, i. e., with only constant, linear and bilinear (quadratic) terms in the unknowns. The Fourier coefficients of the nonlinear terms can in this case be expressed analytically without having to resort to sampling, which corresponds to classical HB. Note that all variables, including the nonlinear forces are represented by 2H + 1 Fourier coefficients. The quadratic form permits the successive computation of X X X k up to high order. Only a single Jacobian matrix has to be factorized to determine all X X X k , k > 0, per expansion point X X X 0 [2]. This contrasts the PreCo where several of such factorizations are required, one per Newton iteration, to get to the next point on the solution branch. The cHB-ANM cannot directly deal with non-smooth nonlinearities. These have to be approximated by appropriate analytic functions (regularization). Moreover, the equations have to be brought into quadratic form, which generally requires introducing auxiliary variables and equations. For details on how the cHB-ANM can be implemented in a computationally efficient and user-friendly way, we refer to [3]. Fig. 1 (a) depicts the frequency response of the Duffing oscillator. For the same H, both methods have the same number of unknowns (Fourier coefficients of q up to order H). However, the results differ as the nonlinear terms are represented in a different way. For the cubic term of the Duffing oscillator, it can be shown that the sampling procedure yields the Fourier coefficients of the nonlinear forces up to order H without aliasing error for N ≥ 4H + 1. For Fig. 1 (a) N was selected accordingly. For 2H cHB + 1 = N AFT , the nonlinear terms are represented with the same number of Fourier coefficients in cHB as samples in the AFT, so that a similar accuracy can be expected. Hence, AFT N AFT = 5 agrees well with cHB H ANM = 2, N AFT = 13 with H ANM = 6 and so on. Fig. 1 (b) shows the results for a single degree of freedom (SDOF) oscillator with elastic dry friction nonlinearity. Without regularization, only AFT can be applied, and a reasonable approximation is achieved with H = 1, N = 45. To regularize, an auxiliary variable describing the force had to be introduced, which is not well-represented with H = 1. Again, for H = 22 = (45 − 1)/2 both methods yield similar results for the regularized model. The remaining deviation to the reference is due to regularization. To compute the results depicted in Fig. 1 (b), AFT-PreCo H = 1, N = 45 (non-smooth) was about 25 times faster than cHB-ANM H = 22. We observed that as the regularization becomes steeper, the number of harmonics in the cHB-ANM has to be increased to continue the whole branch. Sometimes this value seems too high to be useful. In Fig. 1 (c) the computation effort is shown for a modal model (n modes) of a beam with geometric nonlinearity. The results suggest that cHB-ANM becomes more efficient for larger numbers of DOFs. The number of solution points was larger in the case of the cHB-ANM, which overall lead to approximately the same total number of Jacobian factorizations for the entire solution path. Apparently, the AFT calculation of the nonlinear forces and their derivatives is less efficient than the evaluation of the nonlinear terms and calculation of higher-order Taylor series coefficients within the cHB-ANM. Further investigations are in progress to gain a better understanding of this. Fig. 1: Representative results: a) Duffing oscillator; b) SDOF system with elastic dry friction; c) geometrically nonlinear beam with n coupled modal coordinates. A RMS is the root-mean square of the generalized coordinate. More results on different benchmark systems and further analyses of opportunities and limitations of both methods will be published in a journal article.

Conclusions
The results indicate that the AFT-PreCo is better suited for non-smooth nonlinearities such as stick-slip friction or unilateral constraints. The cHB-ANM is highly efficient for low-order polynomial nonlinearities and problems with a large number of DOFs, as in the case of geometrically nonlinear finite element models. It is expected that the cHB-ANM would greatly benefit from selecting a separate, higher truncation order for the auxiliary variables.