Institut für Luftfahrtantriebe

Introduction to
Harmonic Balance
and application to
nonlinear vibrations

M. Krack
Head of the Structural Dynamics group
Institute of Aircraft Propulsion Systems
Department of Aerospace Engineering
University of Stuttgart
1
 Institut für Luftfahrtantriebe

What is Harmonic Balance?

Harmonic Balance is an approximation method for

computing periodic solutions of ordinary
differential equations (ODEs).

The dynamics of many systems (structures, fluids, electrical circuits, …) can be

described by ODEs.
The method is only interesting if we do not know the exact solution  nonlinear ODEs.
Periodic oscillations are often of primary technical relevance.

2
 Institut für Luftfahrtantriebe

Advantages of Harmonic Balance vs. numerical integration

• computational efficiency (usually a few orders of magnitude faster,

because: no integration of long transients, good approximation often already
for small number of variables)

• fewer problems with numerical instability or damping

• consideration of phase-lag boundary conditions

• computation also of physically unstable oscillations

3
 Institut für Luftfahrtantriebe

Outline of talk

• introductory example: Duffing oscillator

• generalization to nonlinear mechanical systems

• implementation in a simple Matlab tool NLvib, application

examples

• limitations, ongoing research

• summary, questions, discussion

4
 Institut für Luftfahrtantriebe

Duffing oscillator

Equation of motion of a single-degree-of-freedom

oscillator with cubic spring (Duffing oscillator),
with damping and harmonic forcing:

Goal: Compute periodic solutions with

Ansatz:

(only one harmonic here)

5
 Institut für Luftfahrtantriebe

Duffing oscillator: Harmonic Balance

Time derivatives of ansatz: With trigonometric identities

Expansion of nonlinear term

Substitute into Duffing equation and collect harmonics

We neglect harmonics with index higher than the ansatz (>1) and balance the
harmonics:
2 algebraic
equations
6 in 2 unknowns
 Institut für Luftfahrtantriebe

Duffing oscillator: Frequency response

Frequency response: We are interested in how the solution evolves with .

We can formulate this as a continuation problem, with as a free parameter:

solve

with respect to

with

in the interval

 We will later apply numerical solution and continuation methods.

 For this simple case, let us develop an analytical solution.
7
 Institut für Luftfahrtantriebe

Duffing oscillator: Frequency response

Transform to polar coordinates

Substitution into algebraic equations

Algebraic manipulations of Eq. (1)-(2)

It is easier to solve Eq. (5) for :

We can have zero, one
or two real-valued
8 solutions .
 Institut für Luftfahrtantriebe

Duffing oscillator: Frequency response

9
 Institut für Luftfahrtantriebe

Duffing oscillator: Frequency response

10
 Institut für Luftfahrtantriebe

Duffing oscillator: Frequency response

11
 Institut für Luftfahrtantriebe

Duffing oscillator: Frequency response

stable

unstable

stable

 For a given excitation frequency, multiple steady-state oscillations can be

12 reached, depending on the initial conditions.
 Institut für Luftfahrtantriebe

Multiple steady states in reality

model

13
 Institut für Luftfahrtantriebe

Duffing oscillator: Summary

Idea of Harmonic Balance

• find periodic solutions of ODEs

• ansatz: truncated Fourier series

• balancing of harmonics  algebraic equation system in Fourier coefficients

To be discussed further

• generalization to multiple harmonics

• systematic derivation of equation system We will focus here

on mechanical
• treatment of generic nonlinearities systems.
• numerical solution

14
 Institut für Luftfahrtantriebe

Outline of talk

• introductory example: Duffing oscillator

• generalization to nonlinear mechanical systems

• implementation in a simple Matlab tool NLvib, application

examples

• limitations, ongoing research

• summary, questions, discussion

15
 Institut für Luftfahrtantriebe

Nonlinear mechanical system

Equations of motion of a time-invariant mechanical system with periodic forcing

with

and

Goal: Compute periodic solutions with

Ansatz:
mathematically
with equivalent
representations of
We will use this
representation.
truncated Fourier
series

16
 Institut für Luftfahrtantriebe

Nonlinear mechanical system

Time derivatives of ansatz:

residual due to Fourier approximation

Substitute into equations of motion of possibly nonsmooth function

with

n(2H+1) algebraic with

Harmonic Balance:
equations
(truncating to order H and setting in n(2H+1)
17 the Fourier coefficients of the unknowns
residual to zero)
 Institut für Luftfahrtantriebe

Harmonic Balance as a Galerkin method

Weighted residual method weight function

Galerkin method: take ansatz functions as weight functions

In our case, the ansatz functions are . We thus obtain:

with

since
18
 Institut für Luftfahrtantriebe

Harmonic Balance vs. numerical integration

numerical integration Harmonic Balance

successive forward time stepping until solve algebraic equation system
transient approaches periodic state in Fourier coefficients

ansatz
periodic solution

20
 Institut für Luftfahrtantriebe

Nonlinear mechanical system

As both the force equilibrium and the unknowns

Interpretation of Harmonic Balance are Fourier coefficients, Harmonic Balance
is a frequency-domain method.
dynamic stiffness matrix

dynamic force
equilibrium in the
frequency domain

linear internal forces nonlinear external

internal forces forces

The true challenge is the calculation of the nonlinear force harmonics

!

formal definition:

21
 Institut für Luftfahrtantriebe

Treatment of nonlinear forces

Opportunities
• polynomial forces: closed formulation using convolution theorem

• piecewise polynomial (incl. piecewise linear) forces: determine transition times, then as above [PETR03,KRAC13]

• generic nonlinear forces: Alternating-Frequency-Time Scheme (AFT)

Alternating-Frequency-Time Scheme (“sampling”)

Number of samples per period:

• theoretical lower limit given by Nyquist-Shannon
theorem

• generic nonlinear, perhaps even non-smooth forces:

generously oversample if you can afford it
22
 Institut für Luftfahrtantriebe

Solution of the algebraic equation system

Task for Harmonic Balance, we have
solve
with respect to

with
Solvers
• pseudo-time solver
• Newton-like solver
• secant solver
• …
Idea of Newton method
Linearization of residual

Iteration step

23
 Institut für Luftfahrtantriebe

Solution of the algebraic equation system

Task for Harmonic Balance, we have
solve
with respect to

with
Solvers
• pseudo-time solver
• Newton-like solver
• secant solver
• …
Idea of Newton method
Linearization of residual

Iteration step

24
 Institut für Luftfahrtantriebe

Solution of the algebraic equation system

Task for Harmonic Balance, we have
solve
with respect to

with
Solvers
• pseudo-time solver
• Newton-like solver
• secant solver
• …
Idea of Newton method
Linearization of residual

Iteration step

25
 Institut für Luftfahrtantriebe

Solution of the algebraic equation system

Task for Harmonic Balance, we have
solve
with respect to

with
Solvers
• pseudo-time solver
• Newton-like solver
• secant solver
• …
Idea of Newton method
Linearization of residual

Iteration step

26
 Institut für Luftfahrtantriebe

Solution of the algebraic equation system

Task for Harmonic Balance, we have
solve
with respect to

with
Solvers
• pseudo-time solver
• Newton-like solver
• secant solver
• …
Idea of Newton method
Linearization of residual

Iteration step

27
 Institut für Luftfahrtantriebe

Solution of the algebraic equation system

Task for Harmonic Balance, we have
solve
with respect to

with
Solvers
• pseudo-time solver
• Newton-like solver
• secant solver
• …
Idea of Newton method
Linearization of residual

Iteration step

28
 Institut für Luftfahrtantriebe

Solution of the algebraic equation system

Task for Harmonic Balance, we have
solve
with respect to

with
Solvers
• pseudo-time solver
• Newton-like solver
• secant solver
• …
Idea of Newton method
Linearization of residual

Iteration step

29
 Institut für Luftfahrtantriebe

Solution of the algebraic equation system

Task for Harmonic Balance, we have
solve
with respect to

with
Solvers
• pseudo-time solver
• Newton-like solver
• secant solver
• …
Idea of Newton method
Linearization of residual

Iteration step

30
 Institut für Luftfahrtantriebe

Solution of the algebraic equation system

Task for Harmonic Balance, we have
solve
with respect to

with
Solvers
• pseudo-time solver
• Newton-like solver
• secant solver
• …
Idea of Newton method
Linearization of residual

Iteration step
 fast convergence near solution
 adjustments required for global convergence
31  analytical gradients greatly reduce computation time
 Institut für Luftfahrtantriebe

Computation of a solution branch (under variation of a free parameter)

Example: Frequency response analysis
 This is a job for a continuation method!
solve
with respect to Numerical path continuation: generate a
sequence of suitably spaced solution
with points within the given parameter range, and
go around turning points (if any).
in the interval

similar analyses
• analysis of self-excited limit cycles
• nonlinear modal analysis
• tracking of resonances
• tracking of bifurcation points

32
 Institut für Luftfahrtantriebe

Ingredients of a continuation method

Predictor
popular variants:
• tangent
• secant
• power series expansion

Parametrization: Quo vadis?

• avoid returning to same solution point or reversing direction on path
• in most cases: additional equation, free parameter as additional unknown
• popular variants: arc length; local; orthogonal

Corrector: apply solver!

Step length control:

• as small as necessary to ensure convergence and not overlook important
characteristics of the solution
• as large as possible to avoid spurious computational effort
33 • empirical rules depending on solver and tolerances, upper and lower bounds
 Institut für Luftfahrtantriebe

Outline of talk

• introductory example: Duffing oscillator

• generalization to nonlinear mechanical systems

• implementation in a simple Matlab tool NLvib, application

examples

• limitations, ongoing research

• summary, questions, discussion

35
 Institut für Luftfahrtantriebe

NLvib – a Matlab tool for nonlinear vibration analysis

Source code and documentation available via www.ila.uni-stuttgart.de/nlvib

Features
• Harmonic Balance with Alternating Frequency-Time Scheme
• Shooting, Newmark numerical time step integration
• solver: predictor-corrector continuation with Newton-like corrector (‘fsolve’),
analytical gradients
• nonlinearities
 local generic nonlinear elements
 (distributed) polynomial stiffness nonlinearity
• analysis types
 frequency response
 nonlinear modal analysis

36
 Institut für Luftfahrtantriebe

NLvib: Duffing oscillator

singleDOFoscillator_cubicSpring.m
% Parameters of the Duffing oscillator
mu = 1;
delta = 0.05;
kappa = 1;
gamma = 1;
P = .1;

% Analysis parameters
H = 1; % harmonic order
N = 2^7; % number of time samples per period
Om_s = .5; % start frequency
Om_e = 1.6; % end frequency

% Initial guess (from underlying linear system)

Q = (-Om_s^2*mu+1i*Om_s*delta+kappa)\P;
x0 = [0;real(Q);-imag(Q);zeros(2*(H-1),1)];

% Solve and continue w.r.t. Om

ds = .01; % Path continuation step size
Sopt = struct('jac','none'); % No analytical Jacobian provided here
X = solve_and_continue(x0,...
@(X) HB_residual_singleDOFcubicSpring(X,mu,delta,kappa,gamma,P,H,N),...
Om_s,Om_e,ds,Sopt);

% Determine excitation frequency and amplitude (magnitude of fundamental

% harmonic)
Om = X(end,:);
a37= sqrt(X(2,:).^2 + X(3,:).^2);
 Institut für Luftfahrtantriebe

NLvib: Duffing oscillator

singleDOFoscillator_cubicSpring.m
% Parameters of the Duffing oscillator
mu = 1;
delta = 0.05;
kappa = 1;
gamma = 1;
P = .1;

% Analysis parameters
H = 7; % harmonic order
N = 2^7; % number of time samples per period
Om_s = .5; % start frequency
Om_e = 1.6; % end frequency

% Initial guess (from underlying linear system)

Q = (-Om_s^2*mu+1i*Om_s*delta+kappa)\P;
x0 = [0;real(Q);-imag(Q);zeros(2*(H-1),1)];

% Solve and continue w.r.t. Om

ds = .01; % Path continuation step size
Sopt = struct('jac','none'); % No analytical Jacobian provided here
X = solve_and_continue(x0,...
@(X) HB_residual_singleDOFcubicSpring(X,mu,delta,kappa,gamma,P,H,N),...
Om_s,Om_e,ds,Sopt);

% Determine excitation frequency and amplitude (magnitude of fundamental

% harmonic)
Om = X(end,:);
a38= sqrt(X(2,:).^2 + X(3,:).^2);
 Institut für Luftfahrtantriebe

NLvib: Duffing oscillator

HB_residual_singleDOFcubicSpring.m
function R = HB_residual_singleDOFcubicSpring(X,mu,delta,kappa,gamma,P,H,N)
% Conversion of real-valued to complex-valued harmonics of generalized coordinates q
Q = [X(1);X(2:2:end-1)-1i*X(3:2:end-1)];

% Excitation frequency
Om = X(end);

% P is the fundamental harmonic of the external forcing

Fex = [0;P;zeros(H-1,1)];

% Specify time samples along period and apply inverse discrete Fourier transform
tau = (0:2*pi/N:2*pi-2*pi/N)';
qnl = real(exp(1i*tau*(0:H))*Q);

% Evaluate nonlinear force in the time domain

fnl = gamma*qnl.^3;

% Forward Discrete Fourier Transform, truncation, conversion to half spectrum notation

Fnlc = fft(fnl)/N;
Fnl = [real(Fnlc(1));2*Fnlc(2:H+1)];

% Dynamic force equilibrium (complex-valued)

Rc = ( -((0:H)'*Om).^2 * mu + 1i*(0:H)'*Om * delta + kappa ).*Q+Fnl-Fex;

% Conversion from complex-valued to real-valued residual

R = [real(Rc(1));real(Rc(2:end));-imag(Rc(2:end))];

39
 Institut für Luftfahrtantriebe

NLvib: friction-damped beam

Equation of motion

with the nonlinear force

and the regularized Coulomb dry friction law

computation times:
• HB, whole frequency response: ~1.3 s
40 • numerical integration, single frequency: ~30 s
 Institut für Luftfahrtantriebe

Vibration prediction of bladed disks with friction joints

static FEA CFD

Component Mode Synthesis

OrAgL
Nonlinear Vibration Analysis
• frequency response analysis
• nonlinear modal analysis
measurements • flutter analysis

validation
• resonances • deflection shapes
42 • damping • dynamic stress
 Institut für Luftfahrtantriebe

Frequency response of a bladed disk with shroud contact

vibration level
rotational speed

For comprehensive rig and

engine validation (with a much
more realistic model), attend the
presentation of paper
GT2018-75186.

43
 Institut für Luftfahrtantriebe

Outline of talk

• introductory example: Duffing oscillator

• generalization to nonlinear mechanical systems

• implementation in a simple Matlab tool NLvib, application

examples

• limitations, ongoing research

• summary, questions, discussion

44
 Institut für Luftfahrtantriebe

Harmonic Balance has to main limitations

Limitation 1: Only periodic oscillations  no quasi-periodic, broadband or chaotic ones

multi-frequency HB ansatz

Extensions
• for quasi-periodic oscillations:
multi-frequency HB [SCHI06,KRAC16],
adjusted HB [GUSK12])
• for broadband or chaotic oscillations:
ongoing research

45 quasi-periodic oscillations on an invariant torus

 Institut für Luftfahrtantriebe

Harmonic Balance has to main limitations

Limitation 2: Harmonic base functions  Gibbs phenomenon near discontinuities

Fourier
approximation reference
signal

examples:
• impacts
• stick-slip transitions

time
Extensions
• enrichment by non-smooth base functions (wavelet balance [JONE15,KIM03])
• numerical integration of non-smooth states (mixed-Shooting-HB [SCHR16])
46
 Institut für Luftfahrtantriebe

Alternatives to Harmonic Balance Shooting problem

solve
Periodic oscillations
• Fourier-based alternatives with respect to
 trigonometric collocation with

 time spectral method where are determined

by forward numerical integration
• Shooting methods

• …

General oscillations:
Forward numerical integration

47
 Institut für Luftfahrtantriebe

Alternatives to Harmonic Balance Shooting problem

solve
Periodic oscillations
• Fourier-based alternatives with respect to
 trigonometric collocation with

 time spectral method where are determined

by forward numerical integration
• Shooting methods

• …

General oscillations:
Forward numerical integration

48
 Institut für Luftfahrtantriebe

Alternatives to Harmonic Balance Shooting problem

solve
Periodic oscillations
• Fourier-based alternatives with respect to
 trigonometric collocation with

 time spectral method where are determined

by forward numerical integration
• Shooting methods

• …

General oscillations:
Forward numerical integration

49
 Institut für Luftfahrtantriebe

Alternatives to Harmonic Balance Shooting problem

solve
Periodic oscillations
• Fourier-based alternatives with respect to
 trigonometric collocation with

 time spectral method where are determined

by forward numerical integration
• Shooting methods

• …

General oscillations:
Forward numerical integration

50
 Institut für Luftfahrtantriebe

Ongoing research on Harmonic Balance

• robust and efficient methods for branching behavior

• reliable and efficient methods for stability assessment (local and global!)

• multi-physical problems

• improvements for non-smooth problems

51
 Institut für Luftfahrtantriebe

Outline of talk

• introductory example: Duffing oscillator

• generalization to nonlinear mechanical systems

• implementation in a simple Matlab tool NLvib, application

examples

• limitations, ongoing research

• summary, questions, discussion

52
 Institut für Luftfahrtantriebe

Summary

Harmonic Balance is a numerical method for the efficient computation

of periodic solutions of nonlinear ordinary differential equations.

It can be interpreted as Galerkin method. It yields an algebraic equation

system, which can be solved using e.g. Newton-like methods and

numerical path continuation.

A relatively simple Harmonic Balance code is implemented in the

Matlab tool NLvib, available via www.ila.uni-stuttgart.de/nlvib.

53
 Institut für Luftfahrtantriebe

Further reading
Harmonic Balance, Alternating Frequency-Time Scheme
[CAME89] Cameron, T. M.; Griffin, J. H.: 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.
[CARD94] Cardona, A.; Coune, T.; Lerusse, A.; Geradin, M.: A Multiharmonic Method for Non-Linear Vibration Analysis. Int.
J. Numer. Meth. Engng 37(9):1593-1608, 1994.
[URAB65] Urabe, M.: Galerkin's Procedure for Nonlinear Periodic Systems. Archive for Rational Mechanics and Analysis
20(2):120-152, 1965.

Treatment of particular nonlinearities, Asymptotic Numerical Method

[COCH09] Cochelin, B.; Vergez, C.: A High Order Purely Frequency-Based Harmonic Balance Formulation for Continuation
of Periodic Solutions. Journal of Sound and Vibration 324(1-2):243-262, 2009.
[KRAC13] Krack, M.; Panning-von Scheidt, L.; Wallaschek, J.: A High-Order Harmonic Balance Method for Systems With
Distinct States. Journal of Sound and Vibration 332(21):5476-5488, 2013.
[PETR03] Petrov, E. P.; Ewins, D. J.: Analytical Formulation of Friction Interface Elements for Analysis of Nonlinear Multi-
Harmonic Vibrations of Bladed Disks. Journal of Turbomachinery 125(2):364-371, 2003.

Computation of quasi-periodic oscillations with Harmonic Balance

[GUSK12] Guskov, M.; Thouverez, F.: Harmonic Balance-Based Approach for Quasi-Periodic Motions and Stability
Analysis. Journal of Vibration and Acoustics 134(3):11pp, 2012.
[KRAC16] Krack, M.; Panning-von Scheidt, L.; Wallaschek, J.: On the interaction of multiple traveling wave modes in the
flutter vibrations of friction-damped tuned bladed disks. J. Eng. Gas Turbines Power 139(4):9pp, 2016.
[SCHI06] Schilder, F.; Vogt, W.; Schreiber, S.; Osinga, H. M.: Fourier Methods for Quasi-Periodic Oscillations. Int. J.
Numer. Meth. Engng 67(5):629-671, 2006.

Computation of non-smooth periodic oscillations

[JONE15] Jones, S.; Legrand, M.: Forced vibrations of a turbine blade undergoing regularized unilateral contact conditions
through the wavelet balance method. Int. J. Numer. Meth. Engng 101(5):351-374, 2015.
[KIM03] Kim, W.-J; Perkins, N. C.: Harmonic Balance/Galerkin Method for Non-Smooth Dynamic Systems. Journal of
Sound and Vibration 261(2):213-224, 2003.
[SCHR16] Schreyer, F.; Leine, R. I.: A Mixed Shooting – Harmonic Balance Method for unilaterally constrained mechanical
systems. Archive of Mechanical Engineering 63:298-313, 2016.
54
 Institut für Luftfahrtantriebe

Questions?

Topics for discussion?

55
 Institut für Luftfahrtantriebe

Appendix: Overview of basic examples in NLvib

name n HB Shooting FRF NMA

01_Duffing 1 o - o -
02_twoDOFoscillator_cubicSpring 2 o o o -
03_twoDOFoscillator_unilateralSpring 2 o o o -
04_twoDOFoscillator_cubicSpring_NM 2 o o - o
05_twoDOFoscillator_tanhDryFriction_NM 2 o o - o
06_twoSprings_geometricNonlinearity 2 o - o o
07_multiDOFoscillator_multipleNonlinearities 3 o - o -
08_beam_tanhDryFriction 16 o o o -
09_beam_cubicSpring_NM 38 o - - o
n: number of degrees of freedom
Run times depend on your computing environment, but should not HB: Harmonic Balance
exceed a minute per example for a standard computer (2017). FRF: nonlinear frequency response analysis
NMA: nonlinear modal analysis
56
 Institut für Luftfahrtantriebe

Appendix: Definition of Mechanical Systems in NLvib

For simplicity, the code comes with

Equations of motion harmonic forcing. Note that you can
easily generalize the external force to
multiple harmonics (which is actually a
good exercise to get familiar with the code).
Local nonlinear elements

Matlab syntax force direction

% Define properties
M = ... % n x n matrix
D = ... % n x n matrix
K = ... % n x n matrix
Fex1 = ...% n x 1 vector
w1 = ... % n x 1 vector
... specific properties and values

% Define nonlinear elements

nonlinear_elements{1} = struct(‘type’,...,‘force_direction’,w1,[‘p1’,v1,‘p2’,v2,...]);
nonlinear_elements{2} = ...

% Define mechanical system

mySystem = MechanicalSystem(M,D,K,nonlinear_elements,Fex1);

57
 Institut für Luftfahrtantriebe

Appendix: Local nonlinear elements in NLvib

Some of the already available local nonlinear elements type

 You can easily add new nonlinearities analogous to the already available ones.
Implementation of nonlinear elements in HB_residual.m (analogous in shooting_residual.m)
...
%% Evaluate nonlinear force in time domain
switch lower(nonlinear_elements{nl}.type)
analytical derivative
case 'cubicspring'
fnl = nonlinear_elements{nl}.stiffness*qnl.^3;
dfnl = ...
case 'mynewnonlinearity'
fnl = ... Tip: When you add a new nonlinearity, run the solver with ‘jac’ option
dfnl = ... ‘none’. If everything is working properly, you can later accelerate the code
... by providing (correct!) analytical derivatives. If you encounter
58 convergence problems at that point, your derivatives are wrong.
 Institut für Luftfahrtantriebe

Appendix: A chain of oscillators in NLvib

Structural matrices

Matlab syntax
% Define properties
mi = ... % vector with length n
ki = ... % vector with length n+1
di = ... % vector with length n+1
Fex1 = ...% n x 1 vector
...

% Define nonlinear elements

nonlinear_elements{1} = ...

% Define chain of oscillators

myChain = ChainOfOscillators(mi,di,ki,nonlinear_elements,Fex1);
59
 Institut für Luftfahrtantriebe

Appendix: An FE model of an Euler-Bernoulli beam in NLvib

Coordinates

matrix containing columns of the identity matrix, so that the generalized

coordinates are compatible with the boundary conditions (BCs)

Matlab syntax
% Define properties
...
BCs = ‘clamped-free’; % example with clamping on the left and free end on the right;
% pinned is also possible; arbitrary combinations are allowed
n_nod = ... % positive integer

% Define beam (rectangular cross section)

myBeam = FE_EulerBernoulliBeam(len,height,thickness,E,rho,BCs,n_nod);

% Add external forcing (add_forcing works in an additive way)

inode = ... % node index
dof = ... % degree of freedom specifier (‘rot’ or ‘trans’)
Fex1 = ... % complex-valued scalar Note that you can add non-grounded elements (as
add_forcing(myBeam,inode,dof,Fex1); in the general MechanicalSystem case, but you will
have to set up the force direction manually.
% Add nonlinear attachment (only grounded nonlinear elements for transversal DOF)
inode_nl = ... % see above
dof_nl = ... % see above
60 add_nonlinear_attachment(myBeam,inode,dof,type,[‘p1’,v1,‘p2’,v2,...]);
 Institut für Luftfahrtantriebe

Appendix: An FE model of an elastic rod in NLvib

Coordinates

matrix containing columns of the identity matrix, so that the generalized

coordinates are compatible with the boundary conditions (BCs)

Matlab syntax
% Define properties
...
BCs = ‘pinned-free’; % example with pinning on the left and free end on the right;
% arbitrary combinations are allowed
n_nod = ... % positive integer

% Define rod
myRod = FE_ElasticRod(len,A,E,rho,BCs,n_nod);

% Add external forcing (add_forcing works in an additive way)

inode = ... % node index
Fex1 = ... % complex-valued scalar Note that you can add non-grounded elements (as
in the general MechanicalSystem case, but you will
add_forcing(myRod,inode,Fex1);
have to set up the force direction manually.

% Add nonlinear attachment (only grounded nonlinear elements)

inode_nl = ... % see above
add_nonlinear_attachment(myRod,inode,type,[‘p1’,v1,‘p2’,v2,...]);

61
 Institut für Luftfahrtantriebe

Appendix: A system with polynomial stiffness in NLvib

Equations of motion of MechanicalSystem, but with nonlinear force

Matlab syntax
% Define properties
p = ... % Nz x n vector of nonnegative integers
E = ... % Nz x n vector of real-valued coefficients
...

% Define system
myPolyStiffSys = System_with_PolynomialStiffnessNonlinearity(M,D,K,p,E,Fex1);

Example: system with geometrical nonlinearity

color scheme

62
 Institut für Luftfahrtantriebe

Appendix: Some practical hints on using NLvib

Strongly simplify your problem first and then successively increase

complexity!
• Always analyze the linearized problem first.
 Do the system matrices have the expected dimensions, symmetries, eigenvalues?

 Derive a suitable initial guess for the nonlinear analysis.

 Derive reference values for linear scaling (‘Dscale’).

• Always start the nonlinear HB analysis with H=1.

• Then increase H successively until the results converge (do not waste resources by
setting it unreasonably high).

63
 Institut für Luftfahrtantriebe

Appendix: Some practical hints on using NLvib

What shall I do if I encounter one or more of the following difficulties?

a) Initial guess not within basin of attraction.
– start in ‘more linear’ regime
– improve the initial guess (e.g. from a suitable linearization or numerical integration)
– if analytical gradients are used, validate them (or run with ‘jac’ parameter set to ‘none’)
b) No convergence during continuation.
 ensure suitable scaling variables (‘Dscale’ vector) and residual funcations
 reduce step length parameter ‘ds’
 ensure numerical path continuation is activated (‘flag’ set to ‘on’ (default))
 increase AFT scheme sampling rate
 analytical gradients, cf. above
c) The computation time is very large.
 scaling, cf. above
 increase step length parameter ‘ds’
 use (correct!) analytical gradients
 lower your expectations 

64
 Institut für Luftfahrtantriebe

Appendix: solve_and_continue in NLvib

Problem statement
solve

with respect to

in the interval

initial guess residual function R(X) interval limits nominal step size
Matlab syntax
continuation
[X,Solinfo,Sol] = solve_and_continue(x0,fun_residual,lam_s,lam_e,ds [,Sopt, ... options structure
fun_postprocess,opt_fsolve]); (optional)

array of structures returned by postprocessing fsolve options

columns are postprocessing functions functions F(X) (optional) (optional)
solution points
Solinfo.FCtotal total function evaluation count
Solinfo.ctime total computation time

per solution step

Solinfo.iEx fsolve exit flag
Solinfo.NIT number of iterations
65 Solinfo.FC function count
 Institut für Luftfahrtantriebe

Appendix: solve_and_continue in NLvib

Most common continuation options (Sopt)
.flag flag whether actual continuation is performed or trivial (sequential)
continuation is employed (default: 1)
.predictor tangent or secant predictors can be specified ['tangent'|'secant']
(default: ‘tangent’)
.parametrization parametrization constraint ['arc_length'|'pseudo_arc_length'|'local'|
'orthogonal'] (default: ‘arc_length’)

.dsmin minimum step size (default: ds/5)

.dsmax maximum step size (default: ds*5)
.stepadapt flag whether step size should be automatically adjusted (default: 1;
recommended if Sopt.flag = 1)
.stepmax maximum number of steps before termination
.termination_criterion cell array of functions (X) returning logic scalar 1 for termination

.jac flag whether analytically provided residual derivatives (Jacobian)

should be used (default), or a finite difference approximation should
be computed (.jac = ‘none’)
.Dscale linear scaling to be applied to vector X

66
 Institut für Luftfahrtantriebe

Appendix: solve_and_continue in NLvib

Why should you apply linear scaling?

Example problem with solution

Rescaled problem with solution

where the scaling matrix

attempts to achieve similar orders of magnitude among the new variables

This suggest that one should scale with the (not a priori known) solution. In practice, one
can achieve good results with typical values for the respective variable, derived e.g.
from a solution of a linearized problem.

67
 Institut für Luftfahrtantriebe

Appendix: solve_and_continue in NLvib

With the suggested scaling, the condition number of the Jacobian change:

A high condition number is likely to cause convergence problems within the Newton method in
the presence of numerical imprecisions. To illustrate this, we apply Matlab’s fsolve to both
problems. Numerical imprecisions are
introduced by letting fsolve approximate
the Jacobian using finite differences.

Without scaling, the solver has apparent

difficulties.

68 iteration number
 Institut für Luftfahrtantriebe

Appendix: solve_and_continue in NLvib

From the viewpoint of the solver, the problem is stretched, making it hard to numerically
approximate the true solution in the unscaled variable space.

iteration
path

iteration
path

69
 Institut für Luftfahrtantriebe

Appendix: Frequency response analysis with NLvib

Harmonic Balance formulation

solve

where

with respect to

in the interval

70
 Institut für Luftfahrtantriebe

Appendix: Frequency response analysis with NLvib

Shooting formulation

solve
positive real-valued scalar

with respect to Rationale behind scaling of residual:

achieve similar orders of magnitude for
quite different vibration levels. Otherwise
in the interval the solver might misinterpret e.g. a small
value as a converged residual.

where are determined by forward

numerical integration

71
 Institut für Luftfahrtantriebe

Appendix: Nonlinear modal analysis with NLvib

Harmonic Balance formulation

positive real-valued scalar

solve
amplitude normalization

phase normalization

where

with respect to

in the interval

Rationale behind scaling of residual: achieve similar orders of magnitude of typical values. Otherwise the dynamic
force equilibrium or the normalization conditions would have unreasonably strong weight, which could have a negative
72 influence the convergence of the solver.
 Institut für Luftfahrtantriebe

Appendix: Nonlinear modal analysis with NLvib

Shooting formulation

solve

with respect to

in the interval amplitude

normalization
where are only with
phase
where are determined by forward numerical normalization
integration

positive real-valued scalar

Rationale behind scaling of residual: achieve similar orders of magnitude for quite different vibration levels.
Otherwise the solver might misinterpret e.g. a small value as a converged residual.

73
 Institut für Luftfahrtantriebe

Appendix: Numerical time step integration (for details see textbooks e.g. [GER15])

The purpose of numerical time step integration is to approximate the solution of an ordinary
differential equation, from given initial values , up to a given end time. Most of the
methods are based on finite difference approximations with respect to time and yield a
quadrature formula governing the values of the unknown states at the next time level .

Some methods directly deal with the second-order differential equation of motion (Newmark, Hilbert-
Hughes-Taylor), other methods require the re-formulation to first-order.
For explicit methods, the quadrature formula can be brought into an form where the unknown
states at the next time level are determined by simplify evaluating a function once. Implicit
methods require one to solve an algebraic equation system at each time level.
The method is (numerically) stable if the states remain finite for a finite step size. There are
conditionally stable methods, which are only stable for sufficiently small step size, and
unconditionally stable methods.
The approximation error depends, among others, on the quadrature formula and the step size.
Important accuracy measures are numerical damping (non-physical decrease of energy), numerical
dispersion (depending of error on contributing oscillation frequencies). Some numerical damping of higher
frequencies can be desirable, particularly if their dynamics is not correctly modeled due to
e.g. finite spatial discretization.
74 [GER15] Geradin, M.; Rixen, D.J.: Mechanical Vibrations, Wiley, 2015
 Institut für Luftfahrtantriebe

Appendix: Newmark method implemented in NLvib

Equation of motion evaluated at end of a time level

Time discretization (constant average acceleration Newmark scheme, see e.g. [GER15])
Note that in the actual
implementation, we
introduce a time
normalization, see
next slide.

Substitution of (4) and (5) into (1) gives an implicit equation in displacements at end of time
step,

which is solved using Newton iterations with Cholesky factorization of the Jacobian.
75 [GER15] Geradin, M.; Rixen, D.J.: Mechanical Vibrations, Wiley, 2015
 Institut für Luftfahrtantriebe

Appendix: Time normalization (for implementation of Newmark method) in NLvib

Normalized time

Reformulation of time derivatives

Equations of motion in non-normalized and normalized time

76
 Institut für Luftfahrtantriebe

Appendix: An almost foolproof approach to analytical gradients (as

used in NLvib)

function [R,dR] = my_function(X,param1,param2)

% Define auxiliary variables from input variables X
x1 = X(1);
x2 = X(2);
Om = X(3);

% Initialize derivative of auxiliary variables (‘Seeding’)

dX = eye(length(X));
dx1 = dX(1,:);
dx2 = dX(2,:);
dOm = dX(3,:);

% Operate on auxiliary variables, determine derivatives in each step using elementary

calculus
z = x1*Om^2;
dz = dx1*Om^2 + x1*2*Om*dOm;
R = z/x2 – x2;
dR = dz/x2 – z/x2^2*dx2 – dx2;

77