Fundamental period of the truncated shear beam as a Galerkin eigenvalue problem
solveWedgePeriod.RdSolves the fundamental period \(T_s\) of the truncated inhomogeneous shear beam
directly from the layers of one realization,
and returns it with its controls: the change under mesh doubling and the
period of the analytic route beside it, with their relative difference.
solveColumnPeriod() documents the shared numerical method, and
solveProfilePeriods() applies both to a list of realizations.
Arguments
- layers
One layer table, as for
solveColumnPeriod(), whose thicknesses add up toHs.- Hs
One positive slope height from crest to base, in metres.
- lo
One truncation ratio \(\lambda_o\) in \([0, 1)\); see
buildGeometry().- mo
One inhomogeneity exponent of the analytic route, the
mothatcharacterizeProfile()fitted on the same layers.- VSo
One basal shear-wave velocity of the analytic route in m/s, the
VSoof the same fit.- method
The solution method of the root model of the analytic route:
"exact", or the surrogate"lm"or"nlm"; seebuildRootModel().- elements
One positive whole number of linear elements per layer of the declared mesh; the control doubles it.
Value
A one-row data.table: Ts (s), the numerical fundamental period on
2 * elements elements per layer; TsMeshChange, the relative change
\(1 - T_e/T_{2e}\); TsAnalytic (s), the period of the analytic route;
TsAnalyticDifference, TsAnalytic / Ts - 1; and an, anExact and
anWithinSupport, the root of the analytic route and its diagnostics
from evaluateRoot(). Each value is a conditional evaluation of the
supplied layers and paired inputs; the controls are numerical
diagnostics, not physical dispersion.
Mathematics
Dakoulas and Gazetas (1985) measure the depth \(z\) from the apex where
the extended slopes meet; the crest lies at \(z = h\), the base at
\(z = H\) and the truncation ratio is \(\lambda = h/H\) (hydrated
copy, Fig. 1 on printed page 167 and eq-free-truncation-ratio). The width of
the dam is proportional to \(z\), and the free lateral vibration obeys
(eq-free-motion)
$$\rho \frac{\partial^2 u}{\partial t^2} = \frac{1}{z}
\frac{\partial}{\partial z}\left[G(z)\, z\, \frac{\partial u}{\partial z}\right],$$
with zero displacement at the base and zero shear at the crest
(eq-free-base-condition and eq-free-crest-condition):
\(u(H) = 0\) and \(G z\, \partial u/\partial z = 0\) at \(z = h\).
For the slope of buildGeometry(), of height \(H_s\) from crest to base
and truncation ratio \(\lambda_o\), the extended slopes meet
\(h = \lambda_o H_s / (1 - \lambda_o)\) above the crest and
\(H = H_s / (1 - \lambda_o)\) (see evaluateBeamPeriod()). For harmonic
motion \(u = U(z) e^{i \omega t}\) (eq-free-harmonic), multiplying by
\(z v\) with \(v(H) = 0\) and integrating by parts, the crest term
vanishes by the free condition and
$$\int_h^H G\, z\, U' v'\, dz = \omega^2 \int_h^H \rho\, z\, U v\, dz,$$
a regular Sturm-Liouville problem with weight \(z\) whose smallest
eigenvalue gives \(T_s = 2\pi/\omega_1\). The layers of the realization
are stacked from the crest, \(z = h\) plus the depth below it, in the
order of their midpoint depths, with \(G_k = \rho_k V_k^2\) and
\(\rho_k = \gamma_k / g\) constant in each; their thicknesses must add up
to \(H_s\).
With \(\lambda_o = 0\) the crest is the apex, where the width and the
crest term vanish with \(z\) itself. The problem stays well posed in the
energy space: the stiffness and mass of the first element remain positive,
and the branch \(Y_q\) of the general solution (eq-free-solution) has
infinite energy \(\int G z U'^2 dz\) at \(z = 0\), so the Galerkin
solution converges to the regular branch, the paper's untruncated case.
The mesh starts at the apex, with no offset; the package tests reproduce
the untruncated periods of Table 7 of the paper (printed page 172).
Analytic route
TsAnalytic is the closed form of the paper (eq-natural-frequency-period)
evaluated by evaluateBeamPeriod() with the apex-to-base
height \(H_s/(1 - \lambda_o)\), the basal velocity VSo and the first
root of the characteristic relation at mo and lo from
evaluateRoot() with the solution method method. It is exact for the
power law \(G = G_b (z/H)^m\) measured from the apex; with the exponent
that characterizeProfile() fits from the crest it is an approximation.
TsAnalyticDifference (TsAnalytic / Ts - 1) is returned so that a
consumer can flag a divergence under a threshold of its own; no
threshold is applied here. The package tests reduce the numerical period
to the closed form on layered power-law profiles with \(m\) = 0, 0.5,
0.95 and 1.5 and \(\lambda\) = 0.075, 0.2, 0.351 and 0.5, and reproduce
the shear-beam periods of the five dams of Table 8 of the paper (printed
page 175) within 0.2 percent.
Numerical method
The weak form is discretised with continuous piecewise-linear functions on
a mesh that follows the layers: every layer is divided into elements
equal elements, so every interface is a node, and the node at the bottom
of the last layer is fixed. With a width \(w(z)\), 1 for the column and
the distance to the apex for the wedge of solveWedgePeriod(), the
integrals of an element of length \(L\) between an upper node \(a\)
and a lower node \(b\) are exact, because \(G\) and \(\rho\) are
constant within a layer and \(w\) is linear: the element stiffness is
$$k_e = G \frac{w_a + w_b}{2L},$$
assembled as \(k_e (\delta_a - \delta_b)(\delta_a - \delta_b)^T\), and
the consistent mass has diagonal entries \(\rho L (3 w_a + w_b)/12\) and
\(\rho L (w_a + 3 w_b)/12\) and coupling \(\rho L (w_a + w_b)/12\).
The generalized eigenvalue problem \(K \phi = \omega^2 M \phi\) has
\(K\) and \(M\) symmetric positive definite for positive \(G\) and
\(\rho\) with the base fixed. The stiffness factors as
\(K = U^T U\) with \(U = D^{1/2} C\), where \(D\) is the diagonal of
the \(k_e\) and \(C\) the bidiagonal difference matrix (1 on the
diagonal, \(-1\) above it: element \(e\) joins nodes \(e\) and
\(e + 1\), the base node removed); \(U\) is the Cholesky factor of
\(K\), and its inverse is explicit because \(C^{-1}\) is the upper
triangle of ones. The Cholesky reduction gives the symmetric matrix
$$B = U^{-T} M U^{-1} = D^{-1/2} S D^{-1/2},\qquad
S_{ij} = c_{\min(i, j)}\ (i \neq j),\qquad S_{ii} = c_i - M_{i, i+1},$$
with \(c_i\) the consistent mass from the top to node \(i\) (the
cumulative row sums of \(M\)). The eigenvalues of \(B\) are
\(1/\omega^2\), so the fundamental mode is its largest eigenvalue
\(\mu_1\), computed by the direct symmetric solver
eigen(symmetric = TRUE), and \(T = 2\pi\sqrt{\mu_1}\); no iterative
solver or convergence loop is involved. Every entry of \(S\) is a sum
of positive masses, and the largest eigenvalue of a symmetric matrix is
resolved to relative machine precision.
The discretisation is a Rayleigh-Ritz method: the numerical period never
exceeds the exact period of the layered beam, and it does not decrease
when every element is halved, since the trial spaces are nested. The
period returned is the one with 2 * elements elements per layer, and the
mesh change \(1 - T_{e}/T_{2e} \ge 0\) between the declared and the
doubled mesh is its convergence certificate. Within a layer the solution
is smooth and the interfaces are nodes, so the error decreases with the
square of the element length and is asymptotically a third of the mesh
change. A profile whose stiffness tends to zero at the surface converges
more slowly. No threshold is applied to the mesh change: it is returned
for the consumer.
Examples
Layers <- data.frame(hs = rep(1, 10), zm = seq(0.5, 9.5, by = 1), gsat = 19, VSm = 200)
solveWedgePeriod(Layers, Hs = 10, lo = 0.35, mo = 0, VSo = 200, method = "exact", elements = 4)
#> Ts TsMeshChange TsAnalytic TsAnalyticDifference an anExact
#> <num> <num> <num> <num> <num> <lgcl>
#> 1: 0.1649457 3.376545e-05 0.1649476 1.125902e-05 2.930155 TRUE
#> anWithinSupport
#> <lgcl>
#> 1: TRUE