Dynamic Analyzers

2. Dynamic Analyzers

2.1 Newmark’s Method for structural problems

Equation of equilibrium that governs the dynamic response of a linear system:

$\mathbf{M} \mathbf{\ddot{U}} + \mathbf{C} \mathbf{\dot{U}} + \mathbf{K} \mathbf{U} = \mathbf{R}$

where M, C and K are the mass, damping and stiffness matrices and U is the displacement vector.

Substituting the inertia forces, the damping forces and the elastic forces, the equation of equilibrium can be written as:

$\mathbf{F}_{I}(t) + \mathbf{F}_{D}(t) + \mathbf{F}_{E}(t) = \mathbf{R}$

For the solution of the dynamic system, we use the Newmark direct integration scheme, that makes the following assumptions:

$^{t+\Delta t}{\mathbf{\dot{U}}} = ^{t}{\mathbf{\dot{U}}} + \left[(1-\delta)^{t}{\mathbf{\ddot{U}}} + \delta ^{t+\Delta t}{\mathbf{\ddot{U}}}\right] \Delta t$

$^{t+\Delta t}{\mathbf{U}} = ^{t}{\mathbf{U}} + ^{t}{\mathbf{\dot{U}}}\Delta t + \left[ \left(\frac{1}{2}-\alpha\right) ^{t}{\mathbf{\ddot{U}}} + \alpha ^{t+\Delta t}{\mathbf{\ddot{U}}} \right] \Delta t^{2}$

where $\alpha$ and $\delta$ are parameters that can be determined to obtain integration accuracy and stability. For the linear acceleration method, the parameters become $\alpha=1/6$ and $\delta=1/2$ , while for Newmark’s method we have $\alpha=1/4$ and $\delta=1/2$ . According to above, the equation of equilibrium becomes:

$\mathbf{M} ^{t+\Delta t}{\mathbf{\ddot{U}}}+ \mathbf{C} ^{t+\Delta t}{\mathbf{\dot{U}} }+ \mathbf{K} ^{t+\Delta t}{\mathbf{U}} = ^{t+\Delta t}{\mathbf{R}}$

Using the Trapezoidal rule, we solve for each time step:

$\left( \frac{4}{\Delta t^{2}} \mathbf{M} + \frac{2}{\Delta t} \mathbf{C} + \mathbf{K}\right) ^{t+\Delta t }{\mathbf{U}} = ^{t+\Delta t}{\mathbf{R}} + \mathbf{M} \left(\frac{4}{\Delta t^{2}} ^{t}{\mathbf{U}} + \frac{4}{\Delta t} ^{t}{\mathbf{\dot{U}}} + ^{t}{\mathbf{\ddot{U}}} \right) + \mathbf{C} \left(\frac{2}{\Delta t} ^{t }{\mathbf{U}} + ^{t}{\mathbf{\dot{U}}} \right)$

and then we can calculate $^{t+\Delta t}{\mathbf{\ddot{U}}}$ and $^{t+\Delta t}{\mathbf{\dot{U}}}$ .

Implementation for Linear Dynamic Analysis

The algorithm for the Newmark integration scheme can be summarized in the following step-by-step procedure:

A. Initial calculations:

Form stiffness matrix K, mass matrix M and damping matrix C.
Initialize $^{0}{\mathbf{U}}$ , $^{0}{\mathbf{\dot{U}}}$ , $^{0}{\mathbf{\ddot{U}}}$ .
Select time step $\Delta t$ and parameters $\delta \geq 0.50$ , $\alpha \geq 0.25(0.5+\delta)^{2}$ .
Calculate integration constants:

$\alpha_{0} = \frac{1}{\alpha \Delta t^{2}}$; $\alpha_{1} = \frac{\delta}{\alpha \Delta t}$; $\alpha_{2} = \frac{1}{\alpha \Delta t}$; $\alpha_{3} = \frac{1}{2\alpha}-1$ $\alpha_{4} = \frac{\delta}{\alpha}-1$; $\alpha_{5}=\frac{\Delta t}{2} \left(\frac{\delta}{\alpha}-2\right)$; $\alpha_{6}=\Delta t (1-\delta) $; $\alpha_{7} = \delta \Delta t$

Form effective stiffness matrix

$\mathbf{\hat{K}} = \mathbf{K}+\alpha_{0}\mathbf{M}+\alpha_{1}\mathbf{C}$

Triangularize:

$\mathbf{\hat{K}} =\mathbf{L}\mathbf{D}\mathbf{L}^{T}$

B. For each time step:

Calculate effective loads:

$^{t+\Delta t}{\mathbf{\hat{R}}}= ^{t+\Delta t}{\mathbf{R}}+\mathbf{M}\left(\alpha_{0} ^{t}{\mathbf{U}} + \alpha_{2} ^{t}{\mathbf{\dot{U}}} + \alpha_{3}^{t}{\mathbf{\ddot{U}}} \right) + \mathbf{C}\left( \alpha_{1} ^{t}{\mathbf{U}} + \alpha_{4} ^{t}{\mathbf{\dot{U}}} + \alpha_{5}^{t}{\mathbf{\ddot{U}}} \right)$

Calculate the displacements vector:

$\mathbf{L}\mathbf{D}\mathbf{L}^{T} \left(^{t+\Delta t}{\mathbf{U}}\right) = ^{t+\Delta t}{\mathbf{\hat{R}}}$ $\mathbf{L}\mathbf{D}\mathbf{L}^{T} \left(^{t+\Delta t}{\mathbf{U}}\right) = ^{t+\Delta t}{\mathbf{\hat{R}}}$

Calculate accelerations and velocities:

$^{t+\Delta t}{}{\mathbf{\ddot{U}}}=\alpha_{0} \left(^{t+\Delta t}{}{\mathbf{U}}-^{t}{}{\mathbf{U}} \right) - \alpha_{2} ^{t}{}{\mathbf{\dot{U}}} - \alpha_{3} ^{t}{}{\mathbf{\ddot{U}}}$

$^{t+\Delta t}{}{\mathbf{\dot{U}}} = ^{t}{}{\mathbf{\dot{U}}} + \alpha_{6} ^{t}{}{\mathbf{\ddot{U}}} + \alpha_{7} ^{t+\Delta t}{}{\mathbf{\ddot{U}}}$

Implementation for Nonlinear Dynamic Analysis

The solution of nonlinear dynamic system of equations can be calculated using Newmark’s method combined with the incremental-iterative procedure of the Newton-Raphson method. For illustration here the modified Newton-Raphson iteration method without damping effects is demonstrated as:

$\mathbf{M}^{t+\Delta t}{}{\mathbf{\ddot{U}}}^{(k)} + \mathbf{K}^{t+\Delta t}{}{\mathbf{U}}^{(k)} = ^{t+\Delta t}{}{\mathbf{\hat{R}}} - ^{t+\Delta t}{}{\mathbf{F}}^{(k-1)}$

$^{t+\Delta t}{}{\mathbf{U}}^{(k)} = ^{t+\Delta t}{}{\mathbf{U}}^{(k-1)} + \Delta\mathbf{U}^{(k)}$

Using the trapezoidal rule of time integration according to Newmark’s method, we get the following assumptions:

$^{t+\Delta t}{}{\mathbf{U}} = ^{t}{}{\mathbf{U}} + \frac{\Delta t}{2} \left( ^{t}{}{\mathbf{\dot{U}}} + ^{t+\Delta t}{}{\mathbf{\dot{U}}} \right)$

$^{t+\Delta t}{}{\mathbf{\dot{U}}} = ^{t}{}{\mathbf{\dot{U}}} + \frac{\Delta t}{2} \left( ^{t}{}{\mathbf{\ddot{U}}} + ^{t+\Delta t}{}{\mathbf{\ddot{U}}} \right)$

According to the above we can calculate the accelerations vector, as:

$^{t+\Delta t}{}{\mathbf{\ddot{U}}}^{(k)} = \frac{4}{\Delta t^{2}}\left(^{t+\Delta t}{}{\mathbf{U}}^{(k-1)} - ^{t}{}{\mathbf{U}} + \Delta\mathbf{U}^{(k)} \right) - \frac{4}{\Delta t}^{t}{}{\mathbf{\dot{U}}} - ^{t}{}{\mathbf{\ddot{U}}}$

The solution can be calculated from the solution of the following system:

$\mathbf{\hat{K}} \Delta\mathbf{U}^{(k)} = ^{t+\Delta t}{}{\mathbf{R}} - ^{t+\Delta t}{}{\mathbf{F}}^{(k-1)} - \mathbf{M} \left[\frac{4}{\Delta t^{2}}\left(^{t+\Delta t}{}{\mathbf{U}}^{(k-1)} - ^{t}{}{\mathbf{U}}\right) - \frac{4}{\Delta t} ^{t}{}{\mathbf{\dot{U}}} - ^{t}{}{\mathbf{\ddot{U}}} \right]$

where:

$\mathbf{\hat{K}} = ^{t}{}{\mathbf{K}} + \frac{4}{\Delta t^{2}}\mathbf{M}$

is the effective stiffness matrix.

2.2 Finite Difference Method for heat transfer problems

The transient heat conduction equation is given by:

$\frac{\partial}{\partial x} \left ( k_x(T)\frac{\partial T}{\partial x} \right)+\frac{\partial}{\partial y}\left ( k_y(T)\frac{\partial T}{\partial y}\right )+\frac{\partial}{\partial z}\left ( k_z(T)\frac{\partial T}{\partial z}\right )+Q=\rho c \frac{\partial T}{\partial t}$

where $T(x,y,z)$ is the temperature field, Q is the heat generation per unit volume and $k,\rho,c$ are the material’s conductivity, density and specific heat capacity, respectively.

The boundary conditions of this type of problem are:

$T=T_b \ \ \ \text{on} \ \ \Gamma_b$

and

$k_x(T)\frac{\partial T}{\partial x}l+k_y(T)\frac{\partial T}{\partial y}m+k_z(T)\frac{\partial T}{\partial z}n + q +h(T-T_a) \ \ \ \text{on} \ \ \Gamma_q$

where, $\Gamma_b \bigcup \Gamma_q=\Gamma \ \ \text{and} \ \ \Gamma_b \bigcap \Gamma_q=\emptyset$ with $\Gamma$ representing the whole boundary. Also, l,m,n are direction cosines, $h$ is the heat transfer coefficient, $T_a$ is the atmospheric temperature and $q$ is the boundary heat flux. The initial condition for the problem is

$T=T_0 \ \ \ \text{at} \ \ t=0.0$

It is now possible to solve the system, provided that appropriate spatial and temporal discretizations are available. The Galerkin method can be used for the spatial approximation and the above equation becomes:

$\mathbf{C}\dot{\mathbf{T}}+\mathbf{K}\mathbf{T}=\mathbf{f}$

with $\mathbf{C}$ , $\mathbf{K}$ and $\mathbf{f}$ being the capacity matrix, the conductivity matrix and the load vector. More specifically:

$C_{ij}=\int_{\Omega}\rho cN_iN_jd\Omega$

$K_{ij}=\int_{\Omega} \left [ k_x(T) \frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+ k_y(T) \frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}+k_z(T) \frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z} \right ] d\Omega + \int_{\Gamma} h N_i N_j d\Gamma$

$f_i=\int_{\Omega}N_iQd\Omega-\int_{\Gamma_q}qN_id\Gamma_q+\int_{\Gamma_q}hN_iT_ad\Gamma_q$

or, in matrix form,

$\mathbf{C}=\int_{\Omega} \rho c \mathbf{N}^T \mathbf{N} d\Omega$

$\mathbf{K}=\int_{\Omega}\mathbf{B}^T \mathbf{D} \mathbf{B} d\Omega+\int_{\Gamma} h\mathbf{N}^T \mathbf{N} d\Gamma$

$\mathbf{f}=\int_{\Omega}Q\mathbf{B}^Td\Omega-\int_{\Gamma_q} q\mathbf{N}^T d\Gamma_q+\int_{\Gamma_q} hT_a\mathbf{N}^T d\Gamma_q$

For the solution of the dynamic system, we use a Finite Difference scheme, that makes the following assumptions:

$\mathbf{C}\frac{\mathbf{T}^{n+1}-\mathbf{T}^n}{\Delta t}+\mathbf{K}\left ( \theta\mathbf{T}^{n+1}+(1-\theta)\mathbf{T}^n \right )=\theta \mathbf{f}^{n+1}+(1-\theta)\mathbf{f}^n$

$(\frac{1}{\Delta t}\mathbf{C}+\theta\mathbf{K}) \mathbf{T}^{n+1}=\left ( \frac{1}{\Delta t}\mathbf{C}-(1-\theta) \mathbf{K} \right )\mathbf{T}^n+\left ( \theta \mathbf{f}^{n+1} + (1-\theta)\mathbf{f}^n \right )$

The above equation gives the nodal values of temperature at the $n+1$ time level. By varying the parameter $\theta$ , different transient schemes can be constructed. For instance, $\theta=0$ results in the forward difference method (fully explicit), $\theta=1$ in the backward difference method (fully implicit) and $\theta=0.5$ in the Crank-Nicolson method (semi-implicit).