ExplicitBathe: document -lnvd, correct syntax and availability

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]
+}

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."