From c07bfae91742308eec94e3ec5b203d5edf62dc0c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Fri, 2 Oct 2026 06:57:45 -0300 Subject: [PATCH 1/4] ExplicitDifferenceStatic: document -alpha/-simple, combined law, examples The page described a fixed alpha = 0.59 and only the simple damping law. Document the options `integrator ExplicitDifferenceStatic <-alpha $alpha> <-simple>`, the default combined law acting on the assembled unbalance P - R(u), Tcl and Python availability, and add two verified examples: SDOF free-vibration decay and a pseudo-static hanging chain that relaxes to the static solution. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- .../codeExample/example_pseudostatic.tcl | 99 +++++++++ .../codeExample/example_sdof_decay.tcl | 120 +++++++++++ .../integrator/ExplicitDifferenceStatic.rst | 191 ++++++++++++++++-- 3 files changed, 389 insertions(+), 21 deletions(-) create mode 100644 source/user/manual/analysis/codeExample/example_pseudostatic.tcl create mode 100644 source/user/manual/analysis/codeExample/example_sdof_decay.tcl diff --git a/source/user/manual/analysis/codeExample/example_pseudostatic.tcl b/source/user/manual/analysis/codeExample/example_pseudostatic.tcl new file mode 100644 index 00000000..62155277 --- /dev/null +++ b/source/user/manual/analysis/codeExample/example_pseudostatic.tcl @@ -0,0 +1,99 @@ +# ============================================================================ +# example_pseudostatic.tcl +# +# Demonstrates using ExplicitDifferenceStatic as a pseudo-static / dynamic +# relaxation tool: a chain of n truss "springs" with lumped masses, loaded +# with constant nodal forces (representing, e.g., a soil column's self +# weight), is relaxed with strong Cundall/FLAC local damping until velocities +# vanish and the displacements converge to the static (hand-calculated) +# solution. +# +# Model: node 0 (fixed) -- k -- node 1 -- k -- node 2 -- k -- node 3 -- k -- node 4 +# each free node i=1..n carries lumped mass m and constant load P +# +# Hand calculation (static equilibrium of the chain): +# Tension in spring i (connecting node i-1 to node i) must carry the sum +# of all loads applied at nodes i, i+1, ..., n: +# T_i = (n - i + 1) * P +# Elongation of spring i: delta_i = T_i / k +# Cumulative displacement at node i: +# u_i = sum_{j=1}^{i} delta_j = sum_{j=1}^{i} (n - j + 1) * P / k +# For n=4, P=2, k=1000: +# u_1 = 4*P/k = 0.008 +# u_2 = u_1 + 3*P/k = 0.014 +# u_3 = u_2 + 2*P/k = 0.018 +# u_4 = u_3 + 1*P/k = 0.020 +# ============================================================================ + +wipe +model basic -ndm 1 -ndf 1 + +set n 4 +set k 1000.0 +set m 1.0 +set P 2.0 + +node 0 0.0 +fix 0 1 +for {set i 1} {$i <= $n} {incr i} { + node $i [expr {double($i)}] + mass $i $m +} + +uniaxialMaterial Elastic 1 $k +for {set i 1} {$i <= $n} {incr i} { + element truss $i [expr {$i-1}] $i 1.0 1 +} + +pattern Plain 1 "Constant" { + for {set i 1} {$i <= $n} {incr i} { + load $i $P + } +} + +# Default local damping (alpha=0.59, combined) -- no -alpha/-simple flags. +integrator ExplicitDifferenceStatic + +constraints Plain +numberer Plain +system FullGeneral +test NormUnbalance 1.0e-12 10 +algorithm Linear +analysis Transient + +# Conservative dt: the n-DOF chain's highest mode is faster than the +# single-DOF estimate omega=sqrt(k/m), so apply a generous safety factor. +set omega [expr {sqrt($k/$m)}] +set dt [expr {2.0/$omega*0.2}] + +set nSteps 4000 +for {set j 1} {$j <= $nSteps} {incr j} { + analyze 1 $dt +} + +puts "========================================================================" +puts " ExplicitDifferenceStatic -- pseudo-static relaxation of a loaded chain" +puts "========================================================================" +puts [format "dt = %.6g, nSteps = %d, t_total = %.4g (s)" $dt $nSteps [expr {$nSteps*$dt}]] +puts "" +puts [format "%-6s %16s %16s %16s %10s" "node" "disp (got)" "disp (expected)" "vel (got)" "err %"] + +set maxErr 0.0 +for {set i 1} {$i <= $n} {incr i} { + set u 0.0 + for {set j 1} {$j <= $i} {incr j} { + set u [expr {$u + ($n-$j+1)*$P/$k}] + } + set got [nodeDisp $i 1] + set v [nodeVel $i 1] + set err [expr {abs($got-$u)/$u*100.0}] + if {$err > $maxErr} { set maxErr $err } + puts [format "%-6d %16.8f %16.8f %16.3e %10.4f" $i $got $u $v $err] +} + +puts "" +if {$maxErr < 0.1} { + puts [format "PASS: max displacement error = %.4f %% (< 0.1 %%), velocities ~ 0" $maxErr] +} else { + puts [format "FAIL: max displacement error = %.4f %%" $maxErr] +} diff --git a/source/user/manual/analysis/codeExample/example_sdof_decay.tcl b/source/user/manual/analysis/codeExample/example_sdof_decay.tcl new file mode 100644 index 00000000..218baa75 --- /dev/null +++ b/source/user/manual/analysis/codeExample/example_sdof_decay.tcl @@ -0,0 +1,120 @@ +# ============================================================================ +# example_sdof_decay.tcl +# +# Verifies the local non-viscous (Cundall/FLAC) damping of the +# ExplicitDifferenceStatic integrator using a single-degree-of-freedom (SDOF) +# spring-mass oscillator in free vibration (no applied load). +# +# Model: node 1 (fixed) -- truss (k) -- node 2 (mass m), 1D, ndf 1 +# Initial condition: u(0) = u0, v(0) = 0 (set with setNodeDisp -commit) +# +# Theory (see ExplicitDifferenceStatic.rst): +# For free vibration (no external load) the unbalance is B = -k*u, so +# dB/dt = -k*v, i.e. sign(dB/dt) = -sign(v) at all times (away from the +# velocity-sign deadband). Substituting into the "combined" damping law +# F_d = 0.5*alpha*|B|*(sign(dB/dt) - sign(v)) +# gives F_d = -alpha*|B|*sign(v), i.e. for this particular problem the +# "combined" formulation collapses onto the "simple" one, +# F_d = -alpha*|B|*sign(v). +# A first-order (small alpha) energy balance over one cycle of amplitude A, +# frequency omega, gives an equivalent viscous damping ratio +# xi = alpha / pi +# and hence a logarithmic decrement per cycle of +# delta = 2*pi*xi = 2*alpha +# so successive displacement peaks should decay by the ratio +# A_(n+1) / A_n ~= exp(-2*alpha) +# (this is the "checked quantity" below; it is a small-alpha asymptote, so +# the match degrades as alpha approaches 1). +# ============================================================================ + +proc runSDOF {alpha mode} { + wipe + model basic -ndm 1 -ndf 1 + + set k 1000.0 + set m 1.0 + set u0 1.0 + + node 1 0.0 + node 2 1.0 + fix 1 1 + mass 2 $m + + uniaxialMaterial Elastic 1 $k + element truss 1 1 2 1.0 1 + + set omega [expr {sqrt($k/$m)}] + set T [expr {2*3.14159265358979/$omega}] + + # Initial displacement u0, zero velocity -- pure free vibration, no load. + setNodeDisp 2 1 $u0 -commit + + if {$mode eq "simple"} { + integrator ExplicitDifferenceStatic -alpha $alpha -simple + } else { + integrator ExplicitDifferenceStatic -alpha $alpha + } + + constraints Plain + numberer Plain + system FullGeneral + test NormUnbalance 1.0e-12 10 + algorithm Linear + analysis Transient + + set dt [expr {$T/400.0}] + set nCycles 6 + set nSteps [expr {int($nCycles*$T/$dt)}] + + set prevDisp $u0 + set prevVel 0.0 + set peaks {} + for {set i 1} {$i <= $nSteps} {incr i} { + analyze 1 $dt + set d [nodeDisp 2 1] + set v [nodeVel 2 1] + if {$prevVel >= 0 && $v < 0} { + lappend peaks $prevDisp + } + set prevDisp $d + set prevVel $v + } + return $peaks +} + +puts "========================================================================" +puts " ExplicitDifferenceStatic -- Cundall/FLAC local damping, SDOF free vibration" +puts "========================================================================" +puts [format "%-6s %-10s %14s %14s %14s %10s" "alpha" "mode" "peak1" "peak2" "ratio" "exp(-2a)"] + +foreach alpha {0.0 0.3 0.8} { + foreach mode {combined simple} { + set peaks [runSDOF $alpha $mode] + set p1 [lindex $peaks 0] + set p2 [lindex $peaks 1] + set ratio [expr {$p2/$p1}] + set expected [expr {exp(-2.0*$alpha)}] + puts [format "%-6s %-10s %14.6f %14.6f %14.6f %10.6f" $alpha $mode $p1 $p2 $ratio $expected] + } +} + +puts "" +puts "Check (small-alpha asymptote xi = alpha/pi, decay ratio = exp(-2*alpha)):" +set peaksC [runSDOF 0.1 combined] +set got [expr {[lindex $peaksC 1]/[lindex $peaksC 0]}] +set expected [expr {exp(-2.0*0.1)}] +set err [expr {abs($got-$expected)/$expected*100.0}] +puts [format " alpha=0.1 combined: got ratio = %.6f, expected = %.6f, error = %.3f %%" $got $expected $err] +if {$err < 1.0} { + puts " PASS: within 1% of the theoretical small-alpha decay ratio" +} else { + puts " FAIL: deviates by more than 1% from theory" +} + +puts "" +puts "Note: for this unloaded free-vibration problem, 'combined' and '-simple'" +puts "give (numerically) identical decay, because sign(dB/dt) = -sign(v)" +puts "always holds when the unbalance B is purely elastic (B=-k*u). The two" +puts "forms only differ once the unbalance's time-derivative sign decouples" +puts "from the velocity sign, e.g. under nonlinear material response or" +puts "time-varying external load (see the pseudo-static example)." diff --git a/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst b/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst index 890df40a..cade922b 100644 --- a/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst +++ b/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst @@ -2,56 +2,205 @@ ExplicitDifferenceStatic ------------------------ -.. function:: integrator 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. + +.. function:: integrator ExplicitDifferenceStatic <-alpha $alpha> <-simple> + +.. csv-table:: + :header: "Argument", "Type", "Description", "Default" + :widths: 15, 10, 55, 15 + + "-alpha $alpha", "|float|", "Local non-viscous damping coefficient, :math:`0 \le \alpha < 1`", "0.59" + "-simple", "flag", "Use the simple damping law :math:`F_d=-\alpha|F|\,\mathrm{sign}(v)` instead of the default combined law", "combined (off)" + +Available in both the Tcl and Python interpreters. .. note:: - * Uses leap-frog integration with velocities at half time steps. - * FLAC-style local non-viscous damping (α = 0.59) for pseudo-static analysis. - * Only mass matrix is required (no tangent matrix assembly). - * Velocity sign memory prevents chatter near zero velocity. - * For stability: :math:`\Delta t \leq \frac{2}{\omega_{max}}` - * Suitable for quasi-static and dynamic problems with adaptive damping. + * Leap-frog scheme: velocities are stored at half time steps + :math:`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 :math:`M`); no stiffness tangent is formed, so + a ``Linear`` algorithm is required and only one ``update()`` call per + step is allowed. + * Stability (undamped): :math:`\Delta t \le 2/\omega_{max}`; with Rayleigh + damping ratio :math:`\xi`: :math:`\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 (:math:`\alpha_M,\beta_K,\beta_{K_i},\beta_{K_c}`), + which is not exposed through the Tcl/Python parser options above. + * Passing ``-alpha 0.0`` disables the damping loop + (``alphaLNVD <= 0.0`` short-circuits in ``formUnbalance``), with no + overhead difference from the plain ``ExplicitDifference`` integrator. + Note that omitting ``-alpha`` entirely does *not* give :math:`\alpha=0`; + it falls back to the default 0.59. + * **Before the fix that added the options above**, the local damping force + was computed from ``DOF_Group::getUnbalance()``, which only holds the + *applied* nodal loads, not the assembled unbalance :math:`P-R(u)`. As a + result the damping term scaled with the applied load only and vanished + identically in free vibration (no applied load), making the integrator's + behavior indistinguishable from plain ``ExplicitDifference`` in that + case. The fix applies the damping to the fully assembled unbalance in + ``formUnbalance``, as documented below. + * See also ``integrator ExplicitBathe $p <0|1> <-lnvd $alpha>`` (on a + separate branch), which adds the same "simple" local non-viscous damping + law to Bathe's composite explicit scheme (no "combined" option). Theory ^^^^^^ -The ExplicitDifferenceStatic method is an explicit difference scheme with FLAC-style local non-viscous damping. The method uses a leap-frog approach where: +**Leap-frog update.** Given velocity at :math:`t-\Delta t/2`, acceleration at +:math:`t`, and displacement at :math:`t`: + +.. math:: -* Velocities are defined at half time steps: :math:`v_{n+1/2}` -* Displacements are defined at full time steps: :math:`d_n` + v_{t+\Delta t/2} = v_{t-\Delta t/2} + \Delta t\, a_t -The damping force uses an adaptive FLAC-style approach: + U_{t+\Delta t} = U_t + \Delta t\, v_{t+\Delta t/2} + +The new acceleration is obtained from :math:`M a_{t+\Delta t} = F_{t+\Delta t}` +where :math:`F = P - R(u)` is the fully assembled nodal unbalance (applied +load minus the assembled resisting force; only :math:`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]_ +(also described in the FLAC/FLAC3D theory manuals, Itasca Consulting Group) +removes energy by adding to the assembled unbalance a force proportional to +the *magnitude of the unbalance itself* and opposing the velocity: .. math:: - F_{damping} = -\alpha \cdot |F_{unbalanced}| \cdot sign(v) + F_d = -\alpha\, |F| \,\mathrm{sign}(v) \qquad \text{("simple", } \texttt{-simple}\text{)} + +or, combining the sign of the unbalance's time rate of change with the sign +of velocity ("combined", the default): -where α = 0.59 by default. The method includes velocity sign memory with a deadband (1e-4) to prevent numerical chatter when velocities approach zero. +.. math:: -The central difference equations are: + F_d = \tfrac{1}{2}\alpha\, |F| \,\big(\mathrm{sign}(\dot F) - \mathrm{sign}(v)\big) + +:math:`\mathrm{sign}(v)` is replaced by a held sign :math:`s_v` with a small +deadband (:math:`|v|<10^{-4}` keeps the previous sign) to avoid chatter as the +velocity passes through zero. :math:`\dot F` is approximated by the backward +difference :math:`(F_t - F_{t-\Delta t})/\Delta t` (only the sign is used, so +the :math:`\Delta t` cancels). This damping acts equation-by-equation on the +fully assembled unbalance :math:`F=P-R(u)` and is added in ``formUnbalance`` +after the base ``TransientIntegrator::formUnbalance`` call; it is unrelated +to (and does not replace) any Rayleigh damping. + +**Equivalent viscous damping ratio (small** :math:`\alpha` **).** For an +undamped SDOF oscillator (:math:`F=-kU`) with no applied load, free vibration +gives :math:`\dot F = -kv`, so :math:`\mathrm{sign}(\dot F)=-\mathrm{sign}(v)` +identically and the combined law reduces exactly to the simple law, +:math:`F_d=-\alpha|F|\,\mathrm{sign}(v)`. A first-order energy balance over +one cycle of amplitude :math:`A` then gives a fractional energy loss per +cycle :math:`\Delta E/E = 4\alpha`, which matches a viscously damped +oscillator with damping ratio .. math:: - v_{n+1/2} = v_{n-1/2} + \Delta t \cdot a_n + \xi \approx \frac{\alpha}{\pi} - d_{n+1} = d_n + \Delta t \cdot v_{n+1/2} +and hence a logarithmic decrement per cycle :math:`\delta=2\pi\xi\approx +2\alpha`, i.e. successive displacement peaks decay by -This integrator is particularly effective for quasi-static analysis and problems where adaptive local damping is beneficial, such as rock mechanics and soil-structure interaction problems. +.. math:: -.. admonition:: Example + \frac{A_{n+1}}{A_n} \approx e^{-2\alpha} + +This relation is accurate to a fraction of a percent for :math:`\alpha +\lesssim 0.1` (for :math:`\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-:math:`\alpha` asymptote as +:math:`\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 :math:`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 :math:`\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 +:math:`\alpha` (or the plain ``ExplicitDifference`` integrator, equivalent to +:math:`\alpha=0`) when running genuine dynamic/earthquake analyses where the +damping should not contaminate the computed response. + +.. admonition:: Example 1. **Tcl Code** .. code-block:: tcl + integrator ExplicitDifferenceStatic -alpha 0.8 + integrator ExplicitDifferenceStatic -alpha 0.3 -simple integrator ExplicitDifferenceStatic + ;# default: -alpha 0.59, combined form 2. **Python Code** .. code-block:: python - integrator('ExplicitDifferenceStatic') - -.. [Cundall1987] Cundall, P.A. (1987). "Distinct Element Models of Rock and Soil Structure." In Analytical and Computational Methods in Engineering Rock Mechanics, 129-163. + ops.integrator('ExplicitDifferenceStatic', '-alpha', 0.8, '-simple') + + 3. **SDOF free-vibration decay example** + + A single-DOF spring-mass oscillator released from rest at :math:`u_0=1` + with no applied load. Theory predicts the ratio of successive + displacement peaks :math:`A_{n+1}/A_n \approx e^{-2\alpha}`; running the + script below with :math:`\alpha=0.1` gives a decay ratio of **0.8182** + against the theoretical **0.8187** (0.07% error), confirming the + small-:math:`\alpha` estimate and that the "combined" and "simple" laws + coincide in free vibration (:math:`\mathrm{sign}(\dot F)=-\mathrm{sign}(v)` + identically for a purely elastic unbalance). + + | :download:`example_sdof_decay.tcl <../codeExample/example_sdof_decay.tcl>` **(TCL)**. + + 4. **Pseudo-static relaxation example** + + A 4-DOF chain of springs and lumped masses, each node carrying a constant + load :math:`P`, relaxed with the default damping (:math:`\alpha=0.59`, + combined) until velocities vanish. The hand-calculated static + displacements are :math:`u_1=0.008`, :math:`u_2=0.014`, :math:`u_3=0.018`, + :math:`u_4=0.020`; running the script below reproduces these to machine + precision with residual velocities of order :math:`10^{-17}`. + + | :download:`example_pseudostatic.tcl <../codeExample/example_pseudostatic.tcl>` **(TCL)**. + +Usage notes and limitations +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +* 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). +* ``-alpha`` is rejected (integrator construction fails, returns ``0``) for + values outside :math:`[0,1)`. +* The ``-alpha``/``-simple`` 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 :math:`\dot F` + decouples from the sign of :math:`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). +* Prior to the fix described in the note above, this damping had **no + effect** in free-vibration / unloaded-structure analyses, since it only + ever saw the (zero or constant) applied-load vector, not the element + resisting forces; any model relying on the old behavior for genuinely + unloaded dynamic relaxation was not actually being damped. + +.. [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*, section on local (Cundall) damping. Code Developed by: |jaabell| From 2ff571b3dbcffafa9c0a8fd7ac7d509316326b78 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Fri, 2 Oct 2026 18:21:29 -0300 Subject: [PATCH 2/4] ExplicitDifferenceStatic: fix Sphinx warnings (FLAC citation, Tcl block) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- .../manual/analysis/integrator/ExplicitDifferenceStatic.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst b/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst index cade922b..1d8925cf 100644 --- a/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst +++ b/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst @@ -67,7 +67,7 @@ load minus the assembled resisting force; only :math:`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]_ +proportional to velocity, Cundall's local damping scheme [Cundall1987]_ [FLAC]_ (also described in the FLAC/FLAC3D theory manuals, Itasca Consulting Group) removes energy by adding to the assembled unbalance a force proportional to the *magnitude of the unbalance itself* and opposing the velocity: @@ -143,8 +143,8 @@ damping should not contaminate the computed response. integrator ExplicitDifferenceStatic -alpha 0.8 integrator ExplicitDifferenceStatic -alpha 0.3 -simple + # default: -alpha 0.59, combined form integrator ExplicitDifferenceStatic - ;# default: -alpha 0.59, combined form 2. **Python Code** From d11e735e2f0e207afe8a6133cbf8373d91ecaa85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Fri, 2 Oct 2026 19:18:26 -0300 Subject: [PATCH 3/4] ExplicitDifferenceStatic: document final behaviour (-vEps, acceleration-based damping, non-fatal options) Describe the damping in user terms (equivalent to F_d = -alpha |F| sign(v) on the unbalanced force for the lumped mass the scheme uses), note that it works with the parallel diagonal systems and is not disturbed by printB, add -vEps, state that invalid options are warned and ignored, and that the damping history resets on domain changes. Update the ExplicitBathe see-also to its current syntax. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- .../integrator/ExplicitDifferenceStatic.rst | 117 +++++++++++------- 1 file changed, 74 insertions(+), 43 deletions(-) diff --git a/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst b/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst index 1d8925cf..a53077ee 100644 --- a/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst +++ b/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst @@ -9,7 +9,7 @@ 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. -.. function:: integrator ExplicitDifferenceStatic <-alpha $alpha> <-simple> +.. function:: integrator ExplicitDifferenceStatic <-alpha $alpha> <-simple> <-vEps $vEps> .. csv-table:: :header: "Argument", "Type", "Description", "Default" @@ -17,6 +17,7 @@ dynamics integrator. "-alpha $alpha", "|float|", "Local non-viscous damping coefficient, :math:`0 \le \alpha < 1`", "0.59" "-simple", "flag", "Use the simple damping law :math:`F_d=-\alpha|F|\,\mathrm{sign}(v)` instead of the default combined law", "combined (off)" + "-vEps $vEps", "|float|", "Velocity-sign dead-band :math:`\varepsilon_v \ge 0`, absolute, in model velocity units: while :math:`|v|\le\varepsilon_v` the previously held sign of :math:`v` is reused instead of recomputed", "1e-4" Available in both the Tcl and Python interpreters. @@ -25,29 +26,42 @@ Available in both the Tcl and Python interpreters. :math:`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 :math:`M`); no stiffness tangent is formed, so - a ``Linear`` algorithm is required and only one ``update()`` call per + a ``Linear`` algorithm is required and only one acceleration solve per step is allowed. * Stability (undamped): :math:`\Delta t \le 2/\omega_{max}`; with Rayleigh damping ratio :math:`\xi`: :math:`\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 (:math:`\alpha_M,\beta_K,\beta_{K_i},\beta_{K_c}`), which is not exposed through the Tcl/Python parser options above. - * Passing ``-alpha 0.0`` disables the damping loop - (``alphaLNVD <= 0.0`` short-circuits in ``formUnbalance``), with no + * The damping is applied directly to the solved nodal acceleration at + each step. Because this scheme requires a lumped (diagonal) mass + matrix, :math:`Ma=F` holds equation-by-equation, so damping the + acceleration this way is equivalent to damping the unbalanced force + :math:`F=P-R(u)` as :math:`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 :math:`\alpha=0`; it falls back to the default 0.59. - * **Before the fix that added the options above**, the local damping force - was computed from ``DOF_Group::getUnbalance()``, which only holds the - *applied* nodal loads, not the assembled unbalance :math:`P-R(u)`. As a - result the damping term scaled with the applied load only and vanished - identically in free vibration (no applied load), making the integrator's - behavior indistinguishable from plain ``ExplicitDifference`` in that - case. The fix applies the damping to the fully assembled unbalance in - ``formUnbalance``, as documented below. - * See also ``integrator ExplicitBathe $p <0|1> <-lnvd $alpha>`` (on a - separate branch), which adds the same "simple" local non-viscous damping - law to Bathe's composite explicit scheme (no "combined" option). + * An out-of-range ``-alpha`` (outside :math:`[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). Theory ^^^^^^ @@ -67,30 +81,45 @@ load minus the assembled resisting force; only :math:`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]_ [FLAC]_ -(also described in the FLAC/FLAC3D theory manuals, Itasca Consulting Group) -removes energy by adding to the assembled unbalance a force proportional to -the *magnitude of the unbalance itself* and opposing the velocity: +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 :math:`v` +at :math:`t+\Delta t/2`: .. math:: - F_d = -\alpha\, |F| \,\mathrm{sign}(v) \qquad \text{("simple", } \texttt{-simple}\text{)} + a_d = a - \alpha\, |a| \,\mathrm{sign}(v) \qquad \text{("simple", } \texttt{-simple}\text{)} -or, combining the sign of the unbalance's time rate of change with the sign -of velocity ("combined", the default): +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]_): .. math:: - F_d = \tfrac{1}{2}\alpha\, |F| \,\big(\mathrm{sign}(\dot F) - \mathrm{sign}(v)\big) + a_d = a + \tfrac{1}{2}\alpha\, |a| \,\big(\mathrm{sign}(\Delta a) - \mathrm{sign}(v)\big) -:math:`\mathrm{sign}(v)` is replaced by a held sign :math:`s_v` with a small -deadband (:math:`|v|<10^{-4}` keeps the previous sign) to avoid chatter as the -velocity passes through zero. :math:`\dot F` is approximated by the backward -difference :math:`(F_t - F_{t-\Delta t})/\Delta t` (only the sign is used, so -the :math:`\Delta t` cancels). This damping acts equation-by-equation on the -fully assembled unbalance :math:`F=P-R(u)` and is added in ``formUnbalance`` -after the base ``TransientIntegrator::formUnbalance`` call; it is unrelated -to (and does not replace) any Rayleigh damping. +where :math:`\Delta a = a - a_{\text{prev}}` is the change from the previous +step's (already damped) acceleration. :math:`\mathrm{sign}(v)` is replaced by +a held sign :math:`s_v` with a dead-band set by ``-vEps`` (default +:math:`10^{-4}`): while :math:`|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, :math:`Ma=F` +holds equation-by-equation, with :math:`F=P-R(u)` the assembled unbalance +(applied load minus resisting force). Substituting :math:`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: + +.. math:: + + 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)} + +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** :math:`\alpha` **).** For an undamped SDOF oscillator (:math:`F=-kU`) with no applied load, free vibration @@ -143,7 +172,8 @@ damping should not contaminate the computed response. integrator ExplicitDifferenceStatic -alpha 0.8 integrator ExplicitDifferenceStatic -alpha 0.3 -simple - # default: -alpha 0.59, combined form + integrator ExplicitDifferenceStatic -alpha 0.59 -vEps 1.0e-6 + # default: -alpha 0.59, combined form, vEps 1e-4 integrator ExplicitDifferenceStatic 2. **Python Code** @@ -182,11 +212,13 @@ Usage notes and limitations * 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). -* ``-alpha`` is rejected (integrator construction fails, returns ``0``) for - values outside :math:`[0,1)`. -* The ``-alpha``/``-simple`` 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). +* An out-of-range ``-alpha`` (outside :math:`[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 :math:`\dot F` @@ -194,13 +226,12 @@ Usage notes and limitations 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). -* Prior to the fix described in the note above, this damping had **no - effect** in free-vibration / unloaded-structure analyses, since it only - ever saw the (zero or constant) applied-load vector, not the element - resisting forces; any model relying on the old behavior for genuinely - unloaded dynamic relaxation was not actually being damped. +* 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*, section on local (Cundall) damping. +.. [FLAC] Itasca Consulting Group, Inc. *FLAC/FLAC3D Theory and Background -- Dynamic Analysis: Damping*, sections on local and combined (local/combined) damping. Code Developed by: |jaabell| From d00a97da9be6fe803e8139222556d8b8bb1c9a21 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Antonio=20Abell?= Date: Fri, 2 Oct 2026 19:19:26 -0300 Subject: [PATCH 4/4] ExplicitDifferenceStatic: note lumped-mass assumption of the local damping Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01K524skQQ24A6JCFWNnb7G1 --- .../manual/analysis/integrator/ExplicitDifferenceStatic.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst b/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst index a53077ee..1157eb1d 100644 --- a/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst +++ b/source/user/manual/analysis/integrator/ExplicitDifferenceStatic.rst @@ -209,6 +209,9 @@ damping should not contaminate the computed response. 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 :math:`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).