MultiPhysicsVault/.raw/FiniteElementProcedures/FiniteElementProcedures_083.md

The solution (9.68) can now be written in the form (9.63); i.e., we have

\left[ \begin{array}{c} t + \Delta t \\ x \\ t \\ x \end{array} \right] = \mathbf {A} \left[ \begin{array}{c} t \\ x \\ t - \Delta t \\ x \end{array} \right] + \mathbf {L} ^ {t} r \tag {9.69}

where

\mathbf {A} = \left[ \begin{array}{c c} \frac {2 - \omega^ {2} \Delta t ^ {2}}{1 + \xi \omega \Delta t} & - \frac {1 - \xi \omega \Delta t}{1 + \xi \omega \Delta t} \\ 1 & 0 \end{array} \right] \tag {9.70}

and

\mathbf {L} = \left[ \begin{array}{c} \frac {\Delta t ^ {2}}{1 + \xi \omega \Delta t} \\ 0 \end{array} \right] \tag {9.71}

As we pointed out in Section 9.2.1, the method is usually employed with \xi=0 .

The Houbolt method. In the Houbolt integration scheme the equilibrium equation (9.62) is considered at time t + \Delta t ; i.e., we use

{ } ^ { t + \Delta t } \ddot { x } + 2 \xi \omega { } ^ { t + \Delta t } \dot { x } + \omega { } ^ { 2 } { } ^ { t + \Delta t } x = { } ^ { t + \Delta t } r \tag {9.72}

with \ddot{x}=\frac{1}{\Delta t^{2}}\left(2^{t+\Delta t}x-5^{t}x+4^{t-\Delta t}x-^{t-2\Delta t}x\right) (9.73)

{ } ^ { t + \Delta t } \dot { x } = \frac { 1 } { 6 \Delta t } \left( 1 1 ^ { t + \Delta t } x - 1 8 ^ { t } x + 9 ^ { t - \Delta t } x - 2 ^ { t - 2 \Delta t } x \right) \tag {9.74}

Substituting (9.73) and (9.74) into (9.72), we can establish the relation

\left[ \begin{array}{c} ^ {t + \Delta t} x \\ ^ {t} x \\ ^ {t - \Delta t} x \end{array} \right] = \mathbf {A} \left[ \begin{array}{c} ^ {t} x \\ ^ {t - \Delta t} x \\ ^ {t - 2 \Delta t} x \end{array} \right] + \mathbf {L} ^ {t + \Delta t} r \tag {9.75}

where

\mathbf {A} = \left[ \begin{array}{c c c} \frac {5 \beta}{\omega^ {2} \Delta t ^ {2}} + 6 \kappa & - \left(\frac {4 \beta}{\omega^ {2} \Delta t ^ {2}} + 3 \kappa\right) & \frac {\beta}{\omega^ {2} \Delta t ^ {2}} + \frac {2 \kappa}{3} \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{array} \right] \tag {9.76}

\beta = \left(\frac {2}{\omega^ {2} \Delta t ^ {2}} + \frac {1 1 \xi}{3 \omega \Delta t} + 1\right) ^ {- 1}; \quad \kappa = \frac {\xi \beta}{\omega \Delta t} \tag {9.77}

and

\mathbf {L} = \left[ \begin{array}{c} \frac {\beta}{\omega^ {2}} \\ 0 \\ 0 \end{array} \right] \tag {9.78}

The Newmark method. In the Newmark integration scheme the equilibrium equation (9.62) is considered at time t + \Delta t ; i.e., we use

{ } ^ { t + \Delta t } \ddot { x } + 2 \xi \omega { } ^ { t + \Delta t } \dot { x } + \omega ^ { 2 } { } ^ { t + \Delta t } x = { } ^ { t + \Delta t } r \tag {9.79}

and the following expansions are employed for the velocity and displacement at time t + \Delta t :

{ } ^ { t + \Delta t } \dot { x } = { } ^ { t } \dot { x } + \left[ ( 1 - \delta ) { } ^ { t } \ddot { x } + \delta { } ^ { t + \Delta t } \ddot { x } \right] \Delta t \tag {9.80}

{ } ^ { t + \Delta t } x = { } ^ { t } x + { } ^ { t } \dot { x } \Delta t + \left[ \left( \frac { 1 } { 2 } - \alpha \right) { } ^ { t } \ddot { x } + \alpha { } ^ { t + \Delta t } \ddot { x } \right] \Delta t ^ { 2 } \tag {9.81}

where \delta and \alpha are parameters to be chosen to obtain optimum stability and accuracy. Newmark proposed as an unconditionally stable scheme the constant-average-acceleration method, in which case \delta = \frac{1}{2} and \alpha = \frac{1}{4} .

Substituting for ^{t + \Delta t}\dot{x} and ^{t + \Delta t}x into (9.79), we can solve for ^{t + \Delta t}\ddot{x} and then use (9.80) and (9.81) to calculate ^{t + \Delta t}\dot{x} and ^{t + \Delta t}x . We can thus establish the relation

\left[ \begin{array}{l} ^ {t + \Delta t} \ddot {x} \\ ^ {t + \Delta t} \dot {x} \\ ^ {t + \Delta t} x \end{array} \right] = \mathbf {A} \left[ \begin{array}{l} ^ {t} \ddot {x} \\ ^ {t} \dot {x} \\ ^ {t} x \end{array} \right] + \mathbf {L} ^ {t + \Delta t} r \tag {9.82}

where \mathbf{A}=\left[\begin{array}{ccc}-\left(\frac{1}{2}-\alpha\right)\beta-2(1-\delta)\kappa&\frac{1}{\Delta t}(-\beta-2\kappa)&\frac{1}{\Delta t^{2}}(-\beta)\\\Delta t\left[1-\delta-\left(\frac{1}{2}-\alpha\right)\delta\beta-2(1-\delta)\delta\kappa\right]&1-\beta\delta-2\delta\kappa&\frac{1}{\Delta t}(-\beta\delta)\\\Delta t^{2}\left[\frac{1}{2}-\alpha-\left(\frac{1}{2}-\alpha\right)\alpha\beta-2(1-\delta)\alpha\kappa\right]&\Delta t(1-\alpha\beta-2\alpha\kappa)&(1-\alpha\beta)\end{array}\right] (9.83)

\beta = \left(\frac {1}{\omega^ {2} \Delta t ^ {2}} + \frac {2 \xi \delta}{\omega \Delta t} + \alpha\right) ^ {- 1}; \quad \kappa = \frac {\xi \beta}{\omega \Delta t} \tag {9.84}

and

\mathbf {L} = \left[ \begin{array}{c} \frac {\beta}{\omega^ {2} \Delta t ^ {2}} \\ \frac {\beta \delta}{\omega^ {2} \Delta t} \\ \frac {\alpha \beta}{\omega^ {2}} \end{array} \right] \tag {9.85}

The Bathe method. In the Bathe integration scheme, the equilibrium equation (9.62) is solved at t + \Delta t

{ } ^ { t + \Delta t } \ddot { x } + 2 \xi \omega { } ^ { t + \Delta t } \dot { x } + \omega { } ^ { 2 t + \Delta t } x = { } ^ { t + \Delta t } r \tag {9.86}

with two equal sub-steps: in the first sub-step the trapezoidal rule is used

{ } ^ { t + \Delta t / 2 } \dot { x } = { } ^ { t } \dot { x } + \frac { \Delta t } { 4 } \left( { } ^ { t } \ddot { x } + { } ^ { t + \Delta t / 2 } \ddot { x } \right) \tag {9.87}

^ {t + \Delta t / 2} x = ^ {t} x + \frac {\Delta t}{4} \left(^ {t} \dot {x} + ^ {t + \Delta t / 2} \dot {x}\right) \tag {9.88}

and in the second sub-step, the 3-point Euler backward method is employed, see L. Collatz [A],

{ } ^ { t + \Delta t } \dot { x } = \frac { 1 } { \Delta t } { } ^ { t } x - \frac { 4 } { \Delta t } { } ^ { t + \Delta t / 2 } x + \frac { 3 } { \Delta t } { } ^ { t + \Delta t } x \tag {9.89}

{ } ^ { t + \Delta t } \ddot { x } = \frac { 1 } { \Delta t } { } ^ { t } \dot { x } - \frac { 4 } { \Delta t } { } ^ { t + \Delta t / 2 } \dot { x } + \frac { 3 } { \Delta t } { } ^ { t + \Delta t } \dot { x } \tag {9.90}

Using (9.87) and (9.88) with the equilibrium equation (9.62) at time t + \Delta t/2 and then (9.89) and (9.90) with (9.86), we obtain

\left[ \begin{array}{l} t + \Delta t \\ \ddot {x} \\ t + \Delta t \\ \dot {x} \\ t + \Delta t \\ x \end{array} \right] = \mathbf {A} \left[ \begin{array}{l} t \\ \ddot {x} \\ t \\ \dot {x} \\ t \\ x \end{array} \right] + \mathbf {L} _ {a} ^ {t + \Delta t / 2} r + \mathbf {L} _ {b} ^ {t + \Delta t} r \tag {9.91}

where A, L_{u} , and L_{b} are the integration approximation and load operators, respectively,

\mathbf {A} = \frac {1}{\beta_ {1} \beta_ {2}} \left[ \begin{array}{c c} - 4 \omega \Delta t (2 4 \xi + 7 \omega \Delta t) & \omega (- 2 8 8 \xi + 1 4 \xi \omega^ {2} \Delta t ^ {2} - 1 4 4 \omega \Delta t + 5 \omega^ {3} \Delta t ^ {3} + 4 8 \xi^ {2} \omega \Delta t) \\ - 4 \Delta t (- 1 2 + \omega^ {2} \Delta t ^ {2}) & 1 4 4 - 4 7 \omega^ {2} \Delta t ^ {2} - 8 \xi \omega^ {3} \Delta t ^ {3} - 2 4 \xi \omega \Delta t \\ 4 \Delta t ^ {2} (7 + 2 \xi \omega \Delta t) & \Delta t (1 4 4 - 5 \omega^ {2} \Delta t ^ {2} + 8 0 \xi \omega \Delta t + 1 6 \xi^ {2} \omega^ {2} \Delta t ^ {2}) \end{array} \right.

\left. \begin{array}{c} \omega^ {2} (2 4 \xi \omega \Delta t + 1 9 \omega^ {2} \Delta t ^ {2} - 1 4 4) \\ \omega^ {2} \Delta t (- 9 6 - 2 4 \xi \omega \Delta t + \omega^ {2} \Delta t ^ {2}) \\ - 1 9 \omega^ {2} \Delta t ^ {2} + 1 4 4 + 1 6 8 \xi \omega \Delta t + 4 8 \xi^ {2} \omega^ {2} \Delta t ^ {2} - 2 \xi \omega^ {3} \Delta t ^ {3} \end{array} \right] \tag {9.92}

\mathbf {L} _ {a} = \frac {1}{\beta_ {1} \beta_ {2}} \left[ \begin{array}{l} - 4 \omega \Delta t (2 4 \xi + 7 \omega \Delta t) \\ - 4 \Delta t (- 1 2 + \omega^ {2} \Delta t ^ {2}) \\ 4 \Delta t ^ {2} (7 + 2 \xi \omega \Delta t) \end{array} \right]; \quad \mathbf {L} _ {b} = \frac {1}{\beta_ {2}} \left[ \begin{array}{l} 9 \\ 3 \Delta t \\ \Delta t ^ {2} \end{array} \right] \tag {9.93}

with \beta_{1} = 16 + 8\xi \omega \Delta t + \omega^{2}\Delta t^{2}; \beta_{2} = 9 + 6\xi \omega \Delta t + \omega^{2}\Delta t^{2} (9.94)

Of course, the same analysis can also be performed when using a general intermediate time position, that is, t + \gamma\Delta t instead of t + \Delta t/2 , where \gamma is a parameter, and it turns out that for

\gamma = 2 - \sqrt {2} \tag {9.95}

the amplitude decay is at its maximum and the period elongation is at its minimum (but not much different than at \gamma = 0.5 ), see K.J. Bathe and G. Noh [A]. The idea of sub-stepping to introduce numerical dissipation (for reducing the solution error) has also been used to design an explicit integration scheme in G. Noh and K.J. Bathe [A].

9.4.2 Stability Analysis

The aim in the numerical integration of the finite element system equilibrium equations is to evaluate a good approximation of the actual dynamic response of the structure under consideration. In order to predict the dynamic response of the structure accurately, it would

seem that all system equilibrium equations in (9.58) must be integrated to high precision, and this means that all n equations of the form (9.62) need to be integrated accurately. Since in direct integration the same time step is used for each equation of the form (9.62), \Delta t would have to be selected corresponding to the smallest period in the system, which may mean that the time step is very small indeed. As an estimate of \Delta t thus required, it appears that if the smallest period is T_{n} , \Delta t would have to be about T_{n}/10 (or smaller; see Section 9.4.3). However, we discussed in Section 9.3.2 that in many analyses practically all dynamic response lies in only some modes of vibration and that for this reason only some mode shapes are considered in mode superposition analysis. In addition, it was pointed out that in many analyses there is little justification to include the response predicted in the higher modes because the frequencies and mode shapes of the finite element mesh can be only crude approximations to the “exact” quantities. Therefore, the finite element idealization has to be chosen in such a way that the lowest p frequencies and mode shapes of the structure are predicted accurately, where p is determined by the distribution and frequency content of the loading.

It can therefore be concluded that in many analyses we are interested in integrating only the first p equations of the n equations in (9.61) accurately. This means that we would be able to revise \Delta t to be about T_{p}/10 , i.e., T_{p}/T_{n} times larger than our first estimate. In practical analysis the ratio T_{p}/T_{n} can be very large, say of the order 1000, meaning that the analysis would be much more effective using \Delta t = T_{p}/10 . However, assuming that we select a time step \Delta t of magnitude T_{p}/10 , we realize that in the direct integration also the response in the higher modes is automatically integrated with the same time step. Since we cannot possibly integrate accurately the response in those modes for which \Delta t is larger than half the natural period T, an important question is: What “response” is predicted in the numerical integration of (9.62) when \Delta t/T is large? This is, in essence, the question of stability of an integration scheme. Stability of an integration method means that the physical initial conditions for the equations with a large value \Delta t/T must not be amplified artificially and thus render worthless any accuracy in the integration of the lower mode response. Stability also means that any “initial” conditions at time t given by errors in the displacements, velocities, and accelerations, which may be due to round-off in the computer, do not grow in the integration. Stability is ensured if the time step is small enough to integrate accurately the response in the highest-frequency component. But this may require a very small step, and, as was pointed out earlier, the accurate integration of the high-frequency response predicted by the finite element assemblage is in many cases not justified and therefore not necessary.

The stability of an integration method is therefore determined by examining the behavior of the numerical solution for arbitrary initial conditions. Therefore, we consider the integration of (9.62) when no load is satisfied; i.e., r = 0. The solution for prescribed initial conditions only as obtained from (9.64) is, hence,

{ } ^ { t + n \Delta t } \hat { \mathbf { X } } = \mathbf { A } ^ { n } { } ^ { t } \hat { \mathbf { X } } \tag {9.96}

Considering the stability of integration methods, we have procedures that are unconditionally stable and that are only conditionally stable. An integration method is unconditionally stable if the solution for any initial conditions does not grow without bound for any time step \Delta t , in particular when \Delta t/T is large. The method is only conditionally stable if the above only holds provided that \Delta t/T is smaller than or equal to a certain value, usually called the stability limit.

For the stability analysis we use the spectral decomposition of A given by A = PJP^{-1} , where P is the matrix of eigenvectors of A, and J is the Jordan canonical form of A with eigenvalues \lambda_{i} of A on its diagonal. We considered in Section 2.5 the case of A being symmetric, in which case J = \Lambda [see (2.108)], P = V , and P^{-1} = V^{T} . However, the approximation operator A is in general a nonsymmetric matrix, and we therefore must use the more general decomposition A = PJP^{-1} , in which J is not necessarily a diagonal matrix but may exhibit unit elements on the superdiagonal line (corresponding to multiple eigenvalues), (see, for example, J. H. Wilkinson [A]).

Of course, using the above spectral decomposition, we have

\mathbf {A} ^ {n} = \mathbf {P} \mathbf {J} ^ {n} \mathbf {P} ^ {- 1} \tag {9.97}

and with this expression we can determine the stability of the time integration scheme.

Let \rho(\mathbf{A}) be the spectral radius \mathbf{A} defined as

\rho (\mathbf {A}) = \max _ {i = 1, 2, \dots} \left| \lambda_ {i} \right| \tag {9.98}

where the absolute value sign requires the evaluation of the absolute magnitude of \lambda_{i} in the complex plane. Then our stability criterion is that

1. If all eigenvalues are distinct, we must have \rho(\mathbf{A}) \leq 1 , whereas
2. If A contains multiple eigenvalues, we require that all such eigenvalues be (in modulus) smaller than 1.

If this stability criterion is fulfilled, we have J^{n} and hence A^{n} bounded for n \to \infty . Furthermore, if \rho(\mathbf{A}) < 1 , we have J^{n} \to 0 and hence A^{n} \to 0 , and the decrease in A^{n} is more rapid when \rho(\mathbf{A}) is small.

Since the stability of an integration method depends only on the eigenvalues of the approximation operator, it may be convenient to apply a similarity transformation on A before evaluating the eigenvalues. For example, in the case of the Newmark method we apply the similarity transformation D^{-1}AD , where D is a diagonal matrix with d_{ii} = (\Delta t)^{i} . As might be expected, we then find that the spectral radius and therefore the stability of the integration method depends only on the time ratio \Delta t/T , the damping ratio \xi , and the integration parameters used. Therefore, for any given \Delta t/T and \xi , it is possible in the Newmark method and the Bathe method to vary the parameters \alpha , \delta , and \gamma respectively, to obtain optimum stability and accuracy characteristics.

Consider as a simple example the stability analysis of the central difference method.

EXAMPLE 9.12: Analyze the central difference method for its integration stability. Consider the case \xi = 0.0 used in (9.8) to (9.12).

We need to calculate the spectral radius of the approximation operator given in (9.70) when \xi = 0 . The eigenvalue problem \mathbf{Au} = \lambda \mathbf{u} to be solved is

\left[ \begin{array}{c c} 2 - \omega^ {2} \Delta t ^ {2} & - 1 \\ 1 & 0 \end{array} \right] \mathbf {u} = \lambda \mathbf {u} \tag {a}

The eigenvalues are the roots of the characteristic polynomial p(\lambda) (see Section 2.5) defined as

p (\lambda) = (2 - \omega^ {2} \Delta t ^ {2} - \lambda) (- \lambda) + 1

\lambda_ {1} = \frac {2 - \omega^ {2} \Delta t ^ {2}}{2} + \sqrt {\frac {\left(2 - \omega^ {2} \Delta t ^ {2}\right) ^ {2}}{4} - 1}

Hence,

\lambda_ {2} = \frac {2 - \omega^ {2} \Delta t ^ {2}}{2} - \sqrt {\frac {\left(2 - \omega^ {2} \Delta t ^ {2}\right) ^ {2}}{4} - 1}

For stability we need that the absolute values of \lambda_{1} and \lambda_{2} be smaller than or equal to 1; i.e., the spectral radius \rho(\mathbf{A}) of the matrix \mathbf{A} in (a) must satisfy the condition \rho(\mathbf{A}) \leq 1 , and this gives the condition \Delta t/T \leq 1/\pi . Hence, the central difference method is stable provided that \Delta t \leq \Delta t_{cr} , where \Delta t_{cr} = T_{n}/\pi . It is interesting to note that this same time step stability limit is also applicable when \xi > 0 (see Exercise 9.14).

By the same procedure as employed in Example 9.12, the Houbolt, Newmark, and Bathe methods can be analyzed for stability using the corresponding approximation operators; Fig. 9.4 shows the stability characteristics. It is noted that the central difference method is only conditionally stable as evaluated in Example 9.12 and that the Newmark, Houbolt, and Bathe methods are unconditionally stable.

line

| Δt / T | Central difference method | Newmark method δ = 1/2, α = 1/4 | δ = 11/20, α = 3/10 | Bathe method | Houbolt method | | ------ | ------------------------- | -------------------------------- | ------------------- | ------------ | -------------- | | 10⁻³ | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | | 10⁻² | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | | 10⁻¹ | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | | 10⁰ | 1.0 | 1.0 | 1.0 | 0.8 | 0.6 | | 10¹ | 1.0 | 1.0 | 0.9 | 0.4 | 0.2 | | 10² | 1.0 | 1.0 | 0.8 | 0.1 | 0.05 | | 10³ | 1.0 | 1.0 | 0.8 | 0.05 | 0.02 | | 10⁴ | 1.0 | 1.0 | 0.8 | 0.02 | 0.01 |

Figure 9.4 Spectral radii of approximation operators, case \xi = 0.0

Fig. 9.5 focuses on the spectral radii of the Newmark method with various parameters, the Houbolt and the Bathe methods and also the Wilson \theta method, see E.L. Wilson, I. Farhoomand and K.J. Bathe [A], for small time steps. The behavior in the range shown is of interest because it is desirable to have a spectral radius very close to 1.0 until about \Delta t / T = 0.3 and then a very rapid decrease to 0.0 thereafter, like displayed in Fig. 9.4 and 9.5 for the Bathe method.

Considering the Newmark method, the two parameters \alpha and \delta can be varied to obtain optimum stability and accuracy. The integration scheme is unconditionally stable provided

that \delta\geq0.5 and \alpha\geq0.25(\delta+0.5)^{2} . The method much used in practice corresponds to \delta=0.5 and \alpha=0.25 (the trapezoidal rule), but we discuss in the next section that the scheme can show undesirable characteristics in practical analyses.

line

Δt / T	ρ(A) - Central difference method	ρ(A) - Newmark method δ=1/2, α=1/4	ρ(A) - Bathe method	ρ(A) - Newton method δ=11/20, α=3/10	ρ(A) - Wilson θ method θ=1.4	ρ(A) - Houbolt method
0.05	1.0	1.0	1.0	1.0	1.0	1.0
0.1	1.0	1.0	1.0	1.0	1.0	1.0
0.15	1.0	1.0	1.0	0.95	0.95	0.95
0.2	1.0	1.0	1.0	0.9	0.9	0.9
0.25	1.0	1.0	1.0	0.85	0.85	0.85
0.3	1.0	1.0	1.0	0.8	0.8	0.8
0.35	1.0	1.0	1.0	0.75	0.75	0.75
0.4	1.0	1.0	1.0	0.7	0.7	0.7

Figure 9.5 Spectral radii for small time steps

Although only some widely used integration schemes have been discussed, the above considerations already show that the analyst has to make a choice as to which method to use. This choice is influenced by the accuracy characteristics of the method, i.e., the accuracy that can be obtained in the integration for a given time step \Delta t .

9.4.3 Accuracy Analysis

The decision as to which integration operator to use in a practical analysis is governed by the cost of solution for a desired accuracy, which in turn is determined by the number of time steps required in the integration. If a conditionally stable algorithm such as the central difference method is employed, the time step size, and hence the number of time steps for a given time range considered, is in practice largely determined by the critical time step \Delta t_{cr} . However, using an unconditionally stable operator, the time step has to be chosen to yield an accurate and effective solution. Because the direct integration of the equilibrium equations in (9.58) is equivalent to integrating simultaneously all n decoupled equations of the form (9.62), we can study the integration accuracy of (9.62) as a function of \Delta t/T , \xi and r. The solution to (9.62) was given in (9.64), and it is this equation that we use to assess the integration errors.

Let us consider for a simple accuracy analysis the solution of the initial value problem defined by

\left. \begin{array}{l} \ddot {x} + \omega^ {2} x = 0 \\ ^ {0} x = 1. 0; \quad^ {0} \dot {x} = 0. 0; \quad^ {0} \ddot {x} = - \omega^ {2} \end{array} \right\} \tag {9.99}

for which the exact solution is x = \cos \omega t . For a complete analysis we would have to consider

also the initial value problem corresponding to ^{0}x = 0.0 , ^{0}\dot{x} = \omega , and ^{0}\ddot{x} = 0 , with the exact solution x = \sin \omega t , and the solution for a general loading condition. In addition, the influence of the damping parameter \xi would need to be investigated. However, the significant numerical characteristics can be demonstrated by considering only the numerical solution of the problem in (9.99).

The Newmark and Bathe methods can be used directly with the initial conditions given in (9.99). However, in the Houbolt method, the initial conditions are defined only by initial displacements, and in the following study the exact displacement values for \Delta t x and 2\Delta t x obtained using the solution x = \cos \omega t have been employed.

The numerical solution of (9.99) using the different integration methods shows that the errors in the integrations can be measured in terms of period elongation and amplitude decay. Figure 9.6 shows the percentage period elongations and amplitude decays in the implicit integration schemes discussed as a function of \Delta t/T . These relationships have been obtained by evaluating (9.64) and comparing the exact solution x = \cos \omega t with the numerical solutions (see K. J. Bathe and E. L. Wilson [A]). For curves of period elongations and amplitude decays of the Newmark method when other parameters are used and some explicit schemes, see K.J. Bathe and G. Noh [A] and G. Noh and K.J. Bathe [A].

line

Δt / T	Houbolt method	Newmark method (δ = 1/2, α = 1/4)	Bathe method
0.00	0.00	0.00	0.00
0.05	2.00	1.00	0.50
0.10	5.00	3.00	2.00
0.15	10.00	6.00	4.00
0.20	15.00	10.00	6.00

line

Δt / T	Houbolt method	Bathe method
0.00	0.0	0.0
0.05	5.0	0.5
0.10	10.0	1.0
0.15	15.0	1.5
0.20	20.0	2.0

Figure 9.6 Percentage period elongations and amplitude decays

The curves in Fig. 9.6 show that, in general, the numerical integrations using any of the methods are accurate when \Delta t/T is smaller than about 0.01. However, when the time step/period ratio is larger, the various integration methods exhibit quite different characteristics. Notably, for a given ratio \Delta t/T , the Bathe method introduces less amplitude decay and period elongation than the Houbolt method, and the Newmark constant-average-acceleration method introduces only period elongation and no amplitude decay. Since for a given \Delta t , the period elongation depends on T, the use of the direct time integration schemes introduces numerical dispersion in wave propagation solutions, see G. Noh, S. Ham, and K.J.

Bathe [A].

The characteristics of the integration errors exhibited in Fig. 9.6 are used in the discussion of the simultaneous integration of all n equations of form (9.62) as required in the solution of (9.61) and thus in the solution of (9.58). We observe that all equations for which the time step to period ratio is small are integrated accurately, but that the response in the equations for which \Delta t/T is large is not obtained with any precision.

These considerations lead to the choice of an appropriate time step \Delta t . Using the central difference method, the time step has to be chosen such that \Delta t is smaller than or equal to \Delta t_{cr} evaluated in Example 9.12. Only in the exceptional case in which the loading or initial conditions significantly excite the highest frequencies do we need to use a much smaller time step than \Delta t_{cr} . However, using one of the unconditionally stable schemes, we realize that the time step \Delta t can be much larger and should only be small enough that the response in all modes that significantly contribute to the total structural response is calculated accurately. The other modal response components are not evaluated accurately in practice (in fact, can not be evaluated accurately for a well-chosen mesh, see Section 9.4.4), and are frequently best rapidly damped out to not pollute the total response solution.

As an example, K.J. Bathe and G. Noh [A] considered the 'model problem' in Figure 9.7. Here, a very stiff spring connected to a flexible spring is considered. While very simple, the model problem represents the essence of stiff and flexible parts in a complex structural system, containing many degrees of freedom. The high stiffness spring represents very stiff parts or constraints (introduced by stiff springs or beams) while the flexible spring represents the flexible parts of the complex structural model. Of course, in a mode superposition solution, the detailed response within the stiff parts would naturally not be included.

When using the Newmark and Bathe methods, the displacements are quite accurately predicted but the acceleration of the mass m_{2} shows very large errors when using the trapezoidal rule, which means that the reaction is also solved for very inaccurately. These errors can be controlled using the Newmark method by choosing parameters to introduce numerical damping but the appropriate values may depend on the problem solved. The same also holds when using the generalized alpha method, see J. Chung and G. M. Hulbert [A].

k _ {1} = 1 0 ^ {7}, k _ {2} = 1, m _ {1} = 0, m _ {2} = 1, m _ {3} = 1, \omega_ {\rho} = 1. 2

Figure 9.7 Two spring model problem; using \Delta t = 0.2618 (hence \Delta t/T_{p} = 0.05 , \Delta t/T_{l} = 0.0417 and \Delta t/T_{2} = 131.76 ) the acceleration at node 2 using the trapezoidal rule is grossly in error whereas the Bathe method gives an accurate prediction (for the first half-step only, the parameters \delta = 3/4 , \alpha = 1 are best used)

In the practical direct integration solution of such systems, it is effective to use an implicit scheme that in essence "works like a mode superposition solution", in that it automatically does not include the undesirable response — here the response in possibly artificially stiff components of the finite element model. Considering explicit schemes, the Noh-Bathe method does have this characteristic as well but with a critical time step and hence is mostly only valuable for wave propagation analyses, see G. Noh and K.J. Bathe [A].

9.4.4 Some Practical Considerations

In order to obtain an effective solution of a dynamic response, it is important to choose an appropriate time integration scheme. This choice depends on the finite element idealization, which in turn depends on the actual physical problem to be analyzed. It follows therefore that the selection of an appropriate finite element idealization of a problem and the choice of an effective integration scheme for the response solution are closely related and must be considered together. The finite element model and time integration scheme are chosen differently depending on whether a structural dynamics or a wave propagation problem is solved.

Structural Dynamics

The basic consideration in the selection of an appropriate finite element model of a structural dynamics problem is that only the lowest modes (or only a few intermediate modes) of a physical system are being excited by the load vector. We discussed this consideration already in Section 9.3.2, where we addressed the problem of how many modes need to be included in a mode superposition analysis. Referring to Section 9.3.2, we can conclude that if a Fourier analysis of the dynamic load input shows that only frequencies below \omega_{u} are contained in the loading, then the finite element mesh should at most represent accurately the frequencies to about \omega_{co} = 4\omega_{u} of the actual system. There is no need to represent the higher frequencies of the actual physical system accurately in the finite element system because the dynamic response contribution in these frequencies is negligible; i.e., for values of \hat{\omega}/\omega in Fig. 9.3 smaller than 0.25, an almost static response is measured and this response is directly included in the direct integration step-by-step dynamic response calculations. We should also note that the static displacement response in the higher modes is small if the mode shapes are almost orthogonal to the load vector and/or the frequencies are high, see (9.46) (i.e., when the structure is very stiff in these mode shapes). When the static response is small it may well be expedient to use \omega_{co} closer to \omega_{u} than given by the factor 4 used above. This reduction of \omega_{co} may be considerable, for example, in earthquake response analysis of certain structures.

The complete procedure for the modeling of a structural vibration problem is therefore:

1. Identify the frequencies significantly contained in the loading, using a Fourier analysis if necessary. These frequencies may change as a function of time. Let the highest frequency significantly contained in the loading be \omega_{u} .
2. Choose a finite element mesh that can accurately represent the static response and accurately represents all frequencies up to about \omega_{co} = 4\omega_{u} .