Mass dashpot spring

Free oscillations

using StructuralDynamicsODESolvers, Plots, LinearAlgebra

k  = 2 ; m  = .5 ;  c = 0 ;
u0 = 1 ; v0 = 0 ;

M = m * ones(1, 1)
C = c * ones(1, 1)
K = k * ones(1, 1)
R = zeros(1)

sys = SecondOrderAffineContinuousSystem(M, C, K, R)

U₀ = u0 * ones(1); V₀ = v0 * ones(1);

ivp_free = InitialValueProblem(sys, (U₀, V₀))

NSTEPS = 1000 ;
Δt = 0.005 ;

alg = Bathe(Δt = Δt )
sol = solve(ivp_free, alg, NSTEPS=NSTEPS);

The following command is the same as plot(times(sol), displacements(sol, 1)).

plot(sol, vars=(0, 1))

Forced oscillations

Problem definition

Let us consider now a forcing term $f(t) = A_f \sin(ω_f \cdot t)$

ωN = sqrt(k/m)
ωf = ωN * 2
Af = 10.0
R  = [ [ Af * sin(ωf * Δt * (i-1) ) ] for i in 1:NSTEPS+1];

Second order problem resolution

X   = nothing # state constraints are ignored
B   = ones(1, 1)
sys = SecondOrderConstrainedLinearControlContinuousSystem(M, C, K, B, X, R)

ivp_forced_secOrder = InitialValueProblem(sys, (U₀, V₀))

alg = Bathe(Δt = Δt )
sol_secOrder = solve(ivp_forced_secOrder, alg, NSTEPS=NSTEPS);

First order homogeneization formulation

The problem can be re-formulated as a first order and homogeneous one given by

\[\left\{ \begin{array}{l} \dot{u} = v \\ \dot{v} = -\omega_N^2 u + u_f/m \\ \dot{u_f} = v_f \\ \dot{v_f} = -\omega_f^2 u_f \end{array} \right.\]

The new vector of variables is

\[\textbf{x} = [ u, v, u_f, v_f ]^T\]

K = [     0 1     0 0 ;
      -ωN^2 0     1/m 0 ;
          0 0     0 1 ;
          0 0 -ωf^2 0 ] ;

C = -Diagonal(ones(4))
M = zeros(4,4)
R = zeros(4)

sys = SecondOrderAffineContinuousSystem(M, C, K, R)

U₀ = [u0; v0; 0; ωf*Af ] ;

ivp_forced_firOrder = InitialValueProblem(sys, (U₀, U₀) )

alg = BackwardEuler(Δt = Δt )
sol_firOrderA = solve(ivp_forced_firOrder, alg, NSTEPS=NSTEPS);

NSTEPS = NSTEPS*3 ; alg = BackwardEuler(Δt = Δt/3.0 )
sol_firOrderB = solve(ivp_forced_firOrder, alg, NSTEPS=NSTEPS);

The solution obtained is

plot(sol_secOrder, vars=(0, 1), xlab="time" )
plot!(sol_firOrderA, vars=(0, 1))
plot!(sol_firOrderB, vars=(0, 1))