Explicit Scheme Stability

In direct integration method, at time

t_{n}

the solutions for the prior steps are known and the solution for the time

t_{n + 1} = t_{n} + Δ t

is required next. The equations to relate displacements, velocities and accelerations in a discrete time scale using the central difference time integration algorithm are given in Central Difference Algorithm. They can be rewritten as:(1)

\begin{array}{l} {\dot{u}}^{n + 1} = {\dot{u}}^{n} + \frac{Δ t}{2} ({\ddot{u}}^{n} + {\ddot{u}}^{n + 1}) = {\dot{u}}^{n + \frac{1}{2}} + \frac{Δ t}{2} {\ddot{u}}^{n + 1} \\ u^{n + 1} = u^{n} + Δ t . {\dot{u}}^{n} + \frac{Δ t^{2}}{2} . {\ddot{u}}^{n} = u^{n} + Δ t . {\dot{u}}^{n + \frac{1}{2}} \end{array}

For stability studies, aim to establish a recursive relationship to link the displacements at three consecutive time steps:(2)

[\begin{matrix} u^{n + 1} \\ u^{n} \end{matrix}] = [A] [\begin{matrix} u^{n} \\ u^{n - 1} \end{matrix}] + [L]

Where, $[A]$ is amplification matrix. A spectral analysis of this matrix highlights the stability of the integration scheme.

The numerical algorithm is stable if and only if the radius spectral of $[A]$ is less than unity. In other words, when the module of all eigen values of $[A]$ are smaller than unity, the numerical stability is ensured.

The stability of a numerical scheme can be studied using the general form of the 2x2 matrix

[A]

:(3)

[A] = [\begin{matrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix}]

Then, the equations are developed for the systems with or without damping. ¹

The eigen values of

[A]

are computed from the characteristic polynomial equation:(4)

\begin{array}{l} \det [A - λ I] = 0 \\ λ^{2} - 2 A_{1} λ + A_{2} = 0 \end{array}

Where,

$A_{1} = \frac{1}{2} t r [A] = \frac{1}{2} (A_{11} + A_{22})$
$A_{2} = \det [A] = A_{11} A_{22} - A_{12} A_{21}$

The eigen values are then obtained as:(5)

λ_{1, 2} = A_{1} \pm \sqrt{A_{1}^{2} - A_{2}}

A_{1}^{2} < A_{2}

, eigen values are complex conjugate; if

A_{1}^{2} = A_{2}

, they are real and identical; if

A_{1}^{2} > A_{2}

, they are real and distinct. You intend to define a stability domain in the

(A_{1}, A_{2})

-space, where the spectral radius

ρ ([A]) = \max (| λ_{i} (A_{1}, A_{2}) |) \leq 1

. The boundary of this domain is given by couples

(A_{1}, A_{2})

such as

ρ ([A]) = 1

. Three cases are to be considered:

Roots are real and one of them is equal to 1:
You then have:(6)
$1 - 2 A_{1} + A_{2} = 0$

This yields:(7)
$\begin{array}{l} A_{2} = 2 A_{1} - 1 \\ λ_{1} = 1 \\ λ_{2} = A_{2} \end{array}$

The corresponding part of the boundary of the stability domain is the segment analytically defined by $1 - 2 A_{1} + A_{2} = 0$ and $- 1 \leq A_{2} \leq 1$ .
Roots are real and one of them is equal to -1:
You then have:(8)
$1 + 2 A_{1} + A_{2} = 0$

This yields:(9)
$\begin{array}{l} A_{2} = - 2 A_{1} - 1 \\ λ_{1} = - 1 \\ λ_{2} = - A_{2} \end{array}$

In this case, the corresponding part of the boundary is the segment given by $1 + 2 A_{1} + A_{2} = 0$ and $- 1 \leq A_{2} \leq 1$ .
Roots are complex conjugate:
Their modulus is equal to 1. You then have, using $λ_{1, 2} = e^{\pm i α}$ :(10)
$\begin{array}{l} 0 = e^{2 i α} - 2 A_{1} e^{i α} + A_{2} \\ = (\cos 2 α - 2 A_{1} \cos α + A_{2}) + i (\sin 2 α - 2 A_{1} \sin α) \\ = (2 \cos α (\cos α - A_{1}) + A_{2} - 1) + i (2 \sin α (\cos α - A_{1})) \end{array}$

This yields:(11)
$\begin{array}{l} 2 \cos α (\cos α - A_{1}) + A_{2} - 1 = 0 \\ 2 \sin α (\cos α - A_{1}) = 0 \end{array}$

Since $\sin α \neq 0$ , you obtain:(12)
$\begin{array}{l} A_{1} = \cos α \\ A_{2} = 1 \end{array}$

The corresponding part of the boundary is thus the segment given by $A_{2} = 1$ and $- 1 \leq A_{2} \leq 1$ .

The 3 segments introduced above define a closed contour. Point $A_{1} = A_{2} = 0$ is located inside this contour and in this case, $ρ ([A]) = 0$ . Since $ρ ([A])$ varies continuously with respect to $A_{1}$ and $A_{2}$ , you can conclude that the stability domain corresponds to the interior of the contour. To precisely define the stability domain, you must also have points leading to double eigen value of modulus 1, that is, the intersections between the parabola $A_{1}^{2} = A_{2}$ and the boundary of the domain. This corresponds to Points $(A_{1}, A_{2}) = (\pm 1, 1)$ .

Figure 1. Stability Domain

You can analytically summarize the description of the stability by means of the following two sets of conditions:(13)
$\begin{array}{l} (1) - \frac{(A_{2} + 1)}{2} \leq A_{1} \leq \frac{(A_{2} + 1)}{2}, - 1 \leq A_{2} < 1 \\ (2) - 1 < A_{1} < 1, A_{2} = 1 \end{array}$

Numerical Stability of Undamped Systems

The stability conditions developed in the previous section can be applied to a one degree-of-freedom of a system without damping. The dynamic equilibrium equation at time

t_{n}

is given by:(14)

m {\ddot{u}}^{n} + k u^{n} = f^{n}

Where,

m

and

k

are respectively the nodal mass and stiffness.

f^{n}

is the external force at time

t_{n}

. Rewriting the central difference time integration equations from Equation 1, you obtain:(15)

\begin{array}{l} u^{n + 1} = u^{n} + Δ t . {\dot{u}}^{n + \frac{1}{2}} = u^{n} + Δ t . {\dot{u}}^{n - \frac{1}{2}} + Δ t^{2} . {\ddot{u}}^{n} \\ u^{n} = u^{n - 1} + Δ t . {\dot{u}}^{n - \frac{1}{2}} \end{array}

and:(16)

{\ddot{u}}^{n} = \frac{u^{n + 1} - 2 u^{n} + u^{n - 1}}{Δ t^{2}}

Substituting these equations into Equation 14, it yields:(17)

m \frac{u^{n + 1} - 2 u^{n} + u^{n - 1}}{Δ t^{2}} + k u^{n} = f^{n}

This equation can be written as Equation 2. Then the amplification matrix takes the expression:(18)

[A] = [\begin{matrix} 2 - ω^{2} Δ t^{2} & - 1 \\ 1 & 0 \end{matrix}]

Where, $ω = \sqrt{\frac{k}{m}}$ is the angular frequency of the considered mode.

Comparing with Equation 3, you have

A_{1} = 1 - \frac{ω^{2} Δ t^{2}}{2}

and

A_{2} = 1

. Stability is then given by:(19)

- 1 < 1 - \frac{ω^{2} Δ t^{2}}{2} < 1

The right inequality is always true if

ω

≠ 0. For, the particular case of

ω

= 0, the scheme is unstable. However, the analytical solution for a system with

ω

= 0 leads to an unbounded solution. The left inequality implies:(20)

Δ t < \frac{2}{ω}

Numerical Stability with Viscous Damping: Velocities at Time Steps

The dynamic equilibrium equation at time step

n

is written as:(21)

m {\ddot{u}}^{n} + c {\dot{u}}^{n} + k u^{n} = f^{n}

Using the equations:(22)

\begin{array}{l} u^{n + 1} = u^{n} + Δ t . {\dot{u}}^{n + \frac{1}{2}} \\ u^{n} = u^{n - 1} + Δ t . {\dot{u}}^{n - \frac{1}{2}} \end{array}

Results in:(23)

u^{n + 1} - u^{n - 1} = Δ t ({\dot{u}}^{n + \frac{1}{2}} + {\dot{u}}^{n - \frac{1}{2}})

For the velocity, write the equations:(24)

\begin{array}{l} {\dot{u}}^{n} = {\dot{u}}^{n - \frac{1}{2}} + \frac{Δ t}{2} {\ddot{u}}^{n} \\ {\dot{u}}^{n + \frac{1}{2}} = {\dot{u}}^{n} + \frac{Δ t}{2} {\ddot{u}}^{n} \end{array}

to obtain:(25)

{\dot{u}}^{n} = \frac{1}{2} ({\dot{u}}^{n - \frac{1}{2}} + {\dot{u}}^{n + \frac{1}{2}}) = \frac{1}{2 Δ t} (u^{n + 1} - u^{n - 1})

Substituting these equations into Equation 21, the recurring continuation equation on the displacement is written in the form:(26)

m \frac{u^{n + 1} - 2 u^{n} + u^{n - 1}}{Δ t^{2}} + \frac{c}{2 Δ t} (u^{n + 1} - u^{n - 1}) + k u^{n} = f^{n}

The equation can be rearranged to obtain the expression of the amplification matrix:(27)

[A] = [\begin{matrix} \frac{2 - ω^{2} Δ t^{2}}{1 + \frac{c Δ t}{2 m}} & \frac{- 1 + \frac{c Δ t}{2 m}}{1 + \frac{c Δ t}{2 m}} \\ \begin{array}{l} 1 \end{array} & \begin{array}{l} 0 \end{array} \end{matrix}]

This yields $A_{1} = \frac{1 - \frac{ω^{2} Δ t^{2}}{2}}{1 + \frac{c Δ t}{2 m}}$ and $A_{2} = \frac{1 - \frac{c Δ t}{2 m}}{1 + \frac{c Δ t}{2 m}}$ .

Stability is given by the set of conditions from Equation 13:(28)

- \frac{1}{1 + \frac{c Δ t}{2 m}} \leq \frac{1 - \frac{ω^{2} Δ t^{2}}{2}}{1 + \frac{c Δ t}{2 m}} \leq \frac{1}{1 + \frac{c Δ t}{2 m}}

$- 1 \leq \frac{1 - \frac{c Δ t}{2 m}}{1 + \frac{c Δ t}{2 m}} < 1$

The second expression is always verified for

c > 0

. It is the same for the right inequality of the first expression. The left inequality of the first expression leads to the condition on the time step:(29)

Δ t \leq \frac{2}{ω}

You find the same condition as in the undamped case, which echoes a conclusion given in ¹. You may yet remark that damping has changed the strict inequality into a large inequality, preventing from weak instability due to a double eigen value of modulus unity.

It is important to note that the relation Equation 29 is obtained by using the expression Equation 25 to compute nodal velocities at time steps. However, in an explicit scheme generally the mid-step velocities ${\dot{u}}^{n + \frac{1}{2}}$ and ${\dot{u}}^{n - \frac{1}{2}}$ are used. This will be studied in the next section.

Numerical Stability with Viscous Damping: Velocities at Mid Steps

Considering the case in which damping effects cannot be neglected, you still would like to deal with decoupled equilibrium equations to be able to use essentially the same computational procedure. Except for the case of full modal projection which is a very expensive technique and practically unused, the damping matrix

[C]

is not diagonal, contrary to

[M]

. The computation of the viscous forces with the exact velocity given by the integration algorithm requires the matrix

[M] + \frac{Δ t}{2} [C]

to be inverted, which can harm the numerical performances. You therefore often compute the viscous forces using the velocities at the preceding mid-step, which are explicit. This leads to an equilibrium at step

n

in the form:(30)

m {\ddot{u}}^{n} + c . {\dot{u}}^{n - \frac{1}{2}} + k u^{n} = f^{n}

The integration algorithm immediately yields:(31)

{\dot{u}}^{n - \frac{1}{2}} = \frac{1}{Δ t} (u^{n} - u^{n - 1})

The recurring continuation becomes:(32)

m \frac{u^{n + 1} - 2 u^{n} + u^{n - 1}}{Δ t^{2}} + \frac{c}{Δ t} (u^{n} - u^{n - 1}) + k u^{n} = f^{n}

As above, you obtain the amplification matrix:(33)

[A] = [\begin{matrix} 2 - ω^{2} Δ t^{2} - \frac{c Δ t}{m} & - 1 + \frac{c Δ t}{m} \\ \begin{array}{l} 1 \end{array} & \begin{array}{l} 0 \end{array} \end{matrix}]

You have in this case $A_{1} = 1 - \frac{ω^{2} Δ t^{2}}{2} - \frac{c Δ t}{2 m}$ and $A_{2} = 1 - \frac{c Δ t}{m}$ .

Stability is again given by the set of conditions Equation 13:(34)

\begin{array}{l} - 1 + \frac{c Δ t}{2 m} \leq 1 - \frac{ω^{2} Δ t^{2}}{2} - \frac{c Δ t}{2 m} \leq 1 - \frac{c Δ t}{2 m} \\ - 1 \leq 1 - \frac{c Δ t}{m} < 1 \end{array}

Right inequalities are always verified in both preceding expressions. Left inequalities now lead to two conditions on the time step:(35)

Δ t \leq \frac{- \frac{c}{m} + \sqrt{\frac{c^{2}}{m^{2}} + 4 ω^{2}}}{ω^{2}}, Δ t \leq \frac{2 m}{c}

Therefore, the critical time step depends not only to $ω$ but also to the mass and the damping. However, the critical time step depends only to $ω$ when using the exact velocities to compute the viscous forces as described in the previous section.

Numerical Stability with Rayleigh Damping

The linearized equations of equilibrium governing the dynamic response of a finite element system can derived from the equations of motion given in Equation of Motion for Translational Velocities and Equation of Motion for Angular Velocities:(36)

[M] {\ddot{U}} + [C] {\dot{U}} + [K] {U} = {F}

In the case of direct step-by-step time integration, it is necessary to evaluate the damping matrix

[C]

explicitly. The Rayleigh damping method assumes that the matrix

[C]

is computed by the following equation:(37)

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

Where,

$[C]$: Viscous damping matrix of the system
$[M]$: Mass matrix
$[K]$: Stiffness matrix

As described in the preceding sections, the computation of the viscous forces by using velocities at time steps leads to obtain a non-diagonal matrix

[C]

which should be inverted in the resolution procedure. To avoid the high cost operations, generally the simplifications are made to obtain a diagonal matrix. Substituting the Rayleigh equation into Equation 36 and using the mid-step velocities for

β [K]

terms and at step nodal velocities for

α [M]

terms, the following expression is obtained:(38)

[M] {{\ddot{u}}^{n}} + α [M] {{\dot{u}}^{n}} + β [K] {{\dot{u}}^{n - \frac{1}{2}}} + [K] {u^{n}} = {f^{n}}

Studying the equilibrium of a node to obtain a one dimensional equation of motion, write:(39)

m \ddot{u} + c \dot{u} + k u = f

Where,

$m$: Modal mass
$c$: Associated modal damping
$k$: Nodal stiffness

This leads to the following recurring continuation on the displacement:(40)

m \frac{u^{n + 1} - 2 u^{n} + u^{n - 1}}{Δ t^{2}} + \frac{α m}{2 Δ t} (u^{n + 1} - u^{n - 1}) + \frac{β k}{Δ t} (u^{n} - u^{n - 1}) + k u^{n} = f^{n}

The amplification matrix is then:(41)

[A] = [\begin{matrix} \frac{2 - ω^{2} Δ t^{2} - β ω^{2} Δ t}{1 + \frac{α Δ t}{2}} & \frac{- 1 + \frac{α Δ t}{2} + β ω^{2} Δ t}{1 + \frac{α Δ t}{2}} \\ \begin{array}{l} 1 \end{array} & \begin{array}{l} 0 \end{array} \end{matrix}]

In this case, $A_{1} = \frac{1 - \frac{ω^{2} Δ t^{2}}{2} - \frac{β ω^{2} Δ t}{2}}{1 + \frac{α Δ t}{2}}$ and $A_{2} = \frac{1 - \frac{α Δ t}{2} - β ω^{2} Δ t}{1 + \frac{α Δ t}{2}}$ .

Stability is obtained as before by means of the set of conditions from Equation 13:(42)

\begin{array}{l} \frac{- 1 + \frac{β ω^{2} Δ t}{2}}{1 + \frac{α Δ t}{2}} \leq \frac{1 - \frac{ω^{2} Δ t^{2}}{2} - \frac{β ω^{2} Δ t}{2}}{1 + \frac{α Δ t}{2}} \leq \frac{1 - \frac{β ω^{2} Δ t}{2}}{1 + \frac{α Δ t}{2}} \\ - 1 \leq \frac{1 - \frac{α Δ t}{2} - β ω^{2} Δ t}{1 + \frac{α Δ t}{2}} < 1 \end{array}

This again yields two conditions on the time step, coming from the left inequalities in both expressions:(43)

Δ t \leq \frac{- β ω + \sqrt{β^{2} ω^{2} + 4}}{ω}, Δ t \leq \frac{2}{β ω^{2}}

It is equivalent to consider only the

β [K]

contribution in the damping for the computation of the time step, which appears to be logical since the

α [M]

contribution is used with the exact velocity. It is advantageous to separate the two contributions, restrictions of the time step then becoming lighter. It can be shown that for the complete treatment of the Rayleigh damping using mid-step velocities, the stability conditions can be given by:(44)

Δ t \leq \frac{- α - β ω^{2} + \sqrt{{(α + β ω^{2})}^{2} + 4 ω^{2}}}{ω}, Δ t \leq \frac{2}{α + β ω^{2}}

Example: Critical Time Step for a Mass-Spring System

Consider a free mass-spring system without damping. The governing differential equation can be written as:

$M \ddot{X} + K X = 0$	(a)

The element matrix expressions are given as:

$[M] = m [\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}]$ ; $[K] = k [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]$

The general solution is assumed in the form of:

$X = C o s (ω t - α) [I] Ψ$	(b)

Where, the angular frequency

ω

and the phase angle

α

are common for all

X_{i}

α

and

Ψ_{i}

are the constants of integration to be determined from the initial conditions of the motion and

ω

is a characteristic value (eigen value) of the system. Substituting (b) into (a) yields:

$(K - ω^{2} M) Ψ C o s (ω t - α) = 0$	(c)

The consistency of (c) for

Ψ C o s (ω t - α) \neq 0

requires that:

$\det (K - ω^{2} M) = 0$ => $ω^{2} = \frac{2 k}{m}$	(d)

Assuming the following numerical values $m = 1$ and $k = 10$ , you have $ω = \sqrt{\frac{2 k}{m}} = 4.472136$ .

The critical time step of the system is given by Equation 20:

$Δ t \leq \frac{2}{ω} \Rightarrow Δ t \leq \frac{2}{4.472136} \Rightarrow Δ t \leq 0.4472$

Example: Critical Time Step for Dropping Body

A dropping body is studied in Body Drop Example with analytical and numerical approaches. As shown in Figure 5, the numerical results are closed to the analytical solution if you use a step-by-step time discretization method with a constant time step $Δ t = 0.1$ . This value is less than the critical time step obtained by:

Δ t_{c r} = \frac{2}{ω_{\max}}

which may be computed as:

$ω = \sqrt{\frac{k}{m}} = \sqrt{20} \Rightarrow Δ t_{c r} = \frac{2}{4.472136} \Rightarrow Δ t_{c r} = 0.4472$

Therefore, the used time step in the Body Drop Example ensures the stability of the numerical scheme as it is less than the critical value. Now using the values higher than or equal to the critical time step, the solution diverges towards the infinity as shown in Figure 4.

It is worthwhile to comment that in a general explicit finite element program as Radioss, the critical time step is computed for an entire element (two nodal masses and stiffness for spring element). In the case of dropping body example, the mass-spring system can be compared by analogy to a two-node mass-spring system where the global stiffness is twice softer. The critical time step is then computed using the nodal time step of the entire element (refer to the following sections ¹ for more details on the computation of nodal time step).

¹ Hughes T.J.R., “Analysis of Transient Algorithms with Particular Reference to Stability Behavior”, Computational Methods for Transient Analysis, eds. T. Belytschko and T.J.R. Hugues, 67-155, North Holland, 1983.