Example C4: Nonlinear Normal Modes of a 2DoF Hysteretic Oscillator using the Extended Periodic Motion Concept

The Extended Periodic Motion Concept

Nonlinear normal modes is a very succinct methodology for describing the dynamics of nonlinear systems. It leads to descriptions of the dynamics in terms of amplitude-dependent natural frequency, daming factor, and mode-shapes for each mode. The Extended Periodic Motion Concept is a very popular methodology for capturing these numerically. The idea here is simple - it iteratively looks for what would be a linear self excitation (in terms of a negative viscous factor) that must be supplied to a system so that it has a periodic solution as an invariant set in its unforced state.

Mathematically the system we will be solving is:

\[\mx{M} \ddot{\vc{u}} + \mx{C} \dot{\vc{u}} + \mx{K} \vc{u} + \vc{f_{nl}} = \xi \mx{M} \dot{\vc{u}}\]

where $\xi$ controls the strength of the self-excitation. Since the system is always solved by the trivial solution, an amplitude constraint is often imposed. In order to simplify calculations later on, the first harmonic amplitude of the degrees of freedom is constrained as:

\[\vc{U}_{c1}^T \mx{M} \vc{U}_{c1} + \vc{U}_{s1}^T \mx{M} \vc{U}_{s1} = a^2\]

where $\vc{U}_{c1} = \frac{\Omega}{\pi} \int_0^{2\pi/\Omega} \vc{u} \cos\Omega t dt$ and $\vc{U}_{s1} = \frac{\Omega}{\pi} \int_0^{2\pi/\Omega} \vc{u} \sin\Omega t dt$. (a) is the user-specified amplitude, which serves as the continuation parameter in this context. In addition to this, there is still a phase indeterminacy present in the system. So we choose an arbitrary DoF (which is expected to be non-zero in the mode under consideration), and set its cosine harmonic to be zero (other constraints are also possible). With this, we will have a square system (same number of unknowns as equations) and it can be solved to obtain the amplitude-dependent nonlinear modal characteristics of the system.

The EPMCRESFUN! routine in MDOFUTILS implements the residue function for EPMC that can be used in tandem with CONTINUATION to obtain nonlinear modal characteristics. This example will illustrate this on the frictional 2DoF oscillator under consideration in Examples C2 and C3. The steps followed are:

First the system is setup (along with the nonlinearity, see descriptions in Example C2.
Next we conduct linear modal analysis on the linearized system. This is important to set the initial guess for the EPMC run next.
Now we setup the Harmonic Balance parameters and then test out the EPMC residue routine using the linearized modes estimated above.
Finally we conduct the contuation over the modal amplitude (a) and then visualize the results.

Preamble: Load Packages

using GLMakie
using LinearAlgebra
using SparseArrays
using ForwardDiff
using NonlinearSolve
using DSP

using juliajim.HARMONIC
using juliajim.CONTINUATION
using juliajim.MDOFUTILS

System Setup

MCK matrices and nonlinearity setup: identical to Example C2.

M = collect(1.0I(2));
K = [2. -1.;-1. 2.];
C = 0.01*M+0.001*K;

mdl = MDOFGEN(M, C, K);

# Nonlinearity
kt = 1.0;
fs = 1.0;
fnl = (t,u,up,fp)-> if all(abs.(fp+kt*(u-up)).<fs)
    return fp+kt*(u-up), kt*ones(size(u)), -kt*ones(size(u)), ones(size(u)); else
        return fs*sign.(fp+kt*(u-up)), zeros(size(u)), zeros(size(u)), zeros(size(u));
end
L = [0.0 1.0];

mdl = ADDNL(mdl, :Hyst, fnl, L);

Linear Modal Analysis

Remember that in the limit of small amplitude, the Jenkins element acts like a linear spring. So we add this additional stiffness to the linearized stiffness matrix and conduct the modal analysis (eigen decomposition). The mode shapes, natural frequencies, and damping are saved since these will be used next as the initial guess for EPMC.

Knl0 = L'kt*L;
K0 = mdl.K+Knl0;
D, V = eigen(K0, mdl.M)
W0s = sqrt.(D);
Xis = diag(V'mdl.C*V)

mi = 1;  # Mode of interest

Setup HB

Setting up parameters for Harmonic Balance: identical to Example C2.

h = 0:5;
N = 256;
t = (0:N-1)*2π/N;

Nhc = sum(all(h.==0, dims=2) + 2any(h.!=0, dims=2));

_, _, zinds, rinds, iinds = HINDS(mdl.Ndofs, h)
Fl = zeros(Nhc*mdl.Ndofs);
Fl[rinds[1]] = 1.0;

E, dEdw = HARMONICSTIFFNESS(mdl.M, mdl.C, mdl.K, W0s[mi], h);

U0 = zeros(mdl.Ndofs*Nhc);
U0[iinds[1:mdl.Ndofs]] = V[:, mi];

Test EPMC Residue

The ordering of unknowns that will be sent to the EPMCRESFUN! routine is (\begin{bmatrix} \vc{u}^T&\omega&\xi&\log_{10}(a) \end{bmatrix}^T), where the last parameter is the log of the amplitude (base 10). The log-scaling is done on the amplitude to help with numerical conditioning issues.

In the following we setup everything and test it to ensure that the residue can be used to converge to the solution.

R = zeros(mdl.Ndofs*Nhc+2);
dRdU = zeros(mdl.Ndofs*Nhc+2, mdl.Ndofs*Nhc+2);
dRda = zeros(mdl.Ndofs*Nhc+2);

EPMCRESFUN!([U0; W0s[mi]; Xis[mi]; -2], mdl, Fl, h, N; R=R, dRdUwx=dRdU, dRda=dRda)

A0 = -2.0;

U0 = zeros(mdl.Ndofs*Nhc);
U0[iinds[1:mdl.Ndofs]] = V[:, mi];

fun = NonlinearFunction((r,uwx,p)->EPMCRESFUN!([uwx;p], mdl, Fl, h, N; R=r),
    jac=(J,uwx,p)->EPMCRESFUN!([uwx;p], mdl, Fl, h, N; dRdUwx=J),
    paramjac=(Jp,uwx,p)->EPMCRESFUN!([uwx;p], mdl, Fl, h, N; dRda=Jp));

prob = NonlinearProblem(fun, [U0;W0s[mi];Xis[mi]], A0);
sol = solve(prob, show_trace=Val(true), maxiters=100);


Algorithm: NewtonRaphson(
    descent = NewtonDescent(),
    autodiff = AutoForwardDiff(),
    vjp_autodiff = AutoFiniteDiff(
        fdtype = Val{:forward}(),
        fdjtype = Val{:forward}(),
        fdhtype = Val{:hcentral}(),
        dir = true
    ),
    jvp_autodiff = AutoForwardDiff(),
    concrete_jac = Val{false}()
)

----    	-------------       	-----------
Iter    	f(u) inf-norm       	Step 2-norm
----    	-------------       	-----------
0       	4.47213595e-06      	0.00000000e+00
1       	7.63931298e-08      	2.76393076e-04
2       	1.33226763e-15      	3.81965696e-08
Final   	1.33226763e-15
----------------------

Continuation

Now we are ready to conduct continuation to obtain the amplitude dependent resonance backbone for the system. The invocation of CONTINUATE is exactly as encountered in the earlier examples.

A0 = -1.;
A1 = 3.;
da = 0.05;
cpars = (parm=:arclength, nmax=1000, Dsc=:none, itopt=4, dpbnds=[da/5, 2da]);

sols, its, dss, xis, Dsc = CONTINUATE([U0;W0s[mi];Xis[mi]], fun, [A0, A1], da; cpars...);

uh = zeros(Complex, maximum(h)+1, length(sols), 2);
for i in 1:2
    uh[h.+1, :, i] = hcat([[s.up[zinds[i:2:end]];
                            s.up[rinds[i:2:end]]+1im*s.up[iinds[i:2:end]]]
                          for s in sols]...);
end
Oms = [s.up[end-2] for s in sols];
Zts = [s.up[end-1] for s in sols]./2Oms;
As = [10^s.up[end] for s in sols];


Algorithm: NewtonRaphson(
    descent = NewtonDescent(),
    autodiff = AutoForwardDiff(),
    vjp_autodiff = AutoFiniteDiff(
        fdtype = Val{:forward}(),
        fdjtype = Val{:forward}(),
        fdhtype = Val{:hcentral}(),
        dir = true
    ),
    jvp_autodiff = AutoForwardDiff(),
    concrete_jac = Val{false}()
)

----    	-------------       	-----------
Iter    	f(u) inf-norm       	Step 2-norm
----    	-------------       	-----------
0       	4.47213595e-05      	NaN
1       	7.63931298e-08      	2.76393076e-04
Final   	7.63931298e-08
----------------------
2. -0.95 with step 0.0500 (0.0500) converged in 1 iterations.
3. -0.85 with step 0.1000 (0.1000) converged in 1 iterations.
4. -0.75 with step 0.1000 (0.1000) converged in 1 iterations.
5. -0.65 with step 0.1000 (0.1000) converged in 1 iterations.
6. -0.55 with step 0.1000 (0.1000) converged in 1 iterations.
7. -0.45 with step 0.1000 (0.1000) converged in 1 iterations.
8. -0.35 with step 0.1000 (0.1000) converged in 1 iterations.
9. -0.25 with step 0.1000 (0.1000) converged in 1 iterations.
10. -0.15 with step 0.1000 (0.1000) converged in 1 iterations.
11. -0.05 with step 0.1000 (0.1000) converged in 1 iterations.
12. 0.05 with step 0.1000 (0.1000) converged in 1 iterations.
13. 0.15 with step 0.1000 (0.1000) converged in 1 iterations.
14. 0.25 with step 0.1000 (0.1000) converged in 1 iterations.
15. 0.33 with step 0.1000 (0.1000) converged in 4 iterations.
16. 0.41 with step 0.0709 (0.1000) converged in 3 iterations.
17. 0.49 with step 0.0768 (0.1000) converged in 3 iterations.
18. 0.58 with step 0.0850 (0.1000) converged in 2 iterations.
19. 0.67 with step 0.0908 (0.1000) converged in 2 iterations.
20. 0.76 with step 0.0940 (0.1000) converged in 2 iterations.
21. 0.86 with step 0.0957 (0.1000) converged in 2 iterations.
22. 0.96 with step 0.0969 (0.1000) converged in 2 iterations.
23. 1.06 with step 0.0977 (0.1000) converged in 2 iterations.
24. 1.15 with step 0.0983 (0.1000) converged in 2 iterations.
25. 1.25 with step 0.0988 (0.1000) converged in 2 iterations.
26. 1.35 with step 0.0992 (0.1000) converged in 2 iterations.
27. 1.45 with step 0.0995 (0.1000) converged in 2 iterations.
28. 1.55 with step 0.0997 (0.1000) converged in 2 iterations.
29. 1.65 with step 0.0998 (0.1000) converged in 2 iterations.
30. 1.75 with step 0.0999 (0.1000) converged in 2 iterations.
31. 1.85 with step 0.0999 (0.1000) converged in 2 iterations.
32. 1.95 with step 0.0999 (0.1000) converged in 2 iterations.
33. 2.05 with step 0.1000 (0.1000) converged in 2 iterations.
34. 2.15 with step 0.1000 (0.1000) converged in 2 iterations.
35. 2.25 with step 0.1000 (0.1000) converged in 2 iterations.
36. 2.35 with step 0.1000 (0.1000) converged in 2 iterations.
37. 2.45 with step 0.1000 (0.1000) converged in 2 iterations.
38. 2.55 with step 0.1000 (0.1000) converged in 2 iterations.
39. 2.65 with step 0.1000 (0.1000) converged in 2 iterations.
40. 2.75 with step 0.1000 (0.1000) converged in 1 iterations.
41. 2.85 with step 0.1000 (0.1000) converged in 1 iterations.
42. 2.95 with step 0.1000 (0.1000) converged in 1 iterations.
43. 3.05 with step 0.1000 (0.1000) converged in 1 iterations.

Plotting

Finally we visualize the results in terms of harmonics as well as in terms of natural frequency and (equivalent) damping factor characteristic curves.

his = [1, 3, 5];

set_theme!(theme_latexfonts());
fsz = 20;
fig = Figure(fontsize=fsz, size=(1000, 600));

axs = [];
for i in eachindex(his[his.<=maximum(h)])
    ax = Axis(fig[1, i],
        ylabel=L"$H_%$(his[i])$ Response (m)", xscale=log10, yscale=log10);
    scatterlines!(ax, As, abs.(uh[his[i].+1, :, 1]), label="x1");
    scatterlines!(ax, As, abs.(uh[his[i].+1, :, 2]), label="x2");
    push!(axs, ax)

    ax = Axis(fig[2, i], xlabel=L"Modal Amplitude $a$",
        ylabel=L"$H_%$(his[i])$ Phase (rad)", xscale=log10);
    scatterlines!(ax, As, unwrap(angle.(uh[his[i].+1, :, 1])), label="x1");
    scatterlines!(ax, As, unwrap(angle.(uh[his[i].+1, :, 2])), label="x2");
    push!(axs, ax)
    if i==1
        axislegend(ax)
    end
end
ax = Axis(fig[1, length(his)+1],
    ylabel=L"Natural Frequency $\omega_n$ (rad/s)", xscale=log10);
scatterlines!(ax, As, Oms);
push!(axs, ax)

ax = Axis(fig[2, length(his)+1], xlabel=L"Modal Amplitude $a$",
    ylabel="Damping Factor (%)", xscale=log10);
scatterlines!(ax, As, 100Zts);
push!(axs, ax)

linkxaxes!(axs...)

    fig

This page was generated using Literate.jl.