GitHub

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:

  1. First we specify the system under study in terms of its "MCK" matrices and the nonlinearities it contains.
  2. Next we setup the parameters for Harmonic Balance and setup the harmonic excitation vector.
  3. 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.02625449e+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=:auto, angopt=deg2rad(1));

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      	1.41421356e+00
1       	2.22044605e-16      	5.35141645e-01
Final   	2.22044605e-16
----------------------
2. 0.29 with step 0.1912 (1.9553) converged in 2 iterations. 0.181317 Δα.
3. 0.53 with step 0.2474 (1.4665) converged in 3 iterations. 0.176125 Δα.
4. 0.75 with step 0.2560 (0.9653) converged in 3 iterations. 0.248763 Δα.
5. 0.86 with step 0.1495 (0.6354) converged in 3 iterations. 0.255357 Δα.
6. 0.92 with step 0.0721 (0.4183) converged in 2 iterations. 0.249693 Δα.
7. 0.97 with step 0.0550 (0.3911) converged in 2 iterations. 0.247214 Δα.
8. 1.00 with step 0.0442 (0.3911) converged in 3 iterations. 0.218686 Δα.
9. 1.00 with step 0.0041 (0.3911) converged in 2 iterations. 0.170406 Δα.
10. 1.01 with step 0.0031 (0.3911) converged in 2 iterations. 0.108919 Δα.
11. 1.01 with step 0.0025 (0.3911) converged in 1 iterations. 0.120217 Δα.
12. 1.01 with step 0.0022 (0.3911) converged in 1 iterations. 0.088590 Δα.
13. 1.01 with step 0.0021 (0.3911) converged in 1 iterations. 0.131196 Δα.
14. 1.02 with step 0.0021 (0.3911) converged in 1 iterations. 0.131026 Δα.
15. 1.02 with step 0.0020 (0.3911) converged in 2 iterations. 0.140195 Δα.
16. 1.02 with step 0.0019 (0.3911) converged in 2 iterations. 0.153944 Δα.
17. 1.02 with step 0.0018 (0.3911) converged in 2 iterations. 0.189248 Δα.
18. 1.02 with step 0.0016 (0.3911) converged in 2 iterations. 0.202397 Δα.
19. 1.02 with step 0.0014 (0.3911) converged in 2 iterations. 0.170198 Δα.
20. 1.02 with step 0.0013 (0.3911) converged in 2 iterations. 0.186758 Δα.
21. 1.02 with step 0.0012 (0.3911) converged in 2 iterations. 0.192790 Δα.
22. 1.03 with step 0.0010 (0.3911) converged in 2 iterations. 0.113272 Δα.
23. 1.03 with step 0.0011 (0.3911) converged in 2 iterations. 0.099523 Δα.
24. 1.03 with step 0.0011 (0.3911) converged in 2 iterations. 0.011468 Δα.
25. 1.03 with step 0.0016 (0.4368) converged in 2 iterations. 0.105579 Δα.
26. 1.03 with step 0.0020 (0.3911) converged in 2 iterations. 0.254097 Δα.
27. 1.03 with step 0.0025 (0.3911) converged in 3 iterations. 0.308583 Δα.
28. 1.04 with step 0.0026 (0.3911) converged in 2 iterations. 0.009560 Δα.
29. 1.04 with step 0.0022 (0.4597) converged in 2 iterations. 0.049746 Δα.
30. 1.04 with step 0.0010 (0.3911) converged in 2 iterations. 0.052659 Δα.
31. 1.04 with step 0.0008 (0.3911) converged in 1 iterations. 0.043061 Δα.
32. 1.04 with step 0.0009 (0.4010) converged in 1 iterations. 0.026527 Δα.
33. 1.04 with step 0.0010 (0.4462) converged in 1 iterations. 0.006343 Δα.
34. 1.05 with step 0.0017 (0.6856) converged in 2 iterations. 0.026882 Δα.
35. 1.05 with step 0.0021 (0.6190) converged in 2 iterations. 0.017039 Δα.
36. 1.05 with step 0.0027 (0.6227) converged in 2 iterations. 0.045791 Δα.
37. 1.05 with step 0.0020 (0.5036) converged in 2 iterations. 0.042988 Δα.
38. 1.05 with step 0.0014 (0.4122) converged in 1 iterations. 0.008760 Δα.
39. 1.06 with step 0.0021 (0.5825) converged in 1 iterations. 0.007407 Δα.
40. 1.06 with step 0.0031 (0.8589) converged in 2 iterations. 0.064532 Δα.
41. 1.06 with step 0.0035 (0.6528) converged in 2 iterations. 0.075737 Δα.
42. 1.07 with step 0.0042 (0.4896) converged in 2 iterations. 0.085864 Δα.
43. 1.07 with step 0.0051 (0.3911) converged in 2 iterations. 0.200574 Δα.
44. 1.08 with step 0.0079 (0.3911) converged in 2 iterations. 0.023944 Δα.
45. 1.10 with step 0.0133 (0.3911) converged in 3 iterations. 0.079149 Δα.
46. 1.13 with step 0.0241 (0.3911) converged in 3 iterations. 0.166580 Δα.
47. 1.17 with step 0.0333 (0.3911) converged in 3 iterations. 0.209305 Δα.
48. 1.22 with step 0.0363 (0.3911) converged in 3 iterations. 0.215816 Δα.
49. 1.26 with step 0.0364 (0.3911) converged in 3 iterations. 0.202333 Δα.
50. 1.30 with step 0.0341 (0.3911) converged in 2 iterations. 0.194790 Δα.
51. 1.33 with step 0.0295 (0.3911) converged in 2 iterations. 0.091683 Δα.
52. 1.36 with step 0.0247 (0.3911) converged in 2 iterations. 0.085698 Δα.
53. 1.38 with step 0.0219 (0.3911) converged in 2 iterations. 0.124087 Δα.
54. 1.41 with step 0.0197 (0.3911) converged in 2 iterations. 0.126618 Δα.
55. 1.46 with step 0.0124 (0.3911) converged in 4 iterations. 0.304006 Δα.
56. 1.56 with step 0.0777 (0.3911) converged in 3 iterations. 0.658107 Δα.
57. 1.66 with step 0.0932 (0.3911) converged in 3 iterations. 2.457037 Δα.
58. 1.71 with step 0.0547 (0.3911) converged in 2 iterations. 0.409148 Δα.
59. 1.73 with step 0.0232 (0.3911) converged in 2 iterations. 0.162320 Δα.
60. 1.74 with step 0.0147 (0.3911) converged in 4 iterations. 2.286520 Δα.
61. 1.74 with step 0.0046 (0.3911) converged in 2 iterations. 0.169740 Δα.
62. 1.74 with step 0.0036 (0.3911) converged in 2 iterations. 0.155536 Δα.
63. 1.75 with step 0.0029 (0.3911) converged in 1 iterations. 0.126541 Δα.
64. 1.75 with step 0.0025 (0.3911) converged in 1 iterations. 0.132522 Δα.
65. 1.75 with step 0.0024 (0.3911) converged in 2 iterations. 0.142450 Δα.
66. 1.75 with step 0.0022 (0.3911) converged in 1 iterations. 0.138671 Δα.
67. 1.76 with step 0.0021 (0.3911) converged in 2 iterations. 0.165725 Δα.
68. 1.76 with step 0.0020 (0.3911) converged in 2 iterations. 0.174734 Δα.
69. 1.76 with step 0.0020 (0.3911) converged in 2 iterations. 0.205094 Δα.
70. 1.76 with step 0.0018 (0.3911) converged in 2 iterations. 0.168021 Δα.
71. 1.76 with step 0.0017 (0.3911) converged in 2 iterations. 0.136844 Δα.
72. 1.76 with step 0.0016 (0.3911) converged in 2 iterations. 0.118312 Δα.
73. 1.77 with step 0.0017 (0.3911) converged in 2 iterations. 0.046899 Δα.
74. 1.77 with step 0.0019 (0.3911) converged in 2 iterations. 0.024784 Δα.
75. 1.77 with step 0.0021 (0.3911) converged in 2 iterations. 0.096782 Δα.
76. 1.77 with step 0.0018 (0.3911) converged in 2 iterations. 0.165317 Δα.
77. 1.77 with step 0.0016 (0.3911) converged in 2 iterations. 0.128657 Δα.
78. 1.77 with step 0.0016 (0.3911) converged in 2 iterations. 0.139956 Δα.
79. 1.78 with step 0.0016 (0.3911) converged in 2 iterations. 0.237978 Δα.
80. 1.78 with step 0.0017 (0.3911) converged in 2 iterations. 0.037295 Δα.
81. 1.78 with step 0.0018 (0.3911) converged in 1 iterations. 0.013030 Δα.
82. 1.78 with step 0.0027 (0.5028) converged in 2 iterations. 0.063216 Δα.
83. 1.79 with step 0.0027 (0.3911) converged in 2 iterations. 0.087029 Δα.
84. 1.79 with step 0.0033 (0.3911) converged in 2 iterations. 0.104102 Δα.
85. 1.79 with step 0.0037 (0.3911) converged in 2 iterations. 0.046973 Δα.
86. 1.80 with step 0.0028 (0.3911) converged in 2 iterations. 0.100378 Δα.
87. 1.80 with step 0.0026 (0.3911) converged in 2 iterations. 0.123493 Δα.
88. 1.80 with step 0.0032 (0.3911) converged in 1 iterations. 0.024386 Δα.
89. 1.81 with step 0.0038 (0.4419) converged in 2 iterations. 0.067037 Δα.
90. 1.81 with step 0.0039 (0.3911) converged in 2 iterations. 0.023156 Δα.
91. 1.82 with step 0.0050 (0.3911) converged in 1 iterations. 0.017939 Δα.
92. 1.82 with step 0.0080 (0.4694) converged in 2 iterations. 0.087318 Δα.
93. 1.84 with step 0.0100 (0.3911) converged in 2 iterations. 0.035556 Δα.
94. 1.85 with step 0.0165 (0.3911) converged in 2 iterations. 0.074423 Δα.
95. 1.88 with step 0.0249 (0.3911) converged in 3 iterations. 0.116355 Δα.
96. 1.92 with step 0.0328 (0.3911) converged in 3 iterations. 0.117185 Δα.
97. 1.95 with step 0.0319 (0.3911) converged in 2 iterations. 0.081583 Δα.
98. 1.98 with step 0.0239 (0.3911) converged in 2 iterations. 0.066118 Δα.
99. 1.99 with step 0.0143 (0.3911) converged in 2 iterations. 0.061020 Δα.
100. 2.00 with step 0.0076 (0.3911) converged in 2 iterations. 0.134896 Δα.
101. 2.01 with step 0.0036 (0.3911) converged in 3 iterations. 0.792522 Δα.
102. 2.05 with step 0.0259 (0.3911) converged in 2 iterations. 0.311127 Δα.
103. 2.11 with step 0.0418 (0.3911) converged in 3 iterations. 0.382495 Δα.
104. 2.21 with step 0.0680 (0.3911) converged in 3 iterations. 0.441913 Δα.
105. 2.38 with step 0.1078 (0.3911) converged in 3 iterations. 0.482185 Δα.
106. 2.63 with step 0.1643 (0.3911) converged in 3 iterations. 0.505426 Δα.
107. 2.99 with step 0.2381 (0.3911) converged in 3 iterations. 0.512471 Δα.
108. 3.54 with step 0.3786 (0.3911) converged in 3 iterations. 0.566013 Δα.

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
Example block output

This page was generated using Literate.jl.