Seismic Response Analysis for Ordinary State-Based Peridynamics in a Linear Isotropic Elastic Material

Earthquakes may induce the failure of soil-structures built with a geo-material. It is therefore a crucially important aspect to study how a crack initiates and propagates in these soil-structures to analyze and improve their earthquake resistance. Peridynamics, although a powerful method for such an analysis, has not been applied to seismic response analysis of soil-structures on a viscous boundary surface that can consider an input seismic wave. To apply seismic response analysis to ordinary state-based peridynamics (OSB-PD), this paper proposed a fictitious viscous boundary layer and Rayleigh damping components constructed with time derivative of force density vectors. OSB-PD was semi-verified through the comparison of calculated response acceleration via the proposed method and via the dynamic finite element method (FEM). The discussions were stressed on the effect of Poisson’s ratio, the frequency, and the Rayleigh damping components relative to the calculated response acceleration in both methods. Moreover, the seismic response acceleration at the top of the analysis model was calculated through the proposed method using an actual seismic acceleration time history, and then the result was compared with that calculated by the dynamic FEM. The results of both methods in the seismic response analysis agreed well with each other. From these results, we verified that the fictitious viscous layer can represent incoming/outgoing acceleration wave from the bottom of the analysis domain, and that Rayleigh damping can be expressed using the time derivative of force density vectors.


Introduction
and Rayleigh damping component. Finally, we compare a seismic response analysis with the proposed method and dynamic FEM, using an actual seismic acceleration time history.
2 Formulation of Seismic Response Analysis in a Linearly Isotropic Elastic Material 2.1 Definition of a Semi-Infinite Domain Figure 1 shows that relationship between a seismic wave and a one-dimensional domain, with the behavior of seismic waves particularly in focus. Here a seismic wave first enters from the bottom of the domain and travels upwards. As soon as it reaches the ground surface, the seismic wave is then reflected and travels downwards. Finally, the seismic wave reaches the bottom of the domain and exits from the domain. Note that the seismic wave separates with an incident wave (E-wave) traveling upwards and a reflected wave (F-wave) traveling downwards, which may combine together to form an E + F wave. Moreover, note that an incoming wave indicates that the seismic wave enters from the bottom of the domain, whereas an outgoing wave indicates its exit from the domain.
Consider a two-dimensional semi-infinite elastic body Ω, as shown in Fig. 2. It is necessary to calculate the semi-infinite region in which Ω has two types of fictitious boundaries: a vertical boundary and horizontal boundary, which can express an incoming/outgoing wave. In the vertical boundary, there are many types of the boundaries (e.g., viscous boundary, vertical fixed boundary). To use a viscous vertical boundary in Ω, a viscous force of a one-dimensional analysis on out-side Ω is required. It is not easy to apply the viscous force to a two-dimensional boundary because of the number of particles within δ is lack near the boundary. Therefore, in this paper, for easy implementation, the vertical boundary defines a fictitious vertical displacement boundary Γ v u . In addition, the horizontal boundary defines a fictitious horizontal viscous boundary Γ h vis , in which Γ h vis indicates an input seismic wave surface. Considering a viscous boundary Γ h vis can express an incoming/outgoing wave. In doing so [27], a viscous force is necessary to act on Γ h vis . Moreover, when the pressure-wave velocity is V p ≔ ffiffiffiffiffiffiffiffi ρ=κ p and , a viscous force assumes the Cartesian coordinates (i.e., x and y indicate the horizontal and vertical directions, respectively): where ρ, μ, and κ are respectively, the material density, shear modulus, and bulk modulus and; u x and u y represent the velocity components of Γ h vis . A viscous force is a term to indicate the theoretical traction force on the surface Γ h vis , caused by E + F waves. It is obtained assuming that the material under Γ h vis is a linear elastic material with shear velocity V s .
For simplicity of discussions, as with many seismic response analyses in a two-dimensional model, only horizontal viscous force is considered in this paper, that is, only the horizontal viscous force by u x is used and the vertical displacement u y is constrained.

Peridynamic Equation of Motion for Seismic Responses Analysis
The peridynamic equation of motion introduced by Silling [9]. In this paper, we adopted the OSB-PD by Madenci and Oterkus [11]. The discretized state-based peridynamics form is where t, N H , x i , j, u, V j , and b are a calculation time, number of material points in a horizon, coordinate of a material point i, other material points in a horizon, a displacement, an incremental volume, and a body force, respectively. In the case of OSB-PD utilizing an isotropic linear elastic material, Madenci and Oterkus [22], Madenci and Oterkus [11], Chen et al. [12], and subsequent studies demonstrated a force density vector T x i ; t ð Þ x j À x i as where Fig. 2 Definition of a semi-infinite domain Ω

Journal of Peridynamics and Nonlocal Modeling
where y i and y j are deformed position of particles x i , x j ; s, ω, h, δ, and θ i are stretch between two material points, influence function, thickness of an analysis model, a radius of a horizon, and a dilatation term of material point i, respectively; a, b, and d are material parameters. Considering a seismic wave entering a body Ω from a horizontal viscous boundary Γ h vis , the already-known displacement u g (i.e., caused by a seismic wave) and the relative displacement u r can be assumed. Therefore, the unknown displacement is given by u = u r + u g . By substituting u into Eq. (2), the equation can be divided into already-known quantities and unknown relative quantities: Because the already-known displacement u g is a rigid motion, it does not induce any force density vectors T x i ; t ð Þ x j À x i and T x j ; t À Á x i À x j . Consequently, the unknown displacement is only the relative displacement u r . Accordingly, u r will be described as unknown displacement u in the succeeding discussions.

Rayleigh Damping Components
For material damping, Rayleigh damping components were considered in this paper. Basically, Rayleigh damping has two types: mass-proportional damping and stiffness-proportional damping. In the classical continuum mechanics of an isotropic linear elastic material, these are respectively denoted by the symbols α and β.
The Rayleigh damping components can be assumed to be equivalent to the following equilibrium equation in the classical continuum mechanics: which leads to the following well-known MCK-type finite element discretized equations of motion: Where σ is stress, [ Comparatively, for Eqs. (5) and (6), a mass-proportional term and a stiffness-proportional term for OSB-PD are required to introduce Rayleigh damping components into Eq. (5). Therefore, with reference to Eq. (6), the mass-proportional term can be defined as, In addition, the stiffness-proportional term is defined by a time derivative of force density vectors on the right-hand side of Eq. (5) as follows: d dt Consequently, the following expression is obtained when Eqs. (8) and (9) are added to Eq. (5): To simplify discussions, two equations are defined as follows: Substituting Eq. (11) into Eq. (10) leads to In this paper, we used a corrected surface effect [11] and a short-range force [14] in OSB-PD. Figure 3 shows a definition of an analysis model with three layers. Here, the "Material layer" indicates the subject to the analysis, the "fictitious viscous layer," a viscous boundary layer that can express an incoming and outgoing seismic wave, and the "Fictitious yconstrained Layer," a fictitious layer to constrain displacement in the y-direction.

Fictitious Viscous Boundary Layer
A fictitious viscous boundary layer means that the same material as the fictitious viscous layer infinitely exists under the material layer. As such, using the fictitious viscous layer is possible to express an incoming/outgoing wave. If u g x i ; t ð Þ in Eq. (12) is twice the input wave (i.e., well-known 2E wave), then only the E-wave enters the material layer.
For easy implementation, the layer is treated the same as a distributed pressure acting in materials [11]. From Eq. (1), a body force density vector can be expressed as, where f vis, x indicates a viscous force density of x-direction and Δ v indicates the height of the fictitious viscous layer. When the fictitious viscous layer is of multiple layers, a non-viscous force is generated as the particles do not necessarily move uniformly. Therefore, in this paper, Δ v is treated the same as particle distance, i.e., only one-line fictitious viscous layer is used. Furthermore, in the fictitious viscous layer, a displacement in the y-direction is constrained. To simplify discussions, we assumed that the material parameters (ρ and V s ) of a fictitious viscous layer, are the same as those of the material layer.

Fictitious Y-Constrained Layer Γ v u
A fictitious y-constrained layer indicates the lateral boundaries of an analysis model that represents no vertical displacement at faraway places. Therefore, a displacement in the ydirection is constrained in the fictitious y-constrained layer as follows. Here, note that a displacement in the x-direction is free: where the subscripts x and y refer to the x and y components of these quantities, respectively.

Time Integration
Time integration of the OSB-PD equation of motion in Eq. (12) can be performed by using explicit forward difference technique as, where n indicates the step number, Δt is the time interval, and t n + 1 =t n +Δt. Substituting Eq.
4 Semi-Verification of OSB-PD We carried out simple calculations to semi-verify the proposed method. Firstly, a sine wave was applied onto a two-dimensional material layer through the fictitious viscous layer. Moreover, herein the proposed method was compared with a dynamic FEM. Specifically, the effect of Poisson's ratio, the input seismic wave number, and the Rayleigh damping components, were discussed following the calculation of response acceleration in both methods. Secondly, to semi-verify the proposed method the seismic response acceleration at the top of an analysis model, as calculated through the proposed method using an actual seismic acceleration time history, was compared with that calculated by dynamic FEM.      150 m/s was selected herein. Thus, in the material layer and the fictitious y-constrained layer, Young's modulus E was set to 112 MPa and the material density ρ was set to 1835.5 kg/m 3 . The height of the fictitious viscous layer Δ v was set to Δ and the corresponding width Δ h set to 5Δ. "0 m," "25 m," "50 m," and "100 m" in Fig. 4 show observation points where calculated response accelerations are compared in both methods.
The input seismic wave (2E) u g;x is as below: where f = 1

/T is the frequency [Hz] and T is the period [s] and t end [s] is the analysis end time.
Note that only the E-wave was allowed to enter the material layer through the bottom of the material layer through the fictitious viscous layer. Herein and in the succeeding section, the maximum frequency of the input seismic wave used was 10 Hz. Moreover, the OSB-PD parameter horizon 2δ was set to 2 × 3.015Δ as approximately one-sixth of the wavelength (V s /   f). Also, note that if the analysis model is large and the wavelength is short, enough degrees of freedom are needed to represent a single wavelength, which consequently leads to the number of particles being extremely large. To address this, shear velocity V s was set to approximately 150 m/s, and the time interval Δt[s] was ascertained with the expression, with reference to BB-PD [16]. Figure 5 displays the constructed FEM model via a 4-node isoparametric quadrilateral element. The FEM model has 2193 nodes and 2100 elements. The height of all elements was set to 2 m, which was roughly one-sixth of the wavelength. The model size and the material parameters were the same as in the OSB-PD model. Moreover, the FEM model has yconstrained condition at the lateral edge, and a viscous boundary was used to express that only the E-wave should enter the FEM model through the bottom. Viscous boundary parameters were the same as in the PD model. More practically, the FEM calculation time interval Δt[s] was set to 0.01 s.

Semi-Verification of the Fictitious Viscous Layer
Seismic wave frequency f, analysis end time t end , and Poisson's ratio was set to 0.5 Hz, 5 s, and (0, 0.1, 0.25, 0.4), respectively. A list of all parameters used in the analysis is given in Table 1. The theoretical arrival times, t = (arrival position)/V s , were T1, T2, T3, T4, T5, and T6 for the respective seismic waves arrival positions 25, 50, 100, 100 + 50 (reflection), 100 + 75 (reflection), and 100 m + 100 m (reflection), respectively. Figures 6, 7, 8, and 9 illustrate a comparison between OSB-PD and FEM. In these figures, the expressions "0 m," "25 m,"  Furthermore, from the figures, the arrival times of OSB-PD and FEM were in good agreement with T1, T2, and T3. Similarly, the arrival times of OSB-PD and FEM, which include a reflection wave, agreed well with T4, T5, and T6. Upon arrival of the reflection wave (F-wave), E-wave becomes E + F wave, and an inflection point is created on the response acceleration. Inflection points at T4, T5, and T6 were consistent in both results.
Thus, we clarified that the fictitious viscous layer in Eq. (13) exhibits the same function as the viscous boundary used in FEM.

Semi-Verification of the Fictitious Viscous Layer for Several Frequencies
Unlike in the preceding section where only a specific frequency was examined, herein several frequency sine waves were applied to the analysis model for evaluating the function of the fictitious viscous layer. Under such case, the Poisson's ratio was set to 0.4, while the input seismic wave (2E) considered and other parameters were as in the previous section. For a clearer demonstration of the analysis results, t end was set to 2.5 s. Figures 10, 11, and 12 depict the analysis results at response accelerations of 1, 5, and 10 Hz, respectively. Qualitatively, results of OSB-PD and FEM agree well, as shown in Figs. 10 and 11. As in the previous section, all the OSB-PD results could express E-wave  transforming to 2E-wave at 100 m. However, a numerical vibration of response acceleration after T1 to T6 was confirmed in Fig. 12. Figures 10c, 11c and 12c indicate the overlap of OSB-PD and FEM curves at t = 0.6-1.6 s. Figure 12 c shows a qualitative difference between OSB-PD and FEM and the numerical vibration of response acceleration subsequent to T3. More specifically, the wavelength became short at high frequency of the input seismic wave, and several degrees of freedom were essentially needed to accurately represent the wavelength deformation. Consequently, the FEM result demonstrated that the E-wave could not be expressed into a 2E-wave at 100 m. When the degree of freedom was greater than that of FEM, the OSB-PD results showed relatively less numerical vibration and could then express the E-wave into a 2E-wave at the same position. Although self-evident, a seismic response analysis should still be performed with careful attention to element height or particle distance. The above results confirmed that the fictitious viscous layer can represent wave propagation and incoming/outgoing acceleration wave at several frequencies.

Semi-Verification of Rayleigh Damping Components
Under this section, we examined the effect of the Rayleigh damping components on the layers. To recall, Rayleigh damping components were used for the material layer and the fictitious yconstrained layer. The fictitious viscous layer had no Rayleigh damping components. The Rayleigh damping components α were 0.1, 1, 5, and 100, whereas the values for β were 0.001, 0.0025, and 0.005. Given that Eq. (9) is sensitive to the time interval, the calculation results normally diverge. To resolve this, the time interval should be smaller than in Eq. (18). Thus, the time interval was set to a value of the Eq. (18) divided by 10, for the verification of β. In addition, in this section, seismic wave frequency f was set to 1 Hz.    [16]. On one hand, when the value of α was high, the analysis result was close to the result of the quasi-static analysis. As such, no difference in time was observed, and the same response acceleration was obtained anywhere, as illustrated in Fig. 16. On the other hand, when the value of β was high, the apparent hardness of the materials increased. Thus, it was shown in Figs. 17, 18, and 19 that the arrival times of the numerical result were different from the arrival times T1 to T6. In addition, due to the influence of the Rayleigh damping β, the E-wave did not transform to 2E-wave at the 100 m position, as shown in Figs. 17, 18, and 19. This scenario was true for both analysis methods. From the above results, we showed that the Rayleigh damping can be expressed by Eqs. (8) and (9).

Seismic Acceleration
Accordingly, to semi-verify our proposed method, we compared the response acceleration at the top of model, as calculated by the proposed method using an actual seismic acceleration time history, with that calculated by the dynamic FEM. Figure 20 shows the input seismic wave [28] from the 1995 Hogyo-ken Nanbu Earthquake (HNE), a direct-type (short-period) earthquake. In the event of HNE, the duration time was 30 s, the sampling rate was 100 Hz (i.e., an observation time interval of 0.01 s), and the maximum frequency f max and minimum frequency f min was 50 and 1/30 Hz, respectively.
We employed the same analysis model as in the previous section. Young's modulus E was set to 112 MPa and the Poisson's ratio to 0.4. The time interval was 6.79 × 10 −4 s. Assuming h = 5% damping for all frequencies of the seismic wave, the corresponding Rayleigh damping components were as follows: Figure 21 demonstrates that the results of OSB-PD and FEM consistently and qualitatively agreed with each other. In particular Fig. 21 a and b give the response acceleration at "0" and "100 m," respectively, whereas Table 2 provides the maximum/minimum response acceleration at these positions. Here, an error, mainly, the difference between OSB-PD and FEM divided by the value of FEM, was observed. Nevertheless, the resulting maximum error was 1.7%, which meant that the OSB-PD and FEM results quantitatively agreed well. Overall, our proposed method is applicable to seismic response analysis.

Conclusion
We formulated the OSB-PD, as well as the Rayleigh damping component, for a seismic response analysis on a viscous boundary surface. We semi-verified the method through a comparison of the dynamic response acceleration results from both our method and the dynamic FEM. The main findings of the study are as follows: 1. Using several Poisson's ratio, i.e., several shear moduli, a qualitative agreement between the results of both methods was seen in all cases. Several frequencies were also examined. Under high frequency, a slight numerical vibration of the response acceleration was verified. However, the results for FEM and OSB-PD qualitatively agreed with each other in all cases. Therefore, the proposed viscous boundary layer is suitable to seismic response analysis. 2. Rayleigh damping was investigated in OSB-PD. On one hand, the mass-proportional damping behaved like the adaptive dynamic relaxation. On the other hand, the stiffnessproportional damping was defined by the time derivative of force density vectors. A few examples were examined, yielding a qualitative agreement between the FEM and OSB-PD results in all cases. Thus, Rayleigh damping can be expressed by the proposed method. 3. Seismic response acceleration at the top of the analysis model was calculated with the proposed method using an actual seismic acceleration time history and with the dynamic FEM. Both results agreed well with each other.
As future work, we will consider the viscous boundary layer and Rayleigh damping to clarify the tensile failure and sliding failure phenomenon of soil embankment caused by an earthquake.
Funding Information This work was supported by JSPS KAKENHI, Grant Number JP 19 K04957.

Compliance with Ethical Standards
Conflict of Interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.