김경종
467 lines
31 KiB
Markdown
467 lines
31 KiB
Markdown
<!-- source-page: 781 -->
|
||
|
||
$$
|
||
\left(^ {\prime} \lambda + \frac {1}{2} \Delta \lambda^ {(1)}\right) \mathbf {R} ^ {T} \Delta \mathbf {U} ^ {(1)} = W \tag {8.112}
|
||
$$
|
||
|
||
and $(^{t + \Delta t}\lambda^{(i - 1)} + \frac{1}{2}\Delta \lambda^{(i)})\mathbf{R}^T\Delta \mathbf{U}^{(i)} = 0$ $i = 2,3,\ldots$
|
||
|
||
where W is selected based on the history of iterations in the previous incremental steps. Figure 8.14(b) illustrates this constraint equation. The use of this scheme can be particularly effective near collapse points.
|
||
|
||
To solve the governing equations we may rewrite (8.107) to obtain
|
||
|
||
$$
|
||
{ } ^ { \tau } \mathbf { K } \Delta \overline { { \mathbf { U } } } ^ { ( i ) } = { } ^ { t + \Delta t } \lambda ^ { ( i - 1 ) } \mathbf { R } - { } ^ { t + \Delta t } \mathbf { F } ^ { ( i - 1 ) } \tag {8.113}
|
||
$$
|
||
|
||
$$
|
||
\mathbf {K} \Delta \overline {{\overline {{\mathbf {U}}}}} = \mathbf {R} \tag {8.114}
|
||
$$
|
||
|
||
and hence, $\Delta \mathbf{U}^{(i)} = \Delta \overline{\mathbf{U}}^{(i)} + \Delta \lambda^{(i)}\Delta \overline{\overline{\mathbf{U}}} \tag{8.115}$
|
||
|
||
Employing the spherical constant arc length criterion (8.111), we then use
|
||
|
||
$$
|
||
\boldsymbol {\lambda} ^ {(i)} = \boldsymbol {\lambda} ^ {(i - 1)} + \Delta \boldsymbol {\lambda} ^ {(i)} \tag {8.116}
|
||
$$
|
||
|
||
and $\mathbf{U}^{(i)} = \mathbf{U}^{(i - 1)} + \Delta \overline{\mathbf{U}}^{(i)} + \Delta \lambda^{(i)}\Delta \overline{\overline{\mathbf{U}}} \tag{8.117}$
|
||
|
||
Substituting from (8.116) and (8.117) into (8.111) gives a quadratic equation in $\Delta\lambda^{(i)}$ . We select the appropriate value to proceed with the solution (see Exercise 8.35).
|
||
|
||
Using the constant increment in external work criterion, $\Delta \lambda^{(1)}$ is directly calculated from (8.112) and then the $\Delta \lambda^{(i)}$ values for $i = 2, 3, \ldots$ , are obtained from (8.112) as
|
||
|
||
$$
|
||
\Delta \lambda^ {(i)} = - \frac {\mathbf {R} ^ {T} \Delta \overline {{{\mathbf {U}}}} ^ {(i)}}{\mathbf {R} ^ {T} \Delta \overline {{{\overline {{{\mathbf {U}}}}}}}} \tag {8.118}
|
||
$$
|
||
|
||
The relation (8.112) also admits the solution $t + \Delta t \lambda^{(i)} = -t + \Delta t \lambda^{(i-1)}$ , but this solution corresponds to a load reversal, which we disregard.
|
||
|
||
A complete solution algorithm based on the above load-displacement-constraint procedures must of course also contain a special scheme to start the incremental solution and must have procedures to self-adaptively select $\Delta l$ and/or W. Also, the algorithm should stop iterating when divergence is imminent and then restart itself with new iterative parameters. Complete solution methods with these ingredients are very valuable and are in common use for the analysis of the collapse response of structures.
|
||
|
||
# 8.4.4 Convergence Criteria
|
||
|
||
If an incremental solution strategy based on iterative methods is to be effective, realistic criteria should be used for the termination of the iteration. At the end of each iteration, the solution obtained should be checked to see whether it has converged within preset tolerances or whether the iteration is diverging. If the convergence tolerances are too loose, inaccurate results are obtained, and if the tolerances are too tight, much computational effort is spent to obtain needless accuracy. Similarly, an ineffective divergence check can terminate the iteration when the solution is not actually diverging or force the iteration to search for an unattainable solution. The objective in this section is to discuss briefly some convergence criteria.
|
||
|
||
Since we are seeking the displacement configuration corresponding to time $t + \Delta t$ , it is natural to require that the displacements at the end of each iteration be within a certain tolerance of the true displacement solution. Hence, a realistic convergence criterion is
|
||
|
||
$$
|
||
\frac {\left\| \Delta \mathbf {U} ^ {(i)} \right\| _ {2}}{\left\| ^ {t + \Delta t} \mathbf {U} \right\| _ {2}} \leq \epsilon_ {D} \tag {8.119}
|
||
$$
|
||
|
||
<!-- source-page: 782 -->
|
||
|
||
where $\epsilon_{D}$ is a displacement convergence tolerance. The vector $^{t+\Delta t}U$ is not known and must be approximated. Frequently, it is appropriate to use in (8.119) the last calculated value $^{t+\Delta t}U^{(i)}$ as an approximation to $^{t+\Delta t}U$ and a sufficiently small value $\epsilon_{D}$ . However, in some analyses the actual solution may still be far from the value obtained when convergence is measured using (8.119) with $^{t+\Delta t}U^{(i)}$ . This is the case when the calculated displacements change only little in each iteration but continue to change for many iterations, as may occur, for example, in elastoplastic analysis under loading conditions when the modified Newton-Raphson iteration is used.
|
||
|
||
A second convergence criterion is obtained by measuring the out-of-balance load vector. For example, we may require that the norm of the out-of-balance load vector be within a preset tolerance $\epsilon_{F}$ of the original load increment
|
||
|
||
$$
|
||
\left\| ^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(i)} \right\| _ {2} \leq \epsilon_ {F} \left\| ^ {t + \Delta t} \mathbf {R} - ^ {t} \mathbf {F} \right\| _ {2} \tag {8.120}
|
||
$$
|
||
|
||
A difficulty with this criterion is that the displacement solution does not enter the termination criterion. As an illustration of this difficulty, consider an elastoplastic truss with a very small strain-hardening modulus entering the plastic region. In this case, the out-of-balance loads may be very small while the displacements may still be much in error. Hence, the convergence criteria in (8.119) and (8.120) may have to be used with very small values of $\epsilon_{D}$ and $\epsilon_{F}$ . Also, the expressions must be modified appropriately when quantities of different units are measured (such as displacements, rotations, pressures, and so on).
|
||
|
||
In order to provide some indication of when both the displacements and the forces are near their equilibrium values, a third convergence criterion may be useful in which the increment in internal energy during each iteration (i.e., the amount of work done by the out-of-balance loads on the displacement increments) is compared to the initial internal energy increment. Convergence is assumed to be reached when, with $\epsilon_{E}$ a preset energy tolerance,
|
||
|
||
$$
|
||
\Delta \mathbf {U} ^ {(i) ^ {T}} \left(^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(i - 1)}\right) \leq \epsilon_ {E} \left(\Delta \mathbf {U} ^ {(1) ^ {T}} \left(^ {t + \Delta t} \mathbf {R} - ^ {t} \mathbf {F}\right)\right) \tag {8.121}
|
||
$$
|
||
|
||
Since this convergence criterion contains both the displacements and the forces, it is in practice an attractive measure. Some experiences with these convergence measures are given by K. J. Bathe and A. P. Cimento [A]. An important point is that the convergence tolerances $\epsilon_{D}$ , $\epsilon_{F}$ , and $\epsilon_{E}$ may need to be quite small in some solutions in order to reach a good solution accuracy. In general, use of the full Newton-Raphson method in the incremental solution leads to a higher accuracy in the solution than use of the modified Newton-Raphson method since, if convergence occurs, the solution error diminishes quite rapidly in the last iterations of the full Newton-Raphson method (consider, for example, Exercises 8.40 and 9.31).
|
||
|
||
# 8.4.5 Exercises
|
||
|
||
8.31. Consider the single degree of freedom system shown.
|
||
|
||
(a) Use the full Newton-Raphson iteration method, the initial stress method, and the BFGS method to calculate the response of the system.
|
||
(b) Establish a value of the constant $c$ for which the full Newton-Raphson method will not converge.
|
||
|
||
|
||
|
||
<details>
|
||
<summary>text_image</summary>
|
||
|
||
R, u
|
||
R = 4
|
||
</details>
|
||
|
||
$F =$ force in the spring $= u + cu^3;c = 0.1$
|
||
|
||
<!-- source-page: 783 -->
|
||
|
||
8.32. Consider the single degree of freedom system in Exercise 8.31 but with $F = \sin(u/L)$ , L = 1.0, and R = 0.5. Perform the solution as in Exercise 8.31(a).
|
||
8.33. Consider the system in Exercise 8.31 and solve for the response by line searching only. (Hence do not perform any Newton-type iteration.)
|
||
8.34. Consider the four-node plane stress element shown.
|
||
|
||
(a) Use a computer program to calculate the stiffness coefficient corresponding to the displacement $u$ . [Hint: Use (8.89) in finite difference form.]
|
||
|
||
(b) Also calculate this stiffness coefficient using the formulation given in Chapter 6.
|
||
|
||
|
||
|
||
<details>
|
||
<summary>text_image</summary>
|
||
|
||
0.3
|
||
0.2
|
||
2
|
||
2
|
||
0.1
|
||
0.3
|
||
0.4
|
||
E = 200,000
|
||
v = 0.3
|
||
Total Lagrangian
|
||
formulation
|
||
Unit thickness
|
||
</details>
|
||
|
||
8.35. Derive the quadratic equation for $\Delta\lambda^{(i)}$ used in the spherical constant arc length criterion and calculate the roots of this equation. Discuss the solutions and identify which solution you would use in the practical implementation.
|
||
8.36. Use a computer program to solve for the collapse and postbuckling response of the simple arch structure considered in Example 6.3.
|
||
8.37. Use a computer program to solve for the collapse response of the three element truss shown. Compare your results with an analytical solution. (Hint: Use a large displacement formulation and a load-displacement-constraint solution method. You may also refer to the paper by P. G. Hodge, K. J. Bathe, and E. N. Dvorkin [A].)
|
||
|
||
|
||
|
||
<details>
|
||
<summary>text_image</summary>
|
||
|
||
H = 5
|
||
Bar areas = 1.0
|
||
E = 200,000
|
||
E_T = 0.0
|
||
σ_y = 100
|
||
①
|
||
②
|
||
σ_y
|
||
element ③
|
||
σ_y
|
||
3σ_y
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
①
|
||
②
|
||
③
|
||
</details>
|
||
|
||
<!-- source-page: 784 -->
|
||
|
||
8.38. Use a computer program to calculate the collapse and postcollapse response of the structure shown. Consider various areas $A_{1}$ that you choose.
|
||
|
||
|
||
|
||
<details>
|
||
<summary>text_image</summary>
|
||
|
||
E = 3 × 10⁶
|
||
Area = 10
|
||
Area = A₁
|
||
Δ₁
|
||
Δ₂
|
||
P
|
||
5
|
||
5
|
||
5
|
||
3
|
||
5
|
||
</details>
|
||
|
||
8.39. Use a computer program to calculate the response of the plane stress cantilever shown. Use the von Mises yield condition with isotropic hardening and increase the load P until full collapse of the structure.
|
||
|
||
Compare the solution efficiencies when using the full Newton-Raphson, modified Newton-Raphson, and the BFGS methods and also use a load-displacement-constraint procedure.
|
||
|
||
|
||
|
||
<details>
|
||
<summary>text_image</summary>
|
||
|
||
E = 200,000 MPa
|
||
v = 0.30
|
||
σy = 200 MPa
|
||
ET = 20 MPa
|
||
4 mm
|
||
100 mm
|
||
</details>
|
||
|
||
Width = 1 mm
|
||
|
||
8.40. Use a computer program to solve for the large displacement response of the cantilever beam shown below. Increase P to reach a tip deflection of $\Delta \doteq 10$ in. Compare the solution efficiencies when using the full Newton-Raphson and modified Newton-Raphson methods, with or without line searches, and the BFGS method.
|
||
|
||
|
||
|
||
<details>
|
||
<summary>text_image</summary>
|
||
|
||
0.1 in
|
||
10.0 in
|
||
P
|
||
Δ
|
||
E = 1.2 × 10⁴ lb/in²
|
||
v = 0.2
|
||
</details>
|
||
|
||
Width = 1.0 in
|
||
|
||
<!-- source-page: 785 -->
|
||
|
||
# Solution of Equilibrium Equations in Dynamic Analysis
|
||
|
||
# 9.1 INTRODUCTION
|
||
|
||
In Section 4.2.1 we derived the equations of equilibrium governing the linear dynamic response of a system of finite elements
|
||
|
||
$$
|
||
\mathbf {M} \ddot {\mathbf {U}} + \mathbf {C} \dot {\mathbf {U}} + \mathbf {K} \mathbf {U} = \mathbf {R} \tag {9.1}
|
||
$$
|
||
|
||
where M, C, and K are the mass, damping, and stiffness matrices; R is the vector of externally applied loads; and U, $\dot{U}$ , and $\ddot{U}$ are the displacement, velocity, and acceleration vectors of the finite element assemblage. It should be recalled that (9.1) was derived from considerations of statics at time t; i.e., (9.1) may be written
|
||
|
||
$$
|
||
\mathbf {F} _ {I} (t) + \mathbf {F} _ {D} (t) + \mathbf {F} _ {E} (t) = \mathbf {R} (t) \tag {9.2}
|
||
$$
|
||
|
||
where $\mathbf{F}_{I}(t)$ are the inertia forces, $\mathbf{F}_{I}(t) = \mathbf{M}\ddot{\mathbf{U}}$ , $\mathbf{F}_{D}(t)$ are the damping forces, $\mathbf{F}_{D}(t) = \mathbf{C}\dot{\mathbf{U}}$ , and $\mathbf{F}_{E}(t)$ are the elastic forces, $\mathbf{F}_{E}(t) = \mathbf{K}\mathbf{U}$ , all of them being time-dependent. Therefore, in dynamic analysis, in principle, static equilibrium at time t, which includes the effect of acceleration-dependent inertia forces and velocity-dependent damping forces, is considered. Vice versa, in static analysis the equations of motion in (9.1) are considered, with inertia and damping effects neglected.
|
||
|
||
The choice for a static or dynamic analysis (i.e., for including or neglecting velocity- and acceleration-dependent forces in the analysis) is usually decided by engineering judgment, the objective thereby being to reduce the analysis effort required. However, it should be realized that the assumptions of a static analysis should be justified since otherwise the analysis results are meaningless. Indeed, in nonlinear analysis the assumption of neglecting inertia and damping forces may be so severe that a solution may be difficult or impossible to obtain.
|
||
|
||
Mathematically, (9.1) represents a system of linear differential equations of second order and, in principle, the solution to the equations can be obtained by standard procedures for the solution of differential equations with constant coefficients (see, for example,
|
||
|
||
<!-- source-page: 786 -->
|
||
|
||
L. Collatz [A]). However, the procedures proposed for the solution of general systems of differential equations can become very expensive if the order of the matrices is large—unless specific advantage is taken of the special characteristics of the coefficient matrices K, C, and M. In practical finite element analysis, we are therefore mainly interested in a few effective methods and we will concentrate in the next sections on the presentation of those techniques. The procedures that we will consider are divided into two methods of solution: direct integration and mode superposition. Although the two techniques may at first sight appear to be quite different, in fact, they are closely related, and the choice for one method or the other is determined only by their numerical effectiveness.
|
||
|
||
In the following we consider first (see Sections 9.2 to 9.4) the solution of the linear equilibrium equations (9.1), and then we discuss the solution of the nonlinear equations of finite element systems idealizing structures and solids (see Section 9.5). Finally, we show in Section 9.6 how the basic concepts discussed are also directly applicable to the analysis of heat transfer and fluid flow.
|
||
|
||
# 9.2 DIRECT INTEGRATION METHODS
|
||
|
||
In direct integration the equations in (9.1) are integrated using a numerical step-by-step procedure, the term “direct” meaning that prior to the numerical integration, no transformation of the equations into a different form is carried out. In essence, direct numerical integration is based on two ideas. First, instead of trying to satisfy (9.1) at any time t, it is aimed to satisfy (9.1) only at discrete time intervals $\Delta t$ apart. This means that, basically, (static) equilibrium, which includes the effect of inertia and damping forces, is sought at discrete time points within the interval of solution. Therefore, it appears that all solution techniques employed in static analysis can probably also be used effectively in direct integration. The second idea on which a direct integration method is based is that a variation of displacements, velocities, and accelerations within each time interval $\Delta t$ is assumed. As will be discussed in detail, it is the form of the assumption on the variation of displacements, velocities, and accelerations within each time interval that determines the accuracy, stability, and cost of the solution procedure.
|
||
|
||
In the following, assume that the displacement, velocity, and acceleration vectors at time 0, denoted by ${}^{0}U, {}^{0}\dot{U}$ , and ${}^{0}\ddot{U}$ , respectively, are known and let the solution to (9.1) be required from time 0 to time T. In the solution the time span under consideration, T, is subdivided into n equal time intervals $\Delta t$ (i.e., $\Delta t = T/n$ ), and the integration scheme employed establishes an approximate solution at times $\Delta t, 2\Delta t, 3\Delta t, \ldots, t, t + \Delta t, \ldots, T$ . Since an algorithm calculates the solution at the next required time from the solutions at the previous times considered, we derive the algorithms by assuming that the solutions at times 0, $\Delta t, 2\Delta t, \ldots, t$ are known and that the solution at time $t + \Delta t$ is required next.
|
||
|
||
In the following sections some effective direct integration methods are presented, but only schemes that solve the full dynamic equilibrium equations at specific times and only use the solutions obtained at those times. Other schemes, like the Wilson $\theta$ method (see E. L. Wilson, I. Farhoomand, and K.J. Bathe [A]), the HHT method (see H.M. Hilber, T. J. R. Hughes, and R. L. Taylor [A]) and the generalized $\alpha$ method (see J. Chung and G. M. Hulbert [A]), and variations thereof, have some undesirable properties (see for example J.M. Benítez and F. J. Montáns [A]). We also assume a constant time step size, a restriction, that could of course be easily removed.
|
||
|
||
<!-- source-page: 787 -->
|
||
|
||
# 9.2.1 The Central Difference Method
|
||
|
||
If the equilibrium relation in (9.1) is regarded as a system of ordinary differential equations with constant coefficients, it follows that any convenient finite difference expressions to approximate the accelerations and velocities in terms of displacements can be used. Therefore, a large number of different finite difference expressions could theoretically be employed. However, the solution scheme should be effective, and it follows that only a few methods need to be considered. A simple widely used scheme (although with limitations, see G. Noh and K.J. Bathe [A]) is the central difference method, in which it is assumed that
|
||
|
||
$$
|
||
{ } ^ { t } \ddot { \mathbf { U } } = \frac { 1 } { \Delta t ^ { 2 } } \left( { } ^ { t - \Delta t } \mathbf { U } - 2 { } ^ { t } \mathbf { U } + { } ^ { t + \Delta t } \mathbf { U } \right) \tag {9.3}
|
||
$$
|
||
|
||
The error in the expansion (9.3) is of order $(\Delta t)^{2}$ , and to have the same order of error in the velocity expansion, we can use
|
||
|
||
$$
|
||
^ {t} \dot {\mathbf {U}} = \frac {1}{2 \Delta t} \left(- ^ {t - \Delta t} \mathbf {U} + ^ {t + \Delta t} \mathbf {U}\right) \tag {9.4}
|
||
$$
|
||
|
||
The displacement solution for time $t + \Delta t$ is obtained by considering (9.1) at time t, i.e.,
|
||
|
||
$$
|
||
\mathbf {M} ^ {\prime} \ddot {\mathbf {U}} + \mathbf {C} ^ {\prime} \dot {\mathbf {U}} + \mathbf {K} ^ {\prime} \mathbf {U} = ^ {\prime} \mathbf {R} \tag {9.5}
|
||
$$
|
||
|
||
Substituting the relations for 'Ü and 'Ü in (9.3) and (9.4), respectively, into (9.5), we obtain
|
||
|
||
$$
|
||
\left(\frac {1}{\Delta t ^ {2}} \mathbf {M} + \frac {1}{2 \Delta t} \mathbf {C}\right) ^ {\prime + \Delta t} \mathbf {U} = ^ {\prime} \mathbf {R} - \left(\mathbf {K} - \frac {2}{\Delta t ^ {2}} \mathbf {M}\right) ^ {\prime} \mathbf {U} - \left(\frac {1}{\Delta t ^ {2}} \mathbf {M} - \frac {1}{2 \Delta t} \mathbf {C}\right) ^ {\prime - \Delta t} \mathbf {U} \tag {9.6}
|
||
$$
|
||
|
||
from which we can solve for $t^{+\Delta t}\mathbf{U}$ . It should be noted that the solution of $t^{+\Delta t}\mathbf{U}$ is thus based on using the equilibrium conditions at time $t$ ; i.e., $t^{+\Delta t}\mathbf{U}$ is calculated by using (9.5). For this reason the integration procedure is called an explicit integration method, and it is noted that such integration schemes do not require a factorization of the (effective) stiffness matrix in the step-by-step solution. On the other hand, the Houbolt, Newmark and Bathe methods, considered in the next sections, use the equilibrium conditions at time $t + \Delta t$ and are called implicit integration methods.
|
||
|
||
A second observation is that using the central difference method, the calculation of $t^{+\Delta t}\mathbf{U}$ involves $^t\mathbf{U}$ and $t^{-\Delta t}\mathbf{U}$ . Therefore, to calculate the solution at time $\Delta t$ , a special starting procedure must be used. Since $^0\mathbf{U}$ , $^0\dot{\mathbf{U}}$ , and $^0\ddot{\mathbf{U}}$ are known [note that with $^0\mathbf{U}$ and $^0\dot{\mathbf{U}}$ known, $^0\ddot{\mathbf{U}}$ can be calculated using (9.1) at time 0; see Example 9.1], the relations in (9.3) and (9.4) can be used to obtain $^{-\Delta t}\mathbf{U}$ ; i.e., we have
|
||
|
||
$$
|
||
{ } ^ { - \Delta t } U _ { i } = { } ^ { 0 } U _ { i } - \Delta t { } ^ { 0 } \dot { U } _ { i } + \frac { \Delta t ^ { 2 } } { 2 } { } ^ { 0 } \ddot { U } _ { i } \tag {9.7}
|
||
$$
|
||
|
||
where the subscript i indicates the ith element of the vector considered. Table 9.1 summarizes the time integration scheme as it might be implemented in the computer.
|
||
|
||
We discuss below the fact that the method is only effective when each time step solution can be performed very efficiently (because a small time step size and therefore a large number of time steps generally need to be used). For this reason, the method is largely
|
||
|
||
<!-- source-page: 788 -->
|
||
|
||
TABLE 9.1 Step-by-step solution using central difference method (general mass and damping matrices)
|
||
|
||
# A. Initial calculations:
|
||
|
||
1. Form stiffness matrix K, mass matrix M, and damping matrix C.
|
||
2. Initialize $^{0}U$ , $^{0}\dot{U}$ , and $^{0}\ddot{U}$ .
|
||
3. Select time step $\Delta t, \Delta t \leq \Delta t_{cr}$ , and calculate integration constants:
|
||
|
||
$$
|
||
a _ {0} = \frac {1}{\Delta t ^ {2}}; \quad a _ {1} = \frac {1}{2 \Delta t}; \quad a _ {2} = 2 a _ {0}; \quad a _ {3} = \frac {1}{a _ {2}}
|
||
$$
|
||
|
||
4. Calculate $^{-\Delta t}U = ^{0}U - \Delta t ^{0}\dot{U} + a_{3} ^{0}\ddot{U}$
|
||
|
||
5. Form effective mass matrix $\hat{M} = a_{0}M + a_{1}C$ .
|
||
|
||
6. Triangularize $\tilde{\mathbf{M}}$ : $\tilde{\mathbf{M}} = \mathbf{LDL}^T$ .
|
||
|
||
# B. For each time step:
|
||
|
||
1. Calculate effective loads at time t:
|
||
|
||
$$
|
||
^ {\prime} \hat {\mathbf {R}} = ^ {\prime} \mathbf {R} - (\mathbf {K} - a _ {2} \mathbf {M}) ^ {\prime} \mathbf {U} - (a _ {0} \mathbf {M} - a _ {1} \mathbf {C}) ^ {\prime - \Delta t} \mathbf {U}
|
||
$$
|
||
|
||
2. Solve for displacements at time $t + \Delta t$ :
|
||
|
||
$$
|
||
\mathbf {L D L} ^ {T} t + \Delta t \mathbf {U} = ^ {t} \hat {\mathbf {R}}
|
||
$$
|
||
|
||
3. If required, evaluate accelerations and velocities at time t:
|
||
|
||
$$
|
||
^ {t} \ddot {\mathbf {U}} = a _ {0} (^ {t - \Delta t} \mathbf {U} - 2 ^ {t} \mathbf {U} + ^ {t + \Delta t} \mathbf {U})
|
||
$$
|
||
|
||
$$
|
||
^ {t} \dot {\mathbf {U}} = a _ {1} (- ^ {t - \Delta t} \mathbf {U} + ^ {t + \Delta t} \mathbf {U})
|
||
$$
|
||
|
||
applied only when a lumped mass matrix can be assumed and velocity-dependent damping can be neglected. Then (9.6) reduces to
|
||
|
||
$$
|
||
\left(\frac {1}{\Delta t ^ {2}} \mathbf {M}\right) ^ {t + \Delta t} \mathbf {U} = ^ {t} \hat {\mathbf {R}} \tag {9.8}
|
||
$$
|
||
|
||
where $\hat{\mathbf{R}} = \mathbf{R} - \left(\mathbf{K} - \frac{2}{\Delta t^{2}}\mathbf{M}\right)\mathbf{U} - \left(\frac{1}{\Delta t^{2}}\mathbf{M}\right)^{t - \Delta t}\mathbf{U}$ (9.9)
|
||
|
||
Therefore, if the mass matrix is diagonal, the system of equations in (9.1) can be solved without factorizing a matrix; i.e., only matrix multiplications are required to obtain the right-hand-side effective load vector $\hat{R}$ after which the displacement components are obtained using
|
||
|
||
$$
|
||
^ {t + \Delta t} U _ {i} = ^ {t} \hat {R} _ {i} \left(\frac {\Delta t ^ {2}}{m _ {i i}}\right) \tag {9.10}
|
||
$$
|
||
|
||
where $^{t+\Delta t}U_{i}$ and $^{t}\hat{R}_{i}$ denote the $i$ th components of the vectors $^{t+\Delta t}U$ and $^{t}\hat{R}$ , respectively, $m_{ii}$ is the $i$ th diagonal element of the mass matrix, and it is assumed that $m_{ii} > 0$ .
|
||
|
||
If the stiffness matrix of the element assemblage is not to be triangularized, it is also not necessary to assemble the matrix. We have shown in Section 4.2.1 [see (4.30)] that
|
||
|
||
$$
|
||
\mathbf {K} ^ {t} \mathbf {U} = \sum_ {i} \mathbf {K} ^ {(i)} {} ^ {t} \mathbf {U} = \sum_ {i} ^ {t} \mathbf {F} ^ {(i)} \tag {9.11}
|
||
$$
|
||
|
||
which means that $\mathbf{K}^{\prime}\mathbf{U}$ as required in (9.9) can be evaluated on the element level by summing the contributions from each element to the effective load vector. Hence, $\hat{\mathbf{R}}$ can be
|
||
|
||
<!-- source-page: 789 -->
|
||
|
||
evaluated efficiently using
|
||
|
||
$$
|
||
^ \prime \hat {\mathbf {R}} = ^ {\prime} \mathbf {R} - \sum_ {i} ^ {\prime} \mathbf {F} ^ {(i)} - \frac {1}{\Delta t ^ {2}} \mathbf {M} (^ {i} - ^ {\Delta t} \mathbf {U} - 2 ^ {\prime} \mathbf {U}) \tag {9.12}
|
||
$$
|
||
|
||
where $^{t}F^{(i)}$ is evaluated and added in compacted form (see Section 12.2.3).
|
||
|
||
Another advantage of using the central difference method in the form given in $(9.10)$ now becomes apparent. Since no stiffness matrix of the complete element assemblage needs to be calculated, the solution can essentially be carried out on the element level and relatively little high speed storage is required. Using this approach, systems of very large order have been solved effectively.
|
||
|
||
However, we have already mentioned that for solution a relatively small time step size must generally be used. Actually, an important consideration in the use of the central difference scheme is that the integration method requires that the time step $\Delta t$ be smaller than a critical value, $\Delta t_{cr}$ , which can be calculated from the mass and stiffness properties of the complete element assemblage. More specifically, we will show in Section 9.4.2 that to obtain a valid solution,
|
||
|
||
$$
|
||
\Delta t \leq \Delta t _ {c r} = \frac {T _ {n}}{\pi} \tag {9.13}
|
||
$$
|
||
|
||
where $T_{n}$ is the smallest period of the finite element assemblage with n degrees of freedom. The period $T_{n}$ could be calculated using one of the techniques discussed in Chapter 11, or a lower bound on $T_{n}$ may be evaluated using norms (see Section 2.7). In practice we frequently estimate an appropriate time step $\Delta t$ using the considerations given in Section 9.4.4.
|
||
|
||
In the solution using (9.10), it was assumed that $m_{ii} > 0$ for all i. The relation in (9.13) states this requirement once more because a zero diagonal element in a diagonal mass matrix means that the element assemblage has a zero period (see Section 10.2.4). In general, all diagonal elements of the mass matrix can be assumed to be larger than zero, in which case (9.13) gives a limit on the magnitude of the time step $\Delta t$ that can be used in the integration. In the analysis of some problems (namely, for wave propagation problems) (9.13) does not require an unduly small time step; but, in other cases (namely, for structural dynamics problems) the time step small enough for accuracy of integration should be much larger than $\Delta t_{cr}$ obtained from (9.13).
|
||
|
||
These thoughts point toward the importance of establishing an effective finite element discretization and time step for a dynamic solution. We discuss these issues in Section 9.4, but already mention the following considerations now.
|
||
|
||
Assume that we solve using the central difference method a large system of equilibrium equations. The time step for the integration would be selected using (9.13). Assume that we now change the smallest diagonal element of the mass matrix to become very small and, in fact, nearly zero. As enumerated above, a diagonal element in the mass matrix cannot be exactly zero because $T_{n}$ would then be zero and the integration would not be possible. However, as the diagonal element in the mass matrix approaches zero, the smallest period of the system, and hence $\Delta t_{cr}$ , approaches zero. Therefore, the reduction of one element $m_{ii}$ necessitates a severe reduction in the time step size that can be used in the integration. On the other hand, since the order of the system is large, we can envisage a certain dynamic loading for which we would hardly expect that the response of the element assemblage changes very much when the smallest element $m_{ii}$ is reduced, even to become
|
||
|
||
<!-- source-page: 790 -->
|
||
|
||
zero. Hence, the cost of analysis would in that case be unduly large only because of one very small diagonal element in the mass matrix. The same condition is also reached when an element in the stiffness matrix is changed to become large.
|
||
|
||
Integration schemes that require the use of a time step $\Delta t$ smaller than a critical time step $\Delta t_{cr}$ , such as the central difference method, are said to be conditionally stable. If a time step larger than $\Delta t_{cr}$ is used, the integration is unstable, meaning that, for example, any errors resulting from round-off in the computer grow and make the response calculations worthless in most cases. The concept of stability of integration is very important, and we will discuss it further in Section 9.4. However, at this stage it is useful to consider the following example.
|
||
|
||
EXAMPLE 9.1: Consider a simple system for which the governing equilibrium equations are
|
||
|
||
$$
|
||
\left[ \begin{array}{l l} 2 & 0 \\ 0 & 1 \end{array} \right] \left[ \begin{array}{l} \ddot {U} _ {1} \\ \ddot {U} _ {2} \end{array} \right] + \left[ \begin{array}{c c} 6 & - 2 \\ - 2 & 4 \end{array} \right] \left[ \begin{array}{l} U _ {1} \\ U _ {2} \end{array} \right] = \left[ \begin{array}{l} 0 \\ 1 0 \end{array} \right] \tag {a}
|
||
$$
|
||
|
||
The free-vibration periods of the system are given in Example 9.6, where we find that $T_{1} = 4.45$ , $T_{2} = 2.8$ . Use the central difference method in direct integration with time steps (1) $\Delta t = T_{2} / 10$ and (2) $\Delta t = 10T_{2}$ to calculate the response of the system for 12 steps. Assume that $^0\mathbf{U} = \mathbf{0}$ and $^0\dot{\mathbf{U}} = \mathbf{0}$ .
|
||
|
||
The first step is to calculate ${}^{0}\ddot{U}$ using the equations in (a) at time 0; i.e., we use
|
||
|
||
$$
|
||
\left[ \begin{array}{l l} 2 & 0 \\ 0 & 1 \end{array} \right] ^ {0} \ddot {\mathbf {U}} + \left[ \begin{array}{c c} 6 & - 2 \\ - 2 & 4 \end{array} \right] \left[ \begin{array}{l} 0 \\ 0 \end{array} \right] = \left[ \begin{array}{l} 0 \\ 1 0 \end{array} \right]
|
||
$$
|
||
|
||
Hence,
|
||
|
||
$$
|
||
^ {0} \ddot {\mathbf {U}} = \left[ \begin{array}{c} 0 \\ 1 0 \end{array} \right]
|
||
$$
|
||
|
||
Now we follow the calculations in Table 9.1.
|
||
|
||
Consider case (1), in which $\Delta t = 0.28$ . We then have (listed to three digits)
|
||
|
||
$$
|
||
a _ {0} = \frac {1}{(0 . 2 8) ^ {2}} = 1 2. 8; \quad a _ {1} = \frac {1}{(2) (0 . 2 8)} = 1. 7 9
|
||
$$
|
||
|
||
$$
|
||
a _ {2} = 2 a _ {0} = 2 5. 5; \quad a _ {3} = \frac {1}{a _ {2}} = 0. 0 3 9 2
|
||
$$
|
||
|
||
Hence, $-\Delta t\mathbf{U} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} - 0.28 \begin{bmatrix} 0 \\ 0 \end{bmatrix} + 0.0392 \begin{bmatrix} 0 \\ 10 \end{bmatrix} = \begin{bmatrix} 0 \\ 0.392 \end{bmatrix}$
|
||
|
||
$$
|
||
\hat {\mathbf {M}} = 1 2. 8 \left[ \begin{array}{l l} 2 & 0 \\ 0 & 1 \end{array} \right] + 1. 7 9 \left[ \begin{array}{l l} 0 & 0 \\ 0 & 0 \end{array} \right]
|
||
$$
|
||
|
||
$$
|
||
= \left[ \begin{array}{c c} 2 5. 5 & 0 \\ 0 & 1 2. 8 \end{array} \right]
|
||
$$
|
||
|
||
The effective loads at time t are
|
||
|
||
$$
|
||
^ {\prime} \hat {\mathbf {R}} = \left[ \begin{array}{l} 0 \\ 1 0 \end{array} \right] + \left[ \begin{array}{c c} 4 5. 0 & 2 \\ 2 & 2 1. 5 \end{array} \right] ^ {\prime} \mathbf {U} - \left[ \begin{array}{c c} 2 5. 5 & 0 \\ 0 & 1 2. 8 \end{array} \right] ^ {\prime - \Delta^ {\prime}} \mathbf {U}
|
||
$$
|
||
|
||
Hence, we need to solve the following equations for each time step,
|
||
|
||
$$
|
||
\left[ \begin{array}{c c} 2 5. 5 & 0 \\ 0 & 1 2. 8 \end{array} \right] ^ {t + \Delta t} \mathbf {U} = ^ {t} \hat {\mathbf {R}} \tag {b}
|
||
$$
|