Time integration of the dynamics - Exact scheme

Time integration for second-order (Lagrangian) systems of the form:

\[\begin{split}M\ddot q(t) + Kq(t) + C\dot q(t) = p \\ q(t_0) = q_0 \\ \dot q(t_0) = v_0\end{split}\]

equivalent to

(1)\[\begin{split}Mdv + Kqdt + cv^+dt = dr \\ v(t) = \dot q^+(t) \\ q(t_0) = q_0, \ \dot q(t_0) = v_0\end{split}\]

with

(2)\[q(t) = q_0 + \int_{t_0}^t v^+(t)\]

\(v^+(t)\) stands for right limit of v in t.

with the vectors \(q = [q_k], v = [v_k] k\in[0,ndof-1]\).

Let us integrate the dynamics over a time step, \(\Delta t = t^{i+1} - t^i\).

We consider two different schemes, the classical (in Siconos) ‘Moreau-Jean’ and ‘Modal-Moreau-Jean’ denoted respectively MJ and MMJ in the following.

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]\)

\[\begin{split}\int_{t^i}^{t^{i+1}} f(q,v,t)dt & \approx \Delta t\theta f(q^{i+1}, v^{i+1}, t^{i+1}) + \Delta t(1 -\theta)f(q^{i}, v^{i}, t^{i}) \\ & \approx \Delta t\theta f^{i+1} + \Delta t(1 -\theta)f^{i} \\\end{split}\]

With MJM we consider diagonal stiffness and damping,

\begin{eqnarray} K = diag(\omega_k^2) \\ C = diag(2\sigma_k) \end{eqnarray}

\(\omega_k\) and \(\sigma_k\) being respectively the modal pulsation and the damping parameter (values to be taken from 2.2 and 2.3 in JSV paper).

Bilbao exact scheme writes:

\[\begin{split}\begin{array}{ccc} Kq &\approx \Theta Kq^i + \frac{(\mathcal{I}-\Theta)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 \(\Theta = diag(\theta_k)\) and \(\Sigma^* = diag(\sigma_k^*)\) some diagonal matrices, with

\[\begin{split}\theta_{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} - \theta_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\Theta Kq^i + \frac{\Delta t}{2}(\mathcal{I}-\Theta)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 MMJ:

\[M(v^{i+1}-v^{i}) + \Delta t\Theta Kq^i + \frac{\Delta t}{2}(\mathcal{I}-\Theta)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} - \Theta)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}) = \tilde 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} - \Theta)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}\]

Notes, remarks, questions

  • Quel nom pour “modal” Moreau-Jean? i.e. qui est à la source (ref?) du schéma de Bilbao?
  • Vérif comportement de W quand \(\Delta t \rightarrow 0\)

Non-smooth problem formulation

\[\begin{split}\dot y &= S_c \dot q \\ P &= S_c^T\lambda\end{split}\]
\[\begin{split}\dot y_c^{i+1} &= S_c \dot q^{i+1} \\ P^{i+1} &= S_c^T\lambda^{i+1}\end{split}\]
\[\begin{split}\dot y^{i+1} &= S_cv^{i} - S_cW^{-1}(\Delta tKq^i + 2\Delta t \Sigma^* v^i) + S_cW^{-1}S_c^T\lambda^{i+1} \\ &= q_{LCP} + M_{LCP}\lambda^{i+1}\end{split}\]

with

\[0 \leq \dot y^{i+1} \perp \lambda^{i+1} \geq 0\]