Stress on Truss
In this example we'll plot an animation of the stress on some truss structure using GaphMakie.jl
and NetworkDynamics.jl
using NetworkDynamics
using OrdinaryDiffEq
using Graphs
using GraphMakie
using LinearAlgebra: norm, ⋅
using Printf
using CairoMakie
CairoMakie.activate!()
Definition of the dynamical system
Define dynamic model for NetworkDynamics
function edge_f!(e, v_s, v_d, (K, L), _)
Δd = view(v_d, 1:2) - view(v_s, 1:2)
nd = norm(Δd)
F = K * ( L - nd )
e[1:2] = F .* view(Δd, 1:2) ./ nd
end
function vertex_free!(du, u, edges, (M, γ, g), _)
F = sum(edges) - γ * view(u, 3:4) + [0, -M*g]
du[3:4] .= F ./ M
du[1:2] = view(u, 3:4)
end
function vertex_fixed!(du, u, edges, _, _)
du[1:2] .= 0.0
end
edge = StaticEdge(f = edge_f!, dim=2, coupling=:antisymmetric)
vertex_free = ODEVertex(f = vertex_free!, dim=4, sym=[:x, :y, :v, :w])
vertex_fix = ODEVertex(f = vertex_fixed!, dim=2, sym=[:x, :y])
Set up graph topology and initial positions.
N = 5
dx = 1.0
shift = 0.2
g = SimpleGraph(2*N + 1)
for i in 1:N
add_edge!(g, i, i+N); add_edge!(g, i, i+N)
if i < N
add_edge!(g, i+1, i+N); add_edge!(g, i, i+1); add_edge!(g, i+N, i+N+1)
end
end
add_edge!(g, 2N, 2N+1)
pos0 = zeros(Point2f, 2N + 1)
pos0[1:N] = [Point((i-1)dx,0) for i in 1:N]
pos0[N+1:2*N] = [Point(i*dx + shift, 1) for i in 1:N]
pos0[2N+1] = Point(N*dx + 1, -1)
fixed = [1,4] # set fixed vertices
create networkdynamics object
verts = ODEVertex[vertex_free for i in 1:nv(g)]
for i in fixed
verts[i] = vertex_fix # use the fixed vertex for the fixed
end
nd = network_dynamics(verts, edge, g);
write some auxiliary functions to map between state vector of DGL u
and graph data
x_idxs = [findfirst(isequal(Symbol("x_$i")), nd.syms) for i in 1:nv(g)]
"Set positions `pos` inside dgl state `u`"
function set_pos!(u, pos)
for (i, idx) in enumerate(x_idxs)
u[idx] = pos[i][1]
u[idx+1] = pos[i][2]
end
end
"Extract vector of Points from dgl state `u`"
function to_pos(u)
pos = Vector{Point2f}(undef, nv(g))
for (i, idx) in enumerate(x_idxs)
pos[i] = Point(u[idx], u[idx+1])
end
return pos
end
"Calculate load on edge for given state."
function get_load(u, p, t=0.0)
gd_nd = nd(u, p, t, GetGD) # exposes underlying graph data struct
force = Vector{Float64}(undef, ne(g))
pos = to_pos(u)
for (i,e) in enumerate(edges(g))
edgeval = get_edge(gd_nd, i)
fvec = Point(edgeval[1], edgeval[2])
dir = pos[dst(e)] .- pos[src(e)]
force[i] = sign(fvec ⋅ dir) * norm(fvec)
end
return force
end
Set parameters.
M = [10 for i in 1:nv(g)] # mass of the nodes
M[end] = 200 # heavy mass
gc = [9.81 for i in 1:nv(g)] # gravitational constant
γ = [200.0 for i in 1:nv(g)] # damping parameter
γ[end] = 100.0
L = [norm(pos0[src(e)] - pos0[dst(e)]) for e in edges(g)] # length of edges
K = [0.5e6 for i in 1:ne(g)] # spring constant of edges
# bundle parameters for NetworkDynamics
para = collect.((zip(M, γ, gc),
zip(K, L)))
Set initial state and solve the system
u0 = zeros(length(nd.syms))
set_pos!(u0, pos0)
tspan = (0.0, 12.0)
prob = ODEProblem(nd, u0, tspan, para)
sol = solve(prob, Tsit5());
Plot the solution
fig = Figure(size=(1000,550))
fig[1,1] = title = Label(fig, "Stress on truss", fontsize=30)
title.tellwidth = false
fig[2,1] = ax = Axis(fig)
# calculate some values for colorscaling
(fmin, fmax) = 0.3 .* extrema([(map(u->get_load(u, para), sol)...)...])
nlabels = [" " for i in 1:nv(g)]
nlabels[fixed] .= "Δ"
elabels = ["edge $i" for i in 1:ne(g)]
p = graphplot!(ax, g;
edge_width=4.0,
node_size=sqrt.(M)*3,
nlabels=nlabels,
nlabels_align=(:center,:top),
nlabels_fontsize=30,
elabels=elabels,
elabels_side=Dict(ne(g) => :right),
edge_color=[0.0 for i in 1:ne(g)],
edge_attr=(colorrange=(fmin,fmax),
colormap=:diverging_bkr_55_10_c35_n256))
hidespines!(ax); hidedecorations!(ax); p[:node_pos][]=to_pos(u0); ax.aspect = DataAspect()
limits!(ax, -0.1, pos0[end][1]+0.3, pos0[end][2]-0.5, 1.15)
# draw colorbar
fig[3,1] = cb = Colorbar(fig, p.plots[1].plots[1], label = "Axial force", vertical=false)
T = tspan[2]
fps = 30
trange = range(0.0, sol.t[end], length=Int(T * fps))
record(fig, "truss.mp4", trange; framerate=fps) do t
title.text = @sprintf "Stress on truss (t = %.2f )" t
p[:node_pos][] = to_pos(sol(t))
load = get_load(sol(t), para)
p.elabels = [@sprintf("%.0f", l) for l in load]
p.edge_color[] = load
end
This page was generated using Literate.jl.