Linear Elastic Dynamic - GCMLab/GCMLab-FEM GitHub Wiki

Formulation

The dynamic solver is based on the $\alpha$-Method by Hilber-Huges and Taylor [1]. The discrete form of the balance of linear momentum is expressed as

$$\textbf{M} \textbf{a}_{n+1} + ( 1 + \alpha ) (\textbf{K} \textbf{d}_{n+1} + \textbf{C} \textbf{v}_{n+1}) - \alpha (\textbf{K} \textbf{d}_{n} + \textbf{C} \textbf{v}_{n}) = \textbf{F}(t_{n+1}+\alpha \Delta t),$$

where $\textbf{a}$ is the acceleration, $\textbf{v}$ is the velocity, and $\textbf{d}$ is the displacement at times $n$ or $n+1$, $\Delta t$ is the time step, $\textbf{F}$ is the vector of external forces, $\textbf{M}$ is the mass matrix, $\textbf{K}$ is the stiffness matrix, and $\textbf{C}$ is the damping matrix. The next equations are used to compute $\textbf{d}_{n+1}$ and $\textbf{v}_{n+1}$

$$ \textbf{d}_{n+1} = d_n + \Delta t \textbf{v}_n + \normalsize{{\Delta t^2}\over{2}}\left[ (1-2\beta )\textbf{a}_n + 2\beta \textbf{a}_{n+1}\right] = \textbf{d}_n^* + \Delta t^2 \beta \textbf{a}_{n+1},$$

$$ \textbf{v}_{n+1} = \textbf{v}_n + \Delta t \left[ (1-\gamma )\textbf{a}_n + \gamma \textbf{a}_{n+1}\right] = \textbf{v}_n^* + \cfrac{\gamma}{\Delta t \beta} \left( \textbf{d}_{n+1} -\textbf{d}_n^* \right),$$

in which $$\gamma = \cfrac{1}{2} - \alpha,$$ $$\beta = \cfrac{\left(1- \alpha \right)^2}{4},$$

with $\gamma \in \left[-\cfrac{1}{3}, 0 \right]$ to ensure the method is unconditionally stable and second-order accurate. A second order backwards difference approximation was used to compute $\textbf{d}_n$ and $\textbf{v}_n$ as

$$\textbf{d}_n = \cfrac{1}{2\Delta t} \left(3\textbf{d}_n - 4\textbf{d}_{n-1} + \textbf{d}_{n-2} \right), $$

$$\textbf{a}_n = \cfrac{1}{\Delta t^2} \left(2\textbf{d}_n - 5\textbf{d}_{n-1} + 4\textbf{d}_{n-2} + \textbf{d}_{n-3} \right). $$

Finally, a discrete system to solve for $\textbf{d}_{n+1}$ is obtained

$$ \left[ \cfrac{1}{\Delta t}\textbf{M} + \cfrac{1+\alpha}{\Delta t \beta}\textbf{C} + (1+\alpha)\textbf{K} \right] \textbf{d}_{n+1} = \textbf{F}(t_{n+1}+\alpha \Delta t) + \alpha(\textbf{K} \textbf{d}_{n} + \textbf{C} \textbf{v}_{n}) + \cfrac{1}{\Delta t^2 \beta}\textbf{M}\textbf{d}_n^* - (1+\alpha)\textbf{C} \left( \textbf{v}_n^* - \cfrac{\gamma}{\Delta t \beta}\textbf{d}_n^*\right) $$

References

[1] Hughes, Thomas J. R.. (2000). Finite Element Method - Linear Static and Dynamic Finite Element Analysis - 9.3.3.3 α-Method (Hilber-Hughes-Taylor Method). Dover Publications. Retrieved from https://app.knovel.com/hotlink/pdf/id:kt00B6PNT1/finite-element-method/method-hilber-hughes