김경종
24 KiB
We will assume for simplicity that the mass and damping matrices do not vary with time.
10.3 Modelling of nonlinearities
10.3.1 Introduction
Dynamic loading of structures often causes excursions of stresses well into the inelastic range and the influence of geometry changes on the response is also significant in many cases. Therefore both material and geometric nonlinear effects should be considered.
Although material behaviour under dynamic loading is very complex and experimental information is scarce, for most structural materials, some general statements can be made.
For example, it has frequently been demonstrated that the instantaneous yield stress is significantly influenced by the rate of straining. Also, the value of the elasticity modulus E_{0} is found to be dependent on the strain rate. For structural materials with limited ductility, such as concrete or rock-like materials, the rate of straining can completely change the material response from elasto-plastic behaviour under low rates to brittle elastic behaviour under high rates of straining. For many structural materials there is still an urgent need for a better understanding of the observed phenomena and underlying microscopic behaviour. However, in attempting to perform an analysis of a dynamically-loaded engineering structure, we must look for an idealized material model, where possibly some compromises have to be made. Furthermore, the model parameters should readily be measurable and easily obtained from reliable experimental data.
For transient dynamic analysis, an elasto-viscoplastic model, as developed in earlier chapters, presents a very good approximation of the true behaviour for many structural materials. The predominant phenomenon of variable instantaneous yield stress is adequately modelled.
In the following, we shall develop the algorithm for the elasto-viscoplastic transient dynamic analysis of plane stress, plane strain and axisymmetric problems. The computer program DYNPAK will be documented and explained and finally, some illustrative examples are given.
10.3.2 Material model
Here we adopt the elasto-viscoplastic material model developed in Chapter 8, where the constitutive relationship is given in the form
\begin{array}{l} \dot {\epsilon} _ {n} = [ \dot {\epsilon} _ {e} ] _ {n} + [ \dot {\epsilon} _ {v p} ] _ {n} \\ = [ D ] ^ {- 1} \dot {\sigma} _ {n}: \gamma \langle \Phi_ {n} (F) \rangle \frac {\dot {\epsilon} F}{\dot {\epsilon} \sigma_ {n}} \tag {10.18} \\ \end{array}
where D is the elasticity matrix, \gamma is the fluidity parameter, F is the yield
function and \dot{\epsilon}_{n} , [\dot{\epsilon}_{e}]_{n} and [\dot{\epsilon}_{vp}]_{n} denote the total, elastic and viscoplastic strain rates at time station t_{n} . We also have the relationships
\sigma_ {n} = D [ \epsilon_ {e} ] _ {n}
\epsilon_ {n} = \left[ \epsilon_ {e} \right] _ {n} + \left[ \epsilon_ {v p} \right] _ {n} \tag {10.19}
and
\langle \Phi_ {n} (F) \rangle = 0 \quad \text { if yield has not occurred. }
= 1 \quad \text { if yield has occurred. } \tag {10.20}
Thus we can rewrite the internal resisting forces as
\boldsymbol {p} _ {n} = \int_ {\Omega} [ \boldsymbol {B} ] ^ {T} \boldsymbol {D} \left\{\epsilon_ {n} - \left[ \epsilon_ {v p} \right] _ {n} \right\} d \Omega \tag {10.21}
The temporal discretization of the equations which govern viscoplastic straining is also based on the assumption that the relationship
[ \dot {\epsilon} _ {v p} ] _ {n} = \gamma \langle \Phi_ {n} (F) \rangle \frac {\partial F}{\partial \sigma_ {n}} \tag {10.22}
is known only for discrete time stations \Delta t apart. The simplest, Euler, integration scheme will here be employed, i.e.,
[ \epsilon_ {v p} ] _ {n + 1} = [ \epsilon_ {v p} ] _ {n} + [ \dot {\epsilon} _ {v p} ] _ {n} \Delta t. \tag {10.23}
The stability limit for the time increment \Delta t , which depends on the specific form of the viscoplastic potential employed in the flow rule, has already been discussed in earlier chapters.
When we adopt the central difference scheme and the viscoplastic material model that we have just described, the algorithm at a particular time station t_{n} follows the sequence shown in Fig. 10.1.
10.3.3 Geometric nonlinearity
If we wish to cater for geometrically nonlinear elastic behaviour we can choose either a total or updated Lagrangian coordinate system. Here we choose a total Lagrangian coordinate system which coincides with the initial undeformed position of the body. ^{(3)}
It transpires that, with the central difference scheme, the only changes required to account for geometrically nonlinear effects are
(i) The modification of the strain-displacement matrix B(d_n) ,
and
(ii) The evaluation of the strains using a deformation Jacobian matrix J_{D}(d_{n}) .
flowchart
graph TD
A["d_{n+1} = \left[ M + \frac{\Delta t}{2} C \right"]^{-1} \left\{ (\Delta t)^2["-p(d_n) + f_n"] + 2M d_n - \left[" M - \frac{\Delta t}{2} C \right"] d_{n-1} \right\} ] --> B["p(d_n) = \int_Ω [B(d_n)"]^T \sigma_n dΩ]
B --> C["ε_n = [B(d_n)"]d_n\n["ε_e"]_n = ε_n_-["ε_vp"]_n]
C --> D["σ_n = D[ε_e"]_n]
D --> E{Φ_n(F) > 0}
E -->|Yes| F["[\dot{\epsilon}_vp"]_n = γ<Φ_n(F) > \frac{\partial F}{\partial \sigma_n}]
F --> G["[\epsilon_vp"]_n + [ε_vp]_n + [\dot{\epsilon}_vp]_n Δt]
G --> H["[\epsilon_vp"]_{n+1} = [ε_vp]_n + [\dot{\epsilon}_vp]_n Δt]
H --> I["[\epsilon_vp"]_{n-1} = [ε_vp]_n + [D(d_n)]^T σ_n dΩ]
I --> J["ε_n = [B(d_n)"]d_n\n["ε_e"]_n = ε_n_-["ε_vp"]_n]
J --> K["σ_n = D[ε_e"]_n]
K --> L["No"]
L --> M["[\dot{\epsilon}_vp"]_n = 0]
M --> N["End"]
subgraph Time steps
A
B
C
D
E
F
G
H
I
end
subgraph Time steps
J
K
L
M
end
Fig. 10.1 Algorithm for elasto viscoplastic straining during a time step.
We will now describe briefly the relevant background theory. All vectors and matrices are given explicitly for the plane stress, plane strain and axisymmetric applications in Table 10.1.
If the initial undeformed position of a particle of material is x_{0} and the total displacement vector at time station t_{n} is u_{n} then the coordinates of the particle are
\boldsymbol {x} _ {n} = \boldsymbol {x} _ {0} + \boldsymbol {u} _ {n} \tag {10.24}
In a total Lagrangian formulation we use Green's strains. The matrix of Green's strains is given as
\boldsymbol {E} _ {n} = \frac {1}{2} \left[ \left[ \boldsymbol {J} _ {D} \right] _ {n} ^ {T} \left[ \boldsymbol {J} _ {D} \right] _ {n} - \boldsymbol {I} \right] \tag {10.25}
Table 10.1 Vectors and matrices used in a total Lagrangian formulation
| Variables | Plane stress/strain | Axisymmetric |
| Coordinates of particle in undeformed initial configuration $x = x_0$ | $[x_0, y_0]^T$ | $[r_0, z_0]^T$ |
| Displacements $u_n$ | $[u_n, v_n]^T$ | $[u_n, w_n]^T$ |
| Coordinates of particle in deformed configuration $x_n$ | $[x_n, y_n]^T = [x_0 + u_n, y_0 + v_n]^T$ | $[r_n, z_n]^T = [r_0 + u_n, z_0 + w_n]$ |
| Vector of Green's strains $\epsilon_n$ | $\begin{bmatrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{bmatrix}_n = \begin{bmatrix} \frac{\partial u_n}{\partial x} + \frac{1}{2} \left( \frac{\partial u_n}{\partial x} \right)^2 + \frac{1}{2} \left( \frac{\partial v_n}{\partial x} \right)^2 \\ \frac{\partial v_n}{\partial y} + \frac{1}{2} \left( \frac{\partial u_n}{\partial y} \right)^2 + \frac{1}{2} \left( \frac{\partial v_n}{\partial y} \right)^2 \\ \frac{\partial u_n}{\partial y} + \frac{\partial v_n}{\partial x} + \frac{\partial u_n}{\partial x} \frac{\partial u_n}{\partial y} + \frac{\partial v_n}{\partial x} \frac{\partial v_n}{\partial y} \end{bmatrix}$ | $\begin{bmatrix} \epsilon_r \\ \epsilon_z \\ \gamma_{rz} \\ \epsilon_0 \end{bmatrix}_n = \begin{bmatrix} \frac{\partial u_n}{\partial r} + \frac{1}{2} \left( \frac{\partial u_n}{\partial r} \right)^2 + \frac{1}{2} \left( \frac{\partial w_n}{\partial r} \right)^2 \\ \frac{\partial w_n}{\partial z} + \frac{1}{2} \left( \frac{\partial u_n}{\partial z} \right)^2 + \frac{1}{2} \left( \frac{\partial w_n}{\partial z} \right)^2 \\ \frac{\partial u_n}{\partial z} + \frac{\partial w_n}{\partial r} + \frac{\partial u_n}{\partial r} \frac{\partial u_n}{\partial z} + \frac{\partial w_n}{\partial r} \frac{\partial w_n}{\partial z} \\ \frac{u_n}{r} + \frac{1}{2} \left( \frac{u_n}{r} \right)^2 \end{bmatrix}$ |
| Deformation Jacobian matrix $J_D(u_n) = [J_D]_n$ | $\begin{bmatrix} \frac{\partial x_n}{\partial x} & \frac{\partial x_n}{\partial y} \\ \frac{\partial y_n}{\partial x} & \frac{\partial y_n}{\partial y} \end{bmatrix}$ | $\begin{bmatrix} \frac{\partial r_n}{\partial r} & \frac{\partial r_n}{\partial z} \\ \frac{\partial z_n}{\partial r} & \frac{\partial z_n}{\partial z} \end{bmatrix}$ |
| Matrix of Green's strains $E_n = \frac{1}{2} \{ [J_D]_n^T [J_D]_n - I \}$ | $\begin{bmatrix} \epsilon_{xx} & \epsilon_{xy} \\ \epsilon_{yx} & \epsilon_{yy} \end{bmatrix}_n$ | $\begin{bmatrix} \epsilon_{rr} & \epsilon_{rz} \\ \epsilon_{zr} & \epsilon_{zz} \end{bmatrix}_n$ |
| Linear strains $[\epsilon_L]_n$ | $\left[ \frac{\partial u_n}{\partial x}, \frac{\partial v_n}{\partial y}, \left( \frac{\partial u_n}{\partial y} + \frac{\partial v_n}{\partial x} \right) \right]^T$ | $\left[ \frac{\partial u_n}{\partial r}, \frac{\partial w_n}{\partial r}, \frac{\partial u_n}{\partial z} + \frac{\partial w_n}{\partial r}, \frac{u_n}{r} \right]^T$ |
Table 10.1 (Cont.)
| Variable | Plane stress/strain | Axisymmetric |
| Nonlinear strains $[\epsilon_{NL}]_n = \frac{1}{2}[A_\theta]_n \theta_n$ where $[A_\theta]_n$ is | $\begin{bmatrix} \frac{\partial u_n}{\partial x} & \frac{\partial v_n}{\partial x} & 0 & 0 \\ 0 & 0 & \frac{\partial u_n}{\partial y} & \frac{\partial v_n}{\partial y} \\ \frac{\partial u_n}{\partial y} & \frac{\partial v_n}{\partial y} & \frac{\partial u_n}{\partial x} & \frac{\partial v_n}{\partial x} \end{bmatrix}$ | $\begin{bmatrix} \frac{\partial u_n}{\partial r} & \frac{\partial w_n}{\partial r} & 0 & 0 & 0 \\ 0 & 0 & \frac{\partial u_n}{\partial z} & \frac{\partial w_n}{\partial z} & 0 \\ \frac{\partial u_n}{\partial z} & \frac{\partial w_n}{\partial z} & \frac{\partial u_n}{\partial r} & \frac{\partial w_n}{\partial r} & 0 \\ 0 & 0 & 0 & 0 & \frac{u_n}{r} \end{bmatrix}$ |
| and displacement gradients $\theta_n$ | $\begin{bmatrix} \frac{\partial u_n}{\partial x} & 0 & \frac{\partial u_n}{\partial x} \\ \frac{\partial v_n}{\partial x} & 0 & \frac{\partial v_n}{\partial x} \\ 0 & \frac{\partial u_n}{\partial y} & \frac{\partial u_n}{\partial y} \\ 0 & \frac{\partial v_n}{\partial y} & \frac{\partial v_n}{\partial y} \end{bmatrix}$ | $\begin{bmatrix} \frac{\partial u_n}{\partial r} & 0 & \frac{\partial u_n}{\partial r} & 0 \\ \frac{\partial w_n}{\partial r} & 0 & \frac{\partial w_n}{\partial r} & 0 \\ 0 & \frac{\partial u_n}{\partial z} & \frac{\partial u_n}{\partial z} & 0 \\ 0 & \frac{\partial w_n}{\partial z} & \frac{\partial w_n}{\partial z} & 0 \\ 0 & 0 & 0 & \frac{u_n}{r} \end{bmatrix}$ |
| Elastic Piola-Kirchoff stresses $\sigma_n = D_n \epsilon_n$ | $[\sigma_x, \sigma_y, \tau_{xy}]_n^T$ | $[\sigma_r, \sigma_z, \tau_{rz}, \sigma_\theta]_n^T$ |
where [J_D]_n is the deformation Jacobian matrix at time station t_n .
The Green's strains can be written as
\epsilon_ {n} = \left[ \epsilon_ {L} \right] _ {n} + \left[ \epsilon_ {N L} \right] _ {n} \tag {10.26}
where [\epsilon_{L}]_{n} are the linear strains given earlier in Chapter 6 and [\epsilon_{NL}]_{n} , the nonlinear strain terms are given as
[ \epsilon_ {N L} ] _ {n} = \frac {1}{2} [ A _ {\theta} ] _ {n} \theta_ {n}. \tag {10.27}
For a set of virtual displacements, the corresponding virtual Green's strains are given as
[ \delta \epsilon ] _ {n} = [ \delta \epsilon_ {L} ] _ {n} + [ A _ {\theta} ] _ {n} \delta \theta_ {n}. \tag {10.28}
Thus the virtual work statement of (10.1) can be rewritten as
\begin{array}{l} \int_ {\Omega} \left[ \delta \boldsymbol {\epsilon} _ {n} \right] ^ {T} \boldsymbol {\sigma} _ {n} d \Omega - \int_ {\Omega} \left[ \delta \boldsymbol {u} _ {n} \right] ^ {T} \left[ \boldsymbol {b} _ {n} - \rho \dot {\boldsymbol {u}} _ {n} - c \dot {\boldsymbol {u}} _ {n} \right] d \Omega \\ - \int_ {\Gamma_ {t}} \left[ \delta \boldsymbol {u} _ {n} \right] ^ {T} \boldsymbol {t} _ {n} d \Gamma = 0 \tag {10.29} \\ \end{array}
where \sigma_{n} are the Piola–Kirchhoff stresses.
As mentioned earlier, all relevant terms are given in Table 10.1.
If we adopt the finite element discretization scheme described earlier, then the displacement gradients \theta_{n} are given in terms of the nodal displacements [d_{i}]_{n} by the linear relation
\boldsymbol {\theta} _ {n} = \sum_ {i = 1} ^ {m} \boldsymbol {G} _ {i} [ \boldsymbol {d} _ {i} ] _ {n} \tag {10.30}
where G_{i} contains Cartesian shape function derivatives as indicated in Table 10.2 for the various applications.
Similarly we have
\delta \theta_ {n} = \sum_ {i = 1} ^ {m} G _ {i} [ \delta d _ {i} ] _ {n}. \tag {10.31}
The linear strain-displacement relationship can be expressed as
[ \boldsymbol {\epsilon} _ {L} ] _ {n} = \sum_ {i = 1} ^ {m} [ \boldsymbol {B} _ {L i} ] _ {n} [ \boldsymbol {d} _ {i} ] _ {n} \tag {10.32}
where [B_{Li}]_{n} is the linear strain displacement matrix introduced earlier.
Table 10.2 The nonlinear strain displacement matrix evaluation in a total Lagrangian finite element formulation
| Variable | Plane stress/strain | Axisymmetric |
| Strain displacement matrix associated with node i $[B_i]_n = [B_{Li}]_n + [A_0]_n G_i$ | $\begin{bmatrix} \frac{\partial x_n}{\partial x} & \frac{\partial N_i}{\partial x} & \frac{\partial y_n}{\partial x} & \frac{\partial N_i}{\partial x} \\ \frac{\partial x_n}{\partial y} & \frac{\partial N_i}{\partial y} & \frac{\partial y_n}{\partial y} & \frac{\partial N_i}{\partial y} \\ \left( \frac{\partial x_n}{\partial y} \frac{\partial N_i}{\partial x} + \frac{\partial x_n}{\partial x} \frac{\partial N_i}{\partial y} \right) & \left( \frac{\partial y_n}{\partial y} \frac{\partial N_i}{\partial x} + \frac{\partial y_n}{\partial x} \frac{\partial N_i}{\partial y} \right) \end{bmatrix}$ | $\begin{bmatrix} \frac{\partial r_n}{\partial r} & \frac{\partial N_i}{\partial r} & \frac{\partial z_n}{\partial r} & \frac{\partial N_i}{\partial r} \\ \frac{\partial r_n}{\partial z} & \frac{\partial N_i}{\partial z} & \frac{\partial z_n}{\partial z} & \frac{\partial N_i}{\partial z} \\ \left( \frac{\partial r_n}{\partial z} \frac{\partial N_i}{\partial r} + \frac{\partial r_n}{\partial r} \frac{\partial N_i}{\partial z} \right) & \frac{\partial z_n}{\partial z} \frac{\partial N_i}{\partial r} + \frac{\partial z_n}{\partial r} \frac{\partial N_i}{\partial z} \\ \left( \frac{r_n}{r} \right) \frac{N_i}{r} & 0 \end{bmatrix}$ |
| where $G_i$ is | $\begin{bmatrix} \frac{\partial N_i}{\partial x} & 0 & \frac{\partial N_i}{\partial y} & 0 \\ 0 & \frac{\partial N_i}{\partial x} & 0 & \frac{\partial N_i}{\partial y} \end{bmatrix}$ | $\begin{bmatrix} \frac{\partial N_i}{\partial r} & 0 & \frac{\partial N_i}{\partial z} & 0 & \frac{N_i}{r} \\ 0 & \frac{\partial N_i}{\partial r} & 0 & \frac{\partial N_i}{\partial z} & 0 \end{bmatrix}$ |
Similarly, we have
[ \delta \epsilon_ {N L} ] _ {n} = \sum_ {i = 1} ^ {m} [ B _ {N L i} ] _ {n} [ \delta d _ {i} ] _ {n} \tag {10.33}
The components of the vector of Green's strains \epsilon_{n} can be written as
\epsilon_ {n} = \sum_ {i = 1} ^ {m} \left[ \left[ B _ {L i} \right] _ {n} + \frac {1}{2} \left[ B _ {N L i} \right] _ {n} \right] \left[ d _ {i} \right] _ {n} \tag {10.34}
where the nonlinear strain-displacement matrix [B_{NLt}]_n is given as
[ \boldsymbol {B} _ {N L i} ] _ {n} = [ \boldsymbol {A} _ {\theta} ] _ {n} \boldsymbol {G} _ {i}. \tag {10.35}
Furthermore it can be shown that the virtual strains can be expressed as
\delta \epsilon_ {n} = \sum_ {i = 1} ^ {m} [ B _ {i} ] _ {n} [ \delta d _ {i} ] _ {n} \tag {10.36}
where
[ \pmb {B} _ {i} ] _ {n} = [ \pmb {B} _ {L i} ] _ {n} + [ \pmb {B} _ {N L i} ] _ {n}
is given in Table 10.2 for the various applications.
If we substitute for \delta \epsilon_{n} and \delta d_{n} in (10.29) and note that the result is true for arbitrary virtual displacements, then we obtain an expression which is identical to (10.4). In the present case we only need to remember that [B_i]_n is defined by (10.36).
We again note that contributions to (10.4) from each element can be obtained separately and assembled appropriately.
Note that we now may evaluate [p_{i}]_{n} as
\int_ {\Omega} [ B _ {i} ] _ {n} ^ {T} \sigma_ {n} d \Omega \quad \text { rather than } \quad \int_ {\Omega} [ B _ {i} ] ^ {T} \sigma_ {n} d \Omega
where [B_i]_n is given by (10.36).
10.4 Explicit time integration scheme
10.4.1 Central difference approximation
We can write the equations (10.4) in matrix form so that at time station t_n we have
\boldsymbol {M} \ddot {\boldsymbol {d}} _ {n} + \boldsymbol {C} \dot {\boldsymbol {d}} _ {n} + \boldsymbol {p} _ {n} = \boldsymbol {f} _ {n} \tag {10.37}
- Note that the body force term -\mathbf{M}\ddot{\mathbf{u}}_g , due to seismic excitation, is included in the body forces which are taken into account in \mathbf{f}_n . Note also that \mathbf{M} and \mathbf{C} may be assembled from the element mass matrices \mathbf{M}^{(e)} and damping matrices \mathbf{C}^{(e)} .
where M and C are the global mass and damping matrices respectively, p_{n} is the global vector of internal resisting nodal forces, f_{n} is the vector of consistent nodal forces for the applied body and surfaces traction forces grouped together, \ddot{d}_{n} is the global vector of nodal accelerations and \dot{d}_{n} is the global vector of nodal velocities.
So far, only spatial discretization has been introduced. We now employ a temporal discretization of the dynamic equilibrium equations by approximating the accelerations and velocities using finite difference expressions.
In particular we adopt a central difference approximation ^{(2)} so that the accelerations can be written as
\ddot {\boldsymbol {d}} _ {n} \simeq \boldsymbol {a} _ {n} = \frac {1}{(\Delta t) ^ {2}} \left\{\boldsymbol {d} _ {n + 1} - 2 \boldsymbol {d} _ {n} + \boldsymbol {d} _ {n - 1} \right\} \tag {10.38}
and the velocities are written as
\dot {\boldsymbol {d}} _ {n} \simeq \boldsymbol {v} _ {n} = \frac {1}{2 \Delta t} \left\{\boldsymbol {d} _ {n - 1} - \boldsymbol {d} _ {n - 1} \right\} \tag {10.39}
in which \Delta t is the time step or interval so that we are sampling the displacements at time stations t_n - \Delta t , t_n and t_n + \Delta t . If we substitute (10.38) and (10.39) into (10.37) we obtain
M \left\{\frac {d _ {n + 1} - 2 d _ {n} - d _ {n - 1}}{(\Delta t) ^ {2}} \right\} - C \left\{\frac {d _ {n + 1} - d _ {n - 1}}{2 \Delta t} \right\} - p _ {n} = f _ {n} \tag {10.40}
which can be rearranged to give
\begin{array}{l} \boldsymbol {d} _ {n + 1} = \left[ \boldsymbol {M} + \frac {\Delta t}{2} \boldsymbol {C} \right] ^ {- 1} \\ \times \left\{(\Delta t) ^ {2} \left[ - p _ {n} + f _ {n} \right] + 2 M d _ {n} - \left[ M - \frac {\Delta t}{2} C \right] d _ {n - 1} \right\}. \tag {10.41} \\ \end{array}
Thus we have
\boldsymbol {d} _ {n + 1} = \boldsymbol {g} \left(\boldsymbol {d} _ {n}, \boldsymbol {d} _ {n - 1}\right). \tag {10.42}
In other words the displacements at time station t_{n} = \Delta t are given explicitly in terms of the displacements at time stations t_{n} and t_{n} = \Delta t .
If the mass matrix M and the damping matrix C are diagonal then the solution of (10.41) becomes trivial and we have for plane stress and plane strain applications the following equations:
\begin{array}{l} \left(d _ {u i}\right) _ {n + 1} = \left(m _ {i i} + \frac {\Delta t}{2} C _ {i i}\right) ^ {- 1} \left[ (\Delta t) ^ {2} \left\{- \left(p _ {u i}\right) _ {n} \dots \left(f _ {u i}\right) _ {n} \right\} \right. \\ \left. + 2 m _ {i i} \left(d _ {u i}\right) _ {n} - \left(m _ {i i} - \frac {\Delta t}{2} c _ {i i}\right) \left(d _ {u i}\right) _ {n - 1} \right] \tag {10.43} \\ \end{array}
and
\begin{array}{l} (d _ {v i}) _ {n + 1} = \left(m _ {i i} + \frac {\Delta t}{2} c _ {i i}\right) ^ {- 1} \left[ (\Delta t) ^ {2} \{- (p _ {v i}) _ {n} + (f _ {v i}) _ {n} \} \right. \\ \left. + 2 m _ {i i} \left(d _ {v i}\right) _ {n} - \left(m _ {i i} - \frac {\Delta t}{2} c _ {i i}\right) \left(d _ {v i}\right) _ {n - 1} \right] \tag {10.44} \\ \end{array}
in which at node i, d_{ut} and d_{vi} are the u and v displacement components in the x and y directions, f_{ut} and f_{vt} are the components of the applied nodal forces in the x and y directions, p_{ut} and p_{vt} are the internal resisting nodal forces in the x and y directions and m_{it} and c_{it} are the diagonal terms of the mass and damping matrices. For axisymmetric problems replace v by w.
From (10.43) and (10.44) we see that for each displacement degree of freedom at time t_{n}+\Delta t we have a separate equation involving information regarding the degree of freedom at times t_{n} and t_{n}-\Delta t . No matrix factorisation or sophisticated equation solving is therefore necessary.
10.4.2 Starting algorithm
As we have seen the governing equilibrium equation at time station t_{n} + \Delta t in the central difference method involves information at the two previous time stations t_{n} and t_{n} - \Delta t . A starting algorithm is therefore necessary and from the initial conditions the values d(0 - \Delta t) may be obtained. We have from (10.39) the condition that
\dot {\boldsymbol {d}} (0) \simeq \boldsymbol {v} (0) = \frac {\boldsymbol {d} (0 + \Delta t) - \boldsymbol {d} (0 - \Delta t)}{2 \Delta t} \tag {10.45}
or
\boldsymbol {d} (0 - \Delta t) = - 2 \Delta t v (0) + \boldsymbol {d} (0 + \Delta t).
If this approximation is substituted in (10.43) then we can write the expression
\begin{array}{l} \left(d _ {u i}\right) _ {1} = \left(m _ {i i} + \frac {\Delta t}{2} c _ {i i}\right) ^ {- 1} \left[ (\Delta t) ^ {2} \{- \left(p _ {u i}\right) _ {0} + \left(f _ {u i}\right) _ {0} \} \right. \\ \left. + 2 m _ {i i} \left(d _ {u i}\right) _ {0} - \left(m _ {i i} - \frac {\Delta^ {t}}{2} c _ {i i}\right) \left\{- 2 \Delta t \left(\dot {d} _ {u i}\right) _ {0} + d _ {u i}\right) _ {1} \right\} \Bigg ] \tag {10.46} \\ \end{array}
or
(d _ {u i}) _ {1} = \frac {(\Delta t) ^ {2}}{2 m _ {i i}} \{- (p _ {u i}) _ {0} + (f _ {u i}) _ {0} \} + (\dot {d} _ {u i}) _ {0} + B \Delta t (d _ {u i}) _ {0}
where
B = 1 - \frac {c _ {i i} \Delta t}{2 m _ {i i}}.
10.4.3 Damping
Very limited information is available on damping in linear solid mechanics problems and there is even less data available for damping in nonlinear situations. It is therefore customary to assume that the damping