ModelingToolkit external base excitation

General Usage
1

hello
I’m trying to simulate a TMD system using ModelingToolkit.

This is the code i’m runing:

using ModelingToolkit
using DifferentialEquations
using LinearAlgebra
using Symbolics: scalarize

@variables t
D = Differential(t)

function Mass(; name, m = 1.0, xy = [0.0, 0.0], u = [0.0, 0.0])
    ps = @parameters m=m
    @variables pos(t)[1:2]=xy v(t)[1:2]=u
    eqs = scalarize(D.(pos) .~ v)
    ODESystem(eqs, t, [pos..., v...], ps; name)
end

function Spring(; name, k = 1e4, l = 1.0)
    ps = @parameters k=k l=l
    @variables x(t), dir(t)[1:2]
    ODESystem(Equation[], t, [x, dir...], ps; name)
end

function Damper(; name, c = 2)
    ps = @parameters c=c
    @variables x(t), dir(t)[1:2]
    ODESystem(Equation[], t, [x, dir...], ps; name)
end

function connect_element(element, x, y)
    [element.x ~ norm(scalarize(x .- y))
    scalarize(element.dir .~ scalarize(x .- y))]
end

function damper_force(damper)
    -damper.c .* scalarize(damper.dir) .* D.(damper.x)
end

function spring_force(spring)
    -spring.k .* scalarize(spring.dir) .* (spring.x - spring.l) ./ spring.x
end

# main mass
m1 = 1.0
xy1 = [0.0, 1.0]
k1 = 1e6
l1 = 1.0
c1 = 50.0

@named mass1 = Mass(m = m1, xy = xy1)
@named damper1 = Damper(c = c1)
@named spring1 = Spring(k = k1, l = l1)

# tmd mass
m2 = m1 * 0.05 # 5% of the main mass
xy2 = [0.0, 1.5] 
k2 = (main_mass_fn*2π)^2*m2
l2 = 0.5
c2 = 50.0
@named mass2 = Mass(m = m2, xy = xy2)
@named damper2 = Damper(c = c2)
@named spring2 = Spring(k = k2, l = l2)

base_pos = [0.0, 0.1]
g = [0.0, -9.81]

eqs = [connect_element(damper1, mass1.pos, base_pos)
            connect_element(spring1, mass1.pos, base_pos)
            connect_element(damper2, mass2.pos, mass1.pos)
            connect_element(spring2, mass2.pos, mass1.pos)
            scalarize(D.(mass2.v) .~ damper_force(damper2) / mass2.m .+ spring_force(spring2) / mass2.m .+ g)
            scalarize(D.(mass1.v) .~ damper_force(damper1) / mass1.m .+ spring_force(spring1) / mass1.m .- damper_force(damper2) / mass1.m .- spring_force(spring2) / mass1.m .+ g)]

@named _model = ODESystem(eqs, t, [ damper1.x; damper1.dir; spring1.x; spring1.dir; mass1.pos; 
                                        damper2.x; damper2.dir; spring2.x; spring2.dir; mass2.pos], [])

@named model = compose(_model, mass1, damper1, spring1, mass2, damper2, spring2)
sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0.0, 1.0))
solve(prob, Rosenbrock23())

my issue is, how to replace the static base position, base_pos = [0.0, 0.1],

to a function in time, base_pos(t)=sin(2π*f*t) for example.

thanks

2

Have you tried

base_pos = [0.0, sin(2π*f*t)]

?

More generally, you could make use of ModelingToolkitStandardLibrary.jl: A Standard Library for ModelingToolkit · ModelingToolkitStandardLibrary.jl in order to avoid having to define your own components for a simple system like this. Using the stdlib, you could use a Sine source as input.

You find a relevant tutorial here

In this tutorial, a Step is used as input rather than a Sine. For a system closer to a mass-spring-damper, see this tutorial instead

1 Like
3

thanks a lot, I’ll look into it.