Hyperelasticity: static and dynamics

A short introduction to the structural models implemented in the FEM_Library

$\newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}} \newcommand{\bvec}[1]{\mathbf{#1}} \newcommand{\bx}{\mathbf{x}} \newcommand{\bd}{\mathbf{d}} \newcommand{\bv}{\dot{\mathbf{d}}} \newcommand{\ba}{\ddot{\mathbf{d}}} \newcommand{\bmw}{\mathbf{w}} \newcommand{\bmv}{\mathbf{v}} \newcommand{\Nh}{N_h} \newcommand{\bU}{\mathbf{U}} \newcommand{\bp}{\bvec{p}} \newcommand{\VV}{\mathbf{V}} \newcommand{\WW}{\mathbf{W}} \newcommand{\QQ}{\mathbf{Q}} \renewcommand{\AA}{\mathbf{A}} \newcommand{\BB}{\mathbf{B}} \newcommand{\CC}{\mathbf{B}} \newcommand{\EE}{\mathbf{E}} \newcommand{\FF}{\mathbf{F}} \newcommand{\FFi}{\mathbf{f}^i} \newcommand{\FFe}{\mathbf{f}^e} \newcommand{\GG}{\mathbf{G}} \newcommand{\CC}{\mathbf{C}} \newcommand{\MM}{\mathbf{M}} \newcommand{\XX}{\mathbf{X}} \newcommand{\KK}{\mathbf{K}} \newcommand{\UU}{\mathbf{U}} \newcommand{\II}{\mathbf{I}} \renewcommand{\SS}{\mathbf{S}} \newcommand{\HH}{\mathbf{H}} \newcommand{\DD}{\mathbf{D}} \newcommand{\NN}{\mathbf{N}} \newcommand{\bmg}{\mathbf{g}} \newcommand{\R}{\mathbb{R}} \newcommand{\JJ}{\mathbf{J}} \newcommand{\RR}{\mathbf{R}} \newcommand{\PP}{\mathbf{P}}$

Basics
- Equation of motion
- Equilibrium equation
Constitutive laws
Finite element approximation
- Static
- Dynamics
  - Newmark scheme
  - Generalized-$\alpha$ method
References

Basics

Let $\Omega \subset \mathbb{R}^{d}$ ($d=2,3$) be a bounded domain, referred to as reference configuration of a body $\mathcal{B}$. The body can be identified as a set of material points. If $P$ is a generic material point, $\mathbf{X}(P)$ indicates the position vector of $P$ in the reference configuration with respect to some reference system.

When the body undergoes deformation, every material point $P$ is being displaced to a new deformed location which is denoted by $\mathbf{x}$. The relation between each material point and its respective deformed location is captured by the deformation function $\bm{\phi}: \mathbb{R}^d \rightarrow \mathbb{R}^3$,

$\mathbf{x} = \bm{\phi}(\mathbf{X}).$

The displacement vector $\bd(t)$ of a material point is defined as

$\bd = \mathbf{x} - \mathbf{X}.$

One of the most important kinematic quantities in the framework of continuum mechanics is the deformation gradient

$\FF = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \nabla \bd + \II.$

We denote by $J = \text{det}(\FF)$ the determinant of the deformation gradient.

Equation of motion

The equation of motion for a continuous medium formulated in the reference configuration reads:

$\begin{equation} \rho \frac{\partial^2 \bd}{\partial t^2} - \text{div}\left( \PP \right) = \mathbf{f} \qquad \text{in } \Omega, \end{equation} \label{eq:CSM_Motion}$

where $\mathbf{f}$ is an external force, $\rho$ is the density of the body and $\PP = \PP(\FF)$ is the first Piola-Kirchhoff stress tensor. Note that the first Piola-Kirchhoff stress tensor is related to the Cauchy stress tensor $\bm{\sigma}$ by the following formula: $\PP = J \bm{\sigma} \FF^{-T}$.

The relation between the stress and the deformation tensor is the so-called constitutive law, peculiar to each specific material. Equation \eqref{eq:CSM_Motion} is referred to as governing equation of nonlinear elastodynamics.

Equilibrium equation

In equilibrium, equation \eqref{eq:CSM_Motion} is replaced by

$\begin{equation} - \text{div}\left( \PP \right) = \mathbf{f} \qquad \text{in } \Omega, \end{equation} \label{eq:CSM_Static}$

and this, complemented by a suitable constitutive law $\PP = \PP(\FF)$, is the governing equation of nonlinear elastostatics.

Constitutive laws

The computation of the displacement field of the considered body as function of the external loads requires to complement equations \eqref{eq:CSM_Motion} and \eqref{eq:CSM_Static} by a suitable constitutive law. In this section we present all the constitutive laws for hyperelastic materials implemented in the redbKIT library.

The characteristic of an hyperelastic (also called Green elastic) material is the existence of a strain energy function $\psi = \psi(\FF)$ such that

$\PP(\FF) = \frac{ \partial \psi(\FF) }{ \partial \FF }.$

Linear elasticity

The simplest practical constitutive model is linear elasticity, defined in terms of the strain energy density as

$\psi(\FF) = \mu \bm{\epsilon} : \bm{\epsilon} + \frac{\lambda}{2} (\text{tr}(\bm{\epsilon}))^2,$

where $\bm{\epsilon} = \tfrac{1}{2}(\FF + \FF^T - \II) = \tfrac{1}{2}(\nabla\bd + \nabla \bd^T)$ is the small (or infinitesimal) strain tensor. The coefficients $\mu$ and $\lambda$ are the Lame constants, which are related to the material properties through the Young modulus $E$ and the Poisson ratio $\nu$:

$\mu = \frac{E}{2(1+E)}, \qquad \lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}.$

Its first Piola-Kirchhoff stress tensor is given by

$\PP(\FF) = 2\mu \bm{\epsilon} + \lambda \text{tr}(\bm{\epsilon}) \II.$

St. Venant-Kirchhoff model

The St. Venant-Kirchhoff material is a an hyperelastic nonlinear model for compressible materials characterised by the following strain energy function:

$\psi(\FF) = \mu \EE : \EE + \frac{\lambda}{2} (\text{tr}(\EE))^2,$

where $\EE = \tfrac{1}{2}(\FF^T\FF - \II) = \tfrac{1}{2} (\CC - \II)$ is the Green-Lagrange strain tensor, while $\lambda$ and $\mu$ are the Lame constants. Its first Piola-Kirchhoff stress tensor is given by

$\PP(\FF) = \FF \left( 2\mu \EE + \lambda \text{tr}(\EE) \II \right).$

Nearly-incompressible Neo-Hookean model

The nearly-incompressible Neo-Hookean material is a an hyperelastic nonlinear model whose strain energy function is specified in terms of an isochoric-volumetric splitting:

$\psi(\FF) = \underbrace{\frac{\mu}{2}\left( \bar{I}_{1} - 3 \right)}_{\psi_{\text{iso}}} + \underbrace{\frac{K}{4}\left( (J-1)^2 + (\ln J)^2\right)}_{\psi_{\text{vol}}}$

where $\bar{I}_{1} = J^{-2/3} I_1 = J^{-2/3} \text{tr}(\CC)$, $J = \text{det}(\FF)$, $\mu$ is the shear modulus and $K$ is the bulk modulus.

The corresponding first Piola-Kirchhoff stress tensor is given by

$\PP(\FF) = \mu J^{-2/3} \left( \FF - \frac{1}{3} I_1 \FF^{-T} \right) + \frac{K}{2} J \left( J - 1 + \frac{1}{J} \ln J\right) \FF^{-T}.$

Nearly-incompressible Raghavan-Vorp model

Like the Neo-Hookean model, the nearly-incompressible Raghavan-Vorp material is a an hyperelastic nonlinear model whose strain energy function is specified in terms of an isochoric-volumetric splitting:

$\psi(\FF) = \underbrace{\alpha \left( \bar{I}_{1} - 3 \right) + \beta \left( \bar{I}_{1} - 3 \right)^2 }_{\psi_{\text{iso}}} + \underbrace{\frac{K}{4}\left( (J-1)^2 + (\ln J)^2\right)}_{\psi_{\text{vol}}}.$

This is a special case of the generalized power law neo-Hookean model; also note that this constitutive model reduces to the Neo-Hookean one for $\alpha = \mu/2$ and $\beta = 0$.

The corresponding first Piola-Kirchhoff stress tensor is given by

$\PP(\FF) = 2\left(\alpha + 2\beta \left( \bar{I}_{1} - 3 \right) \right) J^{-2/3} \left( \FF - \frac{1}{3} I_1 \FF^{-T} \right) + \frac{K}{2} J \left( J - 1 + \frac{1}{J} \ln J\right) \FF^{-T}.$

Finite element approximation

We now consider the finite element approximation of the both nonlinear elastostatic and elastodynamics equations. To this end, it is first useful to complement these equations with suitable boundary and initial conditions, and to introduce their weak formulation.

Static

The boundary-value problem reads as follows: find $\bd : \overline{\Omega} \rightarrow \mathbb{R}^d$ such that

$\begin{equation} \left\{ \begin{aligned} - \text{div}\left( \PP(\bd) \right) &= \mathbf{f} &\qquad &\text{in } \Omega \\ \bd &= \mathbf{g}& \qquad &\text{on } \Gamma_D \\ \PP(\bd) \mathbf{n} &= \mathbf{h} & \qquad &\text{on } \Gamma_N \\ \PP(\bd) \mathbf{n} + \alpha \bd &= \mathbf{0} & \qquad &\text{on } \Gamma_R, \end{aligned} \right. \label{eq:CSM_Static2} \end{equation}$

where we have prescribed suitable boundary displacements $\mathbf{g}$ on $\Gamma_D$ and tractions $\mathbf{h}$ on $\Gamma_N$. On $\Gamma_R$, a Robin-type boundary condition is imposed instead; $\alpha$ is called elastic coefficient.

Let us introduce the following functional spaces

$V = \left \{ \bmw \in [H^1(\Omega)]^d \, : \, \bmw|_{\Gamma_D} = \mathbf{g} \right \}$ $V_D = \left \{ \bmw \in [H^1(\Omega)]^d \, : \, \bmw|_{\Gamma_D} = \mathbf{0} \right \}$

The weak formulation of \eqref{eq:CSM_Static2} reads: find $\bd \in V_D$ such that

$\begin{equation} \int_{\Omega} \PP(\bd) : \nabla \bmw \, d\Omega + \int_{\Gamma_R} \alpha \bd \cdot \bmw \, d\Gamma = \int_{\Omega} \mathbf{f} \cdot \bmw \, d\Omega + \int_{\Gamma_N} \mathbf{h} \cdot \bmw \, d\Gamma \qquad \forall \bmw \in V. \label{eq:CSM_StaticWeak} \end{equation}$

Let us introduce a FE partition $\mathcal{T}_h$ of the domain $\Omega$ from which we construct a conforming finite element space

$X_h = \text{span} \{ \bm{\varphi}_i, \,\, i=1, \dots \Nh \} \subset [H^1(\Omega)]^d.$

We also define $V_h = X_h \cap V$ and $V_{Dh} = X_h \cap V_D$. The finite element approximation of \eqref{eq:CSM_StaticWeak} reads: find $\bd_h \in V_{Dh}$ such that

$\begin{equation} \int_{\Omega} \PP(\bd_h) : \nabla \bmw \, d\Omega + \int_{\Gamma_R} \alpha \bd_h \cdot \bmw \, d\Gamma = \int_{\Omega} \mathbf{f} \cdot \bmw \, d\Omega + \int_{\Gamma_N} \mathbf{h} \cdot \bmw \, d\Gamma \qquad \forall \bmw \in V_h. \label{eq:CSM_StaticFE} \end{equation}$

With a slight abuse of notation, we denote by $\bd \in \mathbb{R}^{\Nh}$ the vector of coefficients in the expansion of $\bd_h$ with respect to the FE basis functions $\bm{\varphi}_i$. At the algebraic level, we obtain the following nonlinear system of equations to be solved:

$\RR( \bd ) = \FFi (\bd) - \FFe = \mathbf{0},$

where $\FFi \in \mathbb{R}^{\Nh}$ denotes the vector of internal forces

$\FFi_i (\bd) = \int_{\Omega} \PP(\bd_h) : \nabla \bm{\varphi}_i \, d\Omega + \int_{\Gamma_R} \alpha \bd \cdot \bm{\varphi}_i \, d\Gamma, \qquad i = 1, \dots, N_h,$

and $\FFe \in \mathbb{R}^{\Nh}$ denotes the vector of external forces

$\FFe_i (\bd) = \int_{\Omega} \mathbf{f} \cdot \bm{\varphi}_i \, d\Omega + \int_{\Gamma_N} \mathbf{h} \cdot \bm{\varphi}_i \, d\Gamma, \qquad i = 1, \dots, N_h.$

This nonlinear system of equations is solved by means of Newton’s method: given $\bd^{0} = \bd_0$, for $k=0,1,\dots$ until convergence we seek $\delta \bd$ such that

$\begin{equation} \JJ(\bd^k) \, \delta \bd = - \RR(\bd^{k}), \end{equation}$

and then set $\bd_{n+1}^{k+1} = \bd_{n+1}^{k} + \delta \bd$. Here $\JJ$ denotes the Jacobian matrix

$\JJ(\bmw) = \frac{\partial \FFi}{\partial \bd} \left( \bmw \right).$

The assembly of the Jacobian matrix $\JJ$ and the residual vector $\mathbf{R}$ is performed by the CSM_Assembler class. Once the system has been assembled, Dirichlet boundary conditions should be imposed; a detailed description of this step can be found here.

Dynamics

The initial/boundary-value problem reads as follows: find $\bd : \overline{\Omega} \times [0,T] \rightarrow \mathbb{R}^d$ such that

$\begin{equation} \left\{ \begin{aligned} \rho \frac{\partial^2 \bd}{\partial t^2} - \text{div}\left( \PP(\bd) \right) &= \mathbf{f} &\qquad &\text{in } \Omega \\ \bd &= \mathbf{g}& \qquad &\text{on } \Gamma_D \\ \PP(\bd) \mathbf{n} &= \mathbf{h} & \qquad &\text{on } \Gamma_N \\ \PP(\bd) \mathbf{n} + \alpha \bd &= \mathbf{0} & \qquad &\text{on } \Gamma_R \end{aligned} \right. \label{eq:CSM_Dynamic2} \end{equation}$

The weak formulation of \eqref{eq:CSM_Dynamic2} reads: for all $t \in (0,T]$, find $\bd = \bd(t) \in V_D$ such that

$\begin{multline} \int_{\Omega} \rho \frac{\partial^2 \bd}{\partial t^2} \cdot \bmw \, d\Omega + \int_{\Omega} \PP(\bd) : \nabla \bmw \, d\Omega + \int_{\Gamma_R} \alpha \bd \cdot \bmw \, d\Gamma \\ = \int_{\Omega} \mathbf{f} \cdot \bmw \, d\Omega + \int_{\Gamma_N} \mathbf{h} \cdot \bmw \, d\Gamma \qquad \forall \bmw \in V. \label{eq:CSM_DynamicWeak} \end{multline}$

At the algebraic level, its finite element approximation leads to the following second-order dynamical system: given initial displacement $\bd_0 \in \R^{\Nh}$ and velocity $\bmv_0 \in \R^{\Nh}$, for all $t \in (0,T]$ find $\bd = \bd(t) \in \R^{\Nh}$ such that

$\begin{equation} \left\{ \begin{aligned} &\MM \ddot{\bd} + \FFi (\bd) = \FFe(t)\\ & \bd(0) = \bd_0\\ & \dot{\bd}(0) = \bmv_0\\ \end{aligned} \right. \label{eq:CSMt_ALG} \end{equation}$

where $\MM \in \R^{\Nh \times \Nh}$ denotes the mass matrix

$\MM_{ij} = \int_{\Omega} \rho \bm{\varphi}_j \cdot \bm{\varphi}_i \, d\Omega.$

Newmark scheme

A popular scheme used for the solution of this initial-values problem is the Newmark scheme, which consists of the following equations:

Equilibrium: $$ \begin{equation} \MM \ddot{\bd}_{n+1} + \FFi (\bd_{n+1}) = \FFe_{n+1} \label{eq:CSMt_Newmark1} \end{equation} $$ Recurrences: $$ \begin{equation} \begin{aligned} & \bd_{n+1} = \bd_n + \Delta t \bv_n + \frac{\Delta t^2}{2} \left( 2\beta \ba_{n+1} + (1-2\beta) \ba_n\right)\\ & \bv_{n+1} = \bv_n + \Delta t \left( \gamma \ba_{n+1} + (1-\gamma) \ba_n \right)\\ \end{aligned} \label{eq:CSMt_Newmark1b} \end{equation} $$

where $\bd_n$, $\dot{\bd}_n$ and $\ddot{\bd}_n$ are the approximations of $\bd(t_n)$, $\dot{\bd}(t_n)$, and $\ddot{\bd}(t_n)$, respectively. The parameters $\gamma \in [0,1]$ and $\beta \in [0 , 1/2]$ determine the stability and accuracy of the approximation. In particular, the method is second order accurate if and only if $\gamma = \tfrac{1}{2}$, while it is unconditionally stable (for linear systems) for $2 \beta \geq \gamma \geq \tfrac{1}{2}$.

We can use the expansion for $\bd_{n+1}$ and $\bv_{n+1}$ to express the acceleration $\ba_{n+1}$ in terms of the displacement:

$\ba_{n+1} = \frac{1}{\beta \Delta t^2} \bd_{n+1} - \bm{\xi}_{n},$

with

$\bm{\xi}_{n} = \frac{1}{\beta \Delta t^2} \left( \bd_{n} + \Delta t \bv_{n}\right) + \frac{1 - 2\beta}{2\beta} \ba_{n}.$

The algebraic system \eqref{eq:CSMt_Newmark1} can thus be written as follows:

$\begin{equation} \frac{1}{\beta \Delta t^2} \MM \bd_{n+1} + \FFi (\bd_{n+1}) = \FFe_{n+1} + \MM \bm{\xi}_{n}. \label{eq:CSMt_Newmark2} \end{equation}$

At each timestep, this nonlinear system of equations is solved by means of Newton’s method: given $\bd_{n+1}^{0} = \bd_n$, for $k=0,1,\dots$ until convergence we seek $\delta \bd$ such that

$\begin{equation} \JJ(\bd_{n+1}^k) \, \delta \bd = - \RR(\bd_{n+1}^{k}), \end{equation}$

and then set $\bd_{n+1}^{k+1} = \bd_{n+1}^{k} + \delta \bd$. Here,

$\JJ(\bd_{n+1}^k) = \frac{1}{\beta \Delta t^2} \MM + \frac{\partial \FFi}{\partial \bd} \left( {\bd_{n+1}^{k}} \right)$

and

$\RR(\bd_{n+1}^k) = \frac{1}{\beta \Delta t^2} \MM \bd_{n+1}^k + \FFi (\bd_{n+1}^k) - \FFe_{n+1} - \MM \bm{\xi}_{n}.$

Generalized-$\alpha$ method

Numerical damping cannot be introduced in the Newmark method without degrading the order of accuracy. To this end, the so called HHT-$\alpha$ method and the generalized-$\alpha$ method were introduced. The latter reads as follows:

Equilibrium: $$ \begin{equation} \MM \ddot{\bd}_{n+1-\alpha_m} + \FFi (\bd_{n+1-\alpha_f}) = \FFe_{n+1-\alpha_f} \label{eq:CSMt_GenAlpha1} \end{equation} $$ Time averagings: $$ \begin{equation} \begin{aligned} &t_{n+1-\alpha_f} = (1-\alpha_f) t_{n+1} + \alpha_f t_n\\ &\FFe_{n+1-\alpha_f} = \FFe(t_{n+1-\alpha_f}) \\ &\bd_{n+1-\alpha_f} = (1-\alpha_f) \bd_{n+1} + \alpha_f \bd_n \\ &\ba_{n+1-\alpha_m} = (1-\alpha_m) \ba_{n+1} + \alpha_m \ba_n \\ \end{aligned} \label{eq:CSMt_GenAlpha2} \end{equation} $$ Recurrences: $$ \begin{equation} \begin{aligned} & \bd_{n+1} = \bd_n + \Delta t \bv_n + \frac{\Delta t^2}{2} \left( 2\beta \ba_{n+1} + (1-2\beta) \ba_n\right)\\ & \bv_{n+1} = \bv_n + \Delta t \left( \gamma \ba_{n+1} + (1-\gamma) \ba_n \right)\\ \end{aligned} \label{eq:CSMt_GenAlpha3} \end{equation} $$

The generalized-$\alpha$ method reduced to the Nemark method for $\alpha_f = 0$ and $\alpha_m = 0$. Moreover, it is a second-order accurate scheme for linear problems, provided that $\gamma = \tfrac{1}{2} - \alpha_m + \alpha_f$; it is unconditionally stable if $\alpha_m \leq \alpha_f \leq \tfrac{1}{2}$ and $\beta \geq \tfrac{1}{4} + \tfrac{1}{2}(\alpha_f - \alpha_m)$. A family of second-order and absolutely stable methods controlling the numerical dissipation can be defined in terms of a unique parameter $\rho_{\infty}$:

$\alpha_m = \frac{2\rho_{\infty}-1}{\rho_{\infty}+1}, \qquad \alpha_f = \frac{\rho_{\infty}}{\rho_{\infty}+1}, \qquad \beta = \frac{1}{(1+\rho_{\infty})^2}, \qquad \gamma = \frac{3 - \rho_{\infty}}{2(1+\rho_{\infty})}$

We can use the expansion for $\bd_{n+1}$ and $\bv_{n+1}$ to express the acceleration $\ba_{n+1-\alpha_m}$ in terms of the displacement $\bd_{n+1}$:

$\ba_{n+1-\alpha_m} = \frac{1-\alpha_m}{\beta \Delta t^2} \bd_{n+1} - \bm{\phi}_n,$

with

$\bm{\phi}_{n} = (1-\alpha_m) \bm{\xi_n} - \alpha_m \ba_n = \frac{1-\alpha_m}{\beta \Delta t^2} \left( \bd_{n} + \Delta t \bv_{n}\right) + \frac{1 - \alpha_m - 2\beta}{2\beta} \ba_{n}.$

The algebraic system \eqref{eq:CSMt_GenAlpha1} can thus be written as follows:

$\begin{equation} \frac{1-\alpha_m}{\beta \Delta t^2} \MM \bd_{n+1} + \FFi (\bd_{n+1-\alpha_f}) = \FFe_{n+1-\alpha_f} + \MM \bm{\phi}_{n}. \label{eq:CSMt_GenAlpha4} \end{equation}$

At each timestep, this nonlinear system of equations is solved by means of Newton’s method: given $\bd_{n+1}^{0} = \bd_n$, for $k=0,1,\dots$ until convergence we seek $\delta \bd$ such that

$\begin{equation} \JJ(\bd_{n+1}^k) \, \delta \bd = - \RR(\bd_{n+1}^{k}), \end{equation}$

and then set $\bd_{n+1}^{k+1} = \bd_{n+1}^{k} + \delta \bd$. Here,

$\JJ(\bd_{n+1}^k) = \frac{(1-\alpha_m)}{\beta \Delta t^2} \MM + (1-\alpha_f)\frac{\partial \FFi}{\partial \bd} \left( (1-\alpha_f)\bd_{n+1}^{k} + \alpha_f \bd_n \right)$

and

$\RR(\bd_{n+1}^k) = \frac{(1-\alpha_m)}{\beta \Delta t^2} \MM \bd_{n+1}^k + \FFi \left( (1-\alpha_f)\bd_{n+1}^{k} + \alpha_f \bd_n \right) - \FFe_{n+1-\alpha_f} - \MM \bm{\phi}_{n}.$

References

Continuum mechanics and constitutive laws:

M. E. Gurtin. An introduction to continuum mechanics. Vol. 158. Academic press, 1982.
R. W. Ogden. Non-linear elastic deformations. Courier Corporation, 1997.
G. A. Holzapfel. Nonlinear solid mechanics. Vol. 24. Chichester: Wiley, 2000.
M. L. Raghavan, and D. A. Vorp. Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysm: identification of a finite strain constitutive model and evaluation of its applicability. Journal of biomechanics 33(4):475-482, 2000.
M. de Luca. Mathematical and Numerical Models for Cerebral Aneurysm Wall Mechanics. PhD Thesis, Politecnico di Milano, 2009.

Time integration:

N. M. Newmark. “A method of computation for structural dynamics.” Journal of the Engineering Mechanics Division, ASCE, 67–94, 1959.
H. M. Hilber, T. J. R. Hughes, and R. L. Taylor. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engineering and Structural Dynamics, 6:283-292, 1977.
J. Chung and G. M. Hulbert. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method. Journal of applied mechanics, 60(2):371-375, 1993.