baram2584/MultiPhysicsVault
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add raw source
Tests / Hermetic test suite (push) Has been cancelled
Details
Tests / Skill frontmatter validation (push) Has been cancelled
Details

Browse Source
This commit is contained in:
김경종
2026-07-02 09:18:17 +09:00
parent 702f28a758
commit 6bca119e6c
2652 changed files with 54284 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.raw/ConceptsApplicationsFiniteElementAnalysis/ConceptsApplicationsFiniteElementAnalysis_043.md
+492
View File
@@ -0,0 +1,492 @@
<!-- source-page: 421 -->
by the central difference method (Eq. 13.10-5 or Eq. 13.10-12 with [C] = [0]) requires
$$
\Delta t \leq \frac {L}{c} \tag {13.10-16}
$$
which is called the CFL condition after Courant, Friedrichs, and Lewy [13.12,13.50]. The physical interpretation of this condition is that $\Delta t$ must be small enough that information does not propagate across more than one element per time step. This observation is true only of linear-displacement elements with lumped mass. Note that if a consistent mass matrix is used, the example of Section 13.5 gives $(\omega_{\mathrm{max}})_{e} = 2\sqrt{3} (c/L)$ so that $\Delta t \leq (L/c)/\sqrt{3}$ ; this is more restrictive than the time step criterion for an element with lumped masses, Eq. 13.10-16. This result is typical. Thus, in addition to providing uncoupled equations and generally more accurate results than consistent mass matrices in explicit integration, lumped mass matrices provide for larger stable time steps.
Higher-order elements yield higher maximum frequencies than lower-order elements. For this reason and for reasons discussed in Section 13.14, one may wish to avoid higher-order elements when doing explicit integration. Similarly, one should avoid penalty constraints, as large penalty numbers make $\omega_{max}$ very large.
For plane and solid finite elements, it is usually difficult to analytically calculate $(\omega_{\mathrm{max}})_{e}$ . Rather, it is possible to use Eq. 13.10-16 with L replaced by an effective element diameter, $L_{e}$ . Effective element diameters are shown in Fig. 13.10-1 for some low-order displacement field elements. For higher-order displacement field elements, effective diameters are difficult to calculate. As an alternative to calculation, the highest frequency of an element (or of a global system of equations) can be accurately bounded by Gerschgorin's theorem [13.12], which, for lumped mass matrices, states that
$$
\omega_ {\max} ^ {2} \leq \max \left(k _ {i i} + \sum_ {\substack {j = 1 \\ j \neq i}} ^ {n _ {e}} | k _ {i j} |\right) / m _ {i i} \quad \text {where} \quad i = 1, 2, \dots , n _ {e} \tag{13.10 - 17}
$$
Here $n_{e}$ is the number of d.o.f. per element. Equation 13.10-17 can be applied to the assembled structure if $k_{ii}$ , $k_{ij}$ , and $m_{ii}$ are replaced by structure coefficients $K_{ii}$ , $K_{ij}$ , and $M_{ii}$ , and $n_{e}$ becomes $n_{eq}$ , the number of global d.o.f. Equation
![](images/page-421_827c1b238433158a93ff68c5ff048a2670fbfd60f253916c5f9e1548ffb6bf7f.jpg)
<details>
<summary>text_image</summary>
√2L
L
L
L_e = √2L
</details>
![](images/page-421_1c1ca18d94e2b621223b4edeaa662ff8172489a003bea1ee5238424debdcc5ca.jpg)
<details>
<summary>text_image</summary>
L
L
L
Lₑ = L
</details>
![](images/page-421_6b341949dc98a3e4db254daa47b3c17b117dcca20f20133bc24ec2c587d43fc8.jpg)
<details>
<summary>text_image</summary>
h₁
h₂
Lₑ = min(h₁, h₂)
</details>
Figure 13.10-1. Effective element diameters for some common two-dimensional elements. Results for three-node elements are exact. Results for the four-node element are exact if the element is rectangular ( $h_{1} = h_{2}$ ) and conservative otherwise. Inscribed circles touch at least three sides of the element.
<!-- source-page: 422 -->
13.10-17 fails if $m_{ii} = 0$ (or if $M_{ii} = 0$ ), as may happen with some mass allocation schemes. The difficulty of $m_{ii} = 0$ (but $M_{ii} \neq 0$ ) is addressed in [13.51].
A useful definition of the relative size of time step used in direct integration is the Courant number, $C_{n}$ , defined as
$$
C _ {n} = \frac {\Delta t _ {\text { actual }}}{\Delta t _ {\text { stable }}} \tag {13.10-18}
$$
where $\Delta t_{actual}$ is the actual time step used and $\Delta t_{stable}$ is the largest time step permitted for stable explicit integration by a specified explicit method (e.g., Eqs. 13.10-8 or 13.10-14). For reasons discussed in Section 13.13, we recommend a maximum $C_{n}$ of about 0.95 to 0.98.
Usually $\Delta t_{stable}$ in Eq. 13.10-18 is not precisely known because $\omega_{max}$ for the finite element mesh is not precisely known. One can determine a precise value of $\omega_{max}$ only by a somewhat inconvenient eigenvalue calculation. Instead, practitioners obtain a close upper bound to $\Delta t_{stable}$ by approximating $\omega_{max}$ by the largest unconstrained frequency among all elements, $(\omega_{\mathrm{max}})_{e}$ , as obtained perhaps from Eq. 13.10-15 or (as an upper bound) from Eq. 13.10-17.
Example: Wave Propagation. Consider the uniform steel bar and tip loading shown in Fig. 13.10-2. The bar is undamped, initially at rest, and is modeled by 40 equal-length linear-displacement elements, so that $L = L_{T}/40 = 0.5$ in. The highest element frequency is given by Eq. 13.10-15 as
$$
(\omega_ {\max}) _ {e} = \frac {2}{0 . 5} \sqrt {\frac {E}{\rho}} = 8. 0 5 3 9 (1 0 ^ {5}) \mathrm{rad/sec} \tag {13.10-19}
$$
This value is very close to the maximum frequency of the finite element mesh, which is $\omega_{\max}=8.0523(10^{5})$ rad/sec. Therefore, stable integration by the central difference method requires, according to Eq. 13.10-14 with $\xi=0$ , that $\Delta t\leq2/(\omega_{\max})_{e}=2.483(10^{-6})$ sec. The computational procedure of Table 13.10-1 is implemented in the Fortran program shown in Fig. 13.10-3.
The stress time-history response at the midpoint of element 20 (at x = 9.75 in.) is shown in Fig. 13.10-4 for $\Delta t = 2.4(10^{-6})$ sec ( $C_{n} = 0.966$ ) and an analysis duration equivalent to two wave traversals along the bar (83 time steps). The exact solution is also shown. Many important features of numerical solution of wave propagation problems are displayed. One is that stress at x = 9.75 in. is zero until sufficient time has elapsed for the stress wave to propagate from the instantaneously loaded tip to the bar midpoint. This is a characteristic feature of a class of partial differential equations called hyperbolic, among which the equation of motion is an example. The increase in mean compressive stress from 100 to 200 psi at 0.15 msec is due to wave reflection from the built-in end.
![](images/page-422_fa3d81f96a830ba5e0babf071e4ce4932a05c1f59b5d09a99cf0ed3814c5403b.jpg)
<details>
<summary>text_image</summary>
P(t)
x
L_T
</details>
![](images/page-422_0a4e0677a736dfb01412f2af75bd224070da45738fde5c5c6ec7e3c9d2d04659.jpg)
<details>
<summary>line</summary>
| t | P(t), lb |
|---|---|
| 0 | 100 |
| 1 | 100 |
| 2 | 100 |
| 3 | 100 |
| 4 | 100 |
| 5 | 100 |
| 6 | 100 |
| 7 | 100 |
| 8 | 100 |
| 9 | 100 |
| 10 | 100 |
| 11 | 100 |
| 12 | 100 |
| 13 | 100 |
| 14 | 100 |
| 15 | 100 |
| 16 | 100 |
| 17 | 100 |
| 18 | 100 |
| 19 | 100 |
| 20 | 100 |
| 21 | 100 |
| 22 | 100 |
| 23 | 100 |
| 24 | 100 |
| 25 | 100 |
| 26 | 100 |
| 27 | 100 |
| 28 | 100 |
| 29 | 100 |
| 30 | 100 |
| 31 | 100 |
| 32 | 100 |
| 33 | 100 |
| 34 | 100 |
| 35 | 100 |
| 36 | 100 |
| 37 | 100 |
| 38 | 100 |
| 39 | 100 |
| 40 | 100 |
| 41 | 100 |
| 42 | 100 |
| 43 | 100 |
| 44 | 100 |
| 45 | 100 |
| 46 | 100 |
| 47 | 100 |
| 48 | 100 |
| 49 | 100 |
| 50 | 100 |
| 51 | 100 |
| 52 | 100 |
| 53 | 100 |
| 54 | 100 |
| 55 | 100 |
| 56 | 100 |
| 57 | 100 |
| 58 | 100 |
| 59 | 100 |
| 60 | 100 |
| 61 | 100 |
| 62 | 100 |
| 63 | 100 |
| 64 | 100 |
| 65 | 100 |
| 66 | 100 |
| 67 | 100 |
| 68 | 100 |
| 69 | 100 |
| 70 | 100 |
| 71 | 100 |
| 72 | 100 |
| 73 | 100 |
| 74 | 100 |
| 75 | 100 |
| 76 | 100 |
| 77 | 100 |
| 78 | 100 |
| 79 | 100 |
| 80 | 100 |
| 81 | 100 |
| 82 | 100 |
| 83 | 100 |
| 84 | 100 |
| 85 | 100 |
| 86 | 100 |
| 87 | 100 |
| 88 | 100 |
| 89 | 100 |
| 90 | 100 |
| 91 | 100 |
| 92 | 100 |
| 93 | 100 |
| 94 | 100 |
| 95 | 100 |
| 96 | 100 |
| 97 | 100 |
| 98 | 100 |
| 99 | 100 |
| 100 | 100 |
</details>
Figure 13.10-2. One-dimensional uniform bar with instantaneous tip loading. The bar is initially at rest. $A = 1.0 \text{ in.}^{2}$ , $E = 30(10^{6}) \text{ psi}$ , $\rho = 7.4(10^{-4}) \text{ lb-sec}^{2}/\text{in.}^{4}$ , $L_{T} = 20 \text{ in.}$ Load $P_{0} = 100 \text{ lb}$ is applied at $t = 0$ .
<!-- source-page: 423 -->
```prolog
C---- Program for 1-dimensional wave propagation through NELE equal-
C---- length finite elements using the central difference method
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION X(101),RM(101),D(101),V(101),FEXT(101),FINT(101)
DATA X,RM,FEXT/101*0.,101*0.,101*0./
C
C NELE = number of elements ELELEN = length of element
C CSA = cross-sectional area DENSTY = mass density
C E = elastic modulus DELT = time step
C NSTEP = number of time steps IFREQ = output interval
C
C---- Input mesh data
READ(5,*)NELE,ELELEN,CSA,DENSTY,E,DELT,NSTEP,IFREQ
NUMNOD=NELE+1
C---- Generate nodal coordinates and form lumped mass matrix
DO 10 I=1,NUMNOD
10 X(I)=FLOAT(I-1)*ELELEN
DO 20 K=1,NELE
RM(K)=RM(K)+DENSTY*CSA*ELELEN/2.
20 RM(K+1)=RM(K+1)+DENSTY*CSA*ELELEN/2.
C---- Set at-rest initial conditions and tip-load force
DO 30 I=1,NUMNOD
D(I)=0.
30 V(I)=0.
FEXT(1)=100.
C---- Integration loop
DO 100 N=1,NSTEP
C---- Get internal force for last time step
CALL INTFOR(D,FINT,X,E,CSA,NELE)
C---- Update displacements using the central difference method and
C---- enforce zero-displacement boundary condition at built-in end
DO 40 I=1,NUMNOD
DOLD=D(I)
D(I)=DELT**2*(FEXT(I)-FINT(I))/RM(I)+D(I)+DELT*V(I)
IF(I.EQ,NUMNOD) D(I)=0.
40 V(I)=(D(I)-DOLD)/DELT
C---- Output if desired
SIG20=E*((D(21)-D(20))/ELELEN)
IF(MOD(N,IFREQ).EQ.0)
WRITE(7,1000)N,FLOAT(N)*DELT,D(20),V(20),SIG20
100 CONTINUE
STOP
1000 FORMAT(I5,4E15.4)
END
SUBROUTINE INTFOR(D,FINT,X,E,CSA,NELE)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION D(1),FINT(1),X(1)
NUMNOD=NELE+1
C---- Zero internal force vector
DO 10 I=1,NUMNOD
10 FINT(I)=0.
C---- Loop over elements
DO 20 K=1,NELE
RL=X(K+1)-X(K)
C---- Compute strain and stress
STRAIN=(D(K+1)-D(K))/RL
STRESS=STRAIN*E
C---- Assemble contribution into internal force vector
F1=-STRESS*CSA
F2=+STRESS*CSA
FINT(K)=FINT(K)+F1
FINT(K+1)=FINT(K+1)+F2
20 CONTINUE
RETURN
END
```
Figure 13.10-3. Fortran program for direct integration by the central-difference method using the procedure in Table 13.10-1.
<!-- source-page: 424 -->
![](images/page-424_f7aac51d8ac8b54fc3828a460d567d9ba2f7c306862d302754b6505133d4f11e.jpg)
<details>
<summary>line</summary>
| Time (milliseconds) | Exact (psi) | Central difference (psi) |
| ------------------- | ----------- | ------------------------ |
| 0.00 | 0 | 0 |
| 0.04 | -100 | -100000 |
| 0.06 | -100 | -100000 |
| 0.08 | -100 | -100000 |
| 0.10 | -100 | -100000 |
| 0.12 | -100 | -100000 |
| 0.14 | -100 | -100000 |
| 0.16 | -200 | -200000 |
| 0.18 | -200 | -200000 |
| 0.20 | -200 | -200000 |
</details>
Figure 13.10-4. Stress time history at x = 9.75 in. for a 40-element model of the 20-in. bar shown in Fig. 13.10-2 using $\Delta t = 2.4(10^{-6})$ sec ( $C_{n} = 0.966$ ). Inset shows instability that results from taking $\Delta t$ too large ( $C_{n} = 1.007$ ).
Figure 13.10-4 also shows severe spurious oscillations, or noise, created by the algorithm, as the stress wave passes. Spurious oscillations are high-frequency motions that the mesh cannot accurately resolve and are similar to the Gibbs phenomenon in fitting a finite Fourier series to a piecewise-continuous function. These spurious oscillations are often a nuisance. Stiffness-proportional damping is sometimes employed to attenuate noise although it is usually better practice to employ a different integration method, such as discussed in Section 13.12, which automatically dissipates high-frequency motion.
To demonstrate the instability in explicit integration that results when the maximum stable time step is exceeded, we repeated the analysis using $\Delta t = 2.5(10^{-6})$ sec ( $C_n = 1.007$ ). The results, shown in the inset of Fig. 13.10-4, are typical of instability in linear analysis and demonstrate a solution that increases wildly with each time step until eventually the computer program aborts.
When $\Delta t$ is slightly less than the stability limit, the numerical solution may display spurious “beating” in which the amplitude of response repeatedly grows and decays. To illustrate beating, we repeated the foregoing analysis using $\Delta t = 2.481(10^{-6})$ sec ( $C_{n} = 0.999$ ). Results are shown in Fig. 13.10-5. We see that average results are good but spurious oscillations are severe. For some problems, particularly problems having few d.o.f., beating can be more detrimental to accuracy than shown by this example.
Concluding Remarks. When an explicit method is used for a single ordinary differential equation, accuracy is markedly time-step-size-dependent. It may therefore appear that the stability criterion is of only academic interest since $\Delta t$ must be considerably smaller to achieve satisfactory accuracy. However, this supposition is not true of systems of finite element equations (about 20 equations or more) in which it is typically observed that excellent accuracy can be obtained using a time step size just under the stability limit. The reason is that the equations of motion are a stiff system of ordinary differential equations, which has a broad frequency spectrum. The stability criterion, Eq. 13.10-8, is based upon the highest frequency, or shortest time scale phenomenon that the mesh can possibly repro-
<!-- source-page: 425 -->
![](images/page-425_bbe6f35acb3e5b722f36c8eb2365d694e288c7864bdb4cb03c9eec5c4441d417.jpg)
<details>
<summary>line</summary>
| Time (milliseconds) | Exact | Central difference, C_N = 0.999 |
| ------------------- | ----- | ------------------------------- |
| 0.00 | -100 | -100 |
| 0.02 | -100 | -100 |
| 0.04 | -100 | -100 |
| 0.06 | -100 | -100 |
| 0.08 | -100 | -100 |
| 0.10 | -100 | -100 |
| 0.12 | -100 | -100 |
| 0.14 | -100 | -100 |
| 0.16 | -200 | -200 |
| 0.18 | -200 | -200 |
| 0.20 | -200 | -200 |
</details>
Figure 13.10-5. Stress time history at x = 9.75 in. for a 40-element model of the bar shown in Fig. 13.10-2 using $\Delta t = 2.481(10^{-6})$ sec ( $C_{n} = 0.999$ ). The response demonstrates “beating” from using a time step close to the stability limit.
duce. It is true that motions of the mesh that occur on this time scale are not accurately resolved by taking $\Delta t$ close to the stability limit. Fortunately, these motions contribute little to the structural response, which is dominated by much lower-frequency, long-time-scale phenomena that are accurately resolved. In other words, we do not require high-frequency phenomena to be accurately resolved; all that we ask is that they be stably integrated. The implication is that in explicit methods, a $\Delta t$ satisfying stability criteria is usually satisfactory to guarantee accuracy. Even so, the maximum allowable $\Delta t$ is often smaller than one would like, because so many time steps may be needed to span the duration of an analysis.
# 13.11 IMPLICIT DIRECT INTEGRATION METHODS
Most of the useful implicit methods are unconditionally stable and have no restriction on the time step size other than as required for accuracy. A popular unconditionally stable implicit method is called the trapezoidal rule or the average acceleration method. (It is also known as the Crank–Nicolson method when applied to parabolic partial differential equations, such as the heat conduction equation.) The trapezoidal rule relates displacements, velocities, and accelerations by
$$
\{\mathbf {D} \} _ {n + 1} = \{\mathbf {D} \} _ {n} + \frac {\Delta t}{2} \left(\{\dot {\mathbf {D}} \} _ {n} + \{\dot {\mathbf {D}} \} _ {n + 1}\right) \tag {13.11-1}
$$
$$
\{\dot {\mathbf {D}} \} _ {n + 1} = \{\dot {\mathbf {D}} \} _ {n} + \frac {\Delta t}{2} \left(\{\ddot {\mathbf {D}} \} _ {n} + \{\ddot {\mathbf {D}} \} _ {n + 1}\right) \tag {13.11-2}
$$
<!-- source-page: 426 -->
TABLE 13.11-1. COMPUTATIONAL PROCEDURE FOR DIRECT INTEGRATION BY THE TRAPEZOIDAL-RULE METHOD.
1. Form [K], [C], and [M].
2. Set initial conditions $\{\mathbf{D}\}_{0} = \{\mathbf{D}(t = 0)\}$ and $\{\dot{\mathbf{D}}\}_{0} = \{\dot{\mathbf{D}}(t = 0)\}$ ; use Eq. 13.9-1 to compute $\{\ddot{\mathbf{D}}\}_{0} = [\mathbf{M}]^{-1}(\{\mathbf{R}^{\mathrm{ext}}\}_{0} - [\mathbf{C}]\{\dot{\mathbf{D}}\}_{0} - [\mathbf{K}]\{\mathbf{D}\}_{0})$ ; $n = 0$ .
3. Form effective stiffness matrix, Eq. 13.11-6, and factor it.
4. Form effective load vector, Eq. 13.11-7.
5. Solve $\{\mathbf{K}^{\mathrm{eff}}\} \{\mathbf{D}\}_{n + 1} = \{\mathbf{R}^{\mathrm{eff}}\}_{n + 1}$ for $\{\mathbf{D}\}_{n + 1}$ by forward- and back-substitution.
6. Update velocity $\{\dot{\mathbf{D}}\}_{n+1}$ and acceleration $\{\ddot{\mathbf{D}}\}_{n+1}$ by Eqs. 13.11-3 and 13.11-4.
7. Output if desired; $n \leftarrow n + 1$ , go to Step 4.
Equations 13.11-1 and 13.11-2 can be obtained using Taylor series (this is left as an exercise) in which it is seen that they are second-order accurate. Alternatively, Eqs. 13.11-1 and 13.11-2 can be solved for $\{\dot{D}\}_{n+1}$ and $\{\ddot{D}\}_{n+1}$ to provide
$$
\{\dot {\mathbf {D}} \} _ {n + 1} = \frac {2}{\Delta t} \left(\{\mathbf {D} \} _ {n + 1} - \{\mathbf {D} \} _ {n}\right) - \{\dot {\mathbf {D}} \} _ {n} \tag {13.11-3}
$$
$$
\{\ddot {\mathbf {D}} \} _ {n + 1} = \frac {4}{\Delta t ^ {2}} \left(\{\mathbf {D} \} _ {n + 1} - \{\mathbf {D} \} _ {n}\right) - \frac {4}{\Delta t} \{\dot {\mathbf {D}} \} _ {n} - \{\ddot {\mathbf {D}} \} _ {n} \tag {13.11-4}
$$
Combination of Eqs. 13.11-3 and 13.11-4 with the equation of motion, Eq. 13.9-1 at time $(n + 1)\Delta t$ , yields
$$
[ \mathbf {K} ^ {\text { eff }} ] \{\mathbf {D} \} _ {n + 1} = \{\mathbf {R} ^ {\text { eff }} \} _ {n + 1} \tag {13.11-5}
$$
where the effective stiffness matrix and effective load vector are, respectively,
$$
[ \mathbf {K} ^ {\mathrm{eff}} ] = \frac {4}{\Delta t ^ {2}} [ \mathbf {M} ] + \frac {2}{\Delta t} [ \mathbf {C} ] + [ \mathbf {K} ] \tag {13.11-6}
$$
$$
\begin{array}{l} \{\mathbf {R} ^ {\text { eff }} \} _ {n + 1} = \{\mathbf {R} ^ {\text { ext }} \} _ {n + 1} + [ \mathbf {M} ] \left(\frac {4}{\Delta t ^ {2}} \{\mathbf {D} \} _ {n} + \frac {4}{\Delta t} \{\dot {\mathbf {D}} \} _ {n} + \{\ddot {\mathbf {D}} \} _ {n}\right) \\ + [ \mathbf {C} ] \left(\frac {2}{\Delta t} \{\mathbf {D} \} _ {n} + \{\dot {\mathbf {D}} \} _ {n}\right) \tag {13.11-7} \\ \end{array}
$$
A flow chart for this algorithm is given in Table 13.11-1.
# Remarks.
1. Equation 13.11-5 is a system of coupled linear algebraic equations even if [M] and [C] are diagonal. For linear problems, $[K^{eff}]$ need be formed and factored only once. After the initial expense of factorization, time stepping can be performed for only the cost of forward- and back-substitution.
2. If [M] is positive definite, then $[\mathbf{K}^{\mathrm{eff}}]$ is nonsingular even if [K] permits rigid-body displacements.
3. The method is easily started from initial conditions $\{\mathbf{D}\}_{0}$ and $\{\dot{\mathbf{D}}\}_{0}$ and Eq. 13.10-7.
<!-- source-page: 427 -->
4. For problems having material nonlinearity, [K] and hence $[K^{eff}]$ are functions of $\{D\}_{n+1}$ , and possibly its time derivative, which are unknowns. Accordingly, [K] must be predicted using an estimate for $\{D\}_{n+1}$ . Equation 13.11-5 is then solved for an improved $\{D\}_{n+1}$ ; hence, the prediction of [K] is improved, and so on. For severe nonlinearity, convergence may be difficult and expensive.
5. In terms of computational efficiency, Eq. 13.11-6 shows that there is no merit to using a lumped mass matrix. In Section 13.13 we show that implicit integration with consistent mass matrices is usually more accurate than with lumped mass matrices.
Choice of Time Step $\Delta t$ . Since the method is numerically stable for any $\Delta t$ , time step selection is based on accuracy considerations alone. Compared with explicit methods, the per-time-step cost of an implicit method is high. Thus implicit methods are economically attractive only when $\Delta t$ can be much larger than would be used in an explicit method. Unconditional stability (which emphatically does not imply unconditional accuracy) coupled with the economic need for large $\Delta t$ tempts many analysts into using time steps that are too large. To select a time step that will provide accurate results, one must identify the highest frequency of interest in the loading or response of a structure. Let this frequency be called $\omega_{u}$ . As an approximation, structure modes with frequency higher than about $3\omega_{u}$ participate quasistatically in the response while modes with frequency lower than $3\omega_{u}$ also participate dynamically. With second-order accurate time integration methods (most popular methods are second-order accurate), a minimum of 20 time steps per period of $\omega_{u}$ should provide very good accuracy for modes that participate dynamically in the response; that is, use $\Delta t < (2\pi/\omega_{u})/20 \approx 0.3/\omega_{u}$ , unless a smaller $\Delta t$ is required because of convergence difficulties in nonlinear analysis. As an additional assurance of an accurate solution, analysis should be repeated using a smaller time step than used in the first analysis. In linear problems of a type well suited to implicit integration, a Courant number $^{3}$ $C_{n}$ of about 20 is typical although occasionally it can be as high as 100 and still yield accurate results.
# 13.12 OTHER IMPLICIT AND EXPLICIT METHODS. MIXED METHODS
The central-difference and trapezoidal-rule methods do not provide automatic dissipation of high-frequency numerical noise, as is sometimes desirable. In what follows, we discuss some methods that often have dissipation (also called artificial viscosity or numerical damping).
The Houbolt method [13.55] is obtained by cubic Lagrange interpolation of $\{D\}$ at times $(n - 2)\Delta t$ through $(n + 1)\Delta t$ . Exact differentiation of the interpolant yields
$$
\{\dot {\mathbf {D}} \} _ {n + 1} = \frac {1}{6 \Delta t} (1 1 \{\mathbf {D} \} _ {n + 1} - 1 8 \{\mathbf {D} \} _ {n} + 9 \{\mathbf {D} \} _ {n - 1} - 2 \{\mathbf {D} \} _ {n - 2}) \tag {13.12-1}
$$
$^{3}$ The definition of $C_{n}$ , Eq. 13.10-18, employs $\Delta t_{stable}$ for explicit integration. This definition remains useful in the present context even though most implicit methods are stable for any time step.
<!-- source-page: 428 -->
$$
\{\ddot {\mathbf {D}} \} _ {n + 1} = \frac {1}{\Delta t ^ {2}} (2 \{\mathbf {D} \} _ {n + 1} - 5 \{\mathbf {D} \} _ {n} + 4 \{\mathbf {D} \} _ {n - 1} - \{\mathbf {D} \} _ {n - 2}) \tag {13.12-2}
$$
This method is implicit and unconditionally stable but provides artificial damping that is too high for low-frequency response. The Houbolt method was once common in general-purpose transient codes but has been supplanted by methods with better algorithmic damping properties and now is more of historical interest.
The Newmark family of methods [13.56] is very popular and is given by
$$
\{\mathbf {D} \} _ {n + 1} = \{\mathbf {D} \} _ {n} + \Delta t \{\dot {\mathbf {D}} \} _ {n} + \frac {\Delta t ^ {2}}{2} [ (1 - 2 \beta) \{\ddot {\mathbf {D}} \} _ {n} + 2 \beta \{\ddot {\mathbf {D}} \} _ {n + 1} ] \tag {13.12-3}
$$
$$
\{\dot {\mathbf {D}} \} _ {n + 1} = \{\dot {\mathbf {D}} \} _ {n} + \Delta t [ (1 - \gamma) \{\ddot {\mathbf {D}} \} _ {n} + \gamma \{\ddot {\mathbf {D}} \} _ {n + 1} ] \tag {13.12-4}
$$
where $\beta$ and $\gamma$ are chosen by the analyst to control stability and accuracy. Substitution of Eqs. 13.12-3 and 13.12-4 into Eq. 13.9-1 at time $(n + 1)\Delta t$ yields equations similar to Eq. 13.10-5 for explicit Newmark methods ( $\beta = 0$ ) and to Eqs. 13.11-5, 13.11-6, and 13.11-7 for implicit Newmark methods ( $\beta > 0$ ). It can be shown that the stability of this algorithm is [13.49]:
unconditional stability when
$$
2 \beta \geq \gamma \geq \frac {1}{2} \tag {13.12-5}
$$
conditional stability when
$$
\gamma \geq \frac {1}{2}, \quad \beta < \frac {1}{2}, \quad \text { and } \quad \Delta t \leq \frac {\xi (\gamma - \frac {1}{2}) + \sqrt {\gamma / 2 - \beta + \xi^ {2} (\gamma - \frac {1}{2}) ^ {2}}}{\omega_ {\max} (\gamma / 2 - \beta)} \tag {13.12-6}
$$
The method is unstable for $\gamma < \frac{1}{2}$ . (As a special case, we note that Eq. 13.12-6 yields infinite $\Delta t$ when $\gamma = \frac{1}{2}$ and $\beta = \frac{1}{4}$ —that is, the trapezoidal rule without damping.) According to our definitions of implicit and explicit methods, Newmark's method is implicit unless $\gamma = \beta = 0$ , which is unstable for any $\Delta t$ and therefore cannot be used. However, most analysts refer to Newmark's method with $\beta = 0$ and $\gamma \geq \frac{1}{2}$ as being explicit in which it is noted that, for practical purposes, [C] must be null or diagonal to avoid the solution of simultaneous equations. An explicit version of the Newmark method, called a predictor-corrector algorithm, that permits nondiagonal [C] is described in [13.61]. Implementation of the explicit and implicit Newmark methods are essentially the same as the computational procedures given in Tables 13.10-1 and 13.11-1, respectively.
A variety of useful techniques obtained from the Newmark family is listed in Table 13.12-1. When $\gamma = \frac{1}{2}$ , the methods have no algorithmic damping and are second-order accurate. An exception is the Fox-Goodwin method with $[C] = [0]$ , which is fourth-order accurate. Taking $\gamma > \frac{1}{2}$ introduces artificial damping, but also reduces the accuracy of the Newmark methods to first order. For implicit Newmark methods, taking
$$
\beta = \frac {1}{4} (\gamma + \frac {1}{2}) ^ {2} \tag {13.12-7}
$$
maximizes the high-frequency dissipation for a given value of $\gamma > \frac{1}{2}$ [13.46].
<!-- source-page: 429 -->
TABLE 13.12-1. SUMMARY OF NEWMARK METHODS: $u =$ UNDAMPED, $d =$ DAMPED. $\Omega_{\mathrm{crit}} = \omega_{\max}\Delta t_{\max}$ . STABILITY REQUIRES $\Delta t\leq \Omega_{\mathrm{crit}} / \omega_{\max}$ .
<table><tr><td></td><td>Method</td><td> $\beta$ </td><td> $\gamma$ </td><td> $\Omega_{\text{crit}}$ </td><td>Accuracy</td></tr><tr><td rowspan="4">Implicit</td><td>Artificially damped</td><td> $>\gamma/2$ </td><td> $>1/2$ </td><td> $\infty^{(u,d)}$ </td><td> $O(\Delta t)$ </td></tr><tr><td>Average acceleration (trapezoidal rule)</td><td>1/4</td><td>1/2</td><td> $\infty^{(u,d)}$ </td><td> $O(\Delta t^{2})$ </td></tr><tr><td>Linear acceleration</td><td>1/6</td><td>1/2</td><td> $2\sqrt{3} \approx 3.464^{(u)}$ Eq. (A) $^{(d)}$ </td><td> $O(\Delta t^{2})$ </td></tr><tr><td>Fox-Goodwin (royal road)</td><td>1/12</td><td>1/2</td><td> $\sqrt{6} \approx 2.449^{(u)}$ Eq. (A) $^{(d)}$ </td><td> $O(\Delta t^{4})^{(u)}$ $O(\Delta t^{2})^{(d)}$ </td></tr><tr><td rowspan="2">Explicit</td><td>Central difference [M], [C] diagonal</td><td>0</td><td>1/2</td><td> $2^{(u)}$ Eq. (A) $^{(d)}$ </td><td> $O(\Delta t^{2})$ </td></tr><tr><td>Artificially damped [M], [C] diagonal</td><td>0</td><td> $>1/2$ </td><td>Eq. (A) $^{(u,d)}$ </td><td> $O(\Delta t)$ </td></tr></table>$$
\Omega_ {\mathrm{crit}} = \frac {\xi (\gamma - \frac {1}{2}) + \sqrt {\gamma / 2 - \beta + \xi^ {2} (\gamma - \frac {1}{2}) ^ {2}}}{\gamma / 2 - \beta} \tag {A}
$$
Interestingly, the presence of damping in the explicit Newmark method raises the stability limit. This is in contrast to Eqs. 13.10-8 and 13.10-14 for other forms of the central-difference method in which no change and a decrease in stability limit, respectively, are observed. Thus, if the fraction $\xi$ of critical damping at $\omega_{\mathrm{max}}$ is not known, a conservative $\Delta t$ is obtained by taking $\xi = 0$ in Eq. 13.12-6. However, practical necessity dictates that [C] be diagonal in the explicit Newmark method, thus the lagged central-difference algorithm, Eq. 13.10-12 (or a predictor-corrector algorithm [13.61]) is preferable for problems in which physically realistic spectral damping is to be modeled.
The Newmark linear-acceleration method appears to be ideal for problems such as earthquake shaking response analysis in which piecewise linear-acceleration records are typically used as excitation. Unfortunately, this implicit method is only conditionally stable. The Wilson- $\theta$ method [13.18,13.28,13.46,13.57] is also a linear-acceleration method, but is unconditionally stable.
A disadvantage of the Newmark methods is that algorithmic damping can only be obtained at the expense of reduced accuracy. The $\alpha$ -method, proposed by Hilber, Hughes, and Taylor [13.58], does not have this weakness, and with appropriate choice of parameters retains second-order accuracy and provides effective high-frequency dissipation. The method uses the Newmark formulas, Eqs. 13.12-3 and 13.12-4, with the modified equation of motion
$$
\begin{array}{l} [ \mathbf {M} ] \{\ddot {\mathbf {D}} \} _ {n + 1} + (1 + \alpha) [ \mathbf {C} ] \{\dot {\mathbf {D}} \} _ {n + 1} - \alpha [ \mathbf {C} ] \{\dot {\mathbf {D}} \} _ {n} + (1 + \alpha) [ \mathbf {K} ] \{\mathbf {D} \} _ {n + 1} \\ - \alpha [ \mathbf {K} ] \{\mathbf {D} \} _ {n} = (1 + \alpha) \{\mathbf {R} _ {n + 1} ^ {\mathrm{ext}} - \alpha \{\mathbf {R} _ {n} ^ {\mathrm{ext}} \} \tag {13.12-8} \\ \end{array}
$$
If the parameters are selected so that $-\frac{1}{3} \leq \alpha \leq 0$ , $\gamma = (1 - 2\alpha)/2$ , and $\beta = (1 - \alpha)^{2}/4$ , the method is implicit, unconditionally stable, and second-order accurate [13.46]. When these guidelines are used, with $\alpha = 0$ , the method reduces to the trapezoidal rule, which has no dissipation. Decreasing $\alpha$ increases the amount of numerical damping.
<!-- source-page: 430 -->
Mixed Methods. A current trend in time integration analysis is to create algorithms that combine explicit and implicit methods, so as to capitalize on the strong points of each. Such schemes are called mixed integration methods. Belystschko and Mullen [13.59,13.60] developed a method that uses a nodal partition to separate nodes into implicit and explicit groups. D.o.f. in the explicit and implicit groups are then integrated by explicit and implicit methods respectively. The motivation for such an approach is that structures often contain spatial subdomains having markedly different time scales. For example, in fluid–structure interaction problems, the time scales associated with the fluid are usually much longer than those associated with the structure. By using an implicit method for the structure and an explicit method for the fluid, we exploit the strong points of each integrator. Hughes and Liu [13.61,13.62] developed an implicit–explicit method similar to Refs. [13.59,13.60] but more implementationally attractive. It uses element partitions rather than nodal partitions. Operator-splitting methods are mixed methods in which a nonlinear material constitutive law is split into parts that give rise to time-dependent and time-independent terms that are integrated explicitly and implicitly, respectively [13.63,13.64]. Element-by-element implicit methods have been developed in which a conventional implicit method is used except that Eqs. 13.11-5 are approximately solved using element-level calculations only [13.65,13.66,13.67]. Thus these methods appear to have the good stability characteristics of implicit methods with the low per-time-step cost of explicit methods. At present, however, these methods are not sufficiently robust and success is very problem-dependent.
# 13.13 STABILITY ANALYSIS. ACCURACY OF DIRECT INTEGRATION METHODS
In Sections 13.10 and 13.11 many stability and accuracy properties of direct integration methods are stated. In the present section we substantiate some of these properties and offer further suggestions for use of the methods. References [13.18,13.46,13.49] contain extensive discussions.
Stability. When we examine stability, it is sufficient to consider the homogeneous form of the equation of motion obtained by taking $\{R^{ext}\} = \{0\}$ . The idea is that if a solution procedure is stable with no external loading, then it will also be stable if $\{R^{ext}\}$ is nonzero but bounded. A number of methods for assessing stability are possible. They fall into two broad categories. In the first, called spectral or Fourier stability, one examines the effects of a time integration method on a single equation of motion obtained by modally uncoupling the original structure equations. In the second, called energy stability, one deals with the original structure matrices and establishes the conditions under which a norm of the solution at time $n \Delta t$ can be bounded by a norm of the solution at time zero. Spectral stability usually provides more insight and sometimes more precise results. Energy methods sometimes provide results that are slightly more conservative than spectral methods, but can be applied to complicated problems to which spectral methods may be inapplicable. References on energy stability techniques include [13.46,13.49, 13.61,13.63].