Non conservation of energy in ModelingToolkit with flexible bouncing ball example

General Usage
,
1

Question

Based on the examples in the ModelingToolkit documentation I’m trying to expand a bit and model the impact behavior between two masses, one being connected to a base through a spring.

To make it simple, an example with only two masses colliding yields expected results, but as soon as I add the spring, this leads to a non-conservative behavior.

The only obvious thing I can think about is that the continuous event is only acting at the very specific gap == 0 condition while it should act on the gap <= 0 inequality condition.

What would be the proper way to account for these inequality conditions ?

Side uncertain note: this seemed to behave as expected in previous versions but I can’t remember which.

Environment

[e9467ef8] GLMakie v0.13.13
[961ee093] ModelingToolkit v11.35.2
[16a59e39] ModelingToolkitStandardLibrary v2.29.5
[1dea7af3] OrdinaryDiffEq v7.1.2

on julia 1.12.6 windows.

Two masses without spring

using ModelingToolkit
import ModelingToolkit: connect, t_nounits as t, D_nounits as D
using ModelingToolkitStandardLibrary
using GLMakie
using OrdinaryDiffEq

@named mass1 = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition.Mass(m=1.)
@named mass2 = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition.Mass(m=0.5)
@named fix = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition.Fixed(s_0=0.)
@named spring = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition.Spring(k=10., l=1.)
@variables gap(t) p(t) E(t) l(t)
@parameters e_newton = 1.


@mtkcompile twomass = System(
    [
        gap ~ mass2.s - mass1.s,
        p ~ mass1.v * mass1.m + mass2.v * mass2.m,
        E ~ 1/2 * mass1.m * mass1.v^2 + 1/2 * mass2.m * mass2.v^2,
    ],
    t;
    systems=[mass1, mass2],
    continuous_events=[gap ~ 0] => [
        mass1.v ~ ((mass1.m - 1 * mass2.m) * Pre(mass1.v) + (1 + 1) * mass2.m * Pre(mass2.v)) / (mass1.m + mass2.m)
        mass2.v ~ ((mass2.m - 1 * mass1.m) * Pre(mass2.v) + (1 + 1) * mass1.m * Pre(mass1.v)) / (mass1.m + mass2.m)
    ]
)

u0 = Dict{Num,Any}([
    mass1.s => 1.,
    mass1.v => 0.,
    mass2.s => 2.,
    mass2.v => -1.,
]
)

prob = ODEProblem(twomass, u0, (0.0, 10.0), fully_determined=true)

sol = solve(prob, Tsit5())
begin
    f = Figure()
    ax = Axis(f[1, 1])
    plot!(sol, idxs=[twomass.mass1.s, twomass.mass2.s, twomass.gap])
    axislegend(ax)
    ax = Axis(f[2, 1])
    plot!(sol, idxs=[twomass.mass1.v, twomass.mass2.v])
    axislegend(ax)
    ax = Axis(f[3, 1])
    plot!(sol, idxs=[twomass.E])
    axislegend(ax)
    f
end

which gives

image
image602×482 28.4 KB

Two masses with spring

@mtkcompile twomassspring = System(
    [
        gap ~ mass2.s - mass1.s,
        p ~ mass1.v * mass1.m + mass2.v * mass2.m,
        l ~ spring.flange_a.s - spring.flange_b.s,
        E ~ 1/2 * mass1.m * mass1.v^2 + 1/2 * mass2.m * mass2.v^2 + 1/2 * spring.k * (l - spring.l)^2,
        connect(fix.flange, spring.flange_b),
        connect(mass1.flange, spring.flange_a),
    ],
    t;
    systems=[mass1, mass2, fix, spring],
    continuous_events=[gap ~ 0] => [
        mass1.v ~ ((mass1.m - 1 * mass2.m) * Pre(mass1.v) + (1 + 1) * mass2.m * Pre(mass2.v)) / (mass1.m + mass2.m)
        mass2.v ~ ((mass2.m - 1 * mass1.m) * Pre(mass2.v) + (1 + 1) * mass1.m * Pre(mass1.v)) / (mass1.m + mass2.m)
    ]
)


u0 = Dict{Num,Any}([
    mass1.s => 1.,
    mass1.v => 0.,
    mass2.s => 2.,
    mass2.v => -1.,
]
)

prob = ODEProblem(twomassspring, u0, (0.0, 10.0), fully_determined=true)

sol = solve(prob, Tsit5())
begin
    f = Figure()
    ax = Axis(f[1, 1])
    plot!(sol, idxs=[twomassspring.mass1.s, twomassspring.mass2.s, twomassspring.gap, twomassspring.l])
    axislegend(ax)
    ax = Axis(f[2, 1])
    plot!(sol, idxs=[twomassspring.mass1.v, twomassspring.mass2.v])
    axislegend(ax)
    ax = Axis(f[3, 1])
    plot!(sol, idxs=[twomassspring.E])
    axislegend(ax)
    f
end

which yields

image
image602×482 39.1 KB

Obviously, the energy conservation is not respected, however events are properly catched.