Next Article in Journal
Optimal Sensor Placement and Multimodal Fusion for Human Activity Recognition in Agricultural Tasks
Previous Article in Journal
Weight Assignment Method and Application of Key Parameters in Shale Gas Resource Evaluation
Previous Article in Special Issue
Graph-Analytical Method for Calculating Settlement of a Single Pile Taking into Account Soil Slippage
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Implementation of PMDL and DRM in OpenSees for Soil-Structure Interaction Analysis

1
Department of Civil Engineering, Kutahya Dumlupinar University, Kutahya 43100, Turkey
2
Department of Civil Engineering, Yildiz Technical University, Istanbul 34220, Turkey
*
Author to whom correspondence should be addressed.
Appl. Sci. 2024, 14(18), 8519; https://doi.org/10.3390/app14188519
Submission received: 12 August 2024 / Revised: 14 September 2024 / Accepted: 19 September 2024 / Published: 21 September 2024
(This article belongs to the Special Issue Soil-Structure Interaction in Structural and Geotechnical Engineering)

Abstract

:
It is widely acknowledged that the effects of soil-structure interaction (SSI) can have substantial implications during periods of intense seismic activity; therefore, accurate quantification of these effects is of paramount importance in the design of earthquake-resistant structures. The analysis of SSI is typically conducted using either direct or substructure methods. Both of these approaches involve the use of numerical models with truncated or reduced-order computational domains. To ensure effective truncation, it is crucial to employ boundary representations that are capable of perfectly absorbing outgoing waves and allowing for the consistent application of input motions. At present, such capabilities are not widely available to researchers and practicing engineers. In order to address this issue, this study implemented the Domain Reduction Method (DRM) and Perfectly Matched Discrete Layers (PMDLs) in OpenSees. The accuracy and stability of these implementations were verified through the use of vertical and inclined incident SV waves in a two-dimensional problem. In terms of computational efficiency, PMDLs require a shorter analysis time (e.g., with PMDLs, the analysis concluded in 35 min as compared to 250 min with extended domain method) and less computational power (one processor for PMDLs against 20 processors for the extended domain method) thus offering a balance between accuracy and efficiency. Furthermore, illustrative examples of the aforementioned implemented features are presented, namely the response analysis of single-cell and double-cell tunnels exposed to plane waves inclined at an angle.

1. Introduction

Infinite domains (unbounded domains) are common in various research fields, especially in wave propagation studies. One particularly complex application in this area is soil–structure interaction (SSI). To tackle this, the soil model is split into two domains: a bounded (interior) domain, and an unbounded (exterior) domain (see Figure 1a). Both the structure and the bounded soil can display nonlinear behavior, while the unbounded soil is usually assumed to be linear. Precise solutions are typically needed only within the interior domain. Therefore, the unbounded soil is replaced by special boundary conditions (see Figure 1b) designed to absorb energy entering the medium and, hence, prevent it from reflecting back into the interior. These boundary conditions are known as absorbing boundary conditions (ABCs).
Different ABCs (e.g., viscous boundaries, transmitting boundaries, non-reflecting boundaries, infinite elements, perfectly matched layers (PMLs), perfectly matched discrete layers (PMDLs)) are used to mimic the exterior domain. Viscous, transmitting, and non-reflecting boundaries are utilized in earthquake engineering simulations to absorb seismic waves, preventing them from reflecting back into the model and artificially amplifying the outcomes. These boundaries are crucial for accurately simulating open or semi-infinite domains by reducing artificial reflections and providing more realistic results [1,2,3,4]. Lysmer and Kuhlemeyer [5] introduced the first local ABC, designed to absorb waves traveling in a specific direction, whereas Lysmer and Drake [6] developed a technique to analyze the propagation of Love waves through structures with non-horizontal layers.
Infinite elements are a type of finite element utilized in numerical simulations to represent problems that extend indefinitely, such as those involving wave propagation or electromagnetic fields. These elements mimic infinite domains by incorporating absorbing boundary conditions, which reduce wave reflections and create a seamless transition from the finite computational area to the infinite region. This method accurately represents far-field behavior without needing a large computational domain, improving efficiency in simulations of unbounded or semi-infinite domains [7,8,9,10,11,12].
A particular type of continuous ABC, designated as PML, entails the substitution of the exterior with an absorbing medium (PML region) that exhibits precise impedance matching with the interior [13]. This substitution is terminated at a specific distance. This perfect impedance matching prevents reflections at the interface, while attenuation in the PML region serves to minimize reflections caused by truncation. The popularity of PML can be attributed to its flexibility and generality as a continuous ABC. Nevertheless, recent research indicates that PML may not be as effective as rational ABCs [14,15].
Pettigrew et al. introduced a framework for implementing a new boundary condition that combines a perfectly matching layer with infinite elements, designed for addressing unbounded elastic wave problems in the time domain [16]. Basu and Chopra presented time-harmonic governing equations for PMLs applicable to anti-plane and plane-strain motions in both elastic and visco-elastic media [17]. Zhang and Taciroglu [18] developed two- and three-dimensional viscoelastic PMLs that account for Rayleigh damping effects, which were then implemented in ABAQUS [19]. Josifovski combined PMLs with frequency dependent finite element formulation [20]. Zhang et al. also developed a PML-based toolbox in ABAQUS (2017) for SSI analysis [21]. The key point of PMLs is the use of complex coordinate stretching, which employs complex-valued functions to enhance wave absorption and minimize boundary reflections in numerical simulations of wave propagation. Kausel and de Oliveira Barbosa developed a coordinate stretch formulation for PMLs and applied it solely to two-dimensional wave propagation problems in homogeneous media, specifically for SH and SV wave propagation [22]. Harari and Albocher [23] applied complex coordinate stretching to modify the material matrix, while Esmaeilzadeh Seylabi et al. [24] and Fontara et al. [25] used it to alter the strain-displacement matrix.
The research on ABCs has resulted in the creation of a new class of continuous ABCs, designated as PMDLs. These PMDLs serve to bridge the gap between PML, and rational ABCs as reported by [26,27,28], among others. In essence, PMDLs represent a specifically discretized version of PML, with the consequence that they are equivalent to rational ABCs. This, in particular, is achieved by employing linear finite elements with midpoint integration for discretization [29,30]. It has been demonstrated that, when approximating the Dirichlet-to-Neumann (DtN) map, the error from midpoint integration cancels out the discretization error, resulting in a perfect match even after exterior discretization. This phenomenon has led to the designation of these discrete layers as “perfectly matched discrete layers,” or PMDLs. The sole source of error in the DtN map is the truncation of PMDLs, which diminishes exponentially with the increase in the number of layers. Initially, PMDLs were developed with the objective of efficiently implementing rational ABCs [27], and they have been demonstrated to be mathematically equivalent to rational ABCs [27]. Given their close connection to both PML and rational ABCs, PMDLs can be regarded as a unification of the two, combining the respective advantages of flexibility and efficiency.
Various studies have validated the effectiveness of PMDLs for different situations e.g., static analysis, and scalar and harmonic wave propagation (see, among others, [26,30,31]). Nguyen and Kim [32] introduced a modified version of PMDLs for SSI analysis in the time-domain, whereas Lee et al. employed PMDLs in layered media for nonlinear SSI analysis in the time-domain [33]. Lee [34] proposed a numerical method incorporating PMDLs for layered poroelastic soil and further developed three-dimensional PMDL formulations as well [35]. Building on prior work, Nguyen and Tassoulas applied the Reciprocal Absorbing Boundary Condition (RABC) in conjunction with PMDLs to facilitate transient analysis of in-plane wave motion within a layered half-space [36]. Nguyen et al. applied PMDLs to the analysis of concrete gravity dams [37].
OpenSees, unlike commercial finite element software, is available for free, which boosts its popularity. Whereas several successfully implemented and widely used models for superstructures (e.g., reinforced concrete, steel, masonry, timber) are available in the OpenSees library, the same cannot be said for soil-structure interactions, where the sheer size of the domain and complexity of the problems render the use of currently available models/methods inefficient in terms of the analysis time and computational resources required. As models used to analyze soil-structure interactions become more complex, analysis times increase significantly. A pertinent example is that of performing seismic fragility analysis, which requires a very high number of simulations to account for the inherent uncertainties and variations in both earthquake events and structural responses.
Therefore, the integration of PMDL and Domain Reduction Method (DRM) methods enhances OpenSees as a powerful tool for seismic analysis and in obtaining better results in less time compared to other boundary conditions discussed in the article. Hence, the addition of PMDLs and DRM to OpenSees facilitates earthquake engineering research by providing researchers with an open-source tool to analyze soil-structure interactions.
In this study, we utilize the PMDL formulation proposed by Nguyen and Kim [32] by developing a custom element in OpenSees (3.5) [38] using the Extended Newmark-β implicit time integration scheme for simulating wave propagation in two-dimensional (2-D) half-spaces. The implementation was verified through a comparison of results derived from the PMDL-truncated domain and those obtained from an extended domain with fixed boundaries, resulting in a notable degree of consistency.
Furthermore, direct modelling of SSI problems frequently encounters difficulties associated with defining appropriate input ground motions [39,40]. The DRM is currently regarded as one of the most effective approaches for simulating semi-infinite domains subjected to far-field excitations [41,42,43]. In this investigation, the DRM was initially employed to calculate the effective nodal forces, utilizing the formulations proposed by [42]. Subsequently, a comprehensive verification process was conducted, in which the DRM-generated free-field responses were compared with either analytical solutions or the results obtained from one-dimensional (1-D) site response analysis (SRA). The findings demonstrate exceptional accuracy for vertical and inclined incident SV waves in a two-dimensional domain with a homogeneous soil layer.
Finally, using the combined DRM-PMDL system, we modeled tunnels (i.e., single-cell tunnel and double-cell tunnel) in a 2-D domain to investigate the effects of the incidence angle. We applied inclined SV waves at various incidence angles and plot accelerations of the corner points, maximum bending moments, drift ratios, rocking rotations, spectral amplifications, and maximum deformations of the tunnels.

2. Perfectly Matched Discrete Layers

In order to simulate the extensive domain of soils, we employed PMDLs, a robust wave-absorbing boundary, to truncate it. The PMDL effectively mitigates reflections at the truncated boundary for waves, regardless of their angle of incidence, thus ensuring rapid attenuation of wave energy within its zones. As in a PML, because there is a perfect match between the domain of interest and the PMDL, waves from the domain pass through the interface without reflection and are damped within the PMDL. The small numerical reflections that occur due to the discretization of the domain are absorbed as they propagate back, minimizing interference with the wavefield, particularly reflections from topographic surfaces. This feature makes PMDLs highly effective in simulating realistic boundary conditions in SSI problems, ensuring that the analysis is not compromised by artificial boundary effects. The PMDL efficiently absorbs both propagating and transient waves due to its incorporation of both real and imaginary lengths. In this work, we implemented 2-D PMDLs introduced by [32] which can be utilized for wave-propagation simulations in both homogeneous and non- homogeneous media.

2.1. Implementation

The half-space problem illustrated in Figure 2a, can be represented through a combination of a domain of interest and an absorbing region like PMDLs (see Figure 2b). According to theoretical principles, PMDLs comprise layers with both real and imaginary lengths. Consequently, four distinct types of PMDL elements are identified as shown in Figure 2c. The first type is the regular PMDL element, characterized by real lengths in both the x and z directions. The second type, the PMDLXX element, has imaginary length in the x direction, and real length in the z direction. The third type, the PMDLZZ element, features real dimensions in the x direction and imaginary dimensions in the z direction. Lastly, the PMDLXZ element possesses imaginary length in both the x and z directions. The PMDLs are modeled using these four types of elements. The top surface of the domain of interest (Figure 2b) can be modeled with complex topography, while all other surfaces related to the PMDLs should be represented as flat surfaces. Additionally, the domain of interest is modeled using conventional four-node quadrilateral elements. The real and imaginary lengths of PMDLs in the x and z directions can be calculated using (1)–(2).
L x j R = 4 H 2 j 1   j = 1,2 , 3 , m L x k I = 2 i V x ω c o s θ x k   k = 1,2 , 3 , n θ x k = π ( k 1 ) 2 n
L z j R = 4 H 2 j 1 2 1 ν 1 2 ν 1 2   j = 1,2 , 3 , m L z k I = 2 i V z ω c o s θ z k   k = 1,2 , 3 , n θ z k = π k 1 2 n
In Equations (1) and (2), L x j R and L z j R refer to the lengths of the j-th real layer in the x- and z-directions, respectively. Similarly, L x k I and L z k I represent the lengths of the k-th imaginary layer in the x- and z-directions. The variables V x and V z indicate the reference velocities in the x- and z-directions. The parameters m and n correspond to the number of real and imaginary layers, respectively. H stands for the height of the interior domain, while ω and ν represent the excitation frequency and Poisson’s ratio, respectively. Figure 3a,b illustrates the global and local coordinates for the PMDLs. Additionally, the dynamic stiffness matrix of the element is calculated using (3), as outlined by [44,45].
S = K + i ω C ω 2 M
The stiffness and mass matrices of the PMDL element shown in Figure 3, can be written as follows:
K = 1 1 1 1 b a B r r T D r r B r r + a b B s s T D s s B s s + B r s T D r s B r s + B s r T D s r B s r d r d s
M = 1 1 1 1 a b ρ N T N d r d s
The strain-displacement matrices ( B r r , B r s , B s r , B s s ), elasticity matrices ( D r r , D r s , D s r , D s s ), shape functions ( N i ), and the shape function matrix ( N ) used herein are provided in Appendix A. In this research, the soil is considered to be elastic, and Rayleigh damping is utilized with a 2% damping ratio to represent a typical soil medium, as defined in [46] and described in (6).
C = α M + β K
Here, α and β are the Rayleigh damping parameters that are proportional to mass and stiffness, respectively. In the case of transient waves, the dimensions of the PMDL elements are real. Consequently, the dynamic stiffness matrix S in (3) is formulated using standard finite elements but with the midpoint integration rule. In contrast, for the propagation of waves, the dimensions of the PMDL elements are imaginary. In order to employ the PMDL elements in the time-domain, it is necessary to eliminate the imaginary lengths. From (3) to (5), using the values for the reference wave velocity, denoted V x   , and the equivalent length, represented by 2 a = 2 i V x / ( ω c o s θ x k ) , the dynamic stiffness of a PMDLXX element may be given by the following expressions:
S x = K x + i ω C x ω 2 M x + 1 i ω R x
K x = 1 1 1 1 B r s T D r s B r s + B s r T D s r B s r d r d s + α 1 1 1 1 b V x cos θ x k ρ N T N d r d s + β 1 1 1 1 V x b cos θ x k B s s T D s s B s s d r d s
C x = 1 1 1 1 b cos θ x k V x B r r T D r r B r r d r d s + β 1 1 1 1 B r s T D r s B r s + B s r T D s r B s r d r d s + 1 1 1 1 b V x cos θ x k ρ N T N d r d s
M x = 1 1 1 1 b cos θ x k V x B r r T D r r B r r d r d s
R x = 1 1 1 1 V x b cos θ x k B s s T D s s B s s d r d s
For a PMDLZZ element, using Equations (3)–(5), the values for the reference wave velocity V z and the equivalent length 2 b = 2 i V z / ( ω c o s θ z k ) , the dynamic stiffness of a PMDLZZ element may be given by the following expression:
S z = K z + i ω C z ω 2 M z + 1 i ω R z
K z = 1 1 1 1 B r s T D r s B r s + B s r T D s r B s r d r d s + α 1 1 1 1 a V z cos θ z k ρ N T N d r d s + β 1 1 1 1 V z a cos θ z k B r r T D r r B r r d r d s
C z = 1 1 1 1 a cos θ z k V z B s s T D s s B s s d r d s + β 1 1 1 1 B r s T D r s B r s + B s r T D s r B s r d r d s + 1 1 1 1 a V z cos θ z k ρ N T N d r d s
M z = 1 1 1 1 a cos θ z k V z B s s T D s s B s s d r d s
R z = 1 1 1 1 V z a cos θ z k B r r T D r r B r r d r d s
For a PMDLXZ element, should at least one of the dimensions, 2a or 2b, be real, the methodology for determining the dynamic stiffness is analogous to that employed for edge elements (PMDLXX and PMDLZZ). In the event that both dimensions 2 a = 2 i V x / ( ω c o s θ x k ) and 2 b = 2 i V z / ( ω c o s θ z k ) are imaginary, the dynamic stiffness must be revised as follows:
S x z = K x z + i ω C x z + 1 i ω R x
K x z = 1 1 1 1 V z cos θ x k V x cos θ z k B r r T D r r B r r + V x cos θ z k V z cos θ x k B s s T D s s B s s + B r s T D r s B r s + B s r T D s r B s r d r d s + 1 1 1 1 V x V z cos θ x k cos θ z k ρ N T N d r d s
C x z = β 1 1 1 1 V z cos θ x k V x cos θ z k B r r T D r r B r r + V x cos θ z k V z cos θ x k B s s T D s s B s s + B r s T D r s B r s + B s r T D s r B s r d r d s
R x z = α 1 1 1 1 V x V z cos θ x k cos θ z k ρ N T N d r d s
In Equations (3)–(20), a and b denote the half-width and half-height of the rectangular PMDL element. S, K, M, C, and R represent the dynamic stiffness, stiffness, mass, damping, and auxiliary matrices, respectively, while ρ represents the material density. PMDLs should be modeled with elastic material, while the interior domain (Figure 2b), along with the substructure and superstructure, can be represented using either elastic or nonlinear materials. A key feature of the PMDL elements is the use of midpoint integration to eliminate discretization errors. As shown in Figure 4, the model that includes PMDL elements utilizes four types of integration rules: 2 × 2 integration for finite elements within the interior domain, 1 × 2 integration for PMDL elements in the X direction, 2 × 1 integration for PMDL elements in the Z direction, and 1 × 1 integration for corner PMDL elements. Table 1 illustrates the abscissas and weights of the numerical integration methodologies.
Guddati and Lim [26] presented the equation of motion for the entire domain, inclusive of the PMLDs, as illustrated in the following Equations (21) and (22):
M u ¨ ( t ) + C u ˙ ( t ) + K u ( t ) + R u ¯ ( t ) = F e x t ( t )
u ¯ ( t ) = u ( t ) d t
Here, F e x t ( t ) denotes the external force vector. Numerical integration methods such as Central Difference, Hilber-Hughes-Taylor and Generalized Alpha are the numerical integration methods used to solve the equation of motion in OpenSees. The classical Newmark-β method is one of these methods which, however, is insufficient to integrate (21). As a result, the Extended Newmark-β method is used to solve (21). To address (23), it is necessary to calculate and update the stiffness matrix and force vector at each step, since their relations are defined incrementally in (24)–(25).
K e f f Δ u = F R
K e f f = K + 1 β ¯ Δ t 2 M + γ ¯ β ¯ Δ t C + α ¯ Δ t β ¯ R
F R = F e x t t + Δ t K u t + Δ t 1 β ¯ Δ t 2 M + γ ¯ β ¯ Δ t C + α ¯ Δ t β ¯ R Δ u + 1 β ¯ Δ t M 1 γ ¯ β ¯ C + Δ t 2 β ¯ 2 α ¯ 2 β ¯ R u ˙ t + 1 2 β ¯ 1 M Δ t 1 γ ¯ 2 β ¯ C Δ t 3 β ¯ 3 α ¯ 6 α β ¯ ¯ R u ¨ t u ¯ t + Δ t u t R
where, F R ,   K e f f ,   Δ u ,   Δ t denote the residual force vector, effective stiffness matrix, incremental displacement vector, and time-step, respectively. The time increment values ( α ¯ , β ¯ , γ ¯ ) for the Extended Newmark-β method are provided in Table 2 with the expressions for acceleration ( u ¨ ), velocity ( u ˙ ), and the time integration of displacements ( u ¯ ) presented through (26), (27), and (28), respectively.
u ¨ t + Δ t = 1 β ¯ Δ t 2 Δ u 1 β ¯ Δ t u ˙ t 1 2 β ¯ 1 u ¨ t
u ˙ t + Δ t = γ ¯ β ¯ Δ t Δ u + 1 γ ¯ β ¯ u ˙ t + Δ t 1 γ ¯ 2 β ¯ u ¨ t
u ¯ t + Δ t = u ¯ t + Δ t u t + α ¯ Δ t β ¯ Δ u + Δ t 2 β ¯ 2 α ¯ 2 β ¯ u ˙ t + Δ t 3 β ¯ 3 α ¯ 6 α ¯ β ¯ u ¨ t

2.2. Verification of PMDL

This section provides a thorough evaluation of the verification and stability of the implemented PMDL. During the verification process, the performance of the PMDL is compared to that of the built-in infinite element in ABAQUS. According to the ABAQUS manual, the infinite elements utilize a damping and reduced stiffness matrix. However, they do not perfectly transmit energy out of the mesh, except when plane body waves strike the boundary perpendicularly in medium.
A half-space scenario (Figure 5a) is analyzed, with the relevant model properties detailed in Table 3. The excitation is provided by a vertical point load (see Figure 5b), which is represented by a Ricker pulse (see Figure 6). The amplitude, dominant frequency, and total duration of the Ricker pulse are 100 kN, 5 Hz, and 3 s, respectively. Two configurations are employed for the analysis: (i) PMDL elements (see Figure 5b), and (ii) infinite elements (replacing PMDL with CINPE4 [19]). Furthermore, an extended domain with dimensions of 1500 m by 750 m is constructed for the purpose of establishing a reference solution. The outer surface of the established reference solution is fixed at the horizontal and vertical directions.
Figure 7 illustrates the comparisons of displacements obtained using the PMDL, infinite elements, and the extended domain for selected points. For points P1 and P2, only vertical displacements are shown, while for points P3 and P4, both horizontal and vertical displacements are presented. As depicted, the PMDL boundary performs excellently across all selected points, whereas the infinite elements do not perform as well. The infinite elements fail to perfectly absorb outgoing waves, leading to reflected waves re-entering the domain of interest. This results in discrepancies between the results from the domain using infinite elements and those from the PMDL and extended domain. These discrepancies are particularly evident in the vertical displacement at point P3. Additionally, Figure 8 displays contour plots of normalized displacement magnitudes at various times (t = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 2.0, 3.0 s) using the PMDL boundary, further demonstrating its effectiveness in absorbing all outgoing waves.
A simulation with a total duration of 10 s was performed to assess the long-term stability of the PMDL. Figure 9 shows the displacement histories for the four selected points, whose locations are indicated in Figure 5b. The results reveal that the entire domain falls silent after approximately t = 1 s, with no reflections occurring between the PMDL boundary and the domain of interest.
Another method to demonstrate the efficiency of the PMDL is by comparing the total computational time. For the extended domain using the ABAQUS solver on 20 processors, the total computational time is approximately 250 min. For the domain using infinite elements with the ABAQUS solver on 20 processors, it takes about 12 min. Meanwhile, the PMDL domain using the OpenSees solver on a single processor requires around 35 min. These comparisons indicate that the PMDL is highly efficient in terms of the computational power required for such an analysis.

3. Domain Reduction Method

The DRM is a technique used to simplify the analysis of wave propagation problems by reducing the computational domain [42,43]. It involves truncating the domain and applying ABCs to simulate the effects of the surrounding, unmodeled region. The DRM is useful for managing complex or large-scale problems, such as those involving seismic, acoustic, or electromagnetic waves, by focusing computational resources on the area of interest. The advantages of the DRM include a reduction in computational cost and resources, as well as the ability to conduct detailed analysis of specific regions. Conversely, in addition to the aforementioned benefits, the precision of the DRM is contingent upon the efficacy of the boundary conditions, and the deployment of the DRM can be a challenging process. In general, the DRM is an effective method for solving wave propagation problems while accounting for the influence of unmodeled regions.
In the framework of the DRM, incoming waves from the unmodeled area are converted into corresponding nodal forces using Equation (29). The aforementioned forces are then applied to the single-layer elements situated between the ABCs and the domain of interest, which is illustrated in Figure 10a.
P e f f = P i e f f P b e f f P e f f = 0 M b e Ω + u ¨ e 0 C b e Ω + u ˙ e 0 K b e Ω + u e 0 M e b Ω + u ¨ b 0 + C e b Ω + u ˙ b 0 + K e b Ω + u b 0
The subscripts e ,   i ,   b denote the nodes between the ABCs layer and the DRM layer, nodes within the domain of interest, and nodes between the DRM layer and the domain of interest, respectively, and where u 0 ,   u ˙ 0 ,   a n d   u ¨ 0 represent the free-field displacement, velocity, and acceleration, respectively. To utilize the coupled PMDL-DRM method in OpenSees, the main steps are as follows:
  • Define the physical domain geometry, identify the boundaries for PMDL and DRM application, and generate the mesh.
  • Extract the nodes and elements corresponding to the DRM layer from the mesh.
  • Use the flowchart given in Appendix A (see Figure A1) to compute the mass, stiffness, and damping matrices ( M Ω + , C Ω + ,   K Ω + ) for the DRM layer.
  • Calculate the free-field responses (displacement, velocity, and acceleration) for the nodes in the DRM layer using 1D SRA or analytical methods.
  • Determine the DRM layer forces ( P e f f ) using Equation (29).
  • Create an OpenSees model and input the calculated DRM layer forces into it.
The flowchart illustrating the sequential steps to implement and use the coupled PMDL-DRM approach in OpenSees can be found in Appendix A (see Figure A2).
Calculating the free-field response for inclined incident SV waves (refer to Figure 10c) using 1D SRA is extremely difficult, as analytical solutions derived by [47] provided in (30)–(32) exist only for a homogeneous half-space in the time-domain. In Equations (30)–(32), U s i ,   U s r ,   and   U p r are the amplitudes of the incident SV wave, reflected SV wave, and reflected P wave, respectively. V s and V p denote the velocities of S-waves and P-waves. The variable f = f ( t ) refers to a time-dependent source-time function. u x and u z correspond to the displacement components for inclined incident SV waves. θ s and θ p are the angles of SV wave incidence and reflected P-wave, respectively.
u x ( x , z ) u z ( x , z ) = U s i c o s θ s s i n θ s f x V s s i n θ s + z V s c o s θ s + t + U s r c o s θ s s i n θ s f x V s s i n θ s z V s c o s θ s + t + U p r s i n θ p c o s θ p f x V p s i n θ p z V p c o s θ p + t
U s r = sin 2 θ s sin 2 θ p V p V s 2 cos 2 2 θ s sin 2 θ s sin 2 θ p + V p V s 2 cos 2 2 θ s U s i
U p r = 2 V p V s 2 s i n 2 θ s c o s 2 θ s V s s i n 2 θ s s i n 2 θ p + V p V s 2 c o s 2 2 θ s V p U s i
In this investigation, the 2-D DRM was implemented in OpenSees. To validate the 2-D DRM, we conducted a comparative analysis of the results obtained from the DRM model with the analytical solutions for inclined incident SV waves (case A) and the corresponding results from the 1-D site response analysis (case B). The material properties of the soil column utilized in the 1-D SRA are identical to those employed in the DRM model (see Figure 10d). The soil column and the domain of interest of the DRM model have the same height. The properties of the DRM model are outlined in Table 3. For all verification cases, a Ricker pulse is employed as the incident wave function, with parameters set to A R i c k e r = 10 4   m , f R i c k e r = 2   Hz , t 0 = 0.8 s. Here, A R i c k e r ,   f R i c k e r , a n d     t 0 are amplitude, dominant frequency, and time offset of the Ricker pulse. The total simulation time is 3 s, with an increment of Δ t = 0.001 s. The time histories and Fourier amplitude spectra of the Ricker pulse are illustrated in Figure 11.
Figure 12 and Figure 13 show the displacement histories for the selected points (refer to Figure 10d), acquired from the DRM, 1-D SRA, and the analytical solutions. The DRM results closely align with both the analytical solutions and 1-D SRA.
After successfully implementing and validating the PMDL and DRM individually, the integrated PMDL-DRM system was rigorously tested. The testing process involved three distinct models to evaluate the accuracy and reliability of this integrated system:
1. The Extended Domain Model: this model (illustrated in Figure 14a) represents a detailed and comprehensive approach, serving as a reference or benchmark for comparison.
2. The PMDL-DRM Combined Model: depicted in Figure 14b, this model incorporates both the PMDL and DRM methodologies, blending their strengths to achieve a more refined simulation.
3. The Winkler Model: shown in Figure 14c, this simpler model is often used for foundational analysis and offers a different perspective on soil-structure interaction.
The primary focus of the comparison was on the roof displacement of the frame structure, particularly the horizontal displacement at a specific point labeled as “A” in Figure 14a–c. This comparison was crucial in determining whether the PMDL-DRM system could accurately replicate the behavior of the structure under seismic loading.
For the purpose of this verification, it was assumed that both the soil and the frame structure behaved elastically. The soil properties across all three models were kept consistent to ensure a fair comparison, with the specific values provided in Table 3. Additionally, the Young Modulus of material of the frame structure is 32 × 10 6 kPa and section properties of the frame structure are detailed in Figure 14a–c.
To simulate real-world conditions, a historical earthquake record—the acceleration time-history of the Kobe earthquake—was used as the input motion, as illustrated in Figure 14d. This record provided a realistic and challenging test for the models.
The results of the analysis are depicted in Figure 14e. The roof displacements calculated by the PMDL-DRM system were found to be in excellent agreement with those obtained from both the extended model and the Winkler model. This close alignment across different models reinforces the credibility of the PMDL-DRM system, demonstrating that it can effectively simulate the structural response under seismic conditions, thereby validating its implementation.

4. Application

Following the successful incorporation and validation of both DRM and PMDL in OpenSees, this section illustrates an application that employs the integrated DRM-PMDL system for examining the influence of the angle of incidence ( θ ).
As documented in [47,48], the angle of incidence has been demonstrated to exert a considerable influence on the response of tunnels, particularly those with a vertical orientation. This section presents an examination of a single-cell tunnel and a double-cell tunnel in homogeneous media subjected to SV waves at varying angles of incidence. The numerical models are illustrated in Figure 15. The DRM interface is employed to define the input motion, while the PMDL is utilized to truncate the domain and mitigate the effects of scattering waves. The input motion is identical to that utilized for validating the DRM implementation (see Figure 11). Four-node quadrilateral plane-strain elements are employed for the domain of interest, while two-node frame elements are utilized for the tunnel structures, with the latter having a thickness of 0.50 m. The contact between the soil and the tunnel is assumed to be bonded. The material properties for the soil and the tunnels are given in Table 4.
R α = t a n 1 X 2 X 1 h T t a n 1 Z 1 Z 4 w T R β = t a n 1 Z 1 Z 4 w T
In Equation (33), R α ,   R β refer to the drift ratio and rocking rotation, respectively. The shear and rocking deformation (see Figure 16) of the tunnel structure are depicted in Figure 17 and Figure 18. In Figure 16, w T and h T denote the width and height of the tunnel. The method for calculating the drift ratio (shear deformation) and rocking rotation is outlined in (33). Initially, we compared the drift ratios and rocking rotations between single-cell and double-cell tunnel structures. As shown in Figure 17 and Figure 18, for this specific scenario, the drift ratios and rocking rotations increase as the incidence angle decreases. Additionally, the Fourier amplitudes of these drift ratios and rocking rotations are calculated and presented in Figure 17 and Figure 18. These figures illustrate that the Fourier amplitude of the drift ratios and rocking rotations decreases with an increasing incidence angle. Figure 19 displays the peak drift ratios and rocking rotations according to the angle of incidence. It is evident from Figure 19 that the peak drift ratios and rocking rotations decrease with an increasing angle of incidence in both single-cell and double-cell tunnels. Moreover, for a specific incidence angle, the lowest drift ratios are observed in the single-cell tunnel. From Figure 19, it is clear that at incidence angles of 0° and 10°, the largest rocking rotations occur in the double-cell tunnel, whereas at angles of 20° and 30°, the largest rocking rotations are observed in the single-cell tunnel.
Figure 20 illustrates the maximum spectral amplification (as percentages) at the ground surface, caused by the presence of a subsurface tunnel structure. The said percentages are obtained by first calculating the difference between the spectral values when considered with and without tunnel structures, and subsequently dividing these differences by the corresponding spectral values for the case where there is no tunnel structure. As illustrated in Figure 20, the amplification percentages of spectral values for the single-cell tunnel decrease as the angle of incidence increases. For the double-cell tunnel, the amplification percentages of Sd and PSV remain nearly constant. Conversely, the amplification percentages of PSA for the double-cell tunnel increase with the angle of incidence. At a specific angle of incidence, the highest amplification percentages are observed in the single-cell tunnel.
Figure 21 and Figure 22 present the deformation plots at the moment of maximum bending moment for the single-cell and double-cell tunnels, respectively, under various incidence angles. It can be observed that greater bending deformation occurs at smaller incidence angles, while more rotation is evident at larger incidence angles. It should be noted that all figures show the deformations magnified by a factor of 10,000. Figure 23 and Figure 24 display the maximum bending moment diagrams for the single-cell and double-cell tunnels, respectively. It is observed that with an increasing incidence angle, the bending moment decreases. The red circle indicates the location of the maximum bending moment in both Figure 23 and Figure 24. As seen in Figure 24, when the incidence angle is 30°, the location of the maximum bending moment changes differently from the other angles. This is because the drift along the vertical axis increases as the incidence angle of the SV wave increases.
Finally, the accelerations at two selected points on the tunnel (P1 and P2, see Figure 15) are compared. As shown in Figure 25 and Figure 26, the maximum horizontal and vertical accelerations are observed at an incidence angle of 30° for both the single-cell and double-cell tunnels, respectively. Additionally, Figure 25 and Figure 26 indicate that the vertical accelerations increase with the angle of incidence.

5. Conclusions

In this manuscript, we provide a detailed account of the methodology for integrating the DRM and PMDLs into OpenSees for two-dimensional problems, including a revised version of the DRM to accommodate inclined incident SV waves. A comprehensive verification process was conducted, encompassing problems involving both vertical and inclined incident SV waves in a two-dimensional domain. Furthermore, we provided practical examples examining the effects of seismic wave incidence angles on the response of single-cell and double-cell tunnels. From the comparisons and figures given in preceding sections, the findings of this study can be summarized as below:
  • PMDLs improve simulation accuracy in wave propagation by minimizing boundary reflections.
  • In terms of resource efficiency, PMDLs are far more computationally efficient than extended domain methods (e.g., PMDLs requiring a single processor against Extended Domain Model requiring 20 processors for the same problem), thus balancing accuracy and efficiency.
  • Regarding the analysis time required, PMDLs significantly reduce computational time when compared to traditional methods (e.g., with PMDLs, the analysis concluded in 35 min as compared to 250 min when the same problem was solved with Extended Domain Model).
  • The coupled PMDL-DRM technique, when implemented in OpenSees, simplifies advanced SSI analyses in truncated domains.
Overall, this research contributes to the field by providing a robust framework for simulating complex seismic wave interactions in 2D problems, paving the way for more advanced seismic analysis in civil engineering applications. While our study successfully implemented and validated the DRM and PMDLs, it is limited by the specific material properties and loading conditions used herein. Future research could explore different soil properties, loading conditions, and three-dimensional extensions of the coupled PMDL-DRM to further enhance its applicability.

Author Contributions

Conceptualization, S.U. and Y.A.; methodology, S.U. and Y.A.; software, S.U.; validation, S.U.; investigation, S.U. and Y.A.; writing—original draft preparation, S.U. and Y.A.; writing—review and editing, S.U. and Y.A.; visualization, S.U.; supervision, Y.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix A

The strain-displacement matrices can be written as follows:
B r r = N 1 r 0 N 2 r 0 N 3 r 0 N 4 r 0 0 N 1 r 0 N 2 r 0 N 3 r 0 N 4 r
B s s = N 1 s 0 N 2 s 0 N 3 s 0 N 4 s 0 0 N 1 s 0 N 2 s 0 N 3 s 0 N 4 s
B r s = N 1 r 0 N 2 r 0 N 3 r 0 N 4 r 0 0 N 1 s 0 N 2 s 0 N 3 s 0 N 4 s
B s r = N 1 s 0 N 2 s 0 N 3 s 0 N 4 s 0 0 N 1 r 0 N 2 r 0 N 3 r 0 N 4 r
Elasticity matrices for the plane-strain conditions are given as:
D r r = λ + 2 μ 0 0 μ
D s s = μ 0 0 λ + 2 μ
D s r = 0 λ λ 0
D s r = 0 μ μ 0
λ , μ ,   and   ν refer to the Lamé constants and Poisson’s ratio. The shape functions and the shape function matrix can be expressed as follows:
N 1 = 1 2 1 r 1 s
N 2 = 1 2 1 + r 1 s
N 3 = 1 2 1 + r 1 + s
N 4 = 1 2 1 r 1 + s
N = N 1 0 N 2 0 N 3 0 N 4 0 0 N 1 0 N 2 0 N 3 0 N 4
Figure A1. Flowchart for computing the stiffness, mass and damping matrices for the DRM layer.
Figure A1. Flowchart for computing the stiffness, mass and damping matrices for the DRM layer.
Applsci 14 08519 g0a1
Figure A2. Flowchart for utilizing the coupled PMDL-DRM method in OpenSees.
Figure A2. Flowchart for utilizing the coupled PMDL-DRM method in OpenSees.
Applsci 14 08519 g0a2

References

  1. Givoli, D. Non-Reflecting Boundary Conditions. J. Comput. Phys. 1991, 94, 1–29. [Google Scholar] [CrossRef]
  2. Lysmer, J.; Waas, G. Shear Waves in Plane Infinite Structures. J. Eng. Mech. Div. 1972, 98, 85–105. [Google Scholar] [CrossRef]
  3. Tabatabaiefar, H.R.; Fatahi, B. Idealisation of Soil–Structure System to Determine Inelastic Seismic Response of Mid-Rise Building Frames. Soil. Dyn. Earthq. Eng. 2014, 66, 339–351. [Google Scholar] [CrossRef]
  4. Hokmabadi, A.S.; Fatahi, B. Influence of Foundation Type on Seismic Performance of Buildings Considering Soil–Structure Interaction. Int. J. Struct. Stab. Dyn. 2016, 16, 1550043. [Google Scholar] [CrossRef]
  5. Lysmer, J.; Kuhlemeyer, R.L. Finite Dynamic Model for Infinite Media. J. Eng. Mech. Div. 1969, 95, 859–877. [Google Scholar] [CrossRef]
  6. Lysmer, J.; Drake, L.A. The Propagation of Love Waves across Nonhorizontally Layered Structures. Bull. Seismol. Soc. Am. 1971, 61, 1233–1251. [Google Scholar] [CrossRef]
  7. Bakhtaoui, Y.; Chelghoum, A. Solution for Soil–Structure Interaction with Direct Infinite Element in Time Domain. Indian Geotech. J. 2020, 50, 655–663. [Google Scholar] [CrossRef]
  8. Gharti, H.N.; Tromp, J.; Zampini, S. Spectral-Infinite-Element Simulations of Gravity Anomalies. Geophys. J. Int. 2018, 215, 1098–1117. [Google Scholar] [CrossRef]
  9. Kim, D.K.; Yun, C.B. Earthquake Response Analysis in the Time Domain for 2D Soil–Structure Systems Using Analytical Frequency-Dependent Infinite Elements. Int. J. Numer. Methods Eng. 2003, 58, 1837–1855. [Google Scholar] [CrossRef]
  10. Lin, S.-Y.; Hung, H.-H.; Yang, J.P.; Yang, Y.B. Seismic Analysis of Twin Tunnels by a Finite/Infinite Element Approach. Int. J. Geomech. 2017, 17, 04017060. [Google Scholar] [CrossRef]
  11. Yang, Y.B.; Hung, H.H.; Lin, K.C.; Cheng, K.W. Dynamic Response of Elastic Half-Space with Cavity Subjected to P and SV Waves by Finite/Infinite Element Approach. Int. J. Struct. Stab. Dyn. 2015, 15, 1540009. [Google Scholar] [CrossRef]
  12. Yun, C.B.; Kim, D.K.; Kim, J.M. Analytical Frequency-Dependent Infinite Elements for Soil–Structure Interaction Analysis in Two-Dimensional Medium. Eng. Struct. 2000, 22, 258–271. [Google Scholar] [CrossRef]
  13. Berenger, J.P. A Perfectly Matched Layer for the Absorption of Electromagnetic Waves. J. Comput. Phys. 1994, 114, 185–200. [Google Scholar] [CrossRef]
  14. Hagstrom, T. New Results on Absorbing Layers and Radiation Boundary Conditions; Springer: Berlin/Heidelberg, Germany, 2003; pp. 1–42. [Google Scholar] [CrossRef]
  15. Rabinovich, D.; Givoli, D.; Bécache, E. Comparison of High-Order Absorbing Boundary Conditions and Perfectly Matched Layers in the Frequency Domain. Int. J. Numer. Method. Biomed. Eng. 2010, 26, 1351–1369. [Google Scholar] [CrossRef]
  16. Pettigrew, J.S.; Mulholland, A.J.; Tant, K.M.M. Towards a Combined Perfectly Matching Layer and Infinite Element Formulation for Unbounded Elastic Wave Problems. Math. Mech. Solids 2021, 27, 794–812. [Google Scholar] [CrossRef]
  17. Basu, U.; Chopra, A.K. Perfectly Matched Layers for Transient Elastodynamics of Unbounded Domains. Int. J. Numer. Methods Eng. 2004, 59, 1039–1074. [Google Scholar] [CrossRef]
  18. Zhang, W.; Taciroglu, E. A Novel Rayleigh-Type Viscoelastic Perfectly-Matched-Layer for Wave Propagation Analysis: Formulation, Implementation and Application. Comput. Methods Appl. Mech. Eng. 2021, 383, 113913. [Google Scholar] [CrossRef]
  19. Hibbitt, K.S. Inc. ABAQUS/Standard Analysis User’s Manual; Hibbitt, K.S. Inc.: Kansas City, KS, USA, 2007. [Google Scholar]
  20. Josifovski, J. Analysis of Wave Propagation and Soil–Structure Interaction Using a Perfectly Matched Layer Model. Soil. Dyn. Earthq. Eng. 2016, 81, 1–13. [Google Scholar] [CrossRef]
  21. Zhang, W.; Seylabi, E.E.; Taciroglu, E. An ABAQUS Toolbox for Soil-Structure Interaction Analysis. Comput. Geotech. 2019, 114, 103143. [Google Scholar] [CrossRef]
  22. Kausel, E.; de Oliveira Barbosa, J.M. PMLs: A Direct Approach. Int. J. Numer. Methods Eng. 2012, 90, 343–352. [Google Scholar] [CrossRef]
  23. Harari, I.; Albocher, U. Studies of FE/PML for Exterior Problems of Time-Harmonic Elastic Waves. Comput. Methods Appl. Mech. Eng. 2006, 195, 3854–3879. [Google Scholar] [CrossRef]
  24. Esmaeilzadeh Seylabi, E.; Jeong, C.; Taciroglu, E. On Numerical Computation of Impedance Functions for Rigid Soil-Structure Interfaces Embedded in Heterogeneous Half-Spaces. Comput. Geotech. 2016, 72, 15–27. [Google Scholar] [CrossRef]
  25. Fontara, I.K.; Schepers, W.; Savidis, S.; Rackwitz, F. Finite Element Implementation of Efficient Absorbing Layers for Time Harmonic Elastodynamics of Unbounded Domains. Soil. Dyn. Earthq. Eng. 2018, 114, 625–638. [Google Scholar] [CrossRef]
  26. Guddati, M.N.; Lim, K.W. Continued Fraction Absorbing Boundary Conditions for Convex Polygonal Domains. Int. J. Numer. Methods Eng. 2006, 66, 949–977. [Google Scholar] [CrossRef]
  27. Guddati, M.N.; Tassoulas, J.L. Continued-Fraction Absorbing Boundary Conditions for the Wave Equation. J. Comput. Acoust. 2000, 8, 139–156. [Google Scholar] [CrossRef]
  28. Zahid, M.A.; Guddati, M.N. Padded Continued Fraction Absorbing Boundary Conditions for Dispersive Waves. Comput. Methods Appl. Mech. Eng. 2006, 195, 3797–3819. [Google Scholar] [CrossRef]
  29. Guddati, M. Perfectly Matched Discrete Layers for Modeling Unbounded Domains. J. Acoust. Soc. Am. 2020, 148, 2452. [Google Scholar] [CrossRef]
  30. Savadatti, S.; Guddati, M.N. A Finite Element Alternative to Infinite Elements. Comput. Methods Appl. Mech. Eng. 2010, 199, 2204–2223. [Google Scholar] [CrossRef]
  31. Thirunavukkarasu, S.; Guddati, M.N. Absorbing Boundary Conditions for Time Harmonic Wave Propagation in Discretized Domains. Comput. Methods Appl. Mech. Eng. 2011, 200, 2483–2497. [Google Scholar] [CrossRef]
  32. Van Nguyen, D.; Kim, D. Dynamic Soil-Structure Interaction Analysis in Time Domain Based on a Modified Version of Perfectly Matched Discrete Layers. J. Rock. Mech. Geotech. Eng. 2020, 12, 168–179. [Google Scholar] [CrossRef]
  33. Lee, J.H.; Kim, J.K.; Kim, J.H. Nonlinear Analysis of Soil–Structure Interaction Using Perfectly Matched Discrete Layers. Comput. Struct. 2014, 142, 28–44. [Google Scholar] [CrossRef]
  34. Lee, J.H. Nonlinear Soil-Structure Interaction Analysis in Poroelastic Soil Using Mid-Point Integrated Finite Elements and Perfectly Matched Discrete Layers. Soil. Dyn. Earthq. Eng. 2018, 108, 160–176. [Google Scholar] [CrossRef]
  35. Lee, J.H.; Kim, J.H.; Kim, J.K. Perfectly Matched Discrete Layers for Three-Dimensional Nonlinear Soil–Structure Interaction Analysis. Comput. Struct. 2016, 165, 34–47. [Google Scholar] [CrossRef]
  36. Nguyen, C.T.; Tassoulas, J.L. Reciprocal Absorbing Boundary Condition with Perfectly Matched Discrete Layers for Transient Analysis of SV-P Waves in a Layered Half-Space. Int. J. Solids Struct. 2018, 155, 89–108. [Google Scholar] [CrossRef]
  37. Van Nguyen, D.; Kim, D.; Park, C.; Choi, B. Seismic Soil–Structure Interaction Analysis of Concrete Gravity Dam Using Perfectly Matched Discrete Layers with Analytical Wavelengths. J. Earthq. Eng. 2021, 25, 1657–1678. [Google Scholar] [CrossRef]
  38. McKenna, F. OpenSees: A Framework for Earthquake Engineering Simulation. Comput. Sci. Eng. 2011, 13, 58–66. [Google Scholar] [CrossRef]
  39. Jeremić, B.; Kunnath, S.; Xiong, F. Influence of Soil–Foundation–Structure Interaction on Seismic Response of the I-880 Viaduct. Eng. Struct. 2004, 26, 391–402. [Google Scholar] [CrossRef]
  40. Zhang, Y.; Conte, J.P.; Yang, Z.; Elgamal, A.; Bielak, J.; Acero, G. Two-Dimensional Nonlinear Earthquake Response Analysis of a Bridge-Foundation-Ground System. Earthq. Spectra 2008, 24, 343–386. [Google Scholar] [CrossRef]
  41. Bielak, J.; Christiano, P. On the Effective Seismic Input for Non-Linear Soil-Structure Interaction Systems. Earthq. Eng. Struct. Dyn. 1984, 12, 107–119. [Google Scholar] [CrossRef]
  42. Bielak, J.; Loukakis, K.; Hisada, Y.; Yoshimura, C. Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions, Part I: Theory. Bull. Seismol. Soc. Am. 2003, 93, 817–824. [Google Scholar] [CrossRef]
  43. Yoshimura, C.; Bielak, J.; Hisada, Y.; Fernández, A. Domain Reduction Method for Three-Dimensional Earthquake Modeling in Localized Regions, Part II: Verification and Applications. Bull. Seismol. Soc. Am. 2003, 93, 825–841. [Google Scholar] [CrossRef]
  44. Weaver, W.; Johnston, P.R. Finite Elements for Structural Analysis; Prentice-Hall: Hoboken, NJ, USA, 1984. [Google Scholar]
  45. Cook, R.D. Concepts and Applications of Finite Element Analysis; John Wiley and Sons: Hoboken, NJ, USA, 2007. [Google Scholar]
  46. Chopra, A.K. Dynamics of Structures: Theory and Applications, 2nd ed.; Prentice Hall: Hoboken, NJ, USA, 2001. [Google Scholar]
  47. Poursartip, B.; Fathi, A.; Kallivokas, L.F. Seismic Wave Amplification by Topographic Features: A Parametric Study. Soil. Dyn. Earthq. Eng. 2017, 92, 503–527. [Google Scholar] [CrossRef]
  48. Zhang, W.; Seylabi, E.E.; Taciroglu, E. Effects of Soil Stratigraphy on Dynamic Soil-Structure Interaction Behavior of Large Underground Structures. In Proceedings of the International Society for Soil Mechanics and Geotechnical Engineering, Seoul, Republic of Korea, 17–21 September 2017. [Google Scholar]
Figure 1. (a) Physical problem; (b) Soil domain with special boundary condition.
Figure 1. (a) Physical problem; (b) Soil domain with special boundary condition.
Applsci 14 08519 g001
Figure 2. (a) Half-space problem; (b) Interior domain with PMDLs; (c) Four types of PMDL elements.
Figure 2. (a) Half-space problem; (b) Interior domain with PMDLs; (c) Four types of PMDL elements.
Applsci 14 08519 g002
Figure 3. PMDL element: (a) global coordinates; (b) local coordinates.
Figure 3. PMDL element: (a) global coordinates; (b) local coordinates.
Applsci 14 08519 g003
Figure 4. Finite elements and PMDL elements.
Figure 4. Finite elements and PMDL elements.
Applsci 14 08519 g004
Figure 5. Illustration of (a) the 2-D half-space problem; (b) the PMDL truncated domain.
Figure 5. Illustration of (a) the 2-D half-space problem; (b) the PMDL truncated domain.
Applsci 14 08519 g005
Figure 6. (a) Input force history and (b) Fourier amplitude.
Figure 6. (a) Input force history and (b) Fourier amplitude.
Applsci 14 08519 g006
Figure 7. Comparison of displacements for the selected points; (a) vertical displacement at P1, (b) vertical displacement at P2, (c) horizontal displacement at P3, (d) vertical displacement at P3, (e) horizontal displacement at P4, and (f) vertical displacement at P4.
Figure 7. Comparison of displacements for the selected points; (a) vertical displacement at P1, (b) vertical displacement at P2, (c) horizontal displacement at P3, (d) vertical displacement at P3, (e) horizontal displacement at P4, and (f) vertical displacement at P4.
Applsci 14 08519 g007aApplsci 14 08519 g007b
Figure 8. Contour plots of the normalized displacement magnitudes obtained by using the PMDL boundary.
Figure 8. Contour plots of the normalized displacement magnitudes obtained by using the PMDL boundary.
Applsci 14 08519 g008
Figure 9. Displacement histories of four selected points; (a) P1, (b) P2, (c) P3, and (d) P4.
Figure 9. Displacement histories of four selected points; (a) P1, (b) P2, (c) P3, and (d) P4.
Applsci 14 08519 g009aApplsci 14 08519 g009b
Figure 10. (a) Modelling of a half-space domain using ABCs and DRM; (b) equivalent 1-D SRA for estimation of free-field response for DRM; (c) schematic representation of SV wave propagation in a 2D homogeneous half-space; (d) numerical model constructed for DRM verification problem.
Figure 10. (a) Modelling of a half-space domain using ABCs and DRM; (b) equivalent 1-D SRA for estimation of free-field response for DRM; (c) schematic representation of SV wave propagation in a 2D homogeneous half-space; (d) numerical model constructed for DRM verification problem.
Applsci 14 08519 g010
Figure 11. Plots for the applied Ricker pulse include; (a) displacement, (b) Fourier amplitude, (c) velocity, and (d) acceleration.
Figure 11. Plots for the applied Ricker pulse include; (a) displacement, (b) Fourier amplitude, (c) velocity, and (d) acceleration.
Applsci 14 08519 g011
Figure 12. Displacement histories for case A; (a) horizontal displacement at P1, (b) vertical displacement at P1, (c) horizontal displacement at P2, (d) vertical displacement at P2, (e) horizontal displacement at P3, and (f) vertical displacement at P3.
Figure 12. Displacement histories for case A; (a) horizontal displacement at P1, (b) vertical displacement at P1, (c) horizontal displacement at P2, (d) vertical displacement at P2, (e) horizontal displacement at P3, and (f) vertical displacement at P3.
Applsci 14 08519 g012
Figure 13. Horizontal displacement histories for case B at (a) P1, (b) P2, and (c) P3.
Figure 13. Horizontal displacement histories for case B at (a) P1, (b) P2, and (c) P3.
Applsci 14 08519 g013
Figure 14. (a) Extended domain model, (b) PMDL-DRM combined model, (c) Winkler model, (d) acceleration time-history of the Kobe earthquake, and (e) comparison of the roof displacements of three distinct models.
Figure 14. (a) Extended domain model, (b) PMDL-DRM combined model, (c) Winkler model, (d) acceleration time-history of the Kobe earthquake, and (e) comparison of the roof displacements of three distinct models.
Applsci 14 08519 g014
Figure 15. Illustration of (a) single-cell, and (b) double-cell tunnels for studying the effects of incidence angle.
Figure 15. Illustration of (a) single-cell, and (b) double-cell tunnels for studying the effects of incidence angle.
Applsci 14 08519 g015
Figure 16. Schematic view of deformation components.
Figure 16. Schematic view of deformation components.
Applsci 14 08519 g016
Figure 17. Illustration of drift ratios: (a) single-cell tunnel, (b) its Fourier amplitudes, (c) double-cell tunnel, and (d) its Fourier amplitudes.
Figure 17. Illustration of drift ratios: (a) single-cell tunnel, (b) its Fourier amplitudes, (c) double-cell tunnel, and (d) its Fourier amplitudes.
Applsci 14 08519 g017aApplsci 14 08519 g017b
Figure 18. Illustration of rocking rotations: (a) single-cell tunnel, (b) its Fourier amplitudes, (c) double-cell tunnel, and (d) its Fourier amplitudes.
Figure 18. Illustration of rocking rotations: (a) single-cell tunnel, (b) its Fourier amplitudes, (c) double-cell tunnel, and (d) its Fourier amplitudes.
Applsci 14 08519 g018
Figure 19. Variation in peak: (a) drift ratios and (b) rocking rotations.
Figure 19. Variation in peak: (a) drift ratios and (b) rocking rotations.
Applsci 14 08519 g019
Figure 20. Variation in maximum spectral amplifications; (a) Sd, (b) PSV, (c) PSA.
Figure 20. Variation in maximum spectral amplifications; (a) Sd, (b) PSV, (c) PSA.
Applsci 14 08519 g020
Figure 21. Plot of peak deformations for various incidence angles in a single-cell tunnel; (a) θ = 0 ° , (b) θ = 10 ° , (c) θ = 20 ° , (d) θ = 30 ° .
Figure 21. Plot of peak deformations for various incidence angles in a single-cell tunnel; (a) θ = 0 ° , (b) θ = 10 ° , (c) θ = 20 ° , (d) θ = 30 ° .
Applsci 14 08519 g021
Figure 22. Plot of peak deformations for various incidence angles in a double-cell tunnel; (a) θ = 0 ° , (b) θ = 10 ° , (c) θ = 20 ° , (d) θ = 30 ° .
Figure 22. Plot of peak deformations for various incidence angles in a double-cell tunnel; (a) θ = 0 ° , (b) θ = 10 ° , (c) θ = 20 ° , (d) θ = 30 ° .
Applsci 14 08519 g022
Figure 23. The maximum bending moment diagrams for the single-cell tunnel; (a) θ = 0 ° , (b) θ = 10 ° , (c) θ = 20 ° , (d) θ = 30 ° .
Figure 23. The maximum bending moment diagrams for the single-cell tunnel; (a) θ = 0 ° , (b) θ = 10 ° , (c) θ = 20 ° , (d) θ = 30 ° .
Applsci 14 08519 g023
Figure 24. The maximum bending moment diagrams for the double-cell tunnel; (a) θ = 0 ° , (b) θ = 10 ° , (c) θ = 20 ° , (d) θ = 30 ° .
Figure 24. The maximum bending moment diagrams for the double-cell tunnel; (a) θ = 0 ° , (b) θ = 10 ° , (c) θ = 20 ° , (d) θ = 30 ° .
Applsci 14 08519 g024
Figure 25. Acceleration histories for the control points of the single-cell tunnel; (a) horizontal acceleration at P1, (b) vertical acceleration at P1, (c) horizontal acceleration at P1, and (d) vertical acceleration at P2.
Figure 25. Acceleration histories for the control points of the single-cell tunnel; (a) horizontal acceleration at P1, (b) vertical acceleration at P1, (c) horizontal acceleration at P1, and (d) vertical acceleration at P2.
Applsci 14 08519 g025
Figure 26. Acceleration histories for the control points of the double-cell tunnel; (a) horizontal acceleration at P1, (b) vertical acceleration at P1, (c) horizontal acceleration at P1, and (d) vertical acceleration at P2.
Figure 26. Acceleration histories for the control points of the double-cell tunnel; (a) horizontal acceleration at P1, (b) vertical acceleration at P1, (c) horizontal acceleration at P1, and (d) vertical acceleration at P2.
Applsci 14 08519 g026
Table 1. Abscissas and weights of the numerical integration methodologies.
Table 1. Abscissas and weights of the numerical integration methodologies.
Integration TypeAbscissaWeights
(Local r Direction)
Weights
(Local s Direction)
r1r2r3r4W1W2W1W2
2 × 2 integration−0.5770.577−0.5770.5771.01.01.01.0
1 × 2 integration0.00.0−0.5770.5772.0-1.01.0
2 × 1 integration−0.5770.5770.00.01.01.02.0-
1 × 1 integration0.00.00.00.02.0-2.0-
Table 2. Time increment values for Extended Newmark-β method.
Table 2. Time increment values for Extended Newmark-β method.
Acceleration Scheme α ¯ β ¯ γ ¯
Average1/121/41/2
Linear1/241/61/2
Table 3. Model properties used for verification.
Table 3. Model properties used for verification.
ParameterValueUnit
Vs400m/s
ρ 2000kg/m3
ν0.3-
m9-
n2-
Table 4. Material properties of the soil and tunnels.
Table 4. Material properties of the soil and tunnels.
Young’s
Modulus (kPa)
Shear Modulus (kPa)Poisson’s RatioShear Wave
Velocity (m/s)
Density (kg/m3)
Soil-80,0000.32002000
Tunnels32,000,000-0.2-2500
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Uzun, S.; Ayvaz, Y. Implementation of PMDL and DRM in OpenSees for Soil-Structure Interaction Analysis. Appl. Sci. 2024, 14, 8519. https://doi.org/10.3390/app14188519

AMA Style

Uzun S, Ayvaz Y. Implementation of PMDL and DRM in OpenSees for Soil-Structure Interaction Analysis. Applied Sciences. 2024; 14(18):8519. https://doi.org/10.3390/app14188519

Chicago/Turabian Style

Uzun, Sefa, and Yusuf Ayvaz. 2024. "Implementation of PMDL and DRM in OpenSees for Soil-Structure Interaction Analysis" Applied Sciences 14, no. 18: 8519. https://doi.org/10.3390/app14188519

APA Style

Uzun, S., & Ayvaz, Y. (2024). Implementation of PMDL and DRM in OpenSees for Soil-Structure Interaction Analysis. Applied Sciences, 14(18), 8519. https://doi.org/10.3390/app14188519

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop