김경종
25 KiB
- Perform the direct integration analysis. The time step
\Delta tfor this solution should equal about\frac{1}{20}T_{co}, whereT_{co} = 2\pi/\omega_{co}(or be smaller for stability reasons when the central difference method is used).
Note that if a mode superposition solution were used for this step (as described in Section 9.3), then \omega_{co} would be the highest frequency to be included in the solution. Hence p in (9.49) is equal to the number of modes with frequencies smaller than or equal to \omega_{co} .
When analyzing a structural dynamics problem, in most cases, an implicit unconditionally stable time integration is most effective. Then the time step \Delta t need generally be only \frac{1}{20}T_{co} (and not smaller, unless convergence difficulties are encountered in the iteration in a nonlinear response calculation; see Section 9.5.2). If an implicit time integration is employed, it is frequently effective to use higher-order finite elements, for example, the 9- and 27-node elements (see Section 5.3) in two- and three-dimensional analysis, respectively, and a consistent mass idealization. The higher-order elements are effective in the representation of bending behavior, but generally need to be employed with a consistent load vector so that the midside and the corner nodes are subjected to their appropriate load contributions in the analysis.
The observation that the use of higher-order elements can be effective with implicit time integration in the analysis of structural dynamics problems is consistent with the fact that higher-order elements have generally been found to be efficient in static analysis, and structural dynamics problems can be thought of as “static problems including inertia effects.” If, on the other hand, the finite element idealization consists of many elements, it can be more efficient to use explicit time integration with a lumped mass matrix, in which case no effective stiffness matrix is assembled and triangularized but a much smaller time step \Delta t must generally be employed in the solution.
Wave Propagation
The major difference between a structural dynamics problem and a wave propagation problem can be understood to be that in a wave propagation problem a large number of frequencies are excited in the system. It follows that one way to analyze a wave propagation problem is to use a sufficiently high cutoff frequency \omega_{co} to obtain enough solution accuracy. The difficulties are in identifying the cutoff frequency to be used and in establishing a corresponding finite element model.
Instead of using these considerations to obtain an appropriate finite element mesh for the analysis of a wave propagation problem, it is generally more effective to employ the concepts used in finite difference solutions.
If we assume that the critical wavelength to be represented is L_{w} , the total time for this wave to travel past a point is
t _ {w} = \frac {L _ {w}}{c} \tag {9.100}
where c is the wave speed. Assuming that n time steps are necessary to represent the travel of the wave,
\Delta t = \frac {t _ {w}}{n} \tag {9.101}
and the "effective length" of a finite element should be
L _ {e} = c \Delta t \tag {9.102}
This effective length and corresponding time step must be able to represent the complete wave travel accurately and are chosen differently depending on the kind of element idealization and time integration scheme used.
Although a special case, the effectiveness of using (9.102) becomes apparent when a one-dimensional analysis of a bar is performed with a lumped mass idealization and the central difference method. If a uniform bar free at both ends and subjected to a sudden constant step load is idealized as an assemblage of two-node truss elements each of length c \Delta t , the exact wave propagation response is obtained in the solution of the model. It is also interesting to note that the time step \Delta t given in (9.102) then corresponds to the stability limit T_n / \pi derived in Example 9.12, i.e., \omega_n = 2c / L_e , and the nonzero (highest) frequency of a single unconstrained element is \omega_n . Hence, the most accurate solution is obtained by integrating with a time step equal to the stability limit and the solution is less accurate when a smaller time step is employed! This deterioration in the accuracy of the predicted solution when \Delta t is smaller than \Delta t_{cr} is most pronounced when a relatively coarse spatial discretization is used. The scheme is studied, to some extent, in G. Noh and K.J. Bathe [A].
In more complex two- and three-dimensional analyses, the exact solution is generally not obtained, and L_{e} is chosen depending on whether the central difference method or an implicit method is employed for solution.
If the explicit central difference method is used, a lumped mass matrix should be employed, and in this case low-order finite elements in uniform meshes are probably most effective; i.e. the four- and eight-node elements in Figs. 5.4 and 5.5 are frequently employed in two- and three-dimensional analyses, respectively. Using these elements, we construct a mesh as uniform as possible and L_{e} is equal to the smallest distance between any two of the nodes of the mesh employed. This length determines \Delta t as given by (9.102). If higher-order (parabolic or cubic) continuum elements are used, again a mesh as uniform as possible should be constructed with the same measure L_{e} , but the time step has to be further reduced because the interior nodes are “stiffer” than the corner nodes. Also, if structural (beam, plate, or shell) elements are included in the mesh, the time step size \Delta t may be governed by the flexural modes in these elements so that the distances between nodes do not alone determine \Delta t (see Table 9.5). Since the condition is always that \Delta t \leq T_{n}/\pi , where T_{n} is the smallest period of the mesh, we aim to use, for an effective solution, an inexpensively calculated lower bound on T_{n} . This bound is given by the smallest period T_{n}^{(m)} of any element, considered individually, measured over all elements in the mesh, as we show in the following example.
EXAMPLE 9.13: Let \omega_{n} be the largest frequency of an assembled finite element mesh and let \omega_{n}^{(m)} be the largest frequency of element m . Show that
\omega_ {n} \leq \max _ {(m)} \omega_ {n} ^ {(m)} \tag {a}
where \max_{(m)} \omega_{n}^{(m)} is the largest element frequency of all elements in the mesh. Hence, for the central difference method we can then use the time step
\Delta t = \frac {2}{\max _ {(m)} \omega_ {n} ^ {(m)}} \leq \Delta t _ {c r}
Using the Rayleigh quotient (see Section 2.6), we have with \mathbf{K}^{(m)} and \mathbf{M}^{(m)} defined in (4.19) and (4.25),
(\omega_ {n}) ^ {2} = \frac {\phi_ {n} ^ {T} \left(\sum_ {m} \mathbf {K} ^ {(m)}\right) \phi_ {n}}{\phi_ {n} ^ {T} \left(\sum_ {m} \mathbf {M} ^ {(m)}\right) \phi_ {n}} \tag {b}
Let
\mathcal {O} u ^ {(m)} = \boldsymbol {\phi} _ {n} ^ {T} \mathbf {K} ^ {(m)} \boldsymbol {\phi} _ {n} \quad \text { and } \quad \mathcal {I} ^ {(m)} = \boldsymbol {\phi} _ {n} ^ {T} \mathbf {M} ^ {(m)} \boldsymbol {\phi} _ {n}
then (\omega_{n})^{2} = \frac{\sum_{m} \mathcal{U}^{(m)}}{\sum_{m} \mathcal{F}^{(m)}} (c)
Now consider the Rayleigh quotient for a single element,
\rho^ {(m)} = \frac {\boldsymbol {\Phi} _ {n} ^ {T} \mathbf {K} ^ {(m)} \boldsymbol {\Phi} _ {n}}{\boldsymbol {\Phi} _ {n} ^ {T} \mathbf {M} ^ {(m)} \boldsymbol {\Phi} _ {n}} = \frac {\mathcal {U} ^ {(m)}}{\mathcal {I} ^ {(m)}} \tag {d}
Since \mathbf{M}^{(m)} and \mathbf{K}^{(m)} are of the same size as \mathbf{K} , we could theoretically imagine \mathcal{U}^{(m)} and \mathcal{I}^{(m)} to be zero (but not for all m ). However, in any case we have for each element (see Section 2.6)
\mathcal {U} ^ {(m)} \leq (\omega_ {n} ^ {(m)}) ^ {2} \mathcal {I} ^ {(m)}
and therefore from (c),
\begin{array}{l} (\omega_ {n}) ^ {2} \leq \frac {\sum_ {m} (\omega_ {n} ^ {(m)}) ^ {2} \mathcal {F} ^ {(m)}}{\sum_ {m} \mathcal {F} ^ {(m)}} \\ \leq \left[ \max _ {m} \left(\omega_ {n} ^ {(m)}\right) ^ {2} \right] \frac {\sum_ {m} \mathcal {J} ^ {(m)}}{\sum_ {m} \mathcal {J} ^ {(m)}} \\ \end{array}
which proves (a). Note that in (b) we used the \mathbf{K}^{(m)} and \mathbf{M}^{(m)} matrices of element m defined in (4.19) and (4.25), that is with all boundary conditions (and the actions of the other elements) removed. Of course, the same proof is applicable if some elements are constrained at certain degrees of freedom (applied to the assemblage of elements).
For some elements the smallest period can be established exactly in closed form, whereas for more complex (distorted and curved) elements, a lower bound on T_{n}^{(m)} may have to be employed. Table 9.5 summarizes some results.
The condition that \Delta t \leq (“element length”/wave speed) is referred to as the CFL condition after R. Courant, K. Friedrichs, and H. Lewy [A]. We also use the Courant (or CFL) number \Delta t/\Delta t_{cr} to indicate the size of time step actually used in a dynamic solution.
The choice of the effective length L_{e} , and hence time step \Delta t , is considerably simpler if an implicit unconditionally stable time integration method is used. In this case, L_{e} can be chosen according to (9.100) to (9.102) as L_{e} = L_{w}/n in the direction of the wave travel, and then \Delta t follows from (9.102). Nonuniform meshes and low- or high-order elements can be used, and when high-order elements are employed, a consistent mass matrix is usually appropriate. A study using the Bathe method is given in G. Noh, S. Ham and K.J. Bathe [A].
These considerations have been put forward for linear dynamic analysis but are largely also applicable in nonlinear analysis. An important point in nonlinear analysis is that the periods and wave velocities represented in the finite element system change during its
TABLE 9.5 Central difference method critical time steps for some elements:
\Delta t _ {c r} ^ {(m)} = T _ {n} ^ {(m)} / \pi = 2 / \omega_ {n} ^ {(m)}
Two-node truss element:
\begin{array}{l} \mathbf {K} ^ {(m)} = \frac {A E}{L} \left[ \begin{array}{c c} 1 & - 1 \\ - 1 & 1 \end{array} \right]; \quad \mathbf {M} ^ {(m)} = \frac {\rho L}{2} \left[ \begin{array}{c c} 1 & 0 \\ 0 & 1 \end{array} \right] \\ \Delta t _ {c r} ^ {(m)} = \frac {L}{c}; \\ \end{array}
Two-node beam element (see Example 4.1):
\mathbf {K} ^ {(m)} = \frac {E I}{L} \left[ \begin{array}{c c c c} \frac {1 2}{L ^ {2}} & - \frac {6}{L} & - \frac {1 2}{L ^ {2}} & - \frac {6}{L} \\ & 4 & \frac {6}{L} & 2 \\ \text { Sym. } & & \frac {1 2}{L ^ {2}} & \frac {6}{L} \\ & & & 4 \end{array} \right]
\mathbf {M} ^ {(m)} = \frac {\rho A L}{2 4} \left[ \begin{array}{c c c c} 1 2 & 0 & 0 & 0 \\ & L ^ {2} & 0 & 0 \\ \text { Sym. } & 1 2 & 0 \\ & & & L ^ {2} \end{array} \right]
\Delta t _ {c r} ^ {(m)} = \sqrt {\frac {A}{4 8 I}} \frac {L ^ {2}}{c}
Four-node square plane stress element (see Example 4.6):
\mathbf {K} ^ {(m)} = \frac {E t}{1 - \nu^ {2}} \left[ \begin{array}{c c c} \frac {3 - \nu}{6} & & \text { Elements are function of } \nu \\ & \ddots & \\ \text { Sym. } & & \frac {3 - \nu}{6} \end{array} \right]
\begin{array}{l} \mathbf {M} ^ {(m)} = \frac {\rho L ^ {2} t}{4} \left[ \begin{array}{c c c c c} 1 & & & & \\ & 1 & & \text {Zeros} \\ & & \ddots & & \\ & & & 1 & \\ & & & & 1 \end{array} \right] \\ \Delta t _ {c r} ^ {(m)} = \frac {L}{c} \sqrt {1 - \nu} \end{array}
where E = Young's modulus, \nu = Poisson's ratio, L = length (side length) of element, A = cross-sectional area of element, \rho = mass density, I = flexural moment of inertia, t = thickness of plane stress element, c = one-dimensional wave speed = \sqrt{E / \rho}
response. Therefore, the selection of the time step size must take into account that in structural dynamics problems the significantly excited frequencies change magnitudes and that in a wave propagation problem the value of c in (9.102) is not constant.
We discuss additional considerations pertaining to nonlinear analysis in Section 9.5. Let us now present an example demonstrating the above modeling features.
EXAMPLE 9.14: Consider the bar shown in Fig. E9.14(a) initially at rest and subjected to a concentrated end load. The response of the bar at time 0.01 sec is sought.
text_image
ρ = 7800 kg/m³ P(t) E = 2 × 10⁵ MPa x 5 cm u(x, t) Cross-sectional area 4 cm² 95 cm ρ = 1560 kg/m³ E = 4.432 MPa P(t) (MPa) 4 sin ω̂t ω̂ = 150 rad/sec Time 0.0419 sec
Figure E9.14(a) Problem description
Solve this problem using two-node truss elements
(i) With the mode superposition method and
(ii) By direct integration with the trapezoidal rule and the central difference method.
We note that the bar is made of two materials (giving a stiff and a flexible section), and that the choice of truss elements to model the bar implies the use of a one-dimensional mathematical model (see Example 3.17). Of course, in practice, the actual problem and solution may be of much more complex nature, and we use this simplified problem statement and mathematical model merely to demonstrate the modeling and solution procedures discussed above.
To solve this problem we need to first select a discretization that accurately represents a sufficient number of the exact frequencies and corresponding mode shapes. Using lumped and consistent mass representations, the frequencies listed in Table E9.14 are calculated with discretizations using 20 and 40 equal-length elements. In this analysis we note that the frequencies of the lumped mass models are always below the frequencies of the consistent mass models. Because of the relatively stiff short section of the bar at its top, the twentieth frequency in the 20-element models is considerably higher than the nineteenth frequency (and the thirty-ninth and fortieth frequencies of the 40-element models are much higher than the thirty-eighth frequency).
TABLE E9.14 Predicted frequencies (radians/second)
| Frequency number | Lumped mass matrix assumption | Consistent mass matrix assumption | ||
| 20-element model | 40-element model | 20-element model | 40-element model | |
| 1 | 7.02516E + 01 | 7.02648E + 01 | 7.02770E + 01 | 7.02712E + 01 |
| 2 | 2.19037E + 02 | 2.19393E + 02 | 2.19812E + 02 | 2.19587E + 02 |
| 3 | 3.78880E + 02 | 3.80576E + 02 | 3.82932E + 02 | 3.81591E + 02 |
| 4 | 5.43097E + 02 | 5.47936E + 02 | 5.55239E + 02 | 5.50977E + 02 |
| 5 | 7.06967E + 02 | 7.17610E + 02 | 7.34395E + 02 | 7.24482E + 02 |
| 6 | 8.67764E + 02 | 8.87724E + 02 | 9.20054E + 02 | 9.00834E + 02 |
| $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |
| 19 | 2.12481E + 03 | 2.89023E + 03 | 3.65556E + 03 | 3.47009E + 03 |
| 20 | 1.93925E + 05 | 3.01715E + 03 | 3.25207E + 05 | 3.69646E + 03 |
| $\vdots$ | $\vdots$ | $\vdots$ | ||
| 38 | 4.26046E + 03 | 7.36679E + 03 | ||
| 39 | 2.73219E + 05 | 3.36596E + 05 | ||
| 40 | 3.97280E + 05 | 6.76577E + 05 | ||
The frequency of the load application lies between the first and second frequencies of the model. Using 4 \times \omega as the cutoff frequency, we note that it should be sufficient to include the response of four modes in the mode superposition solution, i.e., to use p = 4 in (9.49). However, for instructive purposes we consider the response corresponding to p = 1, 2, \ldots, 5 . We also note that the 20-element models are predicting the significantly excited frequencies to sufficient accuracy (we compare the frequencies of the 20-element models with those of the 40-element models), and we use these models for the response solution.
For the mode superposition solution we use the consistent mass model and obtain at time 0.01 sec the results shown in Figs. E9.14(b) and (c).
Figures E9.14(b) and (c) illustrate how, as we increase the number of modes included in the response prediction, the predicted response converges. The predicted response using four modes is almost the same as using five modes. However, we also note that the static correction [i.e., use of (9.52)] improves the response prediction very significantly when only one, two, or three modes are used in the superposition solution. The modal solutions have been obtained by numerical integration of the decoupled equation (9.46) with a time step \Delta t = 0.0004 (which is about T_{5}/20 ).
For the direct integration with the trapezoidal rule we also use the consistent mass 20-element model and the same time step as in the mode superposition solution. This time step ensures that the response in the modes \phi_{1} to \phi_{5} is accurately integrated. Figure E9.14(d) shows the calculated response and the excellent comparison with the mode superposition solution.
For the central difference method solution we use the 20-element lumped mass model. The time step needs to be sufficiently small for stability. If we use the highest frequency in the model, we obtain \Delta t_{cr} = 2/\omega_{20} = 1.03 \times 10^{-5} sec, and if we use the formula given in Table 9.5, we obtain \Delta t_{cr} \geq \min_{m=1,\ldots,20} \Delta t_{cr}^{(m)} = 0.98 \times 10^{-5} sec. Therefore, in practice, we would use \Delta t = 0.98 \times 10^{-5} sec, but here we can use \Delta t = 1.0 \times 10^{-5} sec. Note that whereas we need
line
| x-displacement (m) | Five modes | Four modes | Three modes | Two modes | One mode |
|---|---|---|---|---|---|
| 0 | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 |
| 2e-4 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 4e-4 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
Figure E9.14(b) Solution using mode superposition without static correction
line
| x-displacement (m) | Five modes | Four modes | Three modes | Two modes | One mode |
|---|---|---|---|---|---|
| 0 | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 |
| 0 | 0.4 | 0.4 | 0.4 | 0.4 | 0.4 |
| 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
Figure E9.14(c) Solution using mode superposition with static correction
line
| x-displacement (m) | Trapezoidal rule (___) | Five modes (___) | CDM (--------) |
|---|---|---|---|
| 4e-4 | 0.0 | 0.0 | 0.0 |
| 3e-4 | 0.1 | 0.1 | 0.1 |
| 2e-4 | 0.2 | 0.2 | 0.2 |
| 1e-4 | 0.3 | 0.3 | 0.3 |
| 0 | 0.5 | 0.5 | 0.5 |
| -1e-4 | 1.0 | 1.0 | 1.0 |
Figure E9.14(d) Comparison of solutions obtained by mode superposition and direct integration
line
| x-displacement (m) | x-coordinate (m) |
|---|---|
| 0 | 0.002 |
| 0 | 0.0004 |
Figure E9.14(e) Comparison of direct integration solutions using the trapezoidal rule
only 25 time steps with the trapezoidal rule (hence the CFL number \doteq 40 ), we require 1000 time steps with the central difference method. Of course, this is due to the relatively stiff top part of the bar. The solution using the central difference method compares very well with the mode superposition and implicit direct integration (by the trapezoidal rule) solutions [see Fig. E9.14(d)].
Finally, let us investigate the solution accuracy if we were to increase the time step size using the trapezoidal rule of time integration. If we consider the third frequency \omega_{3}=382.93 rad/sec, we have \hat{\omega}/\omega_{3}\doteq0.39 . Figure 9.3 shows that for this frequency and the larger frequencies the increase in response above the static response is less than 50 percent. However, the static response in a mode \phi_{i} decreases proportionally with the factor 1/\omega_{i}^{2} and is in any case included in the direct time integration. If we use \Delta t=0.002 sec, we still integrate the third mode dynamic response with about 5 percent accuracy (see Fig. 9.6) (and of course the first and second mode responses more accurately). Hence, the choice of the time step \Delta t=0.002 sec appears reasonable (the CFL number is then approximately 200).
Figure E9.14(e) indeed shows that the solution with the time step \Delta t = 0.002 sec is not far from the solution with the smaller time step. This result corresponds also to the results given in Fig. E9.14(c), where it is seen that the mode superposition solution using three modes is already quite accurate, provided the static correction is included. Hence, the selection and use of the finer time step \Delta t = 0.0004 sec for the trapezoidal rule integration was conservative.
9.4.5 Exercises
9.13. Consider the Bathe method using as the first substep (2 - \sqrt{2})\Delta t . The second substep is then of course of size (\sqrt{2} - 1)\Delta t . Show explicitly that with this substep, the amplitude decay is maximized and the period elongation is minimized in the Bathe method. Assume that the physical damping \xi = 0 .
9.14. Assume that the central difference method is used to solve the dynamic equilibrium equations including proportional damping [that is, we have \xi > 0 in (9.62)]. Show that the critical time step is still given by 2 / \omega_{n} .
9.15. Consider the spectral radii corresponding to the Houbolt, Newmark ( \delta = \frac{1}{2} , \alpha = \frac{1}{4} ), and Bathe methods. Show that the values given in Fig. 9.4 for \Delta t/T = 10{,}000 are correct.
9.16. Calculate the percentage period elongations and amplitude decays of the Houbolt, Newmark ( \delta = \frac{1}{2} , \alpha = \frac{1}{4} ), and Bathe methods corresponding to the initial value problem in (9.99) for the case \Delta t/T = 0.10 . Thus, show that the values given in Fig. 9.6 are correct.
9.17. Use a computer program to solve for the lowest six frequencies of the cantilever shown. Consider three different mathematical models: a Hermitian beam model, a plane stress model, and a fully three-dimensional model. In each case, choose appropriate finite element discretizations and consider the lumped and consistent mass matrix assumptions. Verify that accurate results have been obtained.
text_image
0.01 m 0.04 m 0.60 m
Young's modulus E = 200,000 MPa
Poisson's ratio v = 0.3
Mass density \rho = 7800\mathrm{kg / m^3}
9.18. Use a computer program to solve for the lowest six frequencies of the curved cantilever shown. Use the isoparametric mixed interpolated beam element with two, three, or four nodes and use the lumped or consistent mass matrix. Verify that you have obtained accurate results.
text_image
E = 200,000 v = 0.3 ρ = 7800 h = 0.1 Depth = 0.1 R = 10 h Consider vibrations in the x, y plane only y x
9.19. Consider the problem of a concentrated load P traveling over a simply supported beam. Use a computer program to solve this problem for various values of velocity v.
text_image
vt P h L
E = 200,000 MPa
v = 0.3
\rho = 7800\mathrm{kg / m^3}
L = 10 m
h = 0.2 m
Depth = 0.02 m