This page describes the discretisation of the dynamics for second-order (Lagrangian) systems of the form
The feature of the proposed scheme is to combine an exact method for the linear (non-contacting) part of the equations of motion with a Moreau-Jean time-stepping
approach to handle impulses and velocity jumps. This is the implementation of what is proposed in [8] for the simulation of musical string instruments.
Some simulation examples are proposed in examples/Mechanics/Music directory of Siconos repository.
Detailed calculations
For comparison purpose, we consider two different schemes, the classical (in Siconos) ‘Moreau-Jean’ and an implementation of what is proposed in [3], denoted respectively MJ and BMJ in the following.
Let us integrate the dynamics over a time step, \(\Delta t = t^{i+1} - t^i\).
Notations :
\[\begin{split}q_k(t^i) = q_k^i \\
v_k(t^{i}) = v_k^i \\\end{split}\]
In the following, we will use k for space (bottom) indices and i for time (top) indices.
MJ is based on a theta-scheme, for \(\theta \in [0,1]\)
MJ we consider diagonal stiffness and damping,
\begin{eqnarray}
M = diag(\mu) \\
K = \mu.diag(\omega_k^2) \\
C = \mu.diag(2\sigma_k)
\end{eqnarray}
\(\omega_k\) and \(\sigma_k\) being respectively the modal pulsation and the damping parameter (see for instance values taken from [8]).
Bilbao exact scheme writes:
\[\begin{split}\begin{array}{ccc}
Kq &\approx \Gamma Kq^i + \frac{(\mathcal{I}-\Gamma)K}{2}(q^{i+1} + q^{i-1}) \\
C\dot q &\approx \frac{1}{\Delta t}\Sigma^*(q^{i+1} - q^{i-1}) \\
\end{array}\end{split}\]
for \(\Gamma = diag(\gamma_k)\) and \(\Sigma^* = diag(\sigma_k^*)\) some diagonal matrices, with
\[\begin{split}\gamma_{k} &= \frac{2}{\omega_k^2\Delta t^2} - \frac{A_k}{1+e_k-A_k}, \\
\sigma^*_{k} &= \left(\frac{1}{\Delta t} + \frac{\omega_k^2\Delta t}{2} - \gamma_k\frac{\omega_k^2\Delta t}{2} \right)\frac{1-e_k}{1+e_k} \\
A_k &= e^{-\sigma_k\Delta t}\left(e^{\sqrt{\sigma_k^2 - \omega_k^2}\Delta t} + e^{-\sqrt{\sigma_k^2 - \omega_k^2}\Delta t}\right) \\
e_k &= e^{-2\sigma_k\Delta t} \\\end{split}\]
\[\begin{split}\begin{array}{c|c|c}
Dynamics & Moreau-Jean & Modal-Moreau-Jean \\
\int_{t^i}^{t^{i+1}} Mdv & \approx M(v^{i+1}-v^{i}) & \approx M(v^{i+1}-v^{i}) \\
\int_{t^i}^{t^{i+1}} Kqdt & \approx \Delta t(\theta Kq^{i+1} + (1 - \theta) Kq^i) & \approx \Delta t\Gamma Kq^i + \frac{\Delta t}{2}(\mathcal{I}-\Gamma)K(q^{i+1} + q^{i-1}) \\
\int_{t^i}^{t^{i+1}} Cvdt & \approx \Delta t(\theta Cv^{i+1} + (1 - \theta) Cv^i) & \approx \Sigma^*(q^{i+1} - q^{i-1})\\
\int_{t^i}^{t^{i+1}} dr & \approx p^{i+1} & \approx p^{i+1} \\
\end{array}\end{split}\]
For MJ, this leads to
\[\begin{split}M(v^{i+1}-v^{i}) + \Delta t(\theta Kq^{i+1} + (1 - \theta) Kq^i) + \Delta t(\theta Cv^{i+1} + (1 - \theta) Cv^i) &= p^{i+1} \\\end{split}\]
using \(q^{i+1} = q^i + \Delta t(\theta v^{i+1} + (1 - \theta) v^i)\), we get
\[\begin{split}[M + \Delta t^2\theta^2 K + \Delta t\theta C] (v^{i+1}-v^{i}) + \Delta tKq^i + (\Delta t^2\theta K + \Delta tC) v^i = p^{i+1} \\\end{split}\]
And for BMJ:
\[M(v^{i+1}-v^{i}) + \Delta t\Gamma Kq^i + \frac{\Delta t}{2}(\mathcal{I}-\Gamma)K(q^{i+1} + q^{i-1}) +\Sigma^*(q^{i+1} - q^{i-1}) = p^{i+1}\]
With \(q^{i+1} = q^{i} + \Delta tv^{i+1}\), we get
\[\begin{split}q^{i+1} - q^{i-1} &= \Delta t(v^{i+1} + v^i) \\
q^{i+1} + q^{i-1} &= 2q^i + \Delta t(v^{i+1} - v^i) \\\end{split}\]
and
\[\begin{split}[M + \frac{\Delta t^2}{2}(\mathcal{I} - \Gamma)K + \Delta t\Sigma^*] (v^{i+1}-v^{i}) + \Delta tKq^i + 2\Delta t \Sigma^* v^i = p^{i+1} \\\end{split}\]
Both discretisations writes
\[\begin{split}W(v^{i+1}-v^{i}) = v_{free}^i + p^{i+1} \\
or \\
v^{i+1} = v^i_{free} + W^{-1}p^{i+1} \\\end{split}\]
with
\[\begin{split}\begin{array}{c|c|c}
& Moreau-Jean & Modal-Moreau-Jean \\
W & = M + \Delta t^2\theta^2 K + \Delta t\theta C & = M + \frac{\Delta t^2}{2}(\mathcal{I} - \Gamma)K + \Delta t\Sigma^*\\
v_{free}^{i} &= v^i - W^{-1}(\Delta tKq^i + (\Delta t^2\theta K + \Delta tC) v^i) & = v^i - W^{-1}(\Delta tKq^i + 2\Delta t \Sigma^* v^i) \\
\end{array}\end{split}\]
Taylor expansions
For some values of the time step it may be necessary to use Taylor expansions of iteration matrix and \(\Delta t\sigma^*\) to avoid convergence problems.
Those terms write:
\[\begin{split}\Delta t\sigma^* & = \Delta t \sigma + \frac{\Delta t^{3} \sigma}{12} \omega^{2} + \Delta t^{5} \left(\frac{\omega^{4} \sigma}{240} - \frac{\omega^{2} \sigma^{3}}{180}\right) \\
& + \Delta t^{7} \left(\frac{\omega^{6} \sigma}{6048} - \frac{\omega^{4} \sigma^{3}}{1512} + \frac{\omega^{2} \sigma^{5}}{1890}\right) + \mathcal{O}\left(\Delta t^{8}\right) \\
\frac{1}{W_{kk}} &= 1 - \Delta t \sigma + \Delta t^{2} \left(- \frac{\omega^{2}}{12} + \frac{2 \sigma^{2}}{3}\right) \\
&+\Delta t^{3} \left(\frac{\omega^{2} \sigma}{12} - \frac{\sigma^{3}}{3}\right) + \Delta t^{4} \left(\frac{\omega^{4}}{360} - \frac{\omega^{2} \sigma^{2}}{20} + \frac{2 \sigma^{4}}{15}\right)\\
&+ \Delta t^{5} \left(- \frac{\omega^{4} \sigma}{360} + \frac{\omega^{2} \sigma^{3}}{45} - \frac{2 \sigma^{5}}{45}\right) \\
&+ \Delta t^{6} \left(- \frac{\omega^{6}}{20160} + \frac{\omega^{4} \sigma^{2}}{630} - \frac{\omega^{2} \sigma^{4}}{126} + \frac{4 \sigma^{6}}{315}\right)\\
&+ \Delta t^{7} \left(\frac{\omega^{6} \sigma}{20160} - \frac{\omega^{4} \sigma^{3}}{1512} + \frac{\omega^{2} \sigma^{5}}{420} - \frac{\sigma^{7}}{315}\right) + \mathcal{O}\left(\Delta t^{8}\right)\\\end{split}\]