Example C2: Forced Response of a 2DoF Oscillator with Hysteretic Nonlinearity

Much like Examples B1 and B2, this builds off of Example C1 to conduct the forced response analysis of a 2DoF oscillator with a hysteretic support.

The dynamical system that will be studied here is:

\[\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.\]

The steps followed in this file are the same as in Example C1, except for the fact that we will now specify that the nonlinearity type is :Hyst:

First we specify the system under study in terms of its "MCK" matrices and the nonlinearities it contains.
Next we setup the parameters for Harmonic Balance and setup the harmonic excitation vector.
And that's it! We're ready to compute the forced responses and 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

We provide the MCK matrices just as before.

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

mdl = MDOFGEN(M, C, K);

For the nonlinearity setup the signature of the nonlinear force is different (read the documentation of ADDNL), now conducting incremental evaluation and making explicit reference to the displacement and force at the previous step(s). This is added to the MDOFGEN model the same way as before, using the ADDNL routine.

# 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);

Setup HB

We setup the HB, test nonlinear force evaluation through NLEVAL!, and test the HB residue.

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);

Uw0 = [E\Fl; Wst];

Evaluate the Nonlinear Forces

FNL = zeros(mdl.Ndofs*Nhc);
dFNLdU = zeros(mdl.Ndofs*Nhc, mdl.Ndofs*Nhc);
dFNLdw = zeros(mdl.Ndofs*Nhc);
NLEVAL!(2Uw0, mdl, h, N; FNL=FNL, dFNLdU=dFNLdU, dFNLdw=dFNLdw)

22-element Vector{Float64}:
  0.0
  0.0
  0.0
  1.0691716930252375
  0.0
 -0.18391770869534557
  0.0
  0.0
  0.0
  0.0
  ⋮
  0.0007156087781473984
  0.0
  0.0
  0.0
  0.0
  0.0
 -0.012273094865006988
  0.0
  0.02612419391440382

Test HB Residue

R = zeros(mdl.Ndofs*Nhc);
dRdU = zeros(mdl.Ndofs*Nhc, mdl.Ndofs*Nhc);
dRdw = zeros(mdl.Ndofs*Nhc);
HBRESFUN!(Uw0, mdl, Fl, h, N; R=R, dRdU=dRdU, dRdw=dRdw)

Famp = 1.0;

fun = NonlinearFunction((r,u,p)->HBRESFUN!([u;p], mdl, Famp*Fl, h, N; R=r),
    jac=(J,u,p)->HBRESFUN!([u;p], mdl, Famp*Fl, h, N; dRdU=J),
    paramjac=(Jp,u,p)->HBRESFUN!([u;p], mdl, Famp*Fl, h, N; dRdw=Jp));

prob = NonlinearProblem(fun, Uw0[1:end-1], Uw0[end]);
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      	1.41421356e+00
1       	5.55111512e-17      	3.41414040e-01
Final   	5.55111512e-17
----------------------

Continuation

Also like before, we do the continuation to compute the full forced response curve.

Note that the code in the last two sections (Setup HB and Continuation) are identical to their counterparts in Example C1. This is what juliajim is all about - the analysis is abstracted enough that the code should look nearly the same except for the setup of the nonlinearity.

Om0 = 0.1;
Om1 = 3;
dOm = 0.2;
cpars = (parm=:arclength, nmax=1000, Dsc=:none);

sols, its, dss, xis, Dsc = CONTINUATE(Uw0[1:end-1], 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 = [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       	7.98855092e-01      	NaN
1       	2.22044605e-16      	5.35141645e-01
Final   	2.22044605e-16
----------------------
2. 0.30 with step 0.1984 (0.1992) converged in 2 iterations.
3. 0.48 with step 0.1907 (0.1992) converged in 3 iterations.
4. 0.61 with step 0.1387 (0.1626) converged in 3 iterations.
5. 0.69 with step 0.0938 (0.1328) converged in 2 iterations.
6. 0.76 with step 0.0759 (0.1328) converged in 2 iterations.
7. 0.81 with step 0.0598 (0.1328) converged in 2 iterations.
8. 0.85 with step 0.0471 (0.1328) converged in 2 iterations.
9. 0.89 with step 0.0375 (0.1328) converged in 2 iterations.
10. 0.91 with step 0.0303 (0.1328) converged in 2 iterations.
11. 0.94 with step 0.0249 (0.1328) converged in 2 iterations.
12. 0.96 with step 0.0208 (0.1328) converged in 2 iterations.
13. 0.97 with step 0.0176 (0.1328) converged in 2 iterations.
14. 0.99 with step 0.0151 (0.1328) converged in 2 iterations.
15. 1.00 with step 0.0131 (0.1328) converged in 2 iterations.
16. 1.01 with step 0.0115 (0.1328) converged in 4 iterations.
17. 1.01 with step 0.0035 (0.0939) converged in 2 iterations.
18. 1.01 with step 0.0028 (0.0939) converged in 2 iterations.
19. 1.01 with step 0.0023 (0.0939) converged in 1 iterations.
20. 1.02 with step 0.0026 (0.1328) converged in 2 iterations.
21. 1.02 with step 0.0020 (0.1328) converged in 2 iterations.
22. 1.02 with step 0.0016 (0.1328) converged in 1 iterations.
23. 1.02 with step 0.0018 (0.1878) converged in 2 iterations.
24. 1.02 with step 0.0013 (0.1878) converged in 1 iterations.
25. 1.02 with step 0.0015 (0.2656) converged in 2 iterations.
26. 1.02 with step 0.0011 (0.2656) converged in 2 iterations.
27. 1.02 with step 0.0008 (0.2656) converged in 1 iterations.
28. 1.03 with step 0.0010 (0.3756) converged in 2 iterations.
29. 1.03 with step 0.0008 (0.3756) converged in 2 iterations.
30. 1.03 with step 0.0007 (0.3756) converged in 2 iterations.
31. 1.03 with step 0.0007 (0.3756) converged in 2 iterations.
32. 1.03 with step 0.0007 (0.3756) converged in 2 iterations.
33. 1.03 with step 0.0007 (0.3756) converged in 2 iterations.
34. 1.03 with step 0.0008 (0.3756) converged in 2 iterations.
35. 1.03 with step 0.0009 (0.3756) converged in 2 iterations.
36. 1.03 with step 0.0010 (0.3756) converged in 2 iterations.
37. 1.03 with step 0.0011 (0.3756) converged in 2 iterations.
38. 1.03 with step 0.0013 (0.3756) converged in 2 iterations.
39. 1.04 with step 0.0015 (0.3756) converged in 2 iterations.
40. 1.04 with step 0.0017 (0.3756) converged in 2 iterations.
41. 1.04 with step 0.0019 (0.3756) converged in 2 iterations.
42. 1.04 with step 0.0024 (0.3756) converged in 2 iterations.
43. 1.05 with step 0.0025 (0.3756) converged in 2 iterations.
44. 1.05 with step 0.0031 (0.3756) converged in 2 iterations.
45. 1.05 with step 0.0036 (0.3756) converged in 2 iterations.
46. 1.06 with step 0.0040 (0.3756) converged in 2 iterations.
47. 1.06 with step 0.0044 (0.3756) converged in 2 iterations.
48. 1.07 with step 0.0052 (0.3756) converged in 2 iterations.
49. 1.07 with step 0.0057 (0.3756) converged in 2 iterations.
50. 1.08 with step 0.0062 (0.3756) converged in 2 iterations.
51. 1.09 with step 0.0070 (0.3756) converged in 2 iterations.
52. 1.10 with step 0.0079 (0.3756) converged in 2 iterations.
53. 1.11 with step 0.0087 (0.3756) converged in 2 iterations.
54. 1.12 with step 0.0098 (0.3756) converged in 2 iterations.
55. 1.13 with step 0.0111 (0.3756) converged in 2 iterations.
56. 1.14 with step 0.0126 (0.3756) converged in 2 iterations.
57. 1.16 with step 0.0143 (0.3756) converged in 2 iterations.
58. 1.17 with step 0.0164 (0.3756) converged in 2 iterations.
59. 1.19 with step 0.0188 (0.3756) converged in 2 iterations.
60. 1.22 with step 0.0218 (0.3756) converged in 2 iterations.
61. 1.24 with step 0.0254 (0.3756) converged in 2 iterations.
62. 1.28 with step 0.0297 (0.3756) converged in 2 iterations.
63. 1.31 with step 0.0344 (0.3756) converged in 2 iterations.
64. 1.36 with step 0.0408 (0.3756) converged in 2 iterations.
65. 1.41 with step 0.0478 (0.3756) converged in 3 iterations.
66. 1.47 with step 0.0442 (0.3067) converged in 3 iterations.
67. 1.54 with step 0.0573 (0.2504) converged in 3 iterations.
68. 1.61 with step 0.0664 (0.2044) converged in 3 iterations.
69. 1.67 with step 0.0627 (0.1669) converged in 3 iterations.
70. 1.71 with step 0.0441 (0.1363) converged in 2 iterations.
71. 1.74 with step 0.0339 (0.1363) converged in 3 iterations.
72. 1.74 with step 0.0095 (0.1113) converged in 2 iterations.
73. 1.75 with step 0.0059 (0.1113) converged in 2 iterations.
74. 1.75 with step 0.0041 (0.1113) converged in 2 iterations.
75. 1.76 with step 0.0030 (0.1113) converged in 2 iterations.
76. 1.76 with step 0.0022 (0.1113) converged in 1 iterations.
77. 1.76 with step 0.0024 (0.1574) converged in 2 iterations.
78. 1.76 with step 0.0018 (0.1574) converged in 2 iterations.
79. 1.76 with step 0.0014 (0.1574) converged in 2 iterations.
80. 1.76 with step 0.0012 (0.1574) converged in 2 iterations.
81. 1.76 with step 0.0010 (0.1574) converged in 1 iterations.
82. 1.77 with step 0.0013 (0.2226) converged in 2 iterations.
83. 1.77 with step 0.0012 (0.2226) converged in 2 iterations.
84. 1.77 with step 0.0012 (0.2226) converged in 2 iterations.
85. 1.77 with step 0.0013 (0.2226) converged in 2 iterations.
86. 1.77 with step 0.0015 (0.2226) converged in 2 iterations.
87. 1.77 with step 0.0016 (0.2226) converged in 2 iterations.
88. 1.77 with step 0.0020 (0.2226) converged in 2 iterations.
89. 1.78 with step 0.0024 (0.2226) converged in 2 iterations.
90. 1.78 with step 0.0029 (0.2226) converged in 2 iterations.
91. 1.78 with step 0.0037 (0.2226) converged in 2 iterations.
92. 1.79 with step 0.0043 (0.2226) converged in 2 iterations.
93. 1.79 with step 0.0051 (0.2226) converged in 2 iterations.
94. 1.80 with step 0.0057 (0.2226) converged in 2 iterations.
95. 1.81 with step 0.0067 (0.2226) converged in 2 iterations.
96. 1.82 with step 0.0076 (0.2226) converged in 2 iterations.
97. 1.82 with step 0.0086 (0.2226) converged in 2 iterations.
98. 1.83 with step 0.0094 (0.2226) converged in 2 iterations.
99. 1.85 with step 0.0104 (0.2226) converged in 2 iterations.
100. 1.86 with step 0.0117 (0.2226) converged in 2 iterations.
101. 1.87 with step 0.0133 (0.2226) converged in 2 iterations.
102. 1.89 with step 0.0143 (0.2226) converged in 2 iterations.
103. 1.90 with step 0.0161 (0.2226) converged in 2 iterations.
104. 1.92 with step 0.0172 (0.2226) converged in 2 iterations.
105. 1.94 with step 0.0179 (0.2226) converged in 2 iterations.
106. 1.96 with step 0.0191 (0.2226) converged in 2 iterations.
107. 1.98 with step 0.0189 (0.2226) converged in 2 iterations.
108. 1.99 with step 0.0177 (0.2226) converged in 2 iterations.
109. 2.01 with step 0.0135 (0.2226) converged in 4 iterations.
110. 2.03 with step 0.0133 (0.1574) converged in 2 iterations.
111. 2.05 with step 0.0172 (0.1574) converged in 2 iterations.
112. 2.07 with step 0.0230 (0.1574) converged in 2 iterations.
113. 2.11 with step 0.0323 (0.1574) converged in 2 iterations.
114. 2.17 with step 0.0475 (0.1574) converged in 2 iterations.
115. 2.26 with step 0.0722 (0.1574) converged in 3 iterations.
116. 2.36 with step 0.0860 (0.1285) converged in 3 iterations.
117. 2.45 with step 0.0858 (0.1049) converged in 2 iterations.
118. 2.54 with step 0.0938 (0.1049) converged in 2 iterations.
119. 2.64 with step 0.0984 (0.1049) converged in 2 iterations.
120. 2.74 with step 0.1009 (0.1049) converged in 2 iterations.
121. 2.85 with step 0.1023 (0.1049) converged in 2 iterations.
122. 2.95 with step 0.1032 (0.1049) converged in 2 iterations.
123. 3.05 with step 0.1037 (0.1049) converged in 2 iterations.

Plotting

Just like before, we visualize the results. Notice, just like in Example B2 that the non-smooth behavior of the nonlinearity is very striking.

his = [1, 3, 5];

fsz = 24;
fig = Figure(fontsize=fsz, size=(1000, 600));

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

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

    fig

This page was generated using Literate.jl.