|
| 1 | +#= |
| 2 | +# Stress on Truss |
| 3 | +In this exampe we'll simulate the time evolution of a truss structure consisting |
| 4 | +of joints and stiff springs. |
| 5 | +This example can be dowloaded as a normal Julia script [here](@__NAME__.jl). #md |
| 6 | +
|
| 7 | + |
| 8 | +
|
| 9 | +The mathematical model is quite simple: the vertices are point masses with |
| 10 | +positions $x$ and $y$ and velocities $v_x$ and $v_y$. |
| 11 | +For simplicity, we add some damping to the nodal motion based on the velocity. |
| 12 | +```math |
| 13 | +\begin{aligned} |
| 14 | +\dot{x} &= v_x\\ |
| 15 | +\dot{y} &= v_y\\ |
| 16 | +\dot{v}_x &= \frac{1}{M}\left(\sum F_x - γ\,v_x\right)\\ |
| 17 | +\dot{v}_y &= \frac{1}{M}\left(\sum F_y - γ\,v_y\right) - g |
| 18 | +\end{aligned} |
| 19 | +``` |
| 20 | +The vertices cannot absorb any torque, so the beams only exert forces in the direction of the beam. |
| 21 | +```math |
| 22 | +|F| = K\cdot(L - |Δd|) |
| 23 | +``` |
| 24 | +where $L$ is the nominal lenth and $|Δd|$ is the actual length of the beam. |
| 25 | +
|
| 26 | +We start by loading the necessary packages. |
| 27 | +=# |
| 28 | +using NetworkDynamics |
| 29 | +using OrdinaryDiffEqTsit5 |
| 30 | +using Graphs |
| 31 | +using GraphMakie |
| 32 | +using LinearAlgebra: norm |
| 33 | +using Printf |
| 34 | +using CairoMakie |
| 35 | +CairoMakie.activate!() |
| 36 | + |
| 37 | +#= |
| 38 | +## Definition of the dynamical system |
| 39 | +
|
| 40 | +We need 3 models: |
| 41 | +- a fixed vertex which cannot change its position, |
| 42 | +- a free vertex which can move, and |
| 43 | +- a beam which connects two vertices. |
| 44 | +=# |
| 45 | + |
| 46 | +function fixed_g(pos, x, p, t) |
| 47 | + pos .= p |
| 48 | +end |
| 49 | +vertex_fix = VertexFunction(g=fixed_g, psym=[:xfix, :yfix], outsym=[:x, :y], ff=NoFeedForward()) |
| 50 | +#= |
| 51 | +Here we need to specify the `ff` keyword manually, because NetworkDynamics cannot distinguish between |
| 52 | +`g(out, x, p, t)` (`NoFeedForwarwd`) and `g(out, in, p, t)` (`PureFeedForward()`) |
| 53 | +and guesses the latter when $\mathrm{dim}(x)=0$. |
| 54 | +=# |
| 55 | + |
| 56 | +function free_f(dx, x, Fsum, (M, γ, g), t) |
| 57 | + v = view(x, 1:2) |
| 58 | + dx[1:2] .= (Fsum .- γ .* v) ./ M |
| 59 | + dx[2] -= g |
| 60 | + dx[3:4] .= v |
| 61 | + nothing |
| 62 | +end |
| 63 | +vertex_free = VertexFunction(f=free_f, g=3:4, sym=[:vx=>0, :vy=>0, :x, :y], |
| 64 | + psym=[:M=>10, :γ=>200, :g=>9.81], insym=[:Fx, :Fy]) |
| 65 | +#= |
| 66 | +For the edge, we want to do something special. Later on, we want to color the edges according to the force they exert. |
| 67 | +Therefore, are interested in the absolut force rather than just the force vector. |
| 68 | +NetworkDynamics allows you to define so called `Observed` functions, which can recover additional states, |
| 69 | +so called observed, after the simulations. We can use this mechanis, to define a "recipe" for calculating the beam force |
| 70 | +based on the inputs (nodal positions) and the beam parameters. |
| 71 | +=# |
| 72 | + |
| 73 | +function edge_g!(F, pos_src, pos_dst, (K, L), t) |
| 74 | + dx = pos_dst[1] - pos_src[1] |
| 75 | + dy = pos_dst[2] - pos_src[2] |
| 76 | + d = sqrt(dx^2 + dy^2) |
| 77 | + Fabs = K * (L - d) |
| 78 | + F[1] = Fabs * dx / d |
| 79 | + F[2] = Fabs * dy / d |
| 80 | + nothing |
| 81 | +end |
| 82 | +function observedf(obsout, u, pos_src, pos_dst, (K, L), t) |
| 83 | + dx = pos_dst[1] .- pos_src[1] |
| 84 | + dy = pos_dst[2] .- pos_src[2] |
| 85 | + d = sqrt(dx^2 + dy^2) |
| 86 | + obsout[1] = K * (L - d) |
| 87 | + nothing |
| 88 | +end |
| 89 | +beam = EdgeFunction(g=AntiSymmetric(edge_g!), psym=[:K=>0.5e6, :L], outsym=[:Fx, :Fy], obsf=observedf, obssym=[:Fabs]) |
| 90 | + |
| 91 | +#= |
| 92 | +With the models define we can set up graph topology and initial positions. |
| 93 | +=# |
| 94 | +N = 5 |
| 95 | +dx = 1.0 |
| 96 | +shift = 0.2 |
| 97 | +g = SimpleGraph(2*N + 1) |
| 98 | +for i in 1:N |
| 99 | + add_edge!(g, i, i+N); add_edge!(g, i, i+N) |
| 100 | + if i < N |
| 101 | + add_edge!(g, i+1, i+N); add_edge!(g, i, i+1); add_edge!(g, i+N, i+N+1) |
| 102 | + end |
| 103 | +end |
| 104 | +add_edge!(g, 2N, 2N+1) |
| 105 | +pos0 = zeros(Point2f, 2N + 1) |
| 106 | +pos0[1:N] = [Point((i-1)dx,0) for i in 1:N] |
| 107 | +pos0[N+1:2*N] = [Point(i*dx + shift, 1) for i in 1:N] |
| 108 | +pos0[2N+1] = Point(N*dx + 1, -1) |
| 109 | +fixed = [1,4] # set fixed vertices |
| 110 | +nothing #hide |
| 111 | + |
| 112 | +#= |
| 113 | +Now can collect the vertex models and construct the Network object. |
| 114 | +=# |
| 115 | +verts = VertexFunction[vertex_free for i in 1:nv(g)] |
| 116 | +for i in fixed |
| 117 | + verts[i] = vertex_fix # use the fixed vertex for the fixed points |
| 118 | +end |
| 119 | +nw = Network(g, verts, beam) |
| 120 | +#= |
| 121 | +In order to simulate the system we need to initialize the state and parameter vectors. |
| 122 | +Some states and parameters are shared between all vertices/edges. Those have been allready set |
| 123 | +in their constructors. The free symbols are |
| 124 | +- `x` and `y` for the position of the free vertices, |
| 125 | +- `xfix` and `yfix` for the position of the fixed vertices, |
| 126 | +- `L` for the nominal length of the beams. |
| 127 | +=# |
| 128 | +s = NWState(nw) |
| 129 | +## set x/y and xfix/yfix |
| 130 | +for i in eachindex(pos0, verts) |
| 131 | + if i in fixed |
| 132 | + s.p.v[i, :xfix] = pos0[i][1] |
| 133 | + s.p.v[i, :yfix] = pos0[i][2] |
| 134 | + else |
| 135 | + s.v[i, :x] = pos0[i][1] |
| 136 | + s.v[i, :y] = pos0[i][2] |
| 137 | + end |
| 138 | +end |
| 139 | +## set L for edges |
| 140 | +for (i,e) in enumerate(edges(g)) |
| 141 | + s.p.e[i, :L] = norm(pos0[src(e)] - pos0[dst(e)]) |
| 142 | +end |
| 143 | +nothing #hide |
| 144 | + |
| 145 | +#= |
| 146 | +Lastly there is a special vertex at the end of the truss which has a higher mass and reduced damping. |
| 147 | +=# |
| 148 | +s.p.v[11, :M] = 200 |
| 149 | +s.p.v[11, :γ] = 100 |
| 150 | +nothing #hide |
| 151 | + |
| 152 | +#= |
| 153 | +No we have everything ready to build the ODEProblem and simulate the system. |
| 154 | +=# |
| 155 | +tspan = (0.0, 12.0) |
| 156 | +prob = ODEProblem(nw, uflat(s), tspan, pflat(s)) |
| 157 | +sol = solve(prob, Tsit5()) |
| 158 | +nothing #hide |
| 159 | + |
| 160 | +#= |
| 161 | +## Plot the solution |
| 162 | +
|
| 163 | +Plotting trajectories of points is kinda boring. So instead we're going to use |
| 164 | +[`GraphMakie.jl`](https://github.com/MakieOrg/GraphMakie.jl) to create a animation of the timeseries. |
| 165 | +=# |
| 166 | +fig = Figure(size=(1000,550)); |
| 167 | +fig[1,1] = title = Label(fig, "Stress on truss", fontsize=30) |
| 168 | +title.tellwidth = false |
| 169 | + |
| 170 | +fig[2,1] = ax = Axis(fig) |
| 171 | +ax.aspect = DataAspect(); |
| 172 | +hidespines!(ax); # no borders |
| 173 | +hidedecorations!(ax); # no grid, axis, ... |
| 174 | +limits!(ax, -0.1, pos0[end][1]+0.3, pos0[end][2]-0.5, 1.15) # axis limits to show full plot |
| 175 | + |
| 176 | +## get the maximum force during the simulation to get the color scale |
| 177 | +## It is only possible to access `:Fabs` directly becaus we've define the observable function for it! |
| 178 | +(fmin, fmax) = 0.3 .* extrema(Iterators.flatten(sol(sol.t, idxs=eidxs(nw, :, :Fabs)))) |
| 179 | +p = graphplot!(ax, g; |
| 180 | + edge_width = 4.0, |
| 181 | + node_size = 3*sqrt.(try s.p.v[i, :M] catch; 10.0 end for i in 1:nv(g)), |
| 182 | + nlabels = [i in fixed ? "Δ" : "" for i in 1:nv(g)], |
| 183 | + nlabels_align = (:center,:top), |
| 184 | + nlabels_fontsize = 30, |
| 185 | + elabels = ["edge $i" for i in 1:ne(g)], |
| 186 | + elabels_side = Dict(ne(g) => :right), |
| 187 | + edge_color = [0.0 for i in 1:ne(g)], |
| 188 | + edge_attr = (colorrange=(fmin,fmax), |
| 189 | + colormap=:diverging_bkr_55_10_c35_n256)) |
| 190 | + |
| 191 | +## draw colorbar |
| 192 | +fig[3,1] = cb = Colorbar(fig, get_edge_plot(p), label = "Axial force", vertical=false) |
| 193 | + |
| 194 | +T = tspan[2] |
| 195 | +fps = 30 |
| 196 | +trange = range(0.0, sol.t[end], length=Int(T * fps)) |
| 197 | +record(fig, "truss.mp4", trange; framerate=fps) do t |
| 198 | + title.text = @sprintf "Stress on truss (t = %.2f )" t |
| 199 | + s_at_t = NWState(sol, t) |
| 200 | + for i in eachindex(pos0) |
| 201 | + p[:node_pos][][i] = (s_at_t.v[i, :x], s_at_t.v[i, :y]) |
| 202 | + end |
| 203 | + notify(p[:node_pos]) |
| 204 | + load = s_at_t.e[:, :Fabs] |
| 205 | + p.edge_color[] = load |
| 206 | + p.elabels = [@sprintf("%.0f", l) for l in load] |
| 207 | + fig |
| 208 | +end |
0 commit comments