Ok so I think I managed to find the missing part. It does not seem to be enough to state a value for the external forcing directly, so I created an additional Constant3 block to handle such forcings. This, with some corrections based on ModelingToolkitStandardLibrary.jl/src/Mechanical/TranslationalPosition/TranslationalPosition.jl at v2.3.4 · SciML/ModelingToolkitStandardLibrary.jl · GitHub, gives
# Example to demonstrate the trajectory of a moving mass in 3D using ModelingToolkit and modified ModelingToolkitStandardLibrary components
using ModelingToolkit, LinearAlgebra, DifferentialEquations
using ModelingToolkitStandardLibrary.Blocks: RealInput
using GeometryBasics, GLMakie
@variables t
D = Differential(t)
@connector Flange3 begin
s(t)[1:3] = zeros(3), [description = "Translational position of the node in [m/s]"]
f(t)[1:3] = zeros(3), [description = "Force entering the node in [N]", connect = Flow]
end
@mtkmodel Force3 begin
@parameters begin
f_0[1:3] = zeros(3)
end
@components begin
node = Flange3()
f = RealInput(; nin=3, u_start=f_0)
end
@equations begin
collect(node.f .~ -f.u)...
end
end
@mtkmodel Fixed3 begin
@parameters begin
s = zeros(3)
end
@components begin
node = Flange3(; s=s)
end
@equations begin
collect(node.s .~ s)...
end
end
@mtkmodel Free3 begin
@components begin
node = MechanicalPort3()
end
@variables begin
f(t)[1:3] = zeros(3)
end
@equations begin
collect(node.f .~ f)...
end
end
@mtkmodel Constant3 begin
@components begin
output = RealOutput(; nout=3)
end
@parameters begin
k[1:3] = zeros(3), [description = "Constant output value of block"]
end
@equations begin
collect(output.u .~ k)...
end
end
@mtkmodel Mass3 begin
# Definition of the model icon
#@icon "..."
# Definition of the model parameters
@parameters begin
m, [description = "Mass of the moving rigid body in [kg]"]
end
@variables begin
(s(t))[1:3] = zeros(3), [description = "Positions in [m] along X, Y and Z axes of the rigid body"]
(v(t))[1:3] = zeros(3), [description = "Velocities in [m/s] along X, Y and Z axes of the rigid body"]
(f(t))[1:3] = zeros(3), [description = "Forces in [N] along X, Y and Z axes entering the rigid body"]
end
@components begin
node = Flange3(; s=s)
end
@equations begin
collect(node.s .~ s)...
collect(node.f .~ f)...
collect(D.(s) .~ v)...
collect(D.(v) .~ f ./ m)...
end
end
@mtkmodel SpringDamper3 begin
@parameters begin
k, [description = "Spring stiffness in [N/m]"]
c, [description = "Damper viscous damping in [N/m]"]
l0, [description = "Spring length at rest in [m]"]
end
@components begin
origin_node = Flange3()
target_node = Flange3()
end
@variables begin
l(t), [description = "Spring length in [m]"]
v(t), [description = "Damper velocity in [m]"]
f(t)[1:3]
dir(t)[1:3], [description = "Spring direction"]
end
@equations begin
D(l) ~ v
l ~ norm(collect(target_node.s .- origin_node.s))
collect(dir .~ collect(target_node.s .- origin_node.s) ./ l)...
collect(f .~ k .* collect(dir) .* (l - l0) + c .* collect(dir) .* v)...
collect(origin_node.f .~ .+f)...
collect(target_node.f .~ .-f)...
end
end
@named mass = Mass3(m=1, s=[0.0, 0.0, 0.0])
@named springdamper = SpringDamper3(k=10, l0=0.5, c=2)
@named fix = Fixed3(s=[1.0; -1.0; 2.0])
@named g = Force3()
@named sig = Constant3(k=[0.0; 0.0; -9.81])
eqs = Equation[
ModelingToolkit.connect(sig.output, g.f)
ModelingToolkit.connect(g.node, mass.node)
ModelingToolkit.connect(springdamper.target_node, mass.node)
ModelingToolkit.connect(springdamper.origin_node, fix.node)
]
@named model = ODESystem(eqs, t, [], []; systems=[mass, fix, springdamper, g, sig])
sys = structural_simplify(model)
At that point, the structural_simplify works and gives as expected 6 degrees of freedom, 3 positions and 3 velocities.
However, if I try to solve that, it fails on some sort of broadcasting.
initial_conditions = reduce(vcat, [collect(mass.s .=> [0.0, 1.0, -0.1]), collect(mass.v .=> [0.5, 1.0, 0.0])])
prob = ODEProblem(sys, initial_conditions, (0.0, 100.0))
sol = solve(prob, Rodas5(), abstol=1e-8, reltol=1e-8, saveat=0:0.1:200)
leads to
ERROR: MethodError: no method matching +(::Vector{Float64}, ::Float64)
For element-wise addition, use broadcasting with dot syntax: array .+ scalar
Closest candidates are:
+(::Any, ::Any, ::Any, ::Any…)
@ Base operators.jl:578
+(::T, ::T) where T<:Union{Float16, Float32, Float64}
@ Base float.jl:408
+(::AbstractAlgebra.MatrixElem{T}, ::Union{AbstractFloat, Integer, Rational}) where T<:AbstractAlgebra.NCRingElement
@ChrisRackauckas I’m using ModelingToolkit v8.73.2, do you think it is the same bug as in Broadcast inside of `@mtkmodel` macro causes precompilation failures when put into a package - #11 by BambOoxX or something else ?