/
tipforce.jl
140 lines (106 loc) · 3.53 KB
/
tipforce.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
# # [Cantilever with a Tip Force](@id tipforce)
#
# This example shows how to predict the behavior of a cantilever beam that is subjected
# to a constant shear force at the tip.
#
# 
#
#-
#md # !!! tip
#md # This example is also available as a Jupyter notebook:
#md # [`tipforce.ipynb`](@__NBVIEWER_ROOT_URL__/examples/tipforce.ipynb).
#-
using GXBeam, LinearAlgebra
L = 1
EI = 1e6
## shear force (applied at end)
λ = 0:0.5:16
p = EI/L^2
P = λ*p
## create points
nelem = 16
x = range(0, L, length=nelem+1)
y = zero(x)
z = zero(x)
points = [[x[i],y[i],z[i]] for i = 1:length(x)]
## index of endpoints of each beam element
start = 1:nelem
stop = 2:nelem+1
## compliance matrix for each beam element
compliance = fill(Diagonal([0, 0, 0, 0, 1/EI, 0]), nelem)
## create assembly of interconnected nonlinear beams
assembly = Assembly(points, start, stop, compliance=compliance)
## pre-initialize system storage
system = StaticSystem(assembly)
## run an analysis for each prescribed tip load
states = Vector{AssemblyState{Float64}}(undef, length(P))
for i = 1:length(P)
## create dictionary of prescribed conditions
prescribed_conditions = Dict(
## fixed left side
1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0,
theta_z=0),
## shear force on right tip
nelem+1 => PrescribedConditions(Fz = P[i])
)
## perform a static analysis
_, states[i], converged = static_analysis!(system, assembly;
prescribed_conditions=prescribed_conditions)
end
#!jl nothing #hide
# ## Comparison with Analytical Results
#
# The analytical solution to this problem has been presented by several authors. Here
# we follow the solution by H. J. Barten in "On the Deflection of a Cantilever Beam",
# after incorporating the corrections they submitted for finding the tip angle.
#
import Elliptic
δ = range(pi/4, pi/2, length=10^5)[2:end-1]
k = @. cos(pi/4)/sin(δ)
λ_a = @. (Elliptic.F(pi/2, k^2) - Elliptic.F(δ, k^2))^2
θ_a = @. 2*(pi/4 - acos(k))
ξ_a = @. sqrt(2*sin(θ_a)/λ_a) .- 1
η_a = @. 1-2/sqrt(λ_a)*(Elliptic.E(pi/2, k^2) - Elliptic.E(δ, k^2))
#!jl nothing #hide
#
# Plotting the results reveals that the analytical and computational solutions show
# excellent agreement.
#
using Plots
#md using Suppressor #hide
pyplot()
#!jl nothing #hide
#-
#md @suppress_err begin #hide
u = [states[i].points[end].u[1] for i = 1:length(P)]
θ = [states[i].points[end].theta[2] for i = 1:length(P)]
w = [states[i].points[end].u[3] for i = 1:length(P)]
## set up the plot
plot(
xlim = (0, 16),
xticks = 0:1:16,
xlabel = "Nondimensional Force \$\\left(\\frac{PL^2}{EI}\\right)\$",
ylim = (0, 1.2),
yticks = 0.0:0.2:1.2,
ylabel = "Nondimensional Tip Displacements",
legend=:bottomright,
grid = false,
overwrite_figure=false
)
plot!([0], [0], color=:black, label="Analytical")
scatter!([0], [0], color=:black, label="GXBeam")
plot!([0], [0], color=1, label="\$ \\theta/(\\pi/2) \$")
plot!([0], [0], color=2, label="Vertical \$\\left(w/L\\right)\$")
plot!([0], [0], color=3, label="Horizontal \$\\left(-u/L\\right)\$")
plot!(λ_a, θ_a*2/pi, color=1, label="")
scatter!(λ, -4*atan.(θ/4)*2/pi, color=1, label="")
plot!(λ_a, η_a, color=2, label="")
scatter!(λ, w/L, color=2, label="")
plot!(λ_a, -ξ_a, color=3, label="")
scatter!(λ, -u/L, color=3, label="")
plot!(show=true) #!nb
#md savefig("../assets/tipforce-displacement.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide
#md # 