newmark_ss
Produces a discrete-time state-space model of the structural equations
\[\begin{split}\mathbf{\ddot{x}} &= \mathbf{M}^{-1}( -\mathbf{C}\,\mathbf{\dot{x}}-\mathbf{K}\,\mathbf{x}+\mathbf{F} ) \\
\mathbf{y} &= \mathbf{x}\end{split}\]
based on the Newmark-\(\beta\) integration scheme. The output state-space model
has form:
\[\begin{split}\mathbf{X}_{n+1} &= \mathbf{A}\,\mathbf{X}_n + \mathbf{B}\,\mathbf{F}_n \\
\mathbf{Y} &= \mathbf{C}\,\mathbf{X} + \mathbf{D}\,\mathbf{F}\end{split}\]
with \(\mathbf{X} = [\mathbf{x}, \mathbf{\dot{x}}]^T\)
Note that as the state-space representation only requires the input force
\(\mathbf{F}\) to be evaluated at time-step \(n\),the \(\mathbf{C}\) and \(\mathbf{D}\) matrices
are, in general, fully populated.
The Newmark-\(\beta\) integration scheme is carried out following the modifications presented by
Geradin [1] that render it unconditionally stable. The displacement and velocities are estimated as:
\[\begin{split}x_{n+1} &= x_n + \Delta t \dot{x}_n + \left(\frac{1}{2}-\theta_2\right)\Delta t^2 \ddot{x}_n + \theta_2\Delta t
\ddot{x}_{n+1} \\
\dot{x}_{n+1} &= \dot{x}_n + (1-\theta_1)\Delta t \ddot{x}_n + \theta_1\Delta t \ddot{x}_{n+1}\end{split}\]
The stencil is unconditionally stable if the tuning parameters \(\theta_1\) and \(\theta_2\) are chosen as:
\[\begin{split}\theta_1 &= \frac{1}{2} + \alpha \\
\theta_2 &= \frac{1}{4} \left(\theta_1 + \frac{1}{2}\right)^2 \\
\theta_2 &= \frac{5}{80} + \frac{1}{4} (\theta_1 + \theta_1^2) \text{TBC SOURCE}\end{split}\]
where \(\alpha>0\) accounts for small positive algorithmic damping.
The following steps describe how to apply the Newmark-beta scheme to a state-space formulation. The original idea
is based on [1].
The equation of a second order system dynamics reads:
\[M\mathbf{\ddot q} + C\mathbf{\dot q} + K\mathbf{q} = F\]
Applying that equation to the time steps \(n\) and \(n+1\), rearranging terms and multiplying by
\(M^{-1}\):
\[\begin{split}\mathbf{\ddot q}_{n} = - M^{-1}C\mathbf{\dot q}_{n} - M^{-1}K\mathbf{q}_{n} + M^{-1}F_{n} \\
\mathbf{\ddot q}_{n+1} = - M^{-1}C\mathbf{\dot q}_{n+1} - M^{-1}K\mathbf{q}_{n+1} + M^{-1}F_{n+1}\end{split}\]
The relations of the Newmark-beta scheme are:
\[\begin{split}\mathbf{q}_{n+1} &= \mathbf{q}_n + \mathbf{\dot q}_n\Delta t +
(\frac{1}{2}-\beta)\mathbf{\ddot q}_n \Delta t^2 + \beta \mathbf{\ddot q}_{n+1} \Delta t^2 + O(\Delta t^3) \\
\mathbf{\dot q}_{n+1} &= \mathbf{\dot q}_n + (1-\gamma)\mathbf{\ddot q}_n \Delta t +
\gamma \mathbf{\ddot q}_{n+1} \Delta t + O(\Delta t^3)\end{split}\]
Substituting the former relation onto the later ones, rearranging terms, and writing it in state-space form:
\[\begin{split}\begin{bmatrix} I + M^{-1}K \Delta t^2\beta \quad \Delta t^2\beta M^{-1}C \\ (\gamma \Delta t M^{-1}K)
\quad (I + \gamma \Delta t M^{-1}C) \end{bmatrix} \begin{Bmatrix} \mathbf{\dot q}_{n+1} \\
\mathbf{\ddot q}_{n+1} \end{Bmatrix} =
\begin{bmatrix} (I - \Delta t^2(1/2-\beta)M^{-1}K \quad (\Delta t - \Delta t^2(1/2-\beta)M^{-1}C \\
(-(1-\gamma)\Delta t M^{-1}K \quad (I - (1-\gamma)\Delta tM^{-1}C \end{bmatrix}
\begin{Bmatrix} \mathbf{q}_{n} \\ \mathbf{\dot q}_{n} \end{Bmatrix} +
\begin{Bmatrix} (\Delta t^2(1/2-\beta) \\ (1-\gamma)\Delta t \end{Bmatrix} M^{-1}F_n+
\begin{Bmatrix} (\Delta t^2\beta) \\ (\gamma \Delta t) \end{Bmatrix}M^{-1}F_{n+1}\end{split}\]
To understand SHARPy code, it is convenient to apply the following change of notation:
\[\begin{split}\textrm{th1} = \gamma \\
\textrm{th2} = \beta \\
\textrm{a0} = \Delta t^2 (1/2 -\beta) \\
\textrm{b0} = \Delta t (1 -\gamma) \\
\textrm{a1} = \Delta t^2 \beta \\
\textrm{b1} = \Delta t \gamma \\\end{split}\]
Finally:
\[\begin{split}A_{ss1} \begin{Bmatrix} \mathbf{\dot q}_{n+1} \\ \mathbf{\ddot q}_{n+1} \end{Bmatrix} =
A_{ss0} \begin{Bmatrix} \mathbf{\dot q}_{n} \\ \mathbf{\ddot q}_{n} \end{Bmatrix} +
\begin{Bmatrix} (\Delta t^2(1/2-\beta) \\ (1-\gamma)\Delta t \end{Bmatrix} M^{-1}F_n+
\begin{Bmatrix} (\Delta t^2\beta) \\ (\gamma \Delta t) \end{Bmatrix}M^{-1}F_{n+1}\end{split}\]
To finally isolate the vector at \(n+1\), instead of inverting the \(A_{ss1}\) matrix, several systems are
solved. Moreover, the output equation is simply \(y=x\).
- param Minv:
Inverse mass matrix \(\mathbf{M^{-1}}\)
- type Minv:
np.array
- param C:
Damping matrix \(\mathbf{C}\)
- type C:
np.array
- param K:
Stiffness matrix \(\mathbf{K}\)
- type K:
np.array
- param dt:
Timestep increment
- type dt:
float
- param num_damp:
Numerical damping. Default 1e-4
- type num_damp:
float
- returns:
the A, B, C, D matrices of the state space packed in a tuple with the predictor and delay term removed.
- rtype:
tuple
References
[1] - Geradin M., Rixen D. - Mechanical Vibrations: Theory and application to structural dynamics