Time-integration of a mechanical system

29.6. Time-integration of a mechanical system#

Reference to:

29.6.1. Example equation: torsion of elastic beam#

Here, the PDE governing the torsional dynamics of a linear homogeneous elastic beam is used as a toy problem

\[I \partial_{tt} \theta - \partial_x \left( GJ \partial_{x} \theta \right) = m \]

supplied with proper boundary and initial conditions. Here the clamped-free beam is studied, with prescribed torsional moment \(M(t)\) ath teh “free” end so that the boundary conditions read

\[\begin{split}\begin{aligned} \theta(x=0, t) & = 0 \\ GJ \partial_x \theta(x=\ell, t) & = M(t) \ . \end{aligned}\end{split}\]

See

29.6.1.1. Exact solution#

29.6.1.2. Finite element model#

Semi-discretized model in space

\[\begin{split}\begin{aligned} \int_{z=0}^{\ell} w(z) m(z,t) \, dz & = \int_{z=0}^{\ell} w(z) \left\{ I \partial_{tt} \theta - \partial_x \left(GJ \partial_x \theta \right) \right\} \, dz = \\ & = \int_{z=0}^{\ell} w(z) I \partial_{tt} \theta(z,t) \, dz + \int_{z=0}^{\ell} \partial_x w(z) GJ \partial_x \theta(z,t) \, dz - \left.\left[ w(x) GJ \partial_x \theta(x,t) \right]\right|_{z=0}^{\ell} \ . \end{aligned}\end{split}\]

With the boundary conditions of the model, the FEM approximation of the problem becomes

\[\begin{aligned} \int_{z=0}^{\ell} w(z) m(z,t) \, dz + w(\ell) M(t) & = \int_{z=0}^{\ell} w(z) I \partial_{tt} \theta(z,t) \, dz + \int_{z=0}^{\ell} \partial_x w(z) GJ \partial_x \theta(z,t) \, dz \ . \end{aligned}\]

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

from scipy.sparse import csr_matrix

def assemble_fem_matrices(L, N, I_coeff, GJ_coeff):
    """
    Assembles Sparse Mass and Stiffness matrices for P1 Lagrange elements
    based on 'Screenshot from 2026-05-14 18-59-39.png'.
    """
    # Grid setup
    nodes = np.linspace(0, L, N + 1)
    h = L / N  # element width
    num_nodes = N + 1

    # Local Element Matrices for P1 Lagrange
    # Mass Matrix Local: (h/6) * [[2, 1], [1, 2]]
    m_local = (I_coeff * h / 6.0) * np.array([[2.0, 1.0], 
                                              [1.0, 2.0]])
    
    # Stiffness Matrix Local: (1/h) * [[1, -1], [-1, 1]]
    k_local = (GJ_coeff / h) * np.array([[1.0, -1.0], 
                                         [-1.0, 1.0]])

    # Global Assembly using COO format for efficiency
    rows = []
    cols = []
    mass_data = []
    stiff_data = []

    for i in range(N):
        # Indices for the current element (nodes i and i+1)
        idx = [i, i + 1]
        
        for r in range(2):
            for c in range(2):
                rows.append(idx[r])
                cols.append(idx[c])
                mass_data.append(m_local[r, c])
                stiff_data.append(k_local[r, c])

    # Convert to CSR (Compressed Sparse Row) format
    M_sparse = csr_matrix((mass_data, (rows, cols)), shape=(num_nodes, num_nodes))
    K_sparse = csr_matrix((stiff_data, (rows, cols)), shape=(num_nodes, num_nodes))

    return M_sparse, K_sparse

def get_partitioned_matrices(L, N, I_coeff, GJ_coeff, dirichlet_indices):
    # 1. Reuse the assembly logic to get global matrices
    # (Assuming the assembly function from the previous step is defined)
    M_global, K_global = assemble_fem_matrices(L, N, I_coeff, GJ_coeff)
    
    num_nodes = N + 1
    
    # 2. Define Node Sets
    # Node 0 is Dirichlet (index 0), all others are Free
    all_indices = np.arange(num_nodes)
    free_indices = np.setdiff1d(all_indices, dir_indices)
    
    # 3. Partitioning: Extract the FF block (Free-Free)
    # This effectively removes the first row and first column
    M_ff = M_global[free_indices, :][:, free_indices]
    K_ff = K_global[free_indices, :][:, free_indices]
    
    return M_ff, K_ff, free_indices


# Example parameters
L = 1.0           # Length of the domain
N = 10            # Number of elements
I_val = 1.0       # Inertia coefficient (from image)
GJ_val = 1.0      # Torsional stiffness (from image)
dir_indices = [0] # Dirichlet indices of the nodes with prescribed displacement
Mff, Kff, free_indices = get_partitioned_matrices(L, N, I_val, GJ_val, dir_indices)

# print(f"Mass Matrix Shape: {Mff.shape}")
# print(f"Mass      Matrix (First 3x3 block):\n{Mff.toarray()[:3, :3]}")
# print(f"Stiffness Matrix (First 3x3 block):\n{Kff.toarray()[:3, :3]}")
# print(f"Sum(Mass Matrix): {np.sum(Mff.toarray())}")

29.6.2. Boundary conditions#

Matrix partition or augmented system.

29.6.2.1. Matrix partition#

\[\begin{split}\begin{aligned} \mathbf{f} & = \mathbf{M} \ddot{\mathbf{x}} + \mathbf{K} \mathbf{x} \\ \end{aligned}\end{split}\]

\[\begin{split}\begin{aligned} \begin{bmatrix} \mathbf{f}_f \\ \mathbf{f}_d \end{bmatrix} & = \begin{bmatrix} \mathbf{M}_{ff} & \mathbf{M}_{fd} \\ \mathbf{M}_{df} & \mathbf{M}_{dd} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{x}}_{f} \\ \ddot{\mathbf{x}}_d \end{bmatrix} + \begin{bmatrix} \mathbf{K}_{ff} & \mathbf{K}_{fd} \\ \mathbf{K}_{df} & \mathbf{K}_{dd} \end{bmatrix} \begin{bmatrix} \mathbf{x}_{f} \\ \mathbf{x}_d \end{bmatrix} \end{aligned}\end{split}\]

Here the value of the kinematic variables \(\mathbf{x}_d\), and the actions \(\mathbf{f}_f\) is known. The unknowns are \(\mathbf{x}_f\) (the displacement of the free d.o.f.s of the system), and \(\mathbf{f}_d\) (the reactions at the constraints). Separating the unknowns from the known quantities,

\[\begin{split}\begin{aligned} \mathbf{M}_{ff} \ddot{\mathbf{x}}_f + \mathbf{K}_{ff} \mathbf{x}_f & = \mathbf{f}_f - \mathbf{M}_{fd} \ddot{\mathbf{x}}_d - \mathbf{K}_{fd} \mathbf{x}_d \\ \mathbf{f}_d & = - \begin{bmatrix} \mathbf{M}_{df} & \mathbf{M}_{dd} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{x}}_f \\ \ddot{\mathbf{x}}_d \end{bmatrix} - \begin{bmatrix} \mathbf{K}_{df} & \mathbf{K}_{dd} \end{bmatrix} \begin{bmatrix} \mathbf{x}_f \\ \mathbf{x}_d \end{bmatrix} \end{aligned}\end{split}\]

This system is decoupled, as the reactions at the constraints are not involved in the dynamical systems for \(\mathbf{x}_f\), and thus they can be computed a posteriori.

29.6.2.2. Augmented system#

See the discussion about Poisson equation - i.e. the static version of the dynamical equation used here, with no inertia contribution - in Math: Numerical methods for PDEs

29.6.2.2.1. Spectral decomposition of the undamped system#

Eigenvalue problem focuses on the free response of the system, so all the forcing terms - both prescribed actions in \(\mathbf{f}_f\) and prescribed kinetic variables in \(\mathbf{x}_d\) - are set equal to zero. Thus, it’s possible to focus on the \(f\)-part of the partitioned system only, \(\mathbf{M}_{ff} \ddot{\mathbf{x}}_f + \mathbf{K}_{ff} \mathbf{x}_f = \dots\). In the following lines, pedices are dropped.

Either generalized eigenvalue problem of the undamped system as a second-order \(n\)-dimensional system

\[\mathbf{M} \ddot{\mathbf{x}} + \mathbf{K} \mathbf{x} = \mathbf{f} \ ,\]

\[s^2 \mathbf{M} \mathbf{x} + \mathbf{K} \mathbf{x} = \mathbf{0}\]

or the generalized eigenvalue problem as a first-order \(2n\)-dimensional system

\[\begin{split}\begin{aligned} \begin{bmatrix} \mathbf{I} & \cdot \\ \cdot & \mathbf{M} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{x}} \\ \dot{\mathbf{v}} \end{bmatrix} + \begin{bmatrix} \cdot & - \mathbf{I} \\ \mathbf{K} & \cdot \end{bmatrix} \begin{bmatrix} \mathbf{x} \\ \mathbf{v} \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{f} \end{bmatrix} \ . \end{aligned}\end{split}\]

\[s \mathbf{A} \mathbf{z} + \mathbf{B} \mathbf{z} = \mathbf{0}\]

Matrix symmetries and routines.

In the second-orded system, matrices of a mechanical system are usually symmetric. It’s possible to use the routine for Hermitian matrices \(\texttt{.eigsh()}\)
In the first-orded system, the augmented matrices are not symmetric. It’s not possible to use the routine for Hermitian matrices. Use the routine for general matrix \(\texttt{.eigs()}\) instead

Solution of the eigenvalue problem, and diagonalization - modal form. Routines returns mass normalized eigenvectors, so that

\[\mathbf{V}^* \mathbf{M} \mathbf{V} = \mathbf{I} \ .\]

As the matrix form of the generalized eigenvalue problem reads \(\mathbf{M} \mathbf{V} \mathbf{S}^2 + \mathbf{K} \mathbf{V} = \mathbf{0}\).

Changing basis from “nodal” to “modal” variables, \(\mathbf{x} = \mathbf{V} \mathbf{q}\), and projecting on the eigenvectors, the dynamical equation becomes

\[\begin{split}\begin{aligned} \mathbf{V}^* \mathbf{f} & = \mathbf{V}^* \left( \mathbf{M} \mathbf{V} \ddot{\mathbf{u}} + \mathbf{K} \mathbf{V} \mathbf{u} \right) \\ \mathbf{V}^* \mathbf{f} & = \ddot{\mathbf{u}} + \text{diag} \left\{ - s_i^2 \right\} \mathbf{u} = \\ & = \ddot{\mathbf{u}} + \text{diag} \left\{ \Omega_i^2 \right\} \mathbf{u} \ . \end{aligned}\end{split}\]

For a mechanical system with definite positive \(\mathbf{M}\) and \(\mathbf{K}\), the eigenvalues of the system are imaginary so that \(s_i^2 = - \Omega^2_i\).

Show code cell source Hide code cell source

import scipy as sp

#> Solve generalized sparse eigenvalue problem with symmetric matrices, s^2 M x + K x = 0
# evals, evecs = eigsh(A, k=6, M=None, which='LM', ...)
# for solving the A v = s M v problem
# k=number of evals, evecs computed. k < n (not possible to evaluate all the evals with this routine,
# but usually sparse matrix is used for large matrices, and one is not interested in the whole spectrum.
# For small matrices, either convert to full matrix or use this routine twice, one looking for SM and one for LM
# which='LM': largest magnitude, 'SM': smallest magnitude
evals_, evecs = sp.sparse.linalg.eigsh(Kff, k=N-1, M=Mff, which='SM')

#> Here evals^2 = - evals_
# evals are imaginary, \mp j \omega
evals = np.sqrt( -np.real(evals_) - 1j * np.imag(evals_) )

# print(f"First 3 Omega (rad/s): {np.abs(np.imag(evals[:3]))}")
# print("First 3 eigenvectors:")
# print(evecs[:,:3])

#> Run the routine twice to get eigenvalues with largest magnitude, and assemble them with the eigenvalues with smallest magnitude
#> If the full decomposition is required, just used the same routing twice: evaluate SM(k=kSM), and then
# evaluate LM(k=kLM), with kSM+kLM = N
#> As the routine is runned twice without sharing any info, apply normalization (if it's meaningful) later
# If kSM=N-1, then kLM=1 (without coincident eigvalues with largest magnitude), the eigenvectors could be already
# normal. By default the routine returns a mass-normalized eigenvectors
evals_lm_, evecs_lm = sp.sparse.linalg.eigsh(Kff, k=1, M=Mff, which='LM')
evals_lm = np.sqrt( -np.real(evals_lm_) - 1j * np.imag(evals_lm_) )

evals_all = np.block([evals, evals_lm])
evecs_all = np.block([evecs, evecs_lm])

# print(evals_all)
# print(evecs_all)

if ( np.max( evecs_all.T @ Mff @ evecs_all - np.eye(N) ) < 1e-9 ):
    print("Mass-normalization of the eigenvectors within the tol\n")
else:
    print("Mass-normalization of the eigenvectors not satisfied, or above the tol\n")

# print(evals_)

Mass-normalization of the eigenvectors within the tol

29.6.3. Response to a torsional moment at the “free” end#

With force at the free end, the forcing term of the discrete system reads

\[\begin{split}\mathbf{f}_f = \begin{bmatrix} 0 \\ \dots \\ 0 \\ M(t) \end{bmatrix} = \begin{bmatrix} 0 \\ \dots \\ 0 \\ 1 \end{bmatrix} M(t) \ ,\end{split}\]

being the last node the free node. In modal variables

\[\ddot{u}_i + 2 \xi_i \Omega_i \dot{u}_i + \Omega_i^2 u_i = \mathbf{v}^*_i \mathbf{f}_f \ .\]

The input-output relation from the input \(M(t)\) and the \(N-1\)-dimensional output \(u_i\) reads

\[u_i(s) = \frac{\mathbf{v}_i^* \mathbf{F}}{s^2 + 2 \xi_i \Omega_i s + \Omega_i^2} M(t) \ ,\]

i.e. the transfer function between the torsional moment at the free end \(M(t)\) and the amplitude of the \(i^{th}\) mode is

\[G_i(s) = \frac{\mathbf{v}_i^* \mathbf{F}}{s^2 + 2 \xi_i \Omega_i s + \Omega_i^2} \ .\]

../../_images/6ac26704f6ff1d214e9dcdc5f74deb6420534e19fe11070dada5e16ccf888288.png

29.6.4. Time simulations#

Here different approaches are evaluated:

Maximum frequency in forcing spectrum : 3.0
Threshold for slow/fast mode splitting: 12.0
Number of slow modes                  : 4

../../_images/92001cbdaf57d63b0678f481f7b521dfe2febcf4d63ae3607a5008f375e12ab2.png

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 6
      2 u0 = np.zeros(N)
      3 v0 = np.zeros(N)
      4 #> Mode acceleration
      5 u0_m_acc = u0 \
----> 6         - sp.sparse.linalg.spsolve(Kff, F) * M_t(tvec[0]) \
      7         + evecs_all[:,:i_slow] @ np.diag( Om_all[:i_slow]**(-2) ) @ evecs_all[:,:i_slow].T @ F * M_t(tvec[0])
      8 
      9 #> Initial conditions for the reduced (resolved, dynamic, slow) system

NameError: name 'tvec' is not defined

Mode acceleration and static recovery of fast modes.

\[\begin{split}\begin{aligned} \begin{bmatrix} \mathbf{u} \\ \mathbf{K} \mathbf{u} \\ \dot{\mathbf{q}}_s \\ \ddot{\mathbf{q}}_s \end{bmatrix} = \begin{bmatrix} \mathbf{U}_s & \mathbf{0} \\ \mathbf{K} \mathbf{U}_s & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \\ -\text{diag}\{ \omega^2_s \} & -\text{diag}\{ 2 \xi_s \omega_s \} \end{bmatrix} \begin{bmatrix} \mathbf{q}_s \\ \dot{\mathbf{q}}_s \end{bmatrix} + \begin{bmatrix} \mathbf{K}^{-1} - \mathbf{U}_s \text{diag}\{ (m_s \omega_s^2)^{-1} \} \mathbf{U}^*_s \\ \mathbf{I} - \mathbf{M} \mathbf{U}_s \text{diag}\{ m_s^{-1} \} \mathbf{U}^*_s \\ \mathbf{0} \\ \text{diag}\{ m_s^{-1} \} \mathbf{U}_s^T \end{bmatrix} \mathbf{f} \ . \end{aligned}\end{split}\]

Time integration of the slow part

First order system

\[\begin{split}\begin{aligned} \frac{d}{dt} \begin{bmatrix} \mathbf{q}_s \\ \dot{\mathbf{q}}_s \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\text{diag}\{ \omega^2_s \} & -\text{diag}\{ 2 \xi_s \omega_s \} \end{bmatrix} \begin{bmatrix} \mathbf{q}_s \\ \dot{\mathbf{q}}_s \end{bmatrix} + \begin{bmatrix} \mathbf{0} \\ \text{diag}\{ m_s^{-1} \} \mathbf{U}_s^T \end{bmatrix} \mathbf{f} \ . \end{aligned}\end{split}\]

\[\dot{\mathbf{z}} = \mathbf{A} \mathbf{z} + \mathbf{B} M(t) = \mathbf{f}(t, \mathbf{z})\]

Second order system

\[\ddot{\mathbf{q}}_s = - \text{diag}\{ \omega_s^2 \} \mathbf{q}_s - \text{diag}\{ 2 \xi \omega_s \} \dot{\mathbf{q}}_s + \text{diag}\{ m_s^{-1} \} \mathbf{U}_s^* \mathbf{f}\]

29.6.4.1. Explicit Euler#

Explicit Euler.

\[\mathbf{z}_{n+1} = \mathbf{z}_n + h \mathbf{f}(t_n, \mathbf{z}_n)\]

Explicit Euler numerical integration is unstable

…

../../_images/c60ccb28baa3d381575b7054829e2f42e50a92a4bae39297987637db16c6f754.png

29.6.4.2. RK4(5)#

Runge-Kutta. Exploiting scipy.integrate.solve_ivp. RK4(5) is a variable time-step method.

RK4(5) and high-frequencies. This method is not suited for systems with high-frequency contents, as the method uses smaller and smaller time-step to keep local truncation error small (and the method stable). Extreme small values of \(dt\) makes the number of iterations explodes and the performance of the method so poor that it becomes useless.

RK4(5) with fast mode truncation and mode acceleration recovery. Using the fast mode truncation and their static recovery, the dynamical part of the system involves only slow modes. Thus, this approach makes RK4(5) method suitable.

../../_images/8d2a171880c02ede2fd8df39d8035fbb4528e571fc9b3f9c0a3a7645e318f1a0.png

../../_images/7d4d686aac0ad39ffbb0e5b324f8aae276eba13e3f6250dbdebc5a2053c9ab64.png