Fundamental period of a layered soil profile as a Galerkin eigenvalue problem
solveColumnPeriod.RdSolves the fundamental period of the soil column on rigid base, the site
period \(T_r\), directly from the layers of one realization, and returns
it with its controls: the change under mesh doubling and the
density-weighted Rayleigh period beside it. solveWedgePeriod() solves
the fundamental period of the truncated shear beam with the same numerical method,
and solveProfilePeriods() applies both to a list of realizations.
Arguments
- layers
One layer table, such as an element of the
layerscolumn ofsimulateSoilProfiles(): a data.frame or data.table with one row per layer and numeric columnshs(thickness, m),zm(midpoint depth, m; distinct, it orders the layers),gsat(unit weight, kN/m3; any unit, only the ratios between layers matter) andVSm(m/s) or, whenVSmis absent,Gm(MPa), converted with \(V_S = \sqrt{9.81 \cdot 1000\, G_m / \gamma_{sat}}\) as inevaluateSoilLayers(). Every value must be finite and positive; other columns are ignored.- elements
One positive whole number of linear elements per layer of the declared mesh; the control doubles it.
Value
A one-row data.table: Tr (s), the numerical site period on
2 * elements elements per layer; TrMeshChange, the relative change
\(1 - T_e/T_{2e}\); TrRayleighWeighted (s), the density-weighted
Rayleigh period; and TrRayleighWeightedDifference,
TrRayleighWeighted / Tr - 1. Each value is
a conditional evaluation of the supplied layers; the two controls are
numerical diagnostics, not physical dispersion.
Mathematics
Dobry, Oweis and Urzua (1976) model the profile as a shear beam on rigid
rock (hydrated copy, eq-beam-wave-equation,
eq-beam-rock-condition and eq-beam-surface-condition):
$$\rho \frac{\partial^2 u}{\partial t^2} = \frac{\partial}{\partial z}
\left(G \frac{\partial u}{\partial z}\right),\qquad u(H, t) = 0,\qquad
\left(G \frac{\partial u}{\partial z}\right)_{z=0} = 0,$$
with \(z\) the depth below the surface, \(H\) the depth of the rigid
base and \(\rho = \gamma / g\) (section Basic Equation). For harmonic motion
\(u = X(z) e^{i \omega t}\) (eq-beam-harmonic-motion) it is the regular
Sturm-Liouville problem \(-(G X')' = \omega^2 \rho X\) on \([0, H]\),
whose smallest eigenvalue gives \(T_r = 2\pi/\omega_1\). Multiplying by a
test function \(v\) with \(v(H) = 0\) and integrating by parts, the
surface term vanishes by the free condition and
$$\int_0^H G X' v'\, dz = \omega^2 \int_0^H \rho X v\, dz.$$
On layer \(k\), of thickness \(h_k\), \(G_k = \rho_k V_k^2\) and
\(\rho_k = \gamma_k / g\) are constant; the layers are stacked from the
surface in the order of their midpoint depths zm, and \(H\) is the sum
of their thicknesses.
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.
Rayleigh bound
TrRayleighWeighted is estimateWeightedRayleighPeriod() on the same
layers. It evaluates the Rayleigh quotient of Dobry, Oweis and Urzua
(eq-rayleigh-frequency),
\(\omega^2 = \int \rho V^2 X'^2 dz / \int \rho X^2 dz\), on the
first-iteration deflection of eq-rayleigh-mode-iteration. Let \(u\) be
the piecewise-linear shape whose values \(U_k\) at the interfaces are
that deflection accumulated from the base with the slopes
\(s_k = \sigma_k / (\gamma_k V_k^2)\), \(\sigma_k\) the overburden at
the midpoint of layer \(k\), and \(X_k = (U_{k-1} + U_k)/2\). Summation
by parts gives exactly \(\int \rho V^2 u'^2 dz = g^{-1} \sum_k \gamma_k
h_k X_k\), while \(\int \rho u^2 dz = g^{-1} \sum_k \gamma_k h_k
(U_{k-1}^2 + U_{k-1} U_k + U_k^2)/3 \ge g^{-1} \sum_k \gamma_k h_k X_k^2\).
The estimate's \(\omega^2\) is therefore at least the Rayleigh quotient
of \(u\), which is at least \(\omega_1^2\) by Rayleigh's principle:
TrRayleighWeighted never exceeds the exact period of the layered
column, and TrRayleighWeightedDifference cannot exceed the
discretisation error of Tr, asymptotically TrMeshChange / 3. With a
uniform unit weight it is the
constant-density Method 7 of estimateRayleighPeriod().
The package tests reproduce the uniform layer \(T = 4H/V\)
(eq-uniform-period), the power-law stratum of Dobry, Oweis and Urzua
(eq-power-period, with the height exponent \((2 - p)/2\)), their
two-layer equation (eq-twolayer-period-equation) and the n-layer
frequency equation of Charara (2019, eq-frequency-equation-condensed).