diff --git a/source/user/manual/analysis/codeExample/example_bathe_lnvd_pseudostatic.tcl b/source/user/manual/analysis/codeExample/example_bathe_lnvd_pseudostatic.tcl new file mode 100644 index 00000000..6efb3e63 --- /dev/null +++ b/source/user/manual/analysis/codeExample/example_bathe_lnvd_pseudostatic.tcl @@ -0,0 +1,101 @@ +# ============================================================================ +# example_pseudostatic.tcl +# +# Demonstrates using ExplicitBathe with -lnvd 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. Same model as the ExplicitDifferenceStatic pseudo-static example. +# +# 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 + } +} + +# p = 0.5 (no intrinsic numerical damping from the Bathe sub-stepping itself) +# with strong local non-viscous damping alpha=0.8 to drive the system to +# static equilibrium. +integrator ExplicitBathe 0.5 0 -lnvd 0.8 + +constraints Plain +numberer Plain +system Diagonal +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 " ExplicitBathe -lnvd 0.8 -- 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_bathe_lnvd_sdof_decay.tcl b/source/user/manual/analysis/codeExample/example_bathe_lnvd_sdof_decay.tcl new file mode 100644 index 00000000..2cf75c8b --- /dev/null +++ b/source/user/manual/analysis/codeExample/example_bathe_lnvd_sdof_decay.tcl @@ -0,0 +1,122 @@ +# ============================================================================ +# example_sdof_decay.tcl +# +# Verifies the local non-viscous (Cundall/FLAC) damping option (-lnvd) of the +# ExplicitBathe integrator using a single-degree-of-freedom (SDOF) spring-mass +# oscillator in free vibration (no applied load). +# +# Model: node 1 (fixed) -- truss (k=100) -- node 2, 1D, ndf 1 +# element truss 1 1 2 1.0 1 -rho 2.0 (nodes 1 unit apart, A=1 -> total +# mass = rho*A*L = 2.0, lumped half to each node -> 1.0 at the free node 2) +# Initial condition: u(0) = u0 = 0.01, v(0) = 0 (set with setNodeDisp -commit) +# +# The Bathe sub-step parameter is fixed at p = 0.5, for which q0 = q1 = 0 and +# q2 = 0.5 (see ExplicitBathe.rst): this is the sub-case of the scheme with no +# *intrinsic* numerical damping of its own, isolating the effect of -lnvd. +# +# Theory (see ExplicitDifferenceStatic.rst for the derivation, which carries +# over unchanged since -lnvd uses the same "simple" Cundall/FLAC law +# F_d = -alpha*|F|*sign(v) applied to the fully assembled unbalance F=P-R(u), +# here evaluated using the predicted sub-step velocity at each of the two +# Bathe sub-steps): +# For free vibration (no external load), F = -k*u, and, to first order in +# alpha, each cycle of amplitude A loses a fraction of energy +# dE/E = 4*alpha, equivalent to a viscous damping ratio xi = alpha/pi, i.e. +# a logarithmic decrement per cycle delta = 2*alpha, so successive +# displacement peaks should decay by the ratio +# A_(n+1) / A_n ~= exp(-2*alpha) (small-alpha asymptote) +# This is checked below for alpha=0.1 (error < 1%); for the larger +# alpha=0.3 and alpha=0.8 used afterwards the method still damps strongly +# but, as for ExplicitDifferenceStatic, increasingly over-damps relative to +# this small-alpha asymptote. +# ============================================================================ + +proc runSDOF {alpha} { + wipe + model basic -ndm 1 -ndf 1 + + set k 100.0 + set u0 0.01 + + node 1 0.0 + node 2 1.0 + fix 1 1 + + uniaxialMaterial Elastic 1 $k + element truss 1 1 2 1.0 1 -rho 2.0 + + # Initial displacement u0, zero velocity -- pure free vibration, no load. + setNodeDisp 2 1 $u0 -commit + + integrator ExplicitBathe 0.5 0 -lnvd $alpha + + constraints Plain + numberer Plain + system Diagonal + test NormUnbalance 1.0e-12 10 + algorithm Linear + analysis Transient + + set dt 0.01 + set nSteps 400 + set prevV 0.0 + set peaks {} + set maxAbsTail 0.0 + set tailStartStep 300 + for {set i 1} {$i <= $nSteps} {incr i} { + analyze 1 $dt + set d [nodeDisp 2 1] + set v [nodeVel 2 1] + if {$prevV >= 0 && $v < 0} { lappend peaks $d } + set prevV $v + set ad [expr {abs($d)}] + if {$i >= $tailStartStep && $ad > $maxAbsTail} { set maxAbsTail $ad } + } + return [list $peaks $maxAbsTail] +} + +puts "============================================================================" +puts " ExplicitBathe -lnvd -- Cundall/FLAC local damping, SDOF free vibration" +puts "============================================================================" + +# --- 1) Quantitative check of the small-alpha asymptote, alpha = 0.1 ------- +set res [runSDOF 0.1] +set peaks [lindex $res 0] +set p1 [lindex $peaks 0] +set p2 [lindex $peaks 1] +set gotRatio [expr {$p2/$p1}] +set expRatio [expr {exp(-2.0*0.1)}] +set err1 [expr {abs($gotRatio-$expRatio)/$expRatio*100.0}] +puts [format "alpha=0.1: peak1=%.6e peak2=%.6e ratio got=%.6f expected=exp(-2a)=%.6f error=%.3f%%" \ + $p1 $p2 $gotRatio $expRatio $err1] +if {$err1 < 1.0} { + puts " PASS: within 1% of the theoretical small-alpha decay ratio" +} else { + puts " FAIL: deviates by more than 1% from theory" +} + +# --- 2) Amplitude after ~3-4 s of free vibration for alpha = 0, 0.3, 0.8 ---- +puts "" +puts [format "%-6s %18s %18s %10s" "alpha" "max|u|, t in \[3,4\]s (got)" "exp(-2*a*N) * u0 (hand-calc)" "ratio"] +set omega [expr {sqrt(100.0/1.0)}] +set T [expr {2*3.14159265358979/$omega}] +foreach alpha {0.0 0.3 0.8} { + set res [runSDOF $alpha] + set maxAbsTail [lindex $res 1] + # Hand-calc: small-alpha asymptote evaluated at the midpoint of the + # [3,4] s tail window (t = 3.5 s), N = t/T cycles elapsed. + set tMid 3.5 + set N [expr {$tMid/$T}] + set handCalc [expr {0.01*exp(-2.0*$alpha*$N)}] + set ratio [expr {$maxAbsTail/$handCalc}] + puts [format "%-6s %18.6e %18.6e %10.4f" $alpha $maxAbsTail $handCalc $ratio] +} + +puts "" +puts "Note: the 'ratio' column above is got/hand-calc using the small-alpha" +puts "asymptote exp(-2*alpha*N); it is close to 1 for alpha=0 (no damping, as" +puts "expected) and increasingly diverges for alpha=0.3 and alpha=0.8, exactly" +puts "as for the ExplicitDifferenceStatic integrator's own -alpha option: the" +puts "asymptote over-predicts the decay rate once alpha is not small, so the" +puts "true amplitude decays somewhat less than the naive estimate -- the" +puts "scheme does not blow up, it simply departs from the first-order theory." diff --git a/source/user/manual/analysis/integrator/ExplicitBathe.rst b/source/user/manual/analysis/integrator/ExplicitBathe.rst index e92bc407..eeabb338 100644 --- a/source/user/manual/analysis/integrator/ExplicitBathe.rst +++ b/source/user/manual/analysis/integrator/ExplicitBathe.rst @@ -2,26 +2,93 @@ ExplicitBathe ------------- -.. function:: integrator ExplicitBathe $p $compute_critical_timestep? -.. csv-table:: - :header: "Argument", "Type", "Description" - :widths: 10, 10, 40 +This command is used to construct an explicit, two-sub-step transient +integrator based on the composite scheme of Noh and Bathe [Noh2013]_, with an +optional Cundall/FLAC-style **local non-viscous damping** term (``-lnvd``) +that can be used to turn the integrator into a dynamic-relaxation / +pseudo-static solver. - $p, |float|, Damping parameter (0 < p < 1, typically 0.54) - $compute_critical_timestep, |int|, Optional flag to compute critical timestep (0 or 1) +.. function:: integrator ExplicitBathe $p <$computeCriticalDt> <-lnvd $alpha> + +``$computeCriticalDt`` and ``-lnvd $alpha`` are both optional and may be +given in either order after ``$p`` (or omitted altogether). + +.. csv-table:: + :header: "Argument", "Type", "Description", "Default" + :widths: 18, 10, 55, 15 + + "$p", "|float|", "Sub-step / numerical-damping parameter, :math:`0 < p < 1` (typically 0.5-0.95; :math:`p=0.5` gives no intrinsic numerical damping)", "required" + "$computeCriticalDt", "|integer|", "Integer flag; any value greater than ``0`` enables it. If enabled, the critical time step is computed once (from the element mass/initial-stiffness eigenproblem) on the first ``newStep``, and every step's :math:`\Delta t/\Delta t_{crit}` ratio is reported", "0" + "-lnvd $alpha", "|float|", "Local non-viscous (Cundall/FLAC) damping coefficient, :math:`0 \le \alpha < 1`, applied to the solved acceleration of each of the two Bathe sub-steps (equivalent to damping the unbalanced force, since the mass is diagonal -- see Theory below)", "0 (off)" + +Available only in the Tcl interpreter (registered in ``SRC/tcl/commands.cpp``); +there is no corresponding entry in the Python/``OpenSeesCommands`` dispatcher, +so this integrator cannot be used from OpenSeesPy. .. note:: - * The method is second-order accurate and explicit. - * Only mass matrix is assembled on RHS (no tangent matrix required). - * Critical time step is approximately 2× larger than central difference method. - * For stability: :math:`\Delta t \leq \frac{2}{\omega_{max}}` - * Typical p values: 0.54-0.95 for optimal numerical damping. + * The method is second-order accurate and explicit, performing two + sub-steps per :math:`\Delta t`: :math:`t \to t+p\Delta t` and + :math:`t+p\Delta t \to t+\Delta t`. + * Only the mass matrix is ever assembled as the tangent + (``formEleTangent``/``formNodTangent`` only add :math:`M`); as with any + explicit scheme this requires a nonsingular (preferably diagonal) mass + matrix on every active DOF, and exactly two ``update()`` calls per step + (a third call is rejected with a warning). + * :math:`p` must satisfy :math:`0 < p < 1` strictly; :math:`p=0` or + :math:`p=1` are rejected at construction. :math:`p=0.5` is the + special case with zero intrinsic numerical damping (:math:`q_0=q_1=0`, + :math:`q_2=0.5`), useful for isolating the effect of ``-lnvd`` from the + scheme's own damping. + * For stability (undamped): :math:`\Delta t \le 2/\omega_{max}`, with a + critical time step that is approximately twice that of plain central + difference for the same mesh. + * Passing a nonzero ``$computeCriticalDt`` (e.g. ``1``) triggers, on the + first time step only, a generalized-eigenvalue (``dggev``) + critical-time-step scan over every element's lumped mass and initial + stiffness, printing the governing (minimum) damped and undamped + :math:`\Delta t_{crit}` and the element tag at which each occurs; + subsequently every ``newStep`` reports the ratio + :math:`\Delta t/\Delta t_{crit,damped}` with an ``[OK]`` / + ``[WARNING: dt > dt_crit!]`` tag. + * An invalid ``-lnvd`` value (:math:`\alpha` outside :math:`[0,1)`, or a + missing value), or any unrecognized option, prints a warning and is + simply ignored -- local damping stays off and the integrator is still + constructed. Only an invalid ``$p`` (outside :math:`(0,1)`) makes + construction fail. Leaving ``-lnvd`` off, or passing ``-lnvd 0.0``, + reproduces the original (undamped-by-this-mechanism) Noh-Bathe scheme. + * This local damping is unrelated to (and not a substitute for) Rayleigh + damping assigned to elements/regions; it is a purely numerical + dissipation mechanism intended for pseudo-static / dynamic-relaxation + analyses, not for genuine dynamic response. + * ``-lnvd`` is applied to the fully assembled acceleration of each + sub-step, including -- correctly -- when the model is partitioned + across processes with one of the parallel diagonal systems of + equations; re-forming or printing the unbalance for diagnostics does + not change the result. + * The equivalence between damping the acceleration and damping the + unbalanced force (see Theory below) relies on a lumped (diagonal) + mass, as is standard for this explicit scheme. If a consistent mass + matrix were instead assembled and solved with a general (non-diagonal) + solver, ``-lnvd`` would still be applied per equation to the resulting + :math:`a=M^{-1}F`, but the per-equation damping term would no longer + correspond to damping the same equation's unbalanced force, and energy + dissipation would no longer be guaranteed. + * See also :ref:`ExplicitDifferenceStatic`, which applies the same kind + of Cundall/FLAC local damping to a plain central-difference (leap-frog) + integrator, but *by default* combines the velocity-sign and + force-rate forms of the law and adds a small velocity dead-band (``-vEps``) to + avoid sign chatter near zero velocity. ``-lnvd`` here always uses only + the simple :math:`F_d=-\alpha|F|\,\mathrm{sign}(v)` law (with + :math:`\mathrm{sign}(0)=0`) and has no dead-band, so it should not be + expected to behave identically to ``ExplicitDifferenceStatic``'s own + damping for the same :math:`\alpha`. Theory ^^^^^^ -The ExplicitBathe method is a two-step explicit integration scheme with built-in numerical damping. The method performs two sub-steps per time step: +The ExplicitBathe method is a two-step explicit integration scheme with +built-in numerical damping. The method performs two sub-steps per time step: 1. First step: :math:`t \rightarrow t + p\Delta t` 2. Second step: :math:`t + p\Delta t \rightarrow t + \Delta t` @@ -36,22 +103,108 @@ The integration coefficients are computed from the damping parameter p: q_0 = -q_1 - q_2 + 0.5 -The method offers enhanced stability compared to standard central difference, with a critical time step approximately twice as large. The parameter p controls numerical damping, with higher values providing more damping but reduced accuracy. +The method offers enhanced stability compared to standard central +difference, with a critical time step approximately twice as large. The +parameter p controls numerical damping, with higher values providing more +damping but reduced accuracy. -.. admonition:: Example +**Local non-viscous (Cundall/FLAC) damping (-lnvd).** When ``alpha > 0``, +after the linear system is solved for the acceleration of *each* of the two +Bathe sub-steps, that solved acceleration :math:`a` is modified +equation-by-equation: - 1. **Tcl Code** +.. math:: - .. code-block:: tcl + a_d = a - \alpha\, |a| \,\mathrm{sign}(v) - integrator ExplicitBathe 0.54 +where :math:`v` is that sub-step's predicted velocity (:math:`v_t+a_t\,p +\Delta t` for the first sub-step, and the corrected +:math:`v_{t+p\Delta t}+a_{t+p\Delta t}(1-p)\Delta t` for the second) -- the +same velocity that was used to set the model's response before that +sub-step's equations were formed and solved. Because the mass used by this +explicit scheme is diagonal (lumped), this is algebraically identical, +equation by equation, to applying Cundall's local damping scheme +[CundallBathe1987]_, as used in FLAC/FLAC3D, directly to the unbalanced force +:math:`F = P - R(u)` (applied load minus assembled resisting force): + +.. math:: + + F_d = -\alpha\, |F| \,\mathrm{sign}(v) - 2. **Python Code** +It removes energy in proportion to the magnitude of the unbalance itself +rather than to velocity, so the fractional damping per cycle is +amplitude-independent. Because the damping acts on the fully assembled +acceleration, it is applied correctly even when the model is partitioned +across processes using a parallel diagonal system of equations; whether the +underlying unbalance is also re-formed or printed for diagnostics has no +effect on the result. - .. code-block:: python +As for :ref:`ExplicitDifferenceStatic`, a first-order (small-:math:`\alpha`) +energy balance for a free-vibration SDOF oscillator (:math:`F=-kU`, so +:math:`F_d=-\alpha|F|\,\mathrm{sign}(v)` is exact there) gives an equivalent +viscous damping ratio :math:`\xi \approx \alpha/\pi` and a logarithmic +decrement per cycle :math:`\delta \approx 2\alpha`, i.e. - integrator('ExplicitBathe', 0.54) +.. math:: + + \frac{A_{n+1}}{A_n} \approx e^{-2\alpha} + +This asymptote is accurate to a fraction of a percent for +:math:`\alpha \lesssim 0.1` and increasingly over-damps relative to the +true response as :math:`\alpha \to 1` (the scheme remains stable; it simply +dissipates energy faster than the first-order estimate predicts). Because it +removes energy proportional to the *size of the unbalance*, regardless of how +far from equilibrium the system starts, ``-lnvd`` with a large :math:`\alpha` +(close to but below 1) is effective as a dynamic-relaxation driver for +pseudo-static problems (gravity application, settling a soil column, +push-over); a small :math:`\alpha` (or omitting ``-lnvd``) should be used for +genuine dynamic/earthquake analyses, where this damping should not +contaminate the computed response. + +.. admonition:: Example + + 1. **Tcl Code** + + .. code-block:: tcl + + integrator ExplicitBathe 0.54 + integrator ExplicitBathe 0.54 1 + # p=0.5 (no intrinsic damping), strong local non-viscous damping + integrator ExplicitBathe 0.5 0 -lnvd 0.8 + + ExplicitBathe is Tcl-only; there is no Python binding to show here (see + the note above). + + 2. **SDOF free-vibration decay example** + + A single-DOF spring-mass oscillator (:math:`k=100`, lumped mass + :math:`m=1` at the free node, via ``element truss 1 1 2 1.0 1 -rho 2.0``) + released from rest at :math:`u_0=0.01`, integrated with :math:`p=0.5` + (no intrinsic numerical damping) so that any amplitude decay is due to + ``-lnvd`` alone. For :math:`\alpha=0.1`, the ratio of the first two + displacement peaks matches the theoretical small-:math:`\alpha` estimate + :math:`e^{-2\alpha}=0.8187` to within 0.04%; running the script with + :math:`\alpha=0, 0.3, 0.8` gives :math:`\max|u|` over :math:`t\in[3,4]\,s` + of about :math:`1.0\times10^{-2}`, :math:`5.1\times10^{-4}` and + :math:`5.9\times10^{-6}` respectively (versus essentially no decay, + moderate decay, and strong decay predicted qualitatively by the theory + above; the exact values increasingly depart from the small-:math:`\alpha` + asymptote as :math:`\alpha` grows, as expected). + + | :download:`example_bathe_lnvd_sdof_decay.tcl <../codeExample/example_bathe_lnvd_sdof_decay.tcl>` **(TCL)**. + + 3. **Pseudo-static relaxation example** + + The same 4-DOF loaded chain used in the :ref:`ExplicitDifferenceStatic` + pseudo-static example, now relaxed with + ``integrator ExplicitBathe 0.5 0 -lnvd 0.8``. 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^{-16}`-:math:`10^{-17}`. + + | :download:`example_bathe_lnvd_pseudostatic.tcl <../codeExample/example_bathe_lnvd_pseudostatic.tcl>` **(TCL)**. .. [Noh2013] Noh, G., & Bathe, K.J. (2013). "An explicit time integration scheme for the analysis of wave propagations." Computers & Structures, 129, 178-193. +.. [CundallBathe1987] Cundall, P.A. (1987). "Distinct Element Models of Rock and Soil Structure." In *Analytical and Computational Methods in Engineering Rock Mechanics*, 129-163. Code Developed by: |jaabell|