# 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