Nonlinear dynamic analysis of dam-foundation interaction under near-fault seismic loading

2.1.1. Mathematical methods for assessing strong ground motion

The adopted approach employs analytical solutions developed by A.S. Bykovtsev [6-11] for complex curvilinear fault systems, including strike-slip and combined slip mechanisms. Synthetic time histories for the Karzhantau fault are shown in Fig. 1. The Pskem Dam in the Tashkent region is used as the case study (computational scheme in Fig. 2). The dam height is 195 m, crest length 1600 m, crest width 12 m, with slope ratios 1:2.4 (upstream) and 1:2 (downstream). The core base width is 130 m, and the reservoir capacity is about 520.000 m³.

Fig. 1Synthetic displacement, velocity, and acceleration records for the Pskem Dam section 1 km from the Karzhantau fault (L= 36 km, W= 15 km, V= 3 km/s; Cp= 6 km/s, Cs= 4 km/s)

Fig. 2Finite-element computational scheme used for seismic analysis of the Pskem Dam

2.1.2. Regional tectonics

The seismic hazard of the dam site is controlled by four major fault systems within 20-40 km: Karzhantau, North-Ugam, Ugam-Maydontol, and Pskem. Their geometric and seismogenic characteristics are described in [12].

To evaluate the mechanics of the North-Ugam fault system, analytical solutions based on invariant J-integrals were applied. These invariants allow estimation of driving forces at the fault tip depending on rupture length, fault orientation, external stress state, and mechanical properties of the surrounding medium.

Fig. 3Geographic location of the Pskem Dam in the Tashkent administrative area

Fig. 4Regional tectonic framework illustrating major faults in the Pskem Dam region (1:500,000)

2.1.3. Computational model formulation

In this study, a plane heterogeneous model consisting of the dam, foundation, and base is accepted [12], [13]. The computational domain includes a deformable volume $\overrightarrow{V}={\overrightarrow{V}}_1+{\overrightarrow{V}}_2+{\overrightarrow{V}}_3+{\overrightarrow{V}}_4$ and a semi-infinite foundation medium Fig. 5, where all components may exhibit viscoelastic behavior, and the physical properties of their components differ. Body forces and surface loads are included to determine displacements and stresses in the heterogeneous system under dynamic loading. The problems under consideration are posed for a finite domain Fig. 5 of volume $\overrightarrow{V}+{\overrightarrow{V}}_5$ (${\overrightarrow{V}}_5$ is the volume isolated from a half-space) and bounded by surfaces on which conditions of non-reflective impact are set ${\mathrm{\Sigma }}_1^{-}+{\mathrm{\Sigma }}_1^{+}+{\mathrm{\Sigma }}_2^{+}$.

The dynamic process is described using the principle of virtual displacements:

Fig. 5Computational model of the heterogeneous system

Constitutive relations connect stress tensor ${\sigma }_{ij}$ and strain tensor ${\epsilon }_{ij}$ [13]:

For an elastic material of the $n$-th element of the system, quantities $\lambda$ and $\mu$ are the Lamé constants; for a viscoelastic material, they are Volterra integral operators [13]:

For viscoelastic materials, the Lamé constants are replaced by Volterra integral operators $\overrightarrow{u}$:

The system also incorporates Cauchy relations Eq. (4) and non-reflecting boundary conditions Eq. (5) to ensure Rayleigh wave propagation across the finite domain boundaries:

Eq. (5) represents the radiation condition imposed at the artificial boundaries of the truncated computational domain, allowing outgoing waves to propagate without reflection. Displacement components, stress and strain tensors, density parameters, and relaxation operators are defined according to the viscoelastic formulation described above. The subscripts denote element numbering and spatial coordinates within the heterogeneous domain.

2.2.1. Free vibration analysis

The natural dynamic characteristics of the system are obtained by solving the generalized eigenvalue problem:

where [$M$] denotes the global mass matrix and [$K$] represents the effective stiffness matrix incorporating wave energy radiation effects at the artificial boundaries of the truncated computational domain. The frequency-domain representation of Eq. (3) is:

with corresponding kernel transformations:

where ${\mathrm{\Gamma }}_{{\lambda }_n}^S$, ${\mathrm{\Gamma }}_{{\lambda }_n}^C$, ${\mathrm{\Gamma }}_{{\mu }_n}^S$, ${\mathrm{\Gamma }}_{{\mu }_n}^C$ – are the Fourier sine and cosine images of kernels ${\mathrm{\Gamma }}_{{\lambda }_n}\left(\tau \right)$, ${\mathrm{\Gamma }}_{{\mu }_n}\left(\tau \right)$.

Consequently, the eigenvalues in Eq. (6) become complex. The real part defines the oscillation frequency, while the imaginary part represents damping and radiation losses. This approach directly evaluates dissipation within the dam-foundation system. The resulting nonlinear system is solved via iterative root searching and matrix factorization.

2.2.2. Harmonic response analysis

For periodic external excitation, the governing equations are rewritten in the frequency domain, leading to the complex amplitude formulation:

where $\mathrm{\Omega }$ is the excitation frequency, $\left\{X\right\}$ represents the response amplitudes, $\left\{f\right\}$ denotes periodic load components, and $\left\{F\right\}$ includes the complete loading vector comprising hydrostatic and body forces. To simulate transient seismic loading typical of near-fault earthquakes, the dynamic equilibrium equations are solved in the time domain:

with initial conditions:

Damping matrix $\left[C\right]$ includes both material dissipation effects and boundary energy radiation contributions. Time integration is performed using the implicit Newmark-$\beta$ method [13], which provides stable solutions for structural dynamic problems.