Example F: Finite Element Beam with Nonlinear Support

This example considers a finite element model of a Beam (Euler-Bernoulli) that has a nonlinear attachment in one end and is clamped in the other. The same code can be used for simulating 3 different types of nonlinear attachments:

  1. A simple cubic spring of the form $\beta u^3$.
  2. A stiff string attachment with reaction force of the form $T \frac{w_T}{\sqrt{\ell_s^2+w_T^2}}$
  3. A non-smooth frictional joint.

Preamble: Load Packages

using GLMakie
using LinearAlgebra
using Random
using NonlinearSolve
using SparseArrays
using Revise

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

Setup System

Ey = 70e9;
rho = 2700.;
w, h = 3*25e-3, 12.5e-3;
ell = 1.0;
ρA = rho*w*h;
EI = Ey*w*h^3/12;

Cubic Spring

β = 500.0;  # Nonlinearity
fnl = (t, u, ud) -> return β.*u.^3, 3β.*u.^2, zeros(size(u));

Stiffened String

Ts = 1e2;
ls = 0.10;
fnl = (t, u, ud) -> return Ts.*u./sqrt.(ls^2 .+u.^2),
    Ts*ls^2 ./sqrt.(ls^2 .+u.^2).^3, zeros(size(u));
typ = :Inst;

Frictional Support

kt = 500;
fs = 1e-2kt;
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
typ = :Hyst;

Ne = 10;  # Number of elements
Nn = Ne+1;  # Number of nodes
Xn = range(0, ell, Nn);

Me(Le) = ρA*Le/420*[156 22Le 54 -13Le;
	            22Le 4Le^2 13Le -3Le^2;
	            54 13Le 156 -22Le;
	            -13Le -3Le^2 -22Le 4Le^2];
Ke(Le) = EI/Le^3*[12 6Le -12 6Le;
		  6Le 4Le^2 -6Le 2Le^2;
		  -12 -6Le 12 -6Le;
		  6Le 2Le^2 -6Le 4Le^2];

M = zeros(2Nn, 2Nn);
K = zeros(2Nn, 2Nn);
for ei in 1:Ne
    is = 2(ei-1)+1;
    ie = 2(ei+1);

    M[is:ie, is:ie] += Me(diff(Xn[ei:ei+1])[1]);
    K[is:ie, is:ie] += Ke(diff(Xn[ei:ei+1])[1]);
end
F = [zeros(2Nn-2); 1.;0.];

Apply Clamped Boundary Condition

Lb = I(2Nn)[:, 3:end];
Mb = Lb'M*Lb;
Kb = Lb'K*Lb;
Fb = Lb'F;

Wn = sqrt.(eigvals(Kb, Mb));
Zts = [0.2e-2, 0.1e-2];
ab = [1 ./2Wn[1:2] Wn[1:2]/2]\Zts;
Cb = ab[1]*Mb + ab[2]*Kb;

Setup Nonlinearity

mdl = MDOFGEN(Mb, Cb, Kb);
mdl = ADDNL(mdl, typ, fnl, Float64.(Lb[end-1:end-1,:]));

Setup HB

h = 0:5;
N = 256;

Nhc = NHC(h);
_, _, zinds, rinds, iinds = HINDS(mdl.Ndofs, h);

Fl = zeros(Nhc*mdl.Ndofs);
Fl[rinds[1:mdl.Ndofs]] = Fb;

Om0 = 10.0;
Om1 = 450.0;
dOm = 10.0;

Linear Initial Guess

E, _ = HARMONICSTIFFNESS(mdl.M, mdl.C, mdl.K, Om0, h);
U0_(Fa) = E\ (Fa*Fl);

Setup nonlinear function

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

Forced Response Continuation

Famp = 8.0  # 0.1, 1.0, 2.0, 4.0, 8.0, 16.0

U0 = U0_(Famp);
prob = NonlinearProblem(fun(Famp), U0, Om0; abstol=1e-6, reltol=1e-6);
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       	1.59758940e+00      	0.00000000e+00
1       	5.32312255e-12      	2.12720141e-03
Final   	5.32312255e-12
----------------------

Continuation

cpars = (parm=:arclength, nmax=2000, Dsc=:auto);

Om0 = 40.0;
Om1 = 100.0;
dOm = 2.5;
sols, _, _, _, _ = CONTINUATE(U0, fun(Famp), [Om0, Om1], dOm; cpars...);

Uh = [Lb[end-1,:]'reshape(u, mdl.Ndofs,:) for u in sols.u];

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       	1.03236539e+00      	0.00000000e+00
1       	7.79998288e-12      	2.41017353e-03
Final   	7.79998288e-12
----------------------
2. 42.50 with step 2.4997 (0.0625) converged in 1 iterations.
3. 45.16 with step 2.6557 (0.0625) converged in 1 iterations.
4. 47.98 with step 2.8215 (0.0625) converged in 1 iterations.
5. 50.97 with step 2.9973 (0.0625) converged in 1 iterations.
6. 54.16 with step 3.1833 (0.0625) converged in 2 iterations.
7. 56.54 with step 2.3891 (0.0442) converged in 2 iterations.
8. 58.30 with step 1.7606 (0.0312) converged in 2 iterations.
9. 59.58 with step 1.2798 (0.0221) converged in 2 iterations.
10. 60.50 with step 0.9210 (0.0156) converged in 2 iterations.
11. 61.23 with step 0.7444 (0.0125) converged in 3 iterations.
12. 61.91 with step 0.7122 (0.0125) converged in 2 iterations.
13. 62.49 with step 0.6362 (0.0125) converged in 3 iterations.
14. 62.91 with step 0.5034 (0.0125) converged in 3 iterations.
15. 63.21 with step 0.3531 (0.0125) converged in 2 iterations.
16. 63.43 with step 0.2497 (0.0125) converged in 2 iterations.
17. 63.58 with step 0.1805 (0.0125) converged in 2 iterations.
18. 63.71 with step 0.1376 (0.0125) converged in 2 iterations.
19. 63.80 with step 0.1075 (0.0125) converged in 2 iterations.
20. 63.88 with step 0.0863 (0.0125) converged in 2 iterations.
21. 63.95 with step 0.0721 (0.0125) converged in 2 iterations.
22. 64.00 with step 0.0611 (0.0125) converged in 2 iterations.
23. 64.05 with step 0.0525 (0.0125) converged in 2 iterations.
24. 64.10 with step 0.0456 (0.0125) converged in 2 iterations.
25. 64.13 with step 0.0401 (0.0125) converged in 2 iterations.
26. 64.17 with step 0.0355 (0.0125) converged in 2 iterations.
27. 64.20 with step 0.0317 (0.0125) converged in 2 iterations.
28. 64.22 with step 0.0285 (0.0125) converged in 1 iterations.
29. 64.25 with step 0.0264 (0.0125) converged in 2 iterations.
30. 64.27 with step 0.0239 (0.0125) converged in 2 iterations.
31. 64.29 with step 0.0218 (0.0125) converged in 1 iterations.
32. 64.31 with step 0.0206 (0.0125) converged in 2 iterations.
33. 64.33 with step 0.0189 (0.0125) converged in 1 iterations.
34. 64.35 with step 0.0174 (0.0125) converged in 1 iterations.
35. 64.36 with step 0.0166 (0.0125) converged in 1 iterations.
36. 64.38 with step 0.0154 (0.0125) converged in 1 iterations.
37. 64.39 with step 0.0148 (0.0125) converged in 2 iterations.
38. 64.41 with step 0.0139 (0.0125) converged in 1 iterations.
39. 64.42 with step 0.0130 (0.0125) converged in 1 iterations.
40. 64.43 with step 0.0126 (0.0125) converged in 1 iterations.
41. 64.44 with step 0.0119 (0.0125) converged in 1 iterations.
42. 64.46 with step 0.0116 (0.0125) converged in 1 iterations.
43. 64.47 with step 0.0110 (0.0125) converged in 1 iterations.
44. 64.48 with step 0.0108 (0.0125) converged in 2 iterations.
45. 64.49 with step 0.0103 (0.0125) converged in 1 iterations.
46. 64.50 with step 0.0102 (0.0125) converged in 2 iterations.
47. 64.51 with step 0.0097 (0.0125) converged in 1 iterations.
48. 64.52 with step 0.0093 (0.0125) converged in 1 iterations.
49. 64.52 with step 0.0093 (0.0125) converged in 2 iterations.
50. 64.53 with step 0.0089 (0.0125) converged in 1 iterations.
51. 64.54 with step 0.0090 (0.0125) converged in 2 iterations.
52. 64.55 with step 0.0087 (0.0125) converged in 1 iterations.
53. 64.56 with step 0.0084 (0.0125) converged in 2 iterations.
54. 64.57 with step 0.0085 (0.0125) converged in 2 iterations.
55. 64.58 with step 0.0083 (0.0125) converged in 1 iterations.
56. 64.58 with step 0.0085 (0.0125) converged in 2 iterations.
57. 64.59 with step 0.0084 (0.0125) converged in 2 iterations.
58. 64.60 with step 0.0082 (0.0125) converged in 1 iterations.
59. 64.61 with step 0.0085 (0.0125) converged in 2 iterations.
60. 64.62 with step 0.0084 (0.0125) converged in 2 iterations.
61. 64.63 with step 0.0084 (0.0125) converged in 1 iterations.
62. 64.63 with step 0.0083 (0.0125) converged in 2 iterations.
63. 64.64 with step 0.0088 (0.0125) converged in 1 iterations.
64. 64.65 with step 0.0088 (0.0125) converged in 2 iterations.
65. 64.66 with step 0.0088 (0.0125) converged in 2 iterations.
66. 64.67 with step 0.0089 (0.0125) converged in 2 iterations.
67. 64.68 with step 0.0094 (0.0125) converged in 2 iterations.
68. 64.69 with step 0.0096 (0.0125) converged in 2 iterations.
69. 64.70 with step 0.0097 (0.0125) converged in 2 iterations.
70. 64.71 with step 0.0098 (0.0125) converged in 2 iterations.
71. 64.72 with step 0.0100 (0.0125) converged in 1 iterations.
72. 64.73 with step 0.0107 (0.0125) converged in 2 iterations.
73. 64.74 with step 0.0110 (0.0125) converged in 2 iterations.
74. 64.75 with step 0.0112 (0.0125) converged in 2 iterations.
75. 64.76 with step 0.0110 (0.0125) converged in 2 iterations.
76. 64.78 with step 0.0112 (0.0125) converged in 3 iterations.
77. 64.79 with step 0.0121 (0.0125) converged in 2 iterations.
78. 64.80 with step 0.0123 (0.0125) converged in 2 iterations.
79. 64.81 with step 0.0126 (0.0125) converged in 2 iterations.
80. 64.83 with step 0.0128 (0.0125) converged in 2 iterations.
81. 64.84 with step 0.0131 (0.0125) converged in 2 iterations.
82. 64.85 with step 0.0133 (0.0125) converged in 1 iterations.
83. 64.87 with step 0.0136 (0.0125) converged in 1 iterations.
84. 64.88 with step 0.0139 (0.0125) converged in 1 iterations.
85. 64.90 with step 0.0151 (0.0125) converged in 2 iterations.
86. 64.91 with step 0.0154 (0.0125) converged in 2 iterations.
87. 64.93 with step 0.0159 (0.0125) converged in 2 iterations.
88. 64.94 with step 0.0172 (0.0125) converged in 2 iterations.
89. 64.96 with step 0.0165 (0.0125) converged in 1 iterations.
90. 64.98 with step 0.0171 (0.0125) converged in 1 iterations.
91. 65.00 with step 0.0185 (0.0125) converged in 2 iterations.
92. 65.02 with step 0.0192 (0.0125) converged in 2 iterations.
93. 65.04 with step 0.0199 (0.0125) converged in 2 iterations.
94. 65.06 with step 0.0216 (0.0125) converged in 2 iterations.
95. 65.08 with step 0.0212 (0.0125) converged in 2 iterations.
96. 65.10 with step 0.0230 (0.0125) converged in 1 iterations.
97. 65.13 with step 0.0241 (0.0125) converged in 1 iterations.
98. 65.16 with step 0.0252 (0.0125) converged in 2 iterations.
99. 65.18 with step 0.0274 (0.0125) converged in 2 iterations.
100. 65.21 with step 0.0288 (0.0125) converged in 2 iterations.
101. 65.24 with step 0.0313 (0.0125) converged in 2 iterations.
102. 65.28 with step 0.0330 (0.0125) converged in 2 iterations.
103. 65.31 with step 0.0343 (0.0125) converged in 1 iterations.
104. 65.35 with step 0.0364 (0.0125) converged in 2 iterations.
105. 65.39 with step 0.0387 (0.0125) converged in 2 iterations.
106. 65.44 with step 0.0423 (0.0125) converged in 2 iterations.
107. 65.49 with step 0.0452 (0.0125) converged in 2 iterations.
108. 65.54 with step 0.0514 (0.0125) converged in 2 iterations.
109. 65.60 with step 0.0551 (0.0125) converged in 2 iterations.
110. 65.66 with step 0.0604 (0.0125) converged in 2 iterations.
111. 65.73 with step 0.0653 (0.0125) converged in 2 iterations.
112. 65.80 with step 0.0708 (0.0125) converged in 2 iterations.
113. 65.89 with step 0.0806 (0.0125) converged in 2 iterations.
114. 65.98 with step 0.0880 (0.0125) converged in 2 iterations.
115. 66.08 with step 0.0965 (0.0125) converged in 2 iterations.
116. 66.20 with step 0.1104 (0.0125) converged in 2 iterations.
117. 66.33 with step 0.1250 (0.0125) converged in 2 iterations.
118. 66.48 with step 0.1389 (0.0125) converged in 2 iterations.
119. 66.65 with step 0.1585 (0.0125) converged in 2 iterations.
120. 66.85 with step 0.1832 (0.0125) converged in 2 iterations.
121. 67.08 with step 0.2141 (0.0125) converged in 2 iterations.
122. 67.35 with step 0.2464 (0.0125) converged in 2 iterations.
123. 67.66 with step 0.2882 (0.0125) converged in 2 iterations.
124. 68.02 with step 0.3365 (0.0125) converged in 2 iterations.
125. 68.44 with step 0.3892 (0.0125) converged in 2 iterations.
126. 68.93 with step 0.4528 (0.0125) converged in 2 iterations.
127. 69.48 with step 0.5193 (0.0125) converged in 2 iterations.
128. 70.09 with step 0.5839 (0.0125) converged in 2 iterations.
129. 70.77 with step 0.6472 (0.0125) converged in 2 iterations.
130. 71.49 with step 0.7019 (0.0125) converged in 2 iterations.
131. 72.26 with step 0.7490 (0.0125) converged in 2 iterations.
132. 73.06 with step 0.7877 (0.0125) converged in 2 iterations.
133. 73.89 with step 0.8203 (0.0125) converged in 2 iterations.
134. 74.74 with step 0.8478 (0.0125) converged in 2 iterations.
135. 75.62 with step 0.8703 (0.0125) converged in 2 iterations.
136. 76.52 with step 0.8907 (0.0125) converged in 2 iterations.
137. 77.43 with step 0.9083 (0.0125) converged in 2 iterations.
138. 78.35 with step 0.9224 (0.0125) converged in 2 iterations.
139. 79.28 with step 0.9289 (0.0125) converged in 3 iterations.
140. 80.24 with step 0.9509 (0.0125) converged in 2 iterations.
141. 81.22 with step 0.9749 (0.0125) converged in 2 iterations.
142. 82.22 with step 0.9951 (0.0125) converged in 2 iterations.
143. 83.23 with step 1.0129 (0.0125) converged in 2 iterations.
144. 84.26 with step 1.0293 (0.0125) converged in 2 iterations.
145. 85.31 with step 1.0448 (0.0125) converged in 2 iterations.
146. 86.37 with step 1.0598 (0.0125) converged in 1 iterations.
147. 87.45 with step 1.0744 (0.0125) converged in 1 iterations.
148. 88.54 with step 1.0889 (0.0125) converged in 1 iterations.
149. 89.64 with step 1.1033 (0.0125) converged in 1 iterations.
150. 90.76 with step 1.1177 (0.0125) converged in 1 iterations.
151. 91.89 with step 1.1321 (0.0125) converged in 1 iterations.
152. 93.04 with step 1.1467 (0.0125) converged in 1 iterations.
153. 94.20 with step 1.1613 (0.0125) converged in 1 iterations.
154. 95.37 with step 1.1761 (0.0125) converged in 1 iterations.
155. 96.56 with step 1.1910 (0.0125) converged in 1 iterations.
156. 97.77 with step 1.2060 (0.0125) converged in 1 iterations.
157. 98.99 with step 1.2212 (0.0125) converged in 1 iterations.
158. 100.23 with step 1.2366 (0.0125) converged in 1 iterations.

Stability Certification

E0, _ = HARMONICSTIFFNESS(zeros(mdl.Ndofs,mdl.Ndofs), -2mdl.M,
    zeros(mdl.Ndofs, mdl.Ndofs), [1.], h);
indsh1 = [rinds[1:mdl.Ndofs]; iinds[1:mdl.Ndofs]];
E0 = collect(E0);
J = zeros(mdl.Ndofs*Nhc, mdl.Ndofs*Nhc);
stab = zeros(length(sols));
for iw in 1:length(sols)
    fun(Famp).jac(J, sols.u[iw], sols.p[iw])
    eVs = eigvals(J[indsh1,indsh1], sols.p[iw]*E0[indsh1,indsh1]);
    stab[iw] = sum(real(eVs).>=0);
end

Plot Forced Response

set_theme!(theme_latexfonts())
fsz = 18;
fig = Figure(fontsize=fsz);

ax = Axis(fig[1, 1], xlabel="Excitation Frequency (rad/s)",
    ylabel="Response (m)", yscale=log10);
scatterlines!(ax, sols.p./(stab.==0), norm.(Uh)/Famp)
scatterlines!(ax, sols.p./(stab.!=0), norm.(Uh)/Famp)

xlims!(ax, Om0, Om1)
    fig
Example block output

This page was generated using Literate.jl.