Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add translational library #85

Closed
wants to merge 15 commits into from
Closed

Conversation

kianmolani
Copy link

This PR is linked to issue #84

@@ -0,0 +1,17 @@
{
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can add .vscode and .DS_Store in global git ignore setting. Check out https://sebastiandedeyne.com/setting-up-a-global-gitignore-file/

@@ -0,0 +1,36 @@
import ModelingToolkitStandardLibrary.Blocks
include("/Users/kianmolani/Documents/Project Files/Codebases/ModelingToolkitStandardLibrary.jl/src/Mechanical/Translational/Translational.jl")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line should be removed.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How can I add the Translational module to the package and be able to import it using use ModelingToolkitStandardLibrary.Mechanical.Translational? Do I have to first merge this PR? I have tried some possible solutions but without success.

import ModelingToolkitStandardLibrary.Blocks
include("/Users/kianmolani/Documents/Project Files/Codebases/ModelingToolkitStandardLibrary.jl/src/Mechanical/Translational/Translational.jl")
using ModelingToolkit, OrdinaryDiffEq, Test
using Plots
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't import Plots in test. CI cannot actually see plots.

export Flange
include("utils.jl")

export Fixed, Inertia, Spring, Damper, IdealGear
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
export Fixed, Inertia, Spring, Damper, IdealGear
export Fixed, Mass, Spring, Damper, IdealGear

module Translational

using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq
using OffsetArrays
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove not needed packages

@ValentinKaisermayer
Copy link
Contributor

Doc strings are missing and actual tests

@baggepinnen
Copy link
Contributor

Have a look at https://github.com/baggepinnen/AutomaticDocstrings.jl to quickly generate docstring stubs, helps a lot when you have a lot of docstrings to generate

Copy link
Author

@kianmolani kianmolani left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will delete the output folder before the PR is merged - I've just included the figures here for your reference (though lmk if there's a better way to share output results and plots). Also, I don't know whether we want to run tests against Modelica simulation data. I certainly think it's useful as an initial sanity check, but don't know if it's appropriate to include as part of continuous testing. If we do decide to continue testing against Modelica simulation data, I also think that the tests themselves can be better implemented than what's been done here.

Copy link
Contributor

@ValentinKaisermayer ValentinKaisermayer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add the connectors in the doc strings as well? And name the states s(t) to be consistent with the rest of the library.

If you want to add images, simply drag them into a comment or an issue.

Base.@doc """
Flange(;name)

1-dim. translational flange of a shaft.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
1-dim. translational flange of a shaft.
1-dim. translational flange.


# States:
- `s`: [m] Absolute position of flange
- `f`: [N] Cut force in the flange
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
- `f`: [N] Cut force in the flange
- `f`: [N] Cut force into the flange

Base.@doc """
Support(;name)

Support/housing of a 1-dim. translational shaft
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Support/housing of a 1-dim. translational shaft
Support/housing 1-dim. translational flange.


# States:
- `s`: [m] Absolute position of the support/housing
- `f`: [N] Reaction force in the support/housing
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
- `f`: [N] Reaction force in the support/housing
- `f`: [N] Cut force into the flange

"""
PartialCompliant(;name, s_rel_start=0.0, f_start=0.0)

Partial model for the compliant connection of two translational 1-dim. shaft flanges.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Partial model for the compliant connection of two translational 1-dim. shaft flanges.
Partial model for the compliant connection of two translational 1-dim. flanges.

"""
PartialCompliantWithRelativeStates(;name, s_rel_start=0.0, v_rel_start=0.0, a_rel_start=0.0, f_start=0.0)

Partial model for the compliant connection of two translational 1-dim. shaft flanges where the relative position and speed are used as preferred states
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Partial model for the compliant connection of two translational 1-dim. shaft flanges where the relative position and speed are used as preferred states
Partial model for the compliant connection of two translational 1-dim. flanges.

@@ -0,0 +1,58 @@
import ModelingToolkitStandardLibrary.Blocks
include(pwd() * "/src/Mechanical/Translational/Translational.jl")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
include(pwd() * "/src/Mechanical/Translational/Translational.jl")
using Modelingtoolkitstandardlibrary.Mechanical.Translational

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As written, this line returns the error Translational not defined. Is this because the Translational module has not been added to the package using the proper procedures (that is, using JL's package manager)? Note that when creating the Translational library, I merely added a new folder to the repo and proceeded to add the requisite lines of code.

@test mass_v_julia[end] ≈ 0 atol = 1e-3 # all energy has dissipated
@test mean(broadcast(abs, (mass_s_julia - mass_s_modelica))) ≈ 0 atol = 1e-5 # jl and modelica s vs. t results are approx equal
@test mean(broadcast(abs, (mass_v_julia - mass_v_modelica))) ≈ 0 atol = 1e-4 # jl and modelica v vs. t results are approx equal
@test mean(broadcast(abs, (mass_a_julia - mass_a_modelica))) ≈ 0 atol = 1e-2 # jl and modelica a vs. t results are approx equal
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you make the test self-sufficient?

Copy link
Author

@kianmolani kianmolani Jun 24, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can't think of a way to do this without importing the Modelica simulation results. Please let me know if you see a better way. What are your thoughts in general about running these sort of tests? Do you think they better serve as one-offs for sanity checking (i.e. not to be repeated again), or for continuous testing? Please see my original comment on this issue: #85 (review)

"""
PartialElementaryTwoFlangesAndSupport2(;name, use_support=false)

Partial model for a component with two translational 1-dim. shaft flanges and a support used for textual modeling, i.e., for elementary models
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Partial model for a component with two translational 1-dim. shaft flanges and a support used for textual modeling, i.e., for elementary models
Partial model for a component with two translational 1-dim. flanges and a support used for textual modeling, i.e., for elementary models

@@ -0,0 +1,21 @@
"""
Library to model 1-dimensional, translational mechanical systems
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Library to model 1-dimensional, translational mechanical systems
Library to model 1-dimensional, translational mechanical components.

@kianmolani
Copy link
Author

Can you add the connectors in the doc strings as well? And name the states s(t) to be consistent with the rest of the library.

If you want to add images, simply drag them into a comment or an issue.

Thank you for the suggestions, Valentin. The connectors in utils.js already have doc strings added (I mirrored the documentation from the same file in the Rotational lib). Are you referring to something else? Also, what do you mean about naming the states s(t)? The position state is already denoted by s.

@ValentinKaisermayer
Copy link
Contributor

I meant a # Connectors: section in the doc strings, listing the connectors of the components. You can have a look at the doc strings in the other modules.

@kianmolani
Copy link
Author

I meant a # Connectors: section in the doc strings, listing the connectors of the components. You can have a look at the doc strings in the other modules.

Okay, that's clear. Will do, and thanks!

@kianmolani kianmolani marked this pull request as ready for review July 11, 2022 12:48
@kianmolani
Copy link
Author

There appear to be issues with our simulated results for a mass-spring-damper model. Julia and Modelica simulation results do not match, and the difference is very apparent. I've attached below screenshots of Julia and Modelica simulation results of our mass position vs. time for underdamped, overdamped, and critically damped systems. For the underdamped system, the difference appears to be accounted for with an offset of about 0.5, as well as a slight phase shift. For the overdamped and critically damped systems, the solutions possess different equilibrium states.

Screen Shot 2022-07-14 at 11 56 20 AM

Screen Shot 2022-07-14 at 11 58 02 AM

Screen Shot 2022-07-14 at 11 59 15 AM

I've ensured that both Modelica and Julia simulations are initialized in the same way. Note that simulation results between Modelica and Julia for a mass-damper model are very close (tested across a range of parameter values).

This last screenshot shows how I've defined my system connections for a mass-spring-damper model:

Screen Shot 2022-07-14 at 9 46 16 AM

Comment on lines 25 to 28
@named fixed = Fixed(s0=s0)
@named mass = Mass(m=m, s_start=s_start, v_start=v_start)
@named damper = Damper(d=d)
@named spring = Spring(c=c, s_rel0=s_rel0)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
@named fixed = Fixed(s0=s0)
@named mass = Mass(m=m, s_start=s_start, v_start=v_start)
@named damper = Damper(d=d)
@named spring = Spring(c=c, s_rel0=s_rel0)
@named fixed = Fixed(; s0)
@named mass = Mass(; m, s_start, v_start)
@named damper = Damper(; d)
@named spring = Spring(; c, s_rel0)

In case you didn't know, if your variable shares the same name as a keyword argument (a common case), you do not need to repeat the name. You do need to have the semi-colon first though to indicate that keyword arguments are coming.

@baggepinnen
Copy link
Contributor

Could the difference between julia and modelica lie in the initialization? Is there a difference already at the first point saved in the solution?

@kianmolani
Copy link
Author

kianmolani commented Jul 15, 2022

Could the difference between julia and modelica lie in the initialization? Is there a difference already at the first point saved in the solution?

I believe that they've been initialized in the same way. Though I will triple-check this (the solutions match for the mass-damper system, tested across a range of param values, so perhaps the culprit here is the spring). I've attached the Modelica simulation settings below, though I don't know if they'd be of use to you (the simulation interval settings are the same, and so is the integration tolerance level).

modelica_sim_setup

@baggepinnen
Copy link
Contributor

Let me see if I understand this correctly, the fixed point is at s = 4.5 and the spring has resting length s_rel0 = 1, then I'd expect the steady-state position to be either 5.5 or 3.5 depending on the direction the spring is pointing in. It thus looks like the Julia solution is correct since it comes to rest at 3.5? I assume there's no gravity or other external forces acting anywhere? Is the length of the spring and all other params the same in modelica?

@kianmolani
Copy link
Author

kianmolani commented Jul 15, 2022

Let me see if I understand this correctly, the fixed point is at s = 4.5 and the spring has resting length s_rel0 = 1, then I'd expect the steady-state position to be either 5.5 or 3.5 depending on the direction the spring is pointing in. It thus looks like the Julia solution is correct since it comes to rest at 3.5? I assume there's no gravity or other external forces acting anywhere? Is the length of the spring and all other params the same in modelica?

No other forces and all other params confirmed to be the same.

The issue might be that Modelica takes in an additional parameter that we do not: mass length (set to 1 in this simulation). Because Modelica also has the parameter s_start (the initial value of the absolute position of the mass), and because this parameter is set to 3, it means that Modelica is tracking the center of this mass. Therefore, the equilibrium state at x = 3 makes sense for the Modelica simulation.

However, again, we do not take in a parameter for mass length, and so I guess that we're treating the mass as a point mass? In any case, initializing s_start to also equal 3 in our simulation (given that we do not have a mass length) might be problematic in that there's a discontinuity in the connections. I've drawn two schematics to illustrate this point:

For Modelica:

image

For Julia:

image

@kianmolani
Copy link
Author

Indeed, I've just confirmed that setting the spring length to 1.5m in our Julia simulation in order to account for this "gap" yields matching simulation results.

@kianmolani
Copy link
Author

Perhaps mass length is a property that we should include as well.

@codecov
Copy link

codecov bot commented Jul 18, 2022

Codecov Report

Merging #85 (f7f188e) into main (db37027) will decrease coverage by 5.56%.
The diff coverage is 0.00%.

@@            Coverage Diff             @@
##             main      #85      +/-   ##
==========================================
- Coverage   68.59%   63.03%   -5.57%     
==========================================
  Files          23       26       +3     
  Lines         952     1036      +84     
==========================================
  Hits          653      653              
- Misses        299      383      +84     
Impacted Files Coverage Δ
src/Mechanical/Translational/components.jl 0.00% <0.00%> (ø)
src/Mechanical/Translational/sources.jl 0.00% <0.00%> (ø)
src/Mechanical/Translational/utils.jl 0.00% <0.00%> (ø)

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

@YingboMa
Copy link
Member

You should remove .vscode and .DS_Store, then the PR is good to go after formatting.

Delete .DS_Store

Delete .DS_Store

Delete .DS_Store

Delete .DS_Store
@kianmolani
Copy link
Author

kianmolani commented Jul 26, 2022

You should remove .vscode and .DS_Store, then the PR is good to go after formatting.

All .vscode and .DS_Store files have now been deleted.

@YingboMa
Copy link
Member

YingboMa commented Aug 8, 2022

Format failure.

Copy link
Contributor

@ValentinKaisermayer ValentinKaisermayer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Other than the docs it looks good. 👍


Let's start by importing what we need.

```@example
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
```@example
```@example translational


Now let's initialize our system state and define our components and connections. Here, `s0` is the fixed offset position of our flange housing, `s_rel0` is our unstretched spring length, and `s_start` and `v_start` are the initial values of the absolute position and linear velocity of our sliding mass, respectively. Finally, `m`, `c`, and `d` are our mass, spring, and damping constants, respectively. All units are in SI units.

```@parameters t
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
```@parameters t
```@example translational
@parameters t


Before solving our ODE system, let us ensure that the condition for underdamped motion is met.

```
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
```
```@example translational


Now let's solve our ODE system and plot our results (note that we are plotting the position of our mass as a function of time).

```@named model = ODESystem(connections, t, systems=[fixed, mass, damper, spring])
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
```@named model = ODESystem(connections, t, systems=[fixed, mass, damper, spring])
```@example translational
@named model = ODESystem(connections, t, systems=[fixed, mass, damper, spring])


We'll now demonstrate a mass-spring-damper system that is overdamped. For such systems, `d > sqrt(4 * m * c)`, and the zeros of our characteristic polynomial are real. Let us re-initialize our `m`, `c`, and `d` constants accordingly, and re-define our components and connections.

```
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
```
```@example translational


Let us ensure that the condition for overdamped motion is met.

```
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
```
```@example translational


Let's solve our ODE system and plot our results. Observe that for overdamped systems, motion is non-oscillatory.

```@named model = ODESystem(connections, t, systems=[fixed, mass, damper, spring])
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
```@named model = ODESystem(connections, t, systems=[fixed, mass, damper, spring])
```@example translational
@named model = ODESystem(connections, t, systems=[fixed, mass, damper, spring])


Our last demonstration will be of a mass-spring-damper system that is critically damped. The condition for critically damped motion is `d = sqrt(4 * m * c)`, where the roots of our characteristic polynomial are both equal to `d/2m`. Again, let us re-initialize and re-define our system state, components, and connections, solve our ODE system, and plot our results. Observe that for critically damped systems, motion is non-oscillatory.

```
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
```
```@example translational

@bradcarman
Copy link
Contributor

Hi, I'm reviewing the code and I noticed that the Flange connector uses absolute position for the state variable. Might I suggest changing this to velocity? The reason for this is so we have less matching parameters to set. Let me show an example...

The below code requires that each of the 3 components have a position parameter set, and s_rel0 is really just s0-s_start, so it can already be seen as redundant.

s0 = 4.5 # fixed offset position of flange housing
s_start = 3 # initial value of absolute position of sliding mass
s_rel0 = 1.5 # unstretched spring length
v_start = 10 # initial value of absolute linear velocity of sliding mass

@named fixed = Fixed(; s0)
@named mass = Mass(; m=1, s_start, v_start)
@named spring = Spring(; c=10, s_rel0)

eqs = [
    connect(mass.flange_a, spring.flange_a)
    connect(spring.flange_b, fixed.flange)
]

@named model = ODESystem(eqs, t, [], []; systems=[fixed, mass, spring])

sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, ImplicitMidpoint(), dt=0.01)

Plotting sol[mass.s]...
image

If I remove the position parameters from Fixed and Spring, you can see we get the incorrect result...

@named fixed = Fixed()
@named mass = Mass(; m, s_start, v_start)
@named spring = Spring(; c)

image

Now, for comparison (using #107) we can see that only the Body position is needed, which to me makes more sense because we really only need to track absolute positions of actual moving rigid bodies. Additionally the Fixed component is defined by simply setting the port velocity to 0.

x0 = 3 # initial value of absolute position of sliding mass
v0 = 10 # initial value of absolute linear velocity of sliding mass

function Fixed(; name)
    @named port = Port()
    eqs = [port.v ~ 0]
    return compose(ODESystem(eqs, t, [], []; name = name), port)
end

@named fixed = Fixed()
@named mass = Body(; m=1, x0, v0)
@named spring = Spring(; k=10)

eqs = [
    connect(mass.port, spring.port_a)
    connect(spring.port_b, fixed.port)
]

@named model = ODESystem(eqs, t, [], []; systems=[fixed, mass, spring])

sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, ImplicitMidpoint(), dt=0.01)

image

This is also how other acausal modeling tools I know of define their thru and across (flow) variables for translational components.

@baggepinnen
Copy link
Contributor

Might I suggest changing this to velocity?

The intention is to follow modelica stdlib, if we're deviating from them, it's probably an oversight and should be fixed :)

@ValentinKaisermayer
Copy link
Contributor

ValentinKaisermayer commented Sep 15, 2022

Although, Simscape uses velocity and angular velocity for its mechanical domains.

In the example you show, I think you do not need to specify s_rel0 for the spring. It should be fully determined by the mass and fixed position.

@bradcarman
Copy link
Contributor

In this example the spring starts at it's natural length (i.e. uncompressed), therefore one should not need to specify s_rel0. By using position rather than velocity, this overcomplicates the model and now one must specify both the position of the Fixed component as well as s_rel0 for the spring.

Additionally, because the Flange is defined with s defaulting to 0, we now introduce a requirement that the model must solve for it's initial states. For a small model, this is not such a big deal, but for the large models I'm dealing with, initialization is something I'd like to avoid for 1) I need my model to start immediately, waiting for initialization is not an option for real time control, 2) I'd like to guarantee my states.

If we run this simple model without initialization and simplification, it can be seen that the result is not good.

image

Full code below...

using ModelingToolkitStandardLibrary.Mechanical.Translational, ModelingToolkit,
      OrdinaryDiffEq, Test
using Setfield

@parameters t
D = Differential(t)

s0 = 4.5 # fixed offset position of flange housing
s_start = 3 # initial value of absolute position of sliding mass
s_rel0 = 1.5 # unstretched spring length
v_start = 10 # initial value of absolute linear velocity of sliding mass

@named fixed = Fixed(; s0)
@named mass = Mass(; m=1, s_start, v_start)
@named spring = Spring(; c=10, s_rel0)

eqs = [
    connect(mass.flange_a, spring.flange_a)
    connect(spring.flange_b, fixed.flange)
]

@named model = ODESystem(eqs, t, [], []; systems=[fixed, mass, spring])

sys = expand_connections(model)
@set! sys.eqs = [(typeof(eq.lhs) != ModelingToolkit.Connection) & !ModelingToolkit._iszero(eq.lhs) & !ModelingToolkit.isdifferential(eq.lhs) ? 0 ~ eq.rhs - eq.lhs : eq for eq in sys.eqs]

prob = ODEProblem(sys, [], (0, 2.0))
sol = solve(prob, ImplicitMidpoint(), dt=0.01, initializealg=NoInit())

using CairoMakie
begin
    f,a,=lines(sol.t, sol[mass.s], label="sol[mass.s]");    
    lines!(a, sol.t, sol[mass.flange_a.s], label="sol[mass.flange_a.s]")
    lines!(a, sol.t, sol[spring.flange_a.s], label="sol[spring.flange_a.s]", linestyle=:dash)
    lines!(a, sol.t, sol[spring.flange_b.s], label="sol[spring.flange_b.s]")
    lines!(a, sol.t, sol[fixed.flange.s], label="sol[fixed.flange.s]", linestyle=:dash)
    a.xlabel="time [s]"; a.ylabel="position [m]"; 
    Legend(f[1,2], a)
    f
end

Additionally we can see that spring.s_rel is not automatically resolved

image

In contrast, by using velocity and force for the translational Connection, we only need to specify the position of the Body/Mass, and no initialization is required (using #107)...

using ModelingToolkitStandardLibrary.Mechanical.Translational, ModelingToolkit,
      OrdinaryDiffEq, Test
using Setfield

@parameters t
D = Differential(t)

x0 = 3 # initial value of absolute position of sliding mass
v0 = 10 # initial value of absolute linear velocity of sliding mass

function Fixed(; name)
    @named port = Port()
    eqs = [port.v ~ 0]
    return compose(ODESystem(eqs, t, [], []; name = name), port)
end

@named fixed = Fixed()
@named mass = Body(; m=1, x0, v0)
@named spring = Spring(; k=10)

eqs = [
    connect(mass.port, spring.port_a)
    connect(spring.port_b, fixed.port)
]

@named model = ODESystem(eqs, t, [], []; systems=[fixed, mass, spring])

sys = expand_connections(model)
@set! sys.eqs = [(typeof(eq.lhs) != ModelingToolkit.Connection) & !ModelingToolkit._iszero(eq.lhs) & !ModelingToolkit.isdifferential(eq.lhs) ? 0 ~ eq.rhs - eq.lhs : eq for eq in sys.eqs]
prob = ODEProblem(sys, [], (0, 2.0))
sol = solve(prob, ImplicitMidpoint(), dt=0.01, initializealg=NoInit())

using CairoMakie
begin
    f,a,=lines(sol.t, sol[mass.x], label="sol[mass.x]");    
    a.xlabel="time [s]"; a.ylabel="position [m]"; 
    Legend(f[1,2], a)
    f
end

image

@bradcarman
Copy link
Contributor

Making my point in another way, one should be able to define the Fixed component such that position is not needed, we simply specify that the change in position is zero...

function Fixed(; name)
    @named flange = Flange()
    eqs = [D(flange.s) ~ 0]
    return compose(ODESystem(eqs, t, [], []; name = name), flange)
end

But, unfortunately this will not work in the end because connect( fixed.flange, component.flange ) will give the equation...

fixed.flange.s ~ component.flange.s

Therefore, we are forced to specify the position for every single Flange in the model. But in theory, it should not be necessary to specify the position for a Fixed component, this is really overly defining the model.

@Firionus
Copy link

@bradcarman How would you express in your approach a spring that has non-zero length in its relaxed state?
That would certainly be a important for someone that wants to integrate with visualization or if someone wants to model something precisely the way they imagine their coordinate system.

@bradcarman
Copy link
Contributor

One final point (then I'll shut up, ha ha), the necessity that one must specify a spring length for the Spring component is silly when one is considering models that really just deal with stiffness. For example, defining a Hard Stop component that defines the stiffness of the boundary. Should the length/thickness of the boundary really need to be defined? Mathematically this information is irrelevant. My recommend spring definition then is...

function Spring(;name, k = 1e3, delta0 = 0.0)
    pars = @parameters k=k delta0=delta0
    vars = @variables begin
        Δx(t) = delta0
        f(t) = k*delta0
    end
    
    @named port_a = Port()
    @named port_b = Port()
    
    eqs = [
        D(Δx) ~ port_a.v - port_b.v
        f ~ k*Δx
        port_a.f ~ +f
        port_b.f ~ -f
    ]
    return compose(ODESystem(eqs, t, vars, pars; name = name), port_a, port_b)
end

Then, we only need to define 2 parameters, the stiffness, and the initial compression state (delta0).

And to answer @Firionus, for visualization, if we are plotting the position states, the only states tracked by the model are actual moving bodies/masses. So if we have 2 bodies connected by a spring, then we can plot both body positions.

@Firionus
Copy link

And to answer @Firionus, for visualization, if we are plotting the position states, the only states tracked by the model are actual moving bodies/masses. So if we have 2 bodies connected by a spring, then we can plot both body positions.

This doesn't hold true if the spring is connected to a fixed element. Also, visualizing springs with a finite relaxed length is much closer to a real mechanical setup than springs that are infinitely small sometimes, don't you think?

One final point (then I'll shut up, ha ha), the necessity that one must specify a spring length for the Spring component is silly when one is considering models that really just deal with stiffness.

If you find it silly you can just leave it at the default of s_rel0 = 0.0 and achieve that. I realize you'll still have to define the position of fixed points but I think this is a small price to pay for gaining better modelling of real life with the absolute position approach.

Ultimately, it comes down to the choice of programming interface. Should a "mechanical component" only keep track of velocities and forces, i.e. its energy exchange, or should it also track its absolute position. The Modelica Standard Library (MSL 3.2.3) goes with the second approach:

connector Flange "One-dimensional translational flange"

  SI.Position s "Absolute position of flange";
  flow SI.Force f "Cut force directed into flange";
  annotation (...);
end Flange;

And since @baggepinnen said:

The intention is to follow modelica stdlib

This PR gets that right, since it follows the MSL:

@connector function Flange(; name)
    sts = @variables begin
        s(t)
        f(t), [connect = Flow]
    end
    ODESystem(Equation[], t, sts, [], name = name, defaults = Dict(s => 0.0, f => 0.0))
end

source

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants