Example C3: Response Constrained Forced Response of a 2DoF Hysteretic Oscillator

It is sometimes the case that it is necessary to conduct forced response simulations where the response amplitude is kept fixed. In fact, experimentally it is significantly easier to achieve response-controlled forced response measurements (in the context of, say, stepped sine testing) than force-controlled measurements. While there are several reasons for this (large displacement effects on the exciter at resonance, low SNR at anti-resonance points, etc.), we will not go into this here.

MDOFUTILS provides HBRESFUN_A!, a routine that returns the harmonic balance residue along with an amplitude constraint. The usage of this is only slightly different from HBRESFUN! (which we have already encountered in Examples C1 and C2 for forced response analysis).

The same system as before will be studied here:

\[\underbrace{\begin{bmatrix} 1&0\\0&1 \end{bmatrix}}_{\mx{M}} \begin{bmatrix} \ddot{x_1}\\ \ddot{x_2} \end{bmatrix} + \left( 0.01 \mx{M} + 0.001 \mx{K}\right) \begin{bmatrix} \dot{x_1}\\ \dot{x_2} \end{bmatrix} + \underbrace{\begin{bmatrix} 2&-1\\-1&2 \end{bmatrix}}_{\mx{K}} \begin{bmatrix} x_1\\x_2 \end{bmatrix} + \begin{bmatrix} 0\\ f_{fr}(x_2) \end{bmatrix} = \begin{bmatrix} 1\\0 \end{bmatrix}\cos\Omega t.\]

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 = [1. 0.;0. 1.];
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. 1.];

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

Setup HB

HB setup: identical to Example C2.

h = 1:2:5;
# 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;

Wst = 0.6;
E, dEdw = HARMONICSTIFFNESS(mdl.M, mdl.C, mdl.K, Wst, h);
U0 = E\Fl;

Test HB Residue

The residue function HBRESFUN_A! expects a vector of the form (\begin{bmatrix} \vc{u}& F& \Omega \end{bmatrix}^T), where (\vc{u}) is the vector of harmonics, (F) is the (scalar) excitation amplitude (also an unknown), and (\Omega) is the (scalar) excitation frequency. Apart from the same arguments as HBRESFUN!, the required response amplitude Amp must also be provided here. By default, this amplitude is interpreted as the first harmonic amplitude of the first DoF. This can be customized by supplying optional keywords atype (one of :H1 or :RMS, specifying type of amplitude measure) and ashape (either DoF index or vector of weights to DoFs); see the documentation. Here we set this up using NonlinearSolve.jl and solve for a single point.

R = zeros(mdl.Ndofs*Nhc+1);
dRdUf = zeros(mdl.Ndofs*Nhc+1, mdl.Ndofs*Nhc+1);
dRdw = zeros(mdl.Ndofs*Nhc+1);
Amp = 1e0;
dof = 1;

HBRESFUN_A!([U0; 1.0; Wst], mdl, Amp, Fl, h, N; R=R, dRdUf=dRdUf, dRdw=dRdw)

fun = NonlinearFunction((r,uf,p)->HBRESFUN_A!([uf;p], mdl, Amp, Fl, h, N; R=r),
    jac=(J,uf,p)->HBRESFUN_A!([uf;p], mdl, Amp, Fl, h, N; dRdUf=J),
    paramjac=(Jp,uf,p)->HBRESFUN_A!([uf;p], mdl, Amp, Fl, h, N; dRdw=Jp));

prob = NonlinearProblem(fun, [U0;1.0], Wst);
sol = solve(prob, show_trace=Val(true));

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       	5.91774639e-01      	8.57728610e-01
1       	8.98242782e-04      	3.38756045e-01
2       	2.01529002e-07      	7.42307628e-04
3       	1.02140518e-14      	1.66618285e-07
Final   	1.02140518e-14
----------------------

Continuation

Now we continue as before with continuation. In fact this piece of code is exactly the same as in Example C2.

Om0 = 0.1;
Om1 = 3;
dOm = 0.1;

cpars = (parm=:arclength, nmax=400, Dsc=:none, itopt=4);
sols, its, dss, xis, Dsc = CONTINUATE([U0;1.0], fun, [Om0, Om1], dOm; 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 = sols.p;
Fs = [s.up[end-1] 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       	7.98855092e-01      	0.00000000e+00
1       	9.56746203e-04      	7.05641605e-01
2       	2.28622091e-07      	9.38517141e-04
Final   	2.28622091e-07
----------------------
2. 0.19 with step 0.0952 (0.0976) converged in 2 iterations.
3. 0.31 with step 0.1267 (0.1380) converged in 2 iterations.
4. 0.46 with step 0.1596 (0.1952) converged in 2 iterations.
5. 0.63 with step 0.1909 (0.2760) converged in 2 iterations.
6. 0.83 with step 0.2203 (0.3904) converged in 2 iterations.
7. 1.02 with step 0.2191 (0.4879) converged in 3 iterations.
8. 1.16 with step 0.1724 (0.4879) converged in 771 iterations.
9. 1.16 with step 0.0159 (0.2440) converged in 3 iterations.
10. 1.17 with step 0.0063 (0.2817) converged in 3 iterations.
11. 1.17 with step 0.0035 (0.3253) converged in 3 iterations.
12. 1.17 with step 0.0025 (0.3756) converged in 3 iterations.
13. 1.18 with step 0.0022 (0.4337) converged in 3 iterations.
14. 1.18 with step 0.0023 (0.4879) converged in 3 iterations.
15. 1.18 with step 0.0031 (0.4879) converged in 3 iterations.
16. 1.20 with step 0.0068 (0.4879) converged in 3 iterations.
17. 1.18 with step 0.0505 (0.4879) converged in 197 iterations.
18. 1.19 with step 0.0034 (0.2440) converged in 3 iterations.
19. 1.21 with step 0.0083 (0.2817) converged in 3 iterations.
20. 1.28 with step 0.0422 (0.3253) converged in 5 iterations.
21. 1.34 with step 0.0673 (0.2910) converged in 2 iterations.
22. 1.41 with step 0.0820 (0.4115) converged in 3 iterations.
23. 1.45 with step 0.0630 (0.4751) converged in 3 iterations.
24. 1.46 with step 0.0155 (0.4879) converged in 3 iterations.
25. 1.46 with step 0.0061 (0.4879) converged in 3 iterations.
26. 1.46 with step 0.0030 (0.4879) converged in 2 iterations.
27. 1.46 with step 0.0018 (0.4879) converged in 2 iterations.
28. 1.47 with step 0.0019 (0.4879) converged in 2 iterations.
29. 1.47 with step 0.0025 (0.4879) converged in 2 iterations.
30. 1.48 with step 0.0044 (0.4879) converged in 3 iterations.
31. 1.47 with step 0.0115 (0.4879) converged in 193 iterations.
32. 1.47 with step 0.0022 (0.2440) converged in 2 iterations.
33. 1.48 with step 0.0047 (0.3450) converged in 3 iterations.
34. 1.47 with step 0.0128 (0.3984) converged in 26 iterations.
35. 1.47 with step 0.0024 (0.1992) converged in 2 iterations.
36. 1.48 with step 0.0049 (0.2817) converged in 3 iterations.
37. 1.51 with step 0.0169 (0.3253) converged in 89 iterations.
38. 1.50 with step -0.0091 (0.1626) converged in 2 iterations.
39. 1.48 with step -0.0155 (0.2300) converged in 17 iterations.
40. 1.48 with step 0.0040 (0.1150) converged in 2 iterations.
41. 1.50 with step 0.0118 (0.1626) converged in 6 iterations.
42. 1.51 with step 0.0091 (0.1328) converged in 2 iterations.
43. 1.52 with step 0.0109 (0.1878) converged in 2 iterations.
44. 1.53 with step 0.0141 (0.2656) converged in 2 iterations.
45. 1.55 with step 0.0197 (0.3756) converged in 2 iterations.
46. 1.58 with step 0.0271 (0.4879) converged in 2 iterations.
47. 1.61 with step 0.0308 (0.4879) converged in 2 iterations.
48. 1.65 with step 0.0349 (0.4879) converged in 2 iterations.
49. 1.69 with step 0.0397 (0.4879) converged in 2 iterations.
50. 1.74 with step 0.0438 (0.4879) converged in 2 iterations.
51. 1.78 with step 0.0443 (0.4879) converged in 3 iterations.
52. 1.81 with step 0.0386 (0.4879) converged in 3 iterations.
53. 1.84 with step 0.0311 (0.4879) converged in 3 iterations.
54. 1.87 with step 0.0256 (0.4879) converged in 3 iterations.
55. 1.89 with step 0.0225 (0.4879) converged in 3 iterations.
56. 1.91 with step 0.0223 (0.4879) converged in 3 iterations.
57. 1.94 with step 0.0255 (0.4879) converged in 3 iterations.
58. 1.98 with step 0.0329 (0.4879) converged in 3 iterations.
59. 2.04 with step 0.0420 (0.4879) converged in 5 iterations.
60. 2.10 with step 0.0552 (0.4364) converged in 2 iterations.
61. 2.17 with step 0.0719 (0.4879) converged in 2 iterations.
62. 2.25 with step 0.0802 (0.4879) converged in 2 iterations.
63. 2.34 with step 0.0851 (0.4879) converged in 2 iterations.
64. 2.43 with step 0.0874 (0.4879) converged in 2 iterations.
65. 2.52 with step 0.0878 (0.4879) converged in 2 iterations.
66. 2.60 with step 0.0872 (0.4879) converged in 2 iterations.
67. 2.69 with step 0.0860 (0.4879) converged in 2 iterations.
68. 2.77 with step 0.0845 (0.4879) converged in 2 iterations.
69. 2.85 with step 0.0829 (0.4879) converged in 2 iterations.
70. 2.93 with step 0.0812 (0.4879) converged in 2 iterations.
71. 3.01 with step 0.0795 (0.4879) converged in 2 iterations.

Plotting

We plot first just the responses, to demonstrate how the amplitude constraint has been effective.

Amplitude Responses

his = [1, 3, 5];

set_theme!(theme_latexfonts())
fsz = 24;
fig1 = Figure(fontsize=fsz, size=(1000, 600));

for i in eachindex(his[his.<=maximum(h)])
    ax = Axis(fig1[1, i],
        ylabel=L"$H_%$(his[i])$ Response (m)", yscale=log10);
    scatterlines!(ax, Oms, abs.(uh[his[i].+1, :, 1]), label="x1");
    scatterlines!(ax, Oms, abs.(uh[his[i].+1, :, 2]), label="x2");
    if i==1
        axislegend(ax)
    end

    ax = Axis(fig1[2, i], xlabel=L"Excitation Frequency $\Omega$",
        ylabel=L"$H_%$(his[i])$ Phase (rad)");
    scatterlines!(ax, Oms, unwrap(angle.(uh[his[i].+1, :, 1])), label="x1");
    scatterlines!(ax, Oms, unwrap(angle.(uh[his[i].+1, :, 2])), label="x2");
end

    fig1
Example block output

It can be observed that the first harmonic of the first DoF is always fixed to the value of Amp ($1.0$ here), while the other DoF's amplitude changes, showing that we have a fixed amplitude forced response here.

Forced Responses

In order to visualize the nonlinear forced response function, we plot the response over force next. Here the resonances are clearly visible, albeit very different from what was seen in Example C2.

his = [1, 3, 5];

set_theme!(theme_latexfonts())
fsz = 24;
fig2 = Figure(fontsize=fsz, size=(1000, 600));

for i in eachindex(his[his.<=maximum(h)])
    ax = Axis(fig2[1, i],
        ylabel=L"$H_%$(his[i])$ Response (m/N)", yscale=log10);
    scatterlines!(ax, Oms, abs.(uh[his[i].+1, :, 1]./Fs), label="x1");
    scatterlines!(ax, Oms, abs.(uh[his[i].+1, :, 2]./Fs), label="x2");

    ax = Axis(fig2[2, i], xlabel=L"Excitation Frequency $\Omega$",
        ylabel=L"$H_%$(his[i])$ Phase (rad)");
    scatterlines!(ax, Oms, unwrap(angle.(uh[his[i].+1, :, 1])), label="x1");
    scatterlines!(ax, Oms, unwrap(angle.(uh[his[i].+1, :, 2])), label="x2");
end

    fig2
Example block output

This page was generated using Literate.jl.