MSD with 2-DoF

The source files for all examples can be found in /examples.

A more challenging yet simple dynamic system than the harmonic oscillator is the two degrees of freedom mass-spring-damper system, as illustrated below:

Mass-Spring-Damper System

The following differential equation governs the motion of the harmonic oscillator:

TODO: update the equation $m\,\frac{d^2 x}{dt^2} + c\,\frac{d^2 x}{dt^2} + k\,x = F(t)$,

where $x$ is the position, $m$ is the mass, $c$ is the damping coefficient, $k$ is the spring stiffness, and $F(t)$ is an external force applied to the system.

To model the 2-DoF system using the bond graph method, we need to connect the inertance (I), compliance (C), and resistance (R) elements to different 1-junctions and 0-junctions, as depicted in the bond graph diagram:

Bond Graph Representation

In this representation, the inertance element (I) corresponds to the mass (m) in the mechanical system, the compliance element (C) corresponds to the spring constant (k), and the resistance element (R) represents the damping coefficient (c). They were considered, for simplicity, equal for both degrees of freedom.

Bond-graph toolkit

Firstly, we need to import the bond-graph toolkit module to model the system and the DifferentialEquations.jl to solve the resultant ODE. Also, we import the independent variable t from the bond-graph toolkit for defining a custom forcing term.

using BondGraph
using BondGraph: t, D
using DifferentialEquations
using ModelingToolkit
using Plots
using Symbolics.Latexify

Building the model

Setting the elements

Before defining a system, we need to define the single port elements of the system. Then, we define the mass, spring, and damper elements as follows:

@named m = Mass(m=1.0);      # (kg)
@named s = Spring(k=1.0);    # (N/m)
@named d = Damper(c=1.0);    # (N⋅m⋅s)

┌ Warning: power contains 0 flow variables, yet 2 regular (non-flow, non-stream, non-input, non-output) variables. This could lead to imbalanced model that are difficult to debug. Consider marking some of the regular variables as input/output variables.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/connectors.jl:40
┌ Warning: power contains 0 flow variables, yet 2 regular (non-flow, non-stream, non-input, non-output) variables. This could lead to imbalanced model that are difficult to debug. Consider marking some of the regular variables as input/output variables.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/connectors.jl:40
┌ Warning: power contains 0 flow variables, yet 2 regular (non-flow, non-stream, non-input, non-output) variables. This could lead to imbalanced model that are difficult to debug. Consider marking some of the regular variables as input/output variables.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/GFEA2/src/systems/connectors.jl:40

The forcing term can be set as an expression or even a Julia function. In this example, we will first intialize a generic forcing term that will be defined before solving the ODE.

@parameters F(t)
@named f = Se(F)

\[ \begin{align} power_{+}e\left( t \right) =& F\left( t \right) \end{align} \]

Setting the junctions and subsystems

Differently from the 1-DoF, where we could define the entire system with only a 1-junction. For the 2-DoF system we are required to define multiple 1-junctions and one 0-junction. Thus, we define them separately and later we will connect them accordingly.

Firstly we define the ground connected mass ($m_1$).

@named m1_j1 = Junction1([-1, m], [-1, d], [-1, s])

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( power, m_{+}power \right) \\ \mathrm{connect}\left( power, d_{+}power \right) \\ \mathrm{connect}\left( power, s_{+}power \right) \\ \frac{\mathrm{d} m_{+}power_{+}f\left( t \right)}{\mathrm{d}t} = \frac{m_{+}power_{+}e\left( t \right)}{m_{+}I} \\ d_{+}power_{+}e\left( t \right) = d_{+}R d_{+}power_{+}f\left( t \right) \\ s_{+}power_{+}e\left( t \right) = \frac{s_{+}q\left( t \right)}{s_{+}C} \\ \frac{\mathrm{d} s_{+}q\left( t \right)}{\mathrm{d}t} = s_{+}power_{+}f\left( t \right) \\ \end{array} \right] \end{equation} \]

Then, we define the $2^{nd}$ DoF mass (($m_2$)) with two 1-junctions and one 0-junction.

@named m2_j1 = Junction1([-1, m], f)
@named m2_j0 = Junction0()
@named m2_sd = Junction1([-1, s], [-1, d])

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( power, s_{+}power \right) \\ \mathrm{connect}\left( power, d_{+}power \right) \\ s_{+}power_{+}e\left( t \right) = \frac{s_{+}q\left( t \right)}{s_{+}C} \\ \frac{\mathrm{d} s_{+}q\left( t \right)}{\mathrm{d}t} = s_{+}power_{+}f\left( t \right) \\ d_{+}power_{+}e\left( t \right) = d_{+}R d_{+}power_{+}f\left( t \right) \\ \end{array} \right] \end{equation} \]

Connect the subsystems

Define the connections between the junctions.

cons = [
    connect(m2_j0.power, m2_j1.power),
    connect(m2_j0.power, m2_sd.power),
    connect(m1_j1.power, m2_j0.power),
]

Build the system

We build the system by first creating a ODESystem with junctions connections and setting the independent variable as t.

@named mdl = ODESystem(cons, t)

Then, we need to add the junctions subsystems to the model

mdl = compose(mdl, m1_j1, m2_j0, m2_j1, m2_sd)

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( m2_{j0_{+}power}, m2_{j1_{+}power} \right) \\ \mathrm{connect}\left( m2_{j0_{+}power}, m2_{sd_{+}power} \right) \\ \mathrm{connect}\left( m1_{j1_{+}power}, m2_{j0_{+}power} \right) \\ \mathrm{connect}\left( power, m_{+}power \right) \\ \mathrm{connect}\left( power, d_{+}power \right) \\ \mathrm{connect}\left( power, s_{+}power \right) \\ \frac{\mathrm{d} \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} = \frac{\mathrm{m1}_{j1_{+}m_{+}power_{+}e}\left( t \right)}{m1_{j1_{+}m_{+}I}} \\ \mathrm{m1}_{j1_{+}d_{+}power_{+}e}\left( t \right) = m1_{j1_{+}d_{+}R} \mathrm{m1}_{j1_{+}d_{+}power_{+}f}\left( t \right) \\ \mathrm{m1}_{j1_{+}s_{+}power_{+}e}\left( t \right) = \frac{\mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{m1_{j1_{+}s_{+}C}} \\ \frac{\mathrm{d} \mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} = \mathrm{m1}_{j1_{+}s_{+}power_{+}f}\left( t \right) \\ \mathrm{connect}\left( power, m_{+}power \right) \\ \mathrm{connect}\left( f_{+}power, power \right) \\ \frac{\mathrm{d} \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} = \frac{\mathrm{m2}_{j1_{+}m_{+}power_{+}e}\left( t \right)}{m2_{j1_{+}m_{+}I}} \\ \mathrm{m2}_{j1_{+}f_{+}power_{+}e}\left( t \right) = \mathrm{m2}_{j1_{+}f_{+}F}\left( t \right) \\ \mathrm{connect}\left( power, s_{+}power \right) \\ \mathrm{connect}\left( power, d_{+}power \right) \\ \mathrm{m2}_{sd_{+}s_{+}power_{+}e}\left( t \right) = \frac{\mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{m2_{sd_{+}s_{+}C}} \\ \frac{\mathrm{d} \mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} = \mathrm{m2}_{sd_{+}s_{+}power_{+}f}\left( t \right) \\ \mathrm{m2}_{sd_{+}d_{+}power_{+}e}\left( t \right) = m2_{sd_{+}d_{+}R} \mathrm{m2}_{sd_{+}d_{+}power_{+}f}\left( t \right) \\ \end{array} \right] \end{equation} \]

Analysing the model

We can visualize the bond-graph connections by generating a graph plot of the model.

generate_graph(mdl)

We can obtain the bond-graph equations using the code below.

equations(expand_connections(mdl))

\[ \begin{align} \frac{\mathrm{d} \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{\mathrm{m1}_{j1_{+}m_{+}power_{+}e}\left( t \right)}{m1_{j1_{+}m_{+}I}} \\ \mathrm{m1}_{j1_{+}d_{+}power_{+}e}\left( t \right) =& m1_{j1_{+}d_{+}R} \mathrm{m1}_{j1_{+}d_{+}power_{+}f}\left( t \right) \\ \mathrm{m1}_{j1_{+}s_{+}power_{+}e}\left( t \right) =& \frac{\mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{m1_{j1_{+}s_{+}C}} \\ \frac{\mathrm{d} \mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} =& \mathrm{m1}_{j1_{+}s_{+}power_{+}f}\left( t \right) \\ \frac{\mathrm{d} \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{\mathrm{m2}_{j1_{+}m_{+}power_{+}e}\left( t \right)}{m2_{j1_{+}m_{+}I}} \\ \mathrm{m2}_{j1_{+}f_{+}power_{+}e}\left( t \right) =& \mathrm{m2}_{j1_{+}f_{+}F}\left( t \right) \\ \mathrm{m2}_{sd_{+}s_{+}power_{+}e}\left( t \right) =& \frac{\mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{m2_{sd_{+}s_{+}C}} \\ \frac{\mathrm{d} \mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} =& \mathrm{m2}_{sd_{+}s_{+}power_{+}f}\left( t \right) \\ \mathrm{m2}_{sd_{+}d_{+}power_{+}e}\left( t \right) =& m2_{sd_{+}d_{+}R} \mathrm{m2}_{sd_{+}d_{+}power_{+}f}\left( t \right) \\ \mathrm{m2}_{j0_{+}power_{+}f1}\left( t \right) =& - \mathrm{m2}_{sd_{+}d_{+}power_{+}f}\left( t \right) \\ \mathrm{m2}_{j0_{+}power_{+}f1}\left( t \right) =& - \mathrm{m2}_{sd_{+}s_{+}power_{+}f}\left( t \right) \\ - \mathrm{m2}_{j0_{+}power_{+}e1}\left( t \right) =& - \mathrm{m2}_{j0_{+}power_{+}e2}\left( t \right) \\ - \mathrm{m2}_{j0_{+}power_{+}e1}\left( t \right) =& \mathrm{m1}_{j1_{+}power_{+}e}\left( t \right) \\ - \mathrm{m1}_{j1_{+}power_{+}f}\left( t \right) =& - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) \\ - \mathrm{m1}_{j1_{+}power_{+}f}\left( t \right) =& - \mathrm{m1}_{j1_{+}s_{+}power_{+}f}\left( t \right) \\ - \mathrm{m1}_{j1_{+}power_{+}f}\left( t \right) =& - \mathrm{m1}_{j1_{+}d_{+}power_{+}f}\left( t \right) \\ 0 =& - \mathrm{m2}_{j0_{+}power_{+}f1}\left( t \right) - \mathrm{m2}_{j0_{+}power_{+}f2}\left( t \right) + \mathrm{m1}_{j1_{+}power_{+}f}\left( t \right) \\ 0 =& - \mathrm{m2}_{j1_{+}m_{+}power_{+}e}\left( t \right) + \mathrm{m2}_{j0_{+}power_{+}e1}\left( t \right) + \mathrm{m2}_{j1_{+}f_{+}power_{+}e}\left( t \right) \\ 0 =& - \mathrm{m1}_{j1_{+}d_{+}power_{+}e}\left( t \right) - \mathrm{m1}_{j1_{+}m_{+}power_{+}e}\left( t \right) - \mathrm{m1}_{j1_{+}power_{+}e}\left( t \right) - \mathrm{m1}_{j1_{+}s_{+}power_{+}e}\left( t \right) \\ 0 =& - \mathrm{m2}_{sd_{+}d_{+}power_{+}e}\left( t \right) - \mathrm{m2}_{sd_{+}s_{+}power_{+}e}\left( t \right) + \mathrm{m2}_{j0_{+}power_{+}e2}\left( t \right) \\ \mathrm{m2}_{j0_{+}power_{+}f2}\left( t \right) =& \mathrm{m2}_{j1_{+}f_{+}power_{+}f}\left( t \right) \\ \mathrm{m2}_{j0_{+}power_{+}f2}\left( t \right) =& - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \end{align} \]

The equations of the model above were not simplified and represent the bond-graphs equations directly derived from the connections. Then, we simplify the DAE above and obtain the ODE. The simplified set of equations is required for using solve from DifferentialEquations.jl

@named sys = simplifysys(mdl)

\[ \begin{align} \frac{\mathrm{d} \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{m2_{sd_{+}d_{+}R} \left( - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \right) + \frac{ - \mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{m1_{j1_{+}s_{+}C}} + \frac{\mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{m2_{sd_{+}s_{+}C}} - m1_{j1_{+}d_{+}R} \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{m1_{j1_{+}m_{+}I}} \\ \frac{\mathrm{d} \mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} =& \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) \\ \frac{\mathrm{d} \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{m2_{sd_{+}d_{+}R} \left( - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \right) + \frac{\mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{m2_{sd_{+}s_{+}C}} + \mathrm{m2}_{j1_{+}f_{+}F}\left( t \right)}{m2_{j1_{+}m_{+}I}} \\ \frac{\mathrm{d} \mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} =& - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \end{align} \]

Print system states

states(sys)

4-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 m1_j1₊m₊power₊f(t)
 m1_j1₊s₊q(t)
 m2_j1₊m₊power₊f(t)
 m2_sd₊s₊q(t)

Print system parameters

parameters(sys)

7-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 m1_j1₊d₊R
 m1_j1₊s₊C
 m2_sd₊d₊R
 m2_sd₊s₊C
 m1_j1₊m₊I
 m2_j1₊f₊F(t)
 m2_j1₊m₊I

Print the system simplified equations

equations(sys)

By comparison, we can see that the equations above are the same as the canonical form of the harmonic oscillator presented at the beginning.

We can generate LaTeX code from the equations

latexify(equations(sys))

\begin{align} \frac{\mathrm{d} \mathrm{m1}{j1{+}m{+}power{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{m2{sd{+}d{+}R} \left( - \mathrm{m1}{j1{+}m{+}power{+}f}\left( t \right) - \mathrm{m2}{j1{+}m{+}power{+}f}\left( t \right) \right) + \frac{ - \mathrm{m1}{j1{+}s{+}q}\left( t \right)}{m1{j1{+}s{+}C}} + \frac{\mathrm{m2}{sd{+}s{+}q}\left( t \right)}{m2{sd{+}s{+}C}} - m1{j1{+}d{+}R} \mathrm{m1}{j1{+}m{+}power{+}f}\left( t \right)}{m1{j1{+}m{+}I}} \ \frac{\mathrm{d} \mathrm{m1}{j1{+}s{+}q}\left( t \right)}{\mathrm{d}t} =& \mathrm{m1}{j1{+}m{+}power{+}f}\left( t \right) \ \frac{\mathrm{d} \mathrm{m2}{j1{+}m{+}power{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{m2{sd{+}d{+}R} \left( - \mathrm{m1}{j1{+}m{+}power{+}f}\left( t \right) - \mathrm{m2}{j1{+}m{+}power{+}f}\left( t \right) \right) + \frac{\mathrm{m2}{sd{+}s{+}q}\left( t \right)}{m2{sd{+}s{+}C}} + \mathrm{m2}{j1{+}f{+}F}\left( t \right)}{m2{j1{+}m{+}I}} \ \frac{\mathrm{d} \mathrm{m2}{sd{+}s{+}q}\left( t \right)}{\mathrm{d}t} =& - \mathrm{m1}{j1{+}m{+}power{+}f}\left( t \right) - \mathrm{m2}{j1{+}m{+}power{+}f}\left( t \right) \end{align}

Simulate the system

Unforced case

Now, we are going to consider only the unforced case. Therefore, $F(t) = 0$. Thus, we update the force function on the system.

sys_unforced = substitute(sys, Dict(m2_j1.f.F => 0.0))

\[ \begin{align} \frac{\mathrm{d} \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{m2_{sd_{+}d_{+}R} \left( - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \right) + \frac{ - \mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{m1_{j1_{+}s_{+}C}} + \frac{\mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{m2_{sd_{+}s_{+}C}} - m1_{j1_{+}d_{+}R} \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{m1_{j1_{+}m_{+}I}} \\ \frac{\mathrm{d} \mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} =& \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) \\ \frac{\mathrm{d} \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{m2_{sd_{+}d_{+}R} \left( - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \right) + \frac{\mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{m2_{sd_{+}s_{+}C}}}{m2_{j1_{+}m_{+}I}} \\ \frac{\mathrm{d} \mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} =& - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \end{align} \]

Define the simulation time

tspan = (0, 10)

(0, 10)

When defining the ODEProblem as in the ModelingToolkit, we can set the initial value of the states.

prob = ODEProblem(sys_unforced, [m2_j1.m.power.f => 0.0, m2_sd.s.q => 1.0], tspan)

ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 10)
u0: 4-element Vector{Float64}:
 0.0
 0.0
 0.0
 1.0

Use the solve to simulate the system.

sol = solve(prob)
plot(sol, xlabel="Time", ylabel="Amplitude")

Forced case

Now, we are going to consider the forced case, where $F(t) = sin(t)$. Thus, we update the force function on the system.

sys_forced = substitute(sys, Dict(m2_j1.f.F => sin(t)))

\[ \begin{align} \frac{\mathrm{d} \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{m2_{sd_{+}d_{+}R} \left( - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \right) + \frac{ - \mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{m1_{j1_{+}s_{+}C}} + \frac{\mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{m2_{sd_{+}s_{+}C}} - m1_{j1_{+}d_{+}R} \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{m1_{j1_{+}m_{+}I}} \\ \frac{\mathrm{d} \mathrm{m1}_{j1_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} =& \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) \\ \frac{\mathrm{d} \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right)}{\mathrm{d}t} =& \frac{m2_{sd_{+}d_{+}R} \left( - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \right) + \frac{\mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{m2_{sd_{+}s_{+}C}} + \sin\left( t \right)}{m2_{j1_{+}m_{+}I}} \\ \frac{\mathrm{d} \mathrm{m2}_{sd_{+}s_{+}q}\left( t \right)}{\mathrm{d}t} =& - \mathrm{m1}_{j1_{+}m_{+}power_{+}f}\left( t \right) - \mathrm{m2}_{j1_{+}m_{+}power_{+}f}\left( t \right) \end{align} \]

Define the simulation time

tspan = (0, 10)

(0, 10)

When defining the ODEProblem as in the ModelingToolkit we can define the initial value of the states.

prob = ODEProblem(sys_forced, [m2_j1.m.power.f => 0.0, m2_sd.s.q => 0.0], tspan)

ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 10)
u0: 4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0

Use the solve to simulate the system.

sol = solve(prob)
plot(sol, xlabel="Time", ylabel="Amplitude")

This page was generated using Literate.jl.