Static Solution by Explicit Time Integration

Explicit algorithms are very useful for modeling a dynamic simulation. However, they cannot model a quasi-static or static simulation as easily. This is due to the fact that in an explicit approach, first the nodal accelerations are found by resolving the equilibrium equation at time $t_{n}$ $t_{n}$ . Other DOFs are then computed by explicit time integration. This procedure implies that the nodal acceleration must exist; however, some numerical methods may be employed for the simulation of a static process.

Slow Dynamic Computation

The loading is applied at a rate sufficiently slow to minimize the dynamic effects. The final solution is obtained by smoothing the curves.

In case of elasto-plastic problems, one must minimize dynamic overshooting because of the irreversibility of the plastic flow.

Dynamic Relaxation or Nodal Damping (/DYREL)

This method was first introduced by Otter ¹ and has been used in several hydrodynamic codes. A nodal damping

C υ

$C υ$ is added to the momentum equation:(1)

M \frac{d ν}{d t} + C υ = f_{e x t} - f^{int}

$M \frac{d ν}{d t} + C υ = f_{e x t} - f^{int}$

The dashpot force is calculated by:(2)

F_{c} = - C V

$F_{c} = - C V$

The internal force is calculated by:(3)

F_{k} = - K X = F_{k t - Δ t / 2} + K V Δ t

$F_{k} = - K X = F_{k t - Δ t / 2} + K V Δ t$

(4)

γ_{0} = \frac{F_{k}}{M}

$γ_{0} = \frac{F_{k}}{M}$

The total acceleration is given by:(5)

γ_{1} = γ_{0} + \frac{F_{c}}{M} = γ_{0} - \frac{C V}{M}

$γ_{1} = γ_{0} + \frac{F_{c}}{M} = γ_{0} - \frac{C V}{M}$

(6)

V_{t + Δ t / 2} = V_{t - Δ t / 2} + γ_{1} Δ t

$V_{t + Δ t / 2} = V_{t - Δ t / 2} + γ_{1} Δ t$

(7)

V = V_{t - Δ t / 2} + \frac{1}{2} γ_{1} Δ t

$V = V_{t - Δ t / 2} + \frac{1}{2} γ_{1} Δ t$

(8)

V_{t + Δ t / 2} = V_{t - Δ t / 2} + γ_{0} Δ t - \frac{C V}{M} Δ t

$V_{t + Δ t / 2} = V_{t - Δ t / 2} + γ_{0} Δ t - \frac{C V}{M} Δ t$

You have:(9)

V_{t + Δ t / 2} = V_{t - Δ t / 2} + γ_{0} Δ t - \frac{C}{M} (V_{t - Δ t / 2} + \frac{1}{2} γ_{1} Δ t) Δ t

$V_{t + Δ t / 2} = V_{t - Δ t / 2} + γ_{0} Δ t - \frac{C}{M} (V_{t - Δ t / 2} + \frac{1}{2} γ_{1} Δ t) Δ t$

(10)

V_{t + Δ t / 2} = (1 - \frac{C Δ t}{M}) V_{t - Δ t / 2} + (γ_{0} - \frac{1}{2} γ_{1} Δ t \frac{C}{M}) Δ t

$V_{t + Δ t / 2} = (1 - \frac{C Δ t}{M}) V_{t - Δ t / 2} + (γ_{0} - \frac{1}{2} γ_{1} Δ t \frac{C}{M}) Δ t$

(11)

V_{t + Δ t / 2} = (1 - \frac{C Δ t}{M}) V_{t - Δ t / 2} + (1 - Δ t \frac{C γ_{1}}{M γ_{0}}) γ_{0} Δ t

$V_{t + Δ t / 2} = (1 - \frac{C Δ t}{M}) V_{t - Δ t / 2} + (1 - Δ t \frac{C γ_{1}}{M γ_{0}}) γ_{0} Δ t$

Approximation

\frac{γ_{1}}{γ_{0}} = 1

$\frac{γ_{1}}{γ_{0}} = 1$ after the variable is changed,

ω = C \frac{Δ t}{2 M}

$ω = C \frac{Δ t}{2 M}$ , you obtain:(12)

V_{t + Δ t / 2} = (1 - 2 ω) V_{t - Δ t / 2} + (1 - ω) γ_{0} Δ t

$V_{t + Δ t / 2} = (1 - 2 ω) V_{t - Δ t / 2} + (1 - ω) γ_{0} Δ t$

(13)

ω = C \frac{Δ t}{2 M} \to C = \frac{2 M ω}{Δ t}

$ω = C \frac{Δ t}{2 M} \to C = \frac{2 M ω}{Δ t}$

Which gives the expression of

C

$C$ as a proportional matrix to

M

$M$ with:(14)

$C = A M$ $C = A M$ and $A = \frac{2 ω}{Δ t}$ $A = \frac{2 ω}{Δ t}$

or(15)

A = \frac{2 β}{Τ}

$A = \frac{2 β}{Τ}$

(16)

ω = β \frac{Δ t}{Τ}

$ω = β \frac{Δ t}{Τ}$

Combining Equation 14 and Equation 15, and you obtain:(17)

C = \frac{2 β M}{Τ}

$C = \frac{2 β M}{Τ}$

Where, $β$ $β$ is the relaxation coefficient whose recommended value is 1. $T$ $T$ is less than or equal to the highest period of the system. These are the input parameters used in /DYREL option.

The explicit time integration scheme is changed to compute the new velocities. The explicit time integration in Dynamic Analysis gives (Central Difference Algorithm, Equation 7):(18)

V_{t + Δ t / 2} = V_{t - Δ t / 2} + γ_{t} Δ t

$V_{t + Δ t / 2} = V_{t - Δ t / 2} + γ_{t} Δ t$

which is now written as: (19)

V_{t + Δ t / 2} = (1 - 2 ω) V_{t - Δ t / 2} + (1 - ω) γ_{t} Δ t

$V_{t + Δ t / 2} = (1 - 2 ω) V_{t - Δ t / 2} + (1 - ω) γ_{t} Δ t$

Where,(20)

ω = β \frac{Δ t}{Τ}

$ω = β \frac{Δ t}{Τ}$

Energy Discrete Relaxation

This empirical methodology consists in setting to zero the nodal velocities each time the Kinetic Energy reaches a maximum.

The loading is applied at a rate sufficiently slow to minimize the dynamic effects. The final solution is obtained by smoothing the curves.

In case of elasto-plastic problems, one must minimize dynamic overshooting because of the irreversibility of the plastic flow.

Rayleigh Damping (/DAMP)

In this method a proportional damping matrix is defined as:(21)

[C] = α [M] + β [K]

$[C] = α [M] + β [K]$

Where, $α$ $α$ and $β$ $β$ are the pre-defined constants. In modal analysis, the use of a proportional damping matrix allows to reduce the global equilibrium equation to n-uncoupled equations by using an orthogonal transformation.

If the global equilibrium equation is expressed as:(22)

[M] {\ddot{X}} + [C] {\dot{X}} + [K] {X} = {F_{t}}

$[M] {\ddot{X}} + [C] {\dot{X}} + [K] {X} = {F_{t}}$

The transformed uncoupled system of equations can be written as:(23)

{[ϕ]}^{T} [M] [ϕ] {\ddot{ξ}} + {[ϕ]}^{T} [C] [ϕ] {\dot{ξ}} + {[ϕ]}^{T} [K] [ϕ] {ξ} = {[ϕ]}^{T} {F_{T}}

${[ϕ]}^{T} [M] [ϕ] {\ddot{ξ}} + {[ϕ]}^{T} [C] [ϕ] {\dot{ξ}} + {[ϕ]}^{T} [K] [ϕ] {ξ} = {[ϕ]}^{T} {F_{T}}$

With (24)

{[φ]}^{T} [C] [φ] = [\begin{matrix} α + β ω_{1}^{2} & 0 & . & . & 0 \\ 0 & α + β ω_{2}^{2} & . & . & 0 \\ . & . & . & . & . \\ . & . & . & . & . \\ 0 & . & . & . & α + β ω_{n}^{2} \end{matrix}]

${[φ]}^{T} [C] [φ] = [\begin{matrix} α + β ω_{1}^{2} & 0 & . & . & 0 \\ 0 & α + β ω_{2}^{2} & . & . & 0 \\ . & . & . & . & . \\ . & . & . & . & . \\ 0 & . & . & . & α + β ω_{n}^{2} \end{matrix}]$

Each uncoupled equation is written as:(25)

{\ddot{ξ}}_{i} + 2 ω_{i} ζ_{i} {\dot{ξ}}_{i} + ω_{i}^{2} ξ_{i} = f_{i}^{t}

${\ddot{ξ}}_{i} + 2 ω_{i} ζ_{i} {\dot{ξ}}_{i} + ω_{i}^{2} ξ_{i} = f_{i}^{t}$

With (26)

2 ζ_{i} ω_{i} = α + β ω_{i}^{2}

$2 ζ_{i} ω_{i} = α + β ω_{i}^{2}$

Where,

$ω_{i}$ $ω_{i}$: The i^th natural frequency of the system
$ζ_{i}$ $ζ_{i}$: The i^th damping ratio

This leads to a system of

n

$n$ equations with two unknown variables

α

$α$ and

β

$β$ . Regarding to the range of the dominant frequencies of system, two frequencies are chosen. Using the pair of the most significant frequencies, two equations with two unknown variables can be resolved to obtain values for

α

$α$ and

β

$β$ . For high frequencies the role of

β

$β$ is more significant. However, for lower frequencies

α

$α$ plays an important role (Figure 1).

The Rayleigh damping method applied to explicit time-integration method leads to the following equations:(27)

M γ^{t} + C v^{t} = F_{e x t}^{t} - F_{int}^{t}

$M γ^{t} + C v^{t} = F_{e x t}^{t} - F_{int}^{t}$

With

C = α M + β K

$C = α M + β K$

(28)

F_{int}^{t} = F_{int}^{t - d t} + K v^{t - \frac{d t}{2}} d t

$F_{int}^{t} = F_{int}^{t - d t} + K v^{t - \frac{d t}{2}} d t$

(29)

v^{t} = v^{t - \frac{d t}{2}} + γ^{t} \frac{d t}{2}

$v^{t} = v^{t - \frac{d t}{2}} + γ^{t} \frac{d t}{2}$

(30)

M γ_{0}^{t} = F_{e x t}^{t} - F_{int}^{t}

$M γ_{0}^{t} = F_{e x t}^{t} - F_{int}^{t}$

Neglecting

F_{e x t}^{t} - F_{e x t}^{t - d t}

$F_{e x t}^{t} - F_{e x t}^{t - d t}$ and

v^{t} - v^{t - \frac{d t}{2}}

$v^{t} - v^{t - \frac{d t}{2}}$ , in

β K v^{t}

$β K v^{t}$ evaluation you have:(31)

K v^{t} = K v^{t - \frac{d t}{2}} = \frac{F_{int}^{t} - F_{int}^{t - d t}}{d t} \approx M \frac{(γ_{0}^{t} - γ_{0}^{t - d t})}{d t}

$K v^{t} = K v^{t - \frac{d t}{2}} = \frac{F_{int}^{t} - F_{int}^{t - d t}}{d t} \approx M \frac{(γ_{0}^{t} - γ_{0}^{t - d t})}{d t}$

And finally:(32)

γ_{0}^{t} = M^{- 1} (F_{e x t}^{t} - F_{int}^{t})

$γ_{0}^{t} = M^{- 1} (F_{e x t}^{t} - F_{int}^{t})$

(33)

γ^{t} = \frac{γ_{0} - α v^{t - \frac{d t}{2}} - \frac{β}{d t} (γ_{0}^{t} - γ_{0}^{t - d t})}{1 + α \frac{d t}{2}}

$γ^{t} = \frac{γ_{0} - α v^{t - \frac{d t}{2}} - \frac{β}{d t} (γ_{0}^{t} - γ_{0}^{t - d t})}{1 + α \frac{d t}{2}}$

(34)

v^{t + \frac{d t}{2}} = v^{t - \frac{d t}{2}} + γ^{t} d t

$v^{t + \frac{d t}{2}} = v^{t - \frac{d t}{2}} + γ^{t} d t$

The three approaches available in Radioss are Dynamic Relaxation (/DYREL), Energy Discrete Relaxation (/KEREL) and Rayleigh Damping (/DAMP). Refer to Example Guide for application examples.

The loading is applied at a rate sufficiently slow to minimize the dynamic effects. The final solution is obtained by smoothing the curves.

In case of elasto-plastic problems, one must minimize dynamic overshooting because of the irreversibility of the plastic flow.

Acceleration Convergence

For every method, an acceleration of the convergence to the static solution is desirable. The constant time step is one of the more usual methods. In fact, in quasi-static analysis, the duration of the study is proportional to the maximum period of the structure. The total number of computation cycles is then proportional to the ratio $\frac{T}{d t}$ $\frac{T}{d t}$ .

Where,

$T$ $T$: Largest period of the structure
$d t$ $d t$: Time step

The number of time steps necessary to reach the static solution is minimal if all the elements have the same time step. An initial given time step $Δ t_{0}$ $Δ t_{0}$ can be obtained by increasing or decreasing the density of each element. The constant nodal time step option ensures a homogenous time step over the structure. However, in usual static problems the change is expected to be small, but one may think of increasing the density of the element which gives the critical time step in such a way that $Δ t = Δ t_{0}$ $Δ t = Δ t_{0}$ .

Otter J.R.H., Mechanical Engineering and Design 3, 1965.