3.2.6.12. ExplicitDifferenceStatic

This command is used to construct an explicit central-difference (leap-frog) transient integrator that additionally applies Cundall/FLAC-style local non-viscous damping to the assembled nodal unbalance, making it usable as a dynamic-relaxation / pseudo-static solver as well as a plain explicit dynamics integrator.

integrator ExplicitDifferenceStatic <-alpha $alpha> <-simple> <-vEps $vEps>

Argument	Type	Description	Default
-alpha $alpha	float	Local non-viscous damping coefficient, $0 \le \alpha < 1$	0.59
-simple	flag	Use the simple damping law $F_d=-\alpha\|F\|\,\mathrm{sign}(v)$ instead of the default combined law	combined (off)
-vEps $vEps	float	Velocity-sign dead-band $\varepsilon_v \ge 0$, absolute, in model velocity units: while $\|v\|\le\varepsilon_v$ the previously held sign of $v$ is reused instead of recomputed	1e-4

Available in both the Tcl and Python interpreters.

Note

Leap-frog scheme: velocities are stored at half time steps $v_{n+1/2}$, displacements and accelerations at full time steps.
Only the mass matrix is ever assembled as the tangent (formEleTangent/ formNodTangent only add $M$); no stiffness tangent is formed, so a Linear algorithm is required and only one acceleration solve per step is allowed.
Stability (undamped): $\Delta t \le 2/\omega_{max}$; with Rayleigh damping ratio $\xi$: $\Delta t \le 2/(\omega_{max}(\sqrt{1+\xi^2}-\xi))$.
The local damping is applied in addition to any Rayleigh damping set via the four-argument C++ constructor ($\alpha_M,\beta_K,\beta_{K_i},\beta_{K_c}$), which is not exposed through the Tcl/Python parser options above.
The damping is applied directly to the solved nodal acceleration at each step. Because this scheme requires a lumped (diagonal) mass matrix, $Ma=F$ holds equation-by-equation, so damping the acceleration this way is equivalent to damping the unbalanced force $F=P-R(u)$ as $F_d=-\alpha|F|\,\mathrm{sign}(v)$ (or the combined law below) – see Theory. This works correctly with the parallel diagonal systems of equations (MPIDiagonal) because the acceleration used is the fully assembled one; commands that merely re-form and print the unbalanced force (e.g. printB) without advancing the state do not perturb the stored damping history and have no effect on subsequent results.
Passing -alpha 0.0 disables the damping step entirely, with no overhead difference from the plain ExplicitDifference integrator. Note that omitting -alpha entirely does not give $\alpha=0$; it falls back to the default 0.59.
An out-of-range -alpha (outside $[0,1)$), a negative -vEps, or an unrecognized option is not fatal: a WARNING is printed and the offending option is ignored, keeping its default value; the integrator is still constructed.
The per-equation damping history (held velocity sign and previous acceleration) is reset to zero whenever the domain changes (new elements/nodes, pattern changes, numbering changes, …).
History: earlier versions of this integrator computed the damping force from the applied nodal loads only, so it vanished identically in free vibration (no applied load); this has since been corrected to use the actual solved acceleration/unbalance as described above.
See also integrator ExplicitBathe $p <$flag> <-lnvd $alpha>, which adds the same “simple” local non-viscous damping law to Bathe’s composite explicit scheme (no “combined” option).

3.2.6.12.1. Theory

Leap-frog update. Given velocity at $t-\Delta t/2$, acceleration at $t$, and displacement at $t$:

\[ \begin{align}\begin{aligned}v_{t+\Delta t/2} = v_{t-\Delta t/2} + \Delta t\, a_t\\U_{t+\Delta t} = U_t + \Delta t\, v_{t+\Delta t/2}\end{aligned}\end{align} \]

The new acceleration is obtained from $M a_{t+\Delta t} = F_{t+\Delta t}$ where $F = P - R(u)$ is the fully assembled nodal unbalance (applied load minus the assembled resisting force; only $M$ is used as the linear-system tangent).

Local non-viscous (Cundall/FLAC) damping. Rather than a viscous term proportional to velocity, Cundall’s local damping scheme [Cundall1987] removes energy by making the solved acceleration at each step depend on its own magnitude and sign and on the sign of the (leap-frog) velocity $v$ at $t+\Delta t/2$:

\[a_d = a - \alpha\, |a| \,\mathrm{sign}(v) \qquad \text{("simple", } \texttt{-simple}\text{)}\]

or, combining the sign of the step-to-step change in acceleration with the sign of velocity (“combined”, the default, following the Itasca FLAC/FLAC3D manual’s local and combined damping formulation [FLAC]):

\[a_d = a + \tfrac{1}{2}\alpha\, |a| \,\big(\mathrm{sign}(\Delta a) - \mathrm{sign}(v)\big)\]

where $\Delta a = a - a_{\text{prev}}$ is the change from the previous step’s (already damped) acceleration. $\mathrm{sign}(v)$ is replaced by a held sign $s_v$ with a dead-band set by -vEps (default $10^{-4}$): while $|v|\le\varepsilon_v$ the previously held sign is reused, to avoid chatter as the velocity passes through zero.

Because this scheme requires a lumped (diagonal) mass matrix, $Ma=F$ holds equation-by-equation, with $F=P-R(u)$ the assembled unbalance (applied load minus resisting force). Substituting $a=F/m$ term by term shows that damping the acceleration as above is exactly equivalent, for each equation, to damping the unbalanced force itself:

\[ \begin{align}\begin{aligned}F_d = -\alpha\, |F| \,\mathrm{sign}(v) \qquad \text{(simple)}\\F_d = \tfrac{1}{2}\alpha\, |F| \,\big(\mathrm{sign}(\Delta F) - \mathrm{sign}(v)\big) \qquad \text{(combined)}\end{aligned}\end{align} \]

This damping is unrelated to (and does not replace) any Rayleigh damping, and – since it operates on the already fully assembled acceleration – is applied consistently whether the system of equations is sequential or one of the parallel diagonal systems (MPIDiagonal).

Equivalent viscous damping ratio (small $\alpha$ ). For an undamped SDOF oscillator ($F=-kU$) with no applied load, free vibration gives $\dot F = -kv$, so $\mathrm{sign}(\dot F)=-\mathrm{sign}(v)$ identically and the combined law reduces exactly to the simple law, $F_d=-\alpha|F|\,\mathrm{sign}(v)$. A first-order energy balance over one cycle of amplitude $A$ then gives a fractional energy loss per cycle $\Delta E/E = 4\alpha$, which matches a viscously damped oscillator with damping ratio

\[\xi \approx \frac{\alpha}{\pi}\]

and hence a logarithmic decrement per cycle $\delta=2\pi\xi\approx 2\alpha$, i.e. successive displacement peaks decay by

\[\frac{A_{n+1}}{A_n} \approx e^{-2\alpha}\]

This relation is accurate to a fraction of a percent for $\alpha \lesssim 0.1$ (for $\alpha=0.1$ the example below gives a decay ratio of 0.8182 against the theoretical 0.8187, i.e. 0.07% error), and increasingly over-damps relative to the small-$\alpha$ asymptote as $\alpha \to 1$ (the scheme does not blow up; it simply dissipates energy faster than the first-order estimate because the waveform is no longer a small perturbation of a pure sinusoid).

Because $F_d\propto|F|$ rather than a fixed Coulomb force, the fractional damping per cycle is independent of amplitude – this is a hysteretic-type damping, not Coulomb friction, which is what makes it suitable as a general-purpose “numerical damping” for driving a system to static equilibrium (dynamic relaxation) regardless of how far from equilibrium it starts.

When to use it. Set $\alpha$ as large as possible (close to but below 1) while keeping the scheme accurate and stable if the sole goal is a static solution (push-over, gravity application, settling a soil column) – this is the classical FLAC “local damping” use case. Use a small $\alpha$ (or the plain ExplicitDifference integrator, equivalent to $\alpha=0$) when running genuine dynamic/earthquake analyses where the damping should not contaminate the computed response.

Example

Tcl Code

integrator ExplicitDifferenceStatic -alpha 0.8
integrator ExplicitDifferenceStatic -alpha 0.3 -simple
integrator ExplicitDifferenceStatic -alpha 0.59 -vEps 1.0e-6
# default: -alpha 0.59, combined form, vEps 1e-4
integrator ExplicitDifferenceStatic

Python Code

ops.integrator('ExplicitDifferenceStatic', '-alpha', 0.8, '-simple')

SDOF free-vibration decay example

A single-DOF spring-mass oscillator released from rest at $u_0=1$ with no applied load. Theory predicts the ratio of successive displacement peaks $A_{n+1}/A_n \approx e^{-2\alpha}$; running the script below with $\alpha=0.1$ gives a decay ratio of 0.8182 against the theoretical 0.8187 (0.07% error), confirming the small-$\alpha$ estimate and that the “combined” and “simple” laws coincide in free vibration ($\mathrm{sign}(\dot F)=-\mathrm{sign}(v)$ identically for a purely elastic unbalance).

example_sdof_decay.tcl (TCL).

Pseudo-static relaxation example

A 4-DOF chain of springs and lumped masses, each node carrying a constant load $P$, relaxed with the default damping ($\alpha=0.59$, combined) until velocities vanish. The hand-calculated static displacements are $u_1=0.008$, $u_2=0.014$, $u_3=0.018$, $u_4=0.020$; running the script below reproduces these to machine precision with residual velocities of order $10^{-17}$.

example_pseudostatic.tcl (TCL).

3.2.6.12.2. Usage notes and limitations

The damping assumes a lumped (diagonal) mass, as this scheme requires. If a consistent mass is used with a full (non-diagonal) system, the damping is applied equation by equation to $M^{-1}F$ and energy dissipation is no longer guaranteed.
Requires a Linear solution algorithm and (as for any explicit scheme) a nonsingular, preferably diagonal, mass matrix for every active DOF (rotational DOFs with zero mass will make the system singular).
An out-of-range -alpha (outside $[0,1)$) or -vEps (negative), or an unrecognized option, is reported with a WARNING and ignored (the previous/default value is kept) – construction never fails because of a bad option.
The -alpha/-simple/-vEps options only control the local non-viscous damping; separate Rayleigh damping is still available through the four-argument C++ constructor (not exposed by the Tcl/Python parsers).
The “combined” form equals the “simple” form whenever the unbalance is a single-valued elastic function of displacement in free vibration (no applied load); the two forms only diverge once the sign of $\dot F$ decouples from the sign of $v$ – e.g. under nonlinear/path-dependent material response, or while a time-varying load is applied (see example_pseudostatic.tcl, where a constant load is used precisely so that any residual oscillation gets damped consistently).
The damping history is reset whenever the domain changes, so re-forming the domain mid-analysis (adding elements/nodes, changing numberer, etc.) restarts the held velocity sign and previous-acceleration memory from zero for every equation.

[Cundall1987]

Cundall, P.A. (1987). “Distinct Element Models of Rock and Soil Structure.” In Analytical and Computational Methods in Engineering Rock Mechanics, 129-163.

[FLAC]

Itasca Consulting Group, Inc. FLAC/FLAC3D Theory and Background – Dynamic Analysis: Damping, sections on local and combined (local/combined) damping.

Code Developed by: José A. Abell