Constrained oscillator

In this example, we model a mechanical oscillator with two springs connected in series. This imposes the constraint that the forces in both springs must be equal. Consequently, the overall displacement and the potential energy is distributed between the springs to satisfy this constraint.

We define a reversible component representing the constraint that the velocities (efforts) of two springs are equal:

sc = ReversibleComponent(
  Dtry(
    :q => Dtry(ReversiblePort(FlowPort(displacement, CVar(:λ)))),
    :q₂ => Dtry(ReversiblePort(FlowPort(displacement, -CVar(:λ)))),
    :λ => Dtry(ReversiblePort(Constraint(-EVar(:q) + EVar(:q₂))))
  )
)

ReversibleComponent
■
├─ q => Quantity(:displacement, :ℝ, true)
│       f = λ
├─ q₂ => Quantity(:displacement, :ℝ, true)
│        f = -1.0 * λ
└─ λ => constraint: 0 = -1.0 * q.e + q₂.e

\[\begin{bmatrix} \mathtt{q.f} \\ \mathtt{q_2.f} \\ 0 \end{bmatrix} \: = \: \begin{bmatrix} 0 & 0 & +1 \\ 0 & 0 & -1 \\ -1 & +1 & 0 \end{bmatrix} \, \begin{bmatrix} \mathtt{q.e} \\ \mathtt{q_2.e} \\ \lambda \end{bmatrix}\]

The constraint variable $\lambda$ is distributing the overall (change of) displacement between the two springs.

Next, we define the constrained oscillator. To not reinvent the wheel, we use hookean_spring, point_mass, and pkc from the ComponentLibrary.

osc_constraint = CompositeSystem(
  Dtry(
    :q => Dtry(Junction(displacement, Position(1, 2))),
    :p => Dtry(Junction(momentum, Position(1, 4))),
    :q₂ => Dtry(Junction(displacement, Position(3, 2))),
  ),
  Dtry(
    :pe₁ => Dtry(
      InnerBox(
        Dtry(
          :q => Dtry(InnerPort(■.q)),
        ),
        hookean_spring(2.0),
        Position(1, 1)
      ),
    ),
    :pe₂ => Dtry(
      InnerBox(
        Dtry(
          :q => Dtry(InnerPort(■.q₂)),
        ),
        hookean_spring(6.0),
        Position(3, 1)
      ),
    ),
    :sc => Dtry(
      InnerBox(
        Dtry(
          :q => Dtry(InnerPort(■.q)),
          :q₂ => Dtry(InnerPort(■.q₂)),
        ),
        sc,
        Position(2, 2)
      ),
    ),
    :pkc => Dtry(
      InnerBox(
        Dtry(
          :q => Dtry(InnerPort(■.q)),
          :p => Dtry(InnerPort(■.p))
        ),
        pkc,
        Position(1, 3)
      ),
    ),
    :ke => Dtry(
      InnerBox(
        Dtry(
          :p => Dtry(InnerPort(■.p)),
        ),
        point_mass(1.0),
        Position(1, 5)
      ),
    ),
  )
)

This yields the following system of equations:

assemble(osc_constraint)

Flows:
ke.p.f = -1.0 * pe₁.q.e
pe₁.q.f = ke.p.e + -1.0 * sc.λ
pe₂.q.f = sc.λ
Efforts:
ke.p.e = ke.p.x * (ke.m^-1.0)
pe₁.q.e = pe₁.k * pe₁.q.x
pe₂.q.e = pe₂.k * pe₂.q.x
Constraints:
0 = -1.0 * pe₁.q.e + pe₂.q.e  enforced by sc.λ

Simulation

We define an initial condition, simulate the system using the midpoint rule, and plot the resulting displacements:

ic_constraint = Dtry(
  :ke => Dtry(
    :p => Dtry(3.0),
  ),
  :pe₁ => Dtry(
    :q => Dtry(0.0),
  ),
  :pe₂ => Dtry(
    :q => Dtry(0.0),
  )
)
h = 0.1
t = 10
sim = simulate(osc_constraint, midpoint_rule, ic_constraint, h, t)

plot_evolution(
  sim,
  XVar(DtryPath(:pe₁), DtryPath(:q)),
  XVar(DtryPath(:pe₂), DtryPath(:q)),
)

Compared to the first spring, the second spring is stiffer and hence contributes less to the overall displacement.

Comparison with unconstrained case

To verify the numerical solution of the differential-algebraic equation, we compare the system with the unconstrained oscillator. The following code assumes that osc and ic are defined as in the previous example.

comparison = CompositeSystem(
  Dtry(
    :p => Dtry(Junction(momentum, Position(1, 2))),
  ),
  Dtry(
    :osc => Dtry(
      InnerBox(
        Dtry(
          :p => Dtry(InnerPort(■.p)),
        ),
        osc,
        Position(1, 1)
      ),
    ),
    :osc_constraint => Dtry(
      InnerBox(
        Dtry{InnerPort}(),
        osc_constraint,
        Position(1, 4)
      ),
    ),
  )
)

We expect that the overall displacement agrees with the unconstrained case, since the stiffness coefficients satisfy the following relation:

\[\frac{1}{\mathtt{osc\_constraint.pe_1.k}} \, + \, \frac{1}{\mathtt{osc\_constraint.pe_2.k}} \: = \: \frac{1}{\mathtt{osc.pe.k}}\]

ic_comparison = Dtry(
  :osc => ic,
  :osc_constraint => ic_constraint,
)

sim = simulate(comparison, midpoint_rule, ic_comparison, h, t)

maximum(abs.(
  evolution(
    sim,
    XVar(DtryPath(:osc_constraint, :pe₁), DtryPath(:q)) +
    XVar(DtryPath(:osc_constraint, :pe₂), DtryPath(:q)) -
    XVar(DtryPath(:osc, :pe), DtryPath(:q))
  )
))

1.7763568394002505e-15

Indeed, the maximum error is close to the limits of machine precision.