Constrained oscillator

In this example, we model a mechanical oscillator with two springs connected in series:

The series connection 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.

To not reinvent the wheel, this time we use hookean_spring, point_mass, and pkc from the ComponentLibrary. Besides these components, which are defined as in the previous example, we use two_springs_series_connection, which represents the constraint that the velocities (efforts) of two springs are equal:

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₂)),
        ),
        two_springs_series_connection,
        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.