Abstract
Local inclusion topography has significant influence on seismic wave propagation, and the propagation characteristics of seismic waves in poroelastic soils are obviously different from those in single-phase media. Based on Biot’s theory, the scattering of plane P1 wave by inclusion in a three-dimensional poroelastic half-space is studied by using the indirect boundary element method (IBEM). The scattering field is constructed by introducing a virtual wave near the interface between inclusion and half-space and the surface of half-space, and the virtual wave density is obtained by establishing boundary integral equation based on the boundary conditions. The effects of the depth, geometric characteristics, boundary permeability, porosity, incident frequency, and incident angle of the inclusion on elastic wave scattering are systematically analyzed. The results show that due to the soil skeleton-pore water coupling effect, when the porosity is n = 0.3, the surface displacement amplitude of dry soil is larger than that of poroelastic soil. When the porosity is n = 0.36, the surface displacement amplitude of poroelastic soil is larger than that of dry soil. The surface displacement amplitude of poroelastic-drained condition is slightly larger than that of undrained condition. With the increase of inclusion depth, the scattering of elastic wave by inclusion decreases gradually. When P1 wave is incident, the surface displacement amplitude at the depth of H = 0.5 can be increased up to three times as much as that at the depth of H = 1.5. As the inclusion becomes narrower and flatter, the scattering of elastic waves by inclusion decreases gradually. When the ratio between height and length is S = 2/5, the surface displacement magnitude can reach up to 9.5.
1. Introduction
The underground inclusion is a common local topography in coastal sedimentary areas and alluvial plain areas. Inclusion has a significant influence on seismic ground motion, and is one of the important characteristics in the exploration of oil and ore bodies. Therefore, it is of great practical engineering significance to study the scattering of elastic waves by inclusion in the field of earthquake engineering, protective engineering, geophysics, etc.
The research for the scattering of elastic waves by inclusion can generally be divided into analytical and numerical methods: by using boundary integral method, Dravinski [1] studied the scattering effects of single and two elliptical inclusions on SH, P, SV, and Rayleigh waves in two-dimensional half-space. Manoogian and Lee [2] given a semi-analytical solution to the scattering problem of SH waves under the influence of arbitrary shape inclusions by using wave function expansion method and weighted residual method. Yuan [3] used the wave function expansion method to give the closed series expression of scattering of plane SH wave from the circular elastic plug region in two-dimensional elastic half-space. By using the finite element method, Liang et al. [4–6] studied the basic rules of the linear and nonlinear effects of the inclusion on the ground motion in the layered ground under the two-dimensional condition. The study showed that the presence of the inclusion had a significant amplification effect on the ground motion. Because of the scattering of incident wave by inclusion, the vertical acceleration will be generated, and there are more short period components of the vertical acceleration, and the nonlinear amplification effect of inclusion on ground motion is generally smaller than that of linear amplification, the layered site is also significantly different from its equivalent uniform site. He and Liang [7] studied the influence of the randomness of inclusion section shape on the extreme value of surface dynamic response in two dimensions; the results showed that the influence increased with the increase of velocity difference between inclusion medium and half-space medium and increased with the decrease of inclusion burial depth. Han et al. [8] studied the effect of inclusion on the seismic response of superstructure in two-dimensional deep soft soil site by using the coupling method of finite element and indirect boundary element method, the results showed that the presence of inclusion may significantly change the dynamic response of superstructure, and the dynamic response of superstructure may be significantly enlarged.
The above studies are all aimed at the analysis of two-dimensional single-phase medium, but in reality, the soil is mostly two-phase medium or three-phase medium, especially in coastal areas, the site soil is mostly poroelastic. According to Biot’s theory [9, 10], the propagation of elastic waves in poroelastic media is affected by the coupling of soil skeleton and pore water, which is closely related to porosity, soil skeleton stiffness, Poisson’s ratio, and the interface permeability. Gurevich et al. [11] adopted Born’s theory to solve the scattering of plane waves by elliptic plugs in full-space poroelastic media. Yuan et al. [12, 13] studied the influence of poroelastic sand inclusion liquefaction on seismic response of buildings and the ground motion by using two-dimensional effective stress finite element analysis method, and obtained some useful conclusions. By using the finite-indirect boundary element coupling method, Han et al. [14] analyzed the effect of inclusion on the seismic response of superstructure in deep poroelastic sites. Liu et al. [15, 16] studied the scattering of seismic waves by 2-D local sites in a fluid saturated half-space. Recently, Kanaun et al., Panji et al., Wu and Ou also studied the wave scattering around an inclusion considering different types of waves [17–19].
It can be found that previous studies on the effects of poroelastic inclusion on elastic waves have been carried out in full or two-dimensional half-space. However, the effect of 3-D inclusion on elastic waves in three-dimensional poroelastic half-space is rarely studied, especially for the incidence of P waves. In this paper, the indirect boundary element method (IBEM) [20–22] is used to study the scattering of inclusion on plane P1 wave in three-dimensional poroelastic half-space, and the effects of poroelastic inclusion burial depth, geometry and boundary permeability conditions, incident frequency, and incident angle are extensively examined. Several important conclusions are obtained which are beneficial for the antiseismic or antiblast design of engineering near the inclusion site.
2. Model and Solution
2.1. Model
As shown in Figure 1, there is an arbitrary inclusion l in a poroelastic half-space E. It is assumed that both the inclusion and half-space are poroelastic, elastic, homogeneous, and isotropic. For the convenience of calculation and analysis, the coordinate system is established with the center of the circular Earth surface as the origin. Let the incident wave be parallel to the x-z plane and incident angle is . The surface boundary is set as , the interface between inclusion and half-space is set as . Assume that the shear modulus of poroelastic half-space is expressed as , density is expressed as , and Poisson’s ratio is expressed as . Assume that the shear modulus of inclusion is expressed as , density is expressed as , Poisson’s ratio is expressed as , and the remaining calculation parameters are given as follows.
(a)
(b)
As shown in Figure 1(b), the geometric shape of the inclusion can be expressed aswhere , , , respectively, represent the maximum horizontal dimension of the inclusion along the x-axis, y-axis, and z-axis, is the depth between the central position of the inclusion and the ground surface, is the radius of the circular surface. For the convenience of analysis, the length of the inclusion along the x-axis is assumed to remain constant, in order to consider the effect of inclusion geometry on wave scattering. The ratio between width to length () and height to length ratio () are analyzed below.
2.2. Biot’s Poroelastic Two-Phase Medium Theory
It is assumed that the displacement of solid soil skeleton in poroelastic medium is , the displacement of the fluid relative to the soil skeleton is (), and the constitutive relation of homogeneous poroelastic media can be expressed aswhere is pore pressure, is the total stress component of poroelastic soil, and e are strain components and volume strain of soil skeleton, respectively, and are the Lame constants of the soil skeleton, is dirac function, and M are parameters that characterize the compressibility of soil particles and pore fluids (, ). The soil motion equation related to and can be expressed as [10]where is pore water pressure, is total soil density, and are the mass density of soil particles and fluid, is porosity, b is the coefficient reflecting the viscous coupling and b = 0 if the internal friction is ignored, m is a parameter of like mass, which is determined by fluid density, porosity, and pore geometry, , and is the solid and liquid coupling mass density.
The wave propagating in poroelastic medium includes two kinds of expansion waves ( wave and wave) and shear wave. It is assumed that the wave numbers of wave, wave, and SV wave in the half-space are , , and , the wave numbers of wave, wave, and SV wave in the inclusion are , , and .
Dimensionless parameters λ, M, ρ, m, b, and dimensionless frequencies δ are introduced and defined as follows:
The displacement, pore water pressure, and stress field of full-space wave, wave, and SV line wave source acting at the origin are shown in the literature [20].
2.3. Wave Field Analysis
The total wave field is decomposed into free field and scattering field. The free field represents the solution when the inclusion is not in the poroelastic half-space, while the scattering field is constructed by applying virtual load on the half-space surface and the cavity surface according to IBEM principle.
Firstly, the free field analysis is carried out. The expression of potential function under incident plane wave is as follows:
The time factor has been omitted. The stress caused by free wave field is , soil skeleton displacement is , the relative displacement of the fluid is , and the pore water pressure is .
Due to the presence of inclusion, scattering waves will be generated in local topography. The scattering wave in this method is obtained by superimposing three concentrated forces (Fx, Fy, Fz) and flow source (Q) on the boundary, the scattering wave fields of half-space E and inclusion L domain can be constructed by applying virtual wave source load on the boundary. The boundary of the half-space E domain is constructed by S1 and S2, and the boundary of the inclusion L domain is constructed by S2. The soil skeleton displacement , traction force , relative fluid displacement and pore water pressure caused by the single-layer potential theory can be expressed as follows:where x is the boundary point, ξ is the wave source point, and are the wave source density under the action of concentrated force and flow source, respectively, , , , and are green’s function of soil skeleton displacement, total stress, fluid relative displacement, and pore water pressure under the action of concentrated force. , , , and are green’s function of soil skeleton displacement, total stress, fluid relative displacement, and pore water pressure under the action of flow source, respectively.
2.4. Boundary Conditions and Solution
When the poroelastic half-space free surface and the inclusion interface are both drained conditions, the boundary conditions of the problem are stress of free surface is zero, pore water pressure is zero. Solid phase displacement, stress, fluid relative displacement, and pore water pressure are continuous at the interface between half-space and inclusion. The boundary conditions can be expressed as follows:(1)Zero-stress conditions at half-space surface S1 and inclusion internal surface S2(2)Continuity condition at the interface between half-space and sediment S2
when the poroelastic half-space free surface and the inclusion interface are both undrained conditions, the boundary conditions of the problem are zero stress on the free surface and zero relative displacement of the fluid. The solid phase displacement and stress are continuous at the interface between half-space and sediment and fluid relative displacement is zero. The equations of boundary conditions can be expressed as follows:(1)Zero-stress conditions at half-space surface S1 and inclusion internal surface S2(2)Continuity condition at the interface between half-space and sediment S2
In summary, according to the single-layer potential theory, IBEM method is used to discretize the whole space, and virtual loads are applied on the boundary discrete points. From the above-mentioned drained and undrained boundary conditions, matrix equations can be constructed, respectively, virtual wave source density can be obtained, scattering wave field can be constructed, and total wave field can be obtained by superposition with free field.
3. Verification of Accuracy
Since the scattering of elastic waves by a three-dimensional poroelastic half-space inclusion has not yet been completely solved analytically, the calculation accuracy can only be investigated by comparing with the existing results in the case of degradation [23, 24]. Firstly, the inclusion is degenerated into a sphere, and then the poroelastic medium is degenerated into a single-phase medium, the calculated parameters of P1 wave are set as porosity n = 0.1, half-space Poisson’s ratio υ = 0.3, material damping ratio ζ = 0.01, inclusion depth d/a = 2, and dimensionless vertical incidence frequency = 1.0. The calculated parameters of Rayleigh waves are set as porosity n = 0.1, , , and half-space Poisson’s ratio υ = 0.3. As can be seen from Figure 2, the consistency of degradation results is good, which verifies the correctness of this method. In Figure 2, ux and uz represent the surface displacement amplitude in horizontal and vertical directions, respectively.
(a)
(b)
4. Example Analyses
For the convenience of analysis, the parameters of poroelastic medium in half-space are selected as follows: porosity n = 0.3 and 0.36, Poisson’s ratio , material damping ratio ζ = 0.001, critical bulk modulus of soil Kcr = 200 MPa, bulk modulus of soil skeleton = 36000 MPa, mass density of the solid frame , fluid bulk modulus = 2000 MPa, and fluid mass density . Poroelastic inclusion parameters are set as porosity n = 0.36, Poisson’s ratio , material damping ratio ζ = 0.001, critical bulk modulus of soil Kcr = 22 MPa, bulk modulus of soil skeleton = 4000 MPa, mass density of the solid frame , fluid bulk modulus = 2000 MPa, and mass density of fluid . The corresponding dimensionless parameters of poroelastic soils are shown in Table 1. In order to discuss the influence of inclusion geometry, the ratio of width to length is donated as D = ay/ax and the ratio of height to length is donated as .
4.1. Displacement Amplitudes Comparison of Dry Soil, Poroelastic-Drained Soil, and Poroelastic-Undrained Soil
Figure 3 shows the comparison of the surface displacement amplitudes of dry soil, poroelastic-drained soil, and poroelastic-undrained soil near the three-dimensional inclusion body by incident P wave. The buried depth of the inclusion is H = 0.5, the ratio between width and length is D = 1/1, the ratio between height and length is S = 1/4, the incident frequency is η = 0.5, 1.0 and 2.0, the incident angle is , the porosity is n = 0.3 and 0.36. In the figure, “dry” means dry soil, “drained” means boundary permeable, and “undrained” means boundary impervious. The abscissa x/a represents the ratio between the distance from each point to the origin on the surface of the x-z section (marked section A in the figure) and Y-Z (marked section B in the figure) and the long axis of the inclusion. The surface displacement amplitude corresponding to the longitudinal coordinates in half-space has been standardized. For the convenience of analysis, only the principal direction displacement is given in the figure.
From Figure 3, it can be seen that under P1 wave incidence, due to the influence of soil skeleton-pore water coupling and boundary permeability conditions of poroelastic soil, the spatial distribution characteristics of surface displacement amplitude in dry soil, saturated-drained soil, and poroelastic-undrained soil are obviously different. When the porosity is n = 0.3, the displacement amplitude of dry soil is generally greater than that of poroelastic soil, and the displacement amplitude of poroelastic-drained soil is generally greater than that of poroelastic-undrained soil. When the porosity n = 0.36, the surface displacement amplitude of poroelastic soil is greater than that of dry soil. For example, when the incident frequency is η = 0.5 and the incident angle is θ = 0°, the displacement amplitudes of principal direction are 3.4, 2.5, and 2.2. When the incident frequency is η = 0.5, the incident angle is θ = 0°, the porosity is n = 0.36, the amplitudes of principal direction are 3.4, 5.1, and 4.9. With the increase of incident angle, the principal direction displacement amplitude of P1 wave decreases gradually in both dry and poroelastic soils. For example, when incident frequency is η = 0.5, the porosity is n = 0.3, and the surface displacement amplitude is 5.1, 4.2, and 2.4 at incident angles of 0°, 30°, and 60 °in poroelastic-drained situation. With the increase of incident frequency, the spatial distribution of P1 wave displacement becomes more complex, when incident frequency is η = 2.0, the surface displacement amplitudes and spatial distribution of dry soil, poroelastic-drained soil, and poroelastic-undrained soil are close to each other.
4.2. Influence of Inclusion Depth on Surface Displacement Amplitude
Figure 4 shows the surface displacement amplitudes comparisons of different buried depths at drained and undrained boundary conditions. The inclusion buried depths is donated as = 0.5, 1.0, and 1.5, the incident frequency is donated as η = 0.5, 1.0, and 1.5, and the incident angle of P1 wave is θ = 0° (i.e., vertical incidence). Because the surface displacement amplitude of cross section x-z and cross section y-z are exactly the same (as can be seen from the previous section), only the surface displacement amplitude of cross section x-z is given here, it can be seen that the surface displacement amplitude generally appears when the burial depth is H = 0.5, which is similar to that of dry soil. Under low frequency incidence, the displacement amplitude decreases gradually with the increase of the burial depth of inclusion. For example, under the drained boundary condition, when the incident frequency is η = 1.0, the porosity is n = 0.36, and the buried depth is H = 0.5, 1.0, and 1.5, the principal direction displacement amplitude are 5.9, 2.6, and 2.0. The displacement amplitude at H = 0.5 is three times larger than that at H = 1.5. Under high frequency incidence, sometimes the displacement amplitudes at the depth H = 1.5 are larger than those at the depth H = 1.0. For example, under undrained boundary condition, when the incident frequency is η = 2.0, the porosity is n = 0.3, the displacement amplitudes of the principal direction corresponding to the buried depth H = 1.0 and 1.5 are 2.4 and 3.0. Therefore, attention should be paid to the high frequency incidence in practical engineering. The change of the depth of inclusion sometimes changes the spatial distribution of P1 wave displacement. When the buried depth is high (H = 1.5) and the incident frequency is low (η = 0.5 and 1.0), the scattering effect of the inclusion on the plane P1 wave is very small and the surface displacement amplitude is close to the free field.
(a)
(b)
4.3. Influence of the Geometry of Inclusion on the Amplitude of Surface Displacement
Figure 5 shows the surface displacement amplitude near the inclusion when P1 wave incidence; in the figure, “drained” represents the boundary permeability and “undrained” represents the boundary impermeability. In order to systematically study the influence of the geometry of inclusion on elastic wave scattering, two cases are discussed. The first is to keep the ratio between height and length S = 1/4 unchanged, the ratio between width and length of inclusion is changed as D = 1/4, 1/2, and 1/4, that is, the inclusion becomes narrower gradually. In the second case, keep the ratio between width and length D = 1/1 unchanged, the ratio between height and length is changed as S = 2/5, 1/4, and 1/10, that is, the inclusion becomes flatter gradually. Incident frequency is donated as η = 0.5, 1.0, and 1.5, and incident angle is θ = 0° (i.e., vertical incidence). Because the surface displacement amplitude of cross section x-z and cross section y-z are exactly the same (as can be seen from the previous section), only the surface displacement amplitude of cross section x-z is given here.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
It can be seen from the figure that the change of inclusion geometry has a great influence on the surface displacement amplitude and displacement spatial distribution. Generally, the scattering effect of the P1 wave decreases with the decrease of the ratio between width and length D (gradually narrowing) while keeping the ratio between height and length S unchanged. When the incident frequency is η = 1.0 and the boundary condition is undrained, the surface maximum displacement amplitudes of principal direction are is 3.3, 3.1, and 2.6, when the ratio between width and length of the inclusion is D = 1/1, 1/2, and 1/4, respectively, comparing with the free field, the maximum displacement amplitudes will increase by 70%, 55%, and 30%, respectively. With the increase of incident frequency, the displacement amplitudes spatial distribution will be different between drained boundary and undrained boundary. For example, when P1 wave incidence is vertical and the incident frequency is η = 1.0, the surface displacement amplitude forms two peaks near the edge of the inclusion under the boundary drained condition, while the surface displacement amplitude forms a peak at the center of the inclusion body under the undrained condition. When the ratio between width and length of the inclusion is D = 1/4, the effect of the inclusion on the scattering of P1 wave at low frequency (η = 0.5) is very small. When the ratio between height and length of inclusion is changed as S = 2/5, S = 1/4, and S = 1/10 with the ratio between width and length D = 1/1 keeping constant (the inclusion flattens gradually), the scattering effect of the inclusion on elastic waves gradually decreases. For example, when the incident frequency is η =1.0 at drained condition, the ratio between height and length of inclusion is S = 2/5, S = 1/4, and S = 1/10, the surface displacement amplitudes are 5.2, 3.9, and 2.9, respectively, which increase by 195%, 95%, and 45% compared with the free field. It is noteworthy that the response of P1 wave to the change of the ratio between height and length is particularly remarkable. When P1 wave incidence is vertical, the incident frequency is η = 1.0 at the undrained condition, the ratio between height and length is S = 2/5, and the displacement amplitude can reach 9.5, which is 2.9 times and 3.4 times of the ratio between height and length of S = 1/4 and 1/10. When the inclusion is flat (the ratio between height and length is S = 1/10) and the incident frequency is low, the scattering of P1 wave by the inclusion is very small.
4.4. Spectrum Analysis of Displacement Amplitude
Figure 6 shows the surface displacement amplitude spectrum near the inclusion when P1 wave incidence is vertical at drained boundary, calculation parameters are denoted as inclusion depth H = 0.5, incident frequency η is 0∼3. The discussion is divided into two cases. First, keep the ratio between height and length S = 1/4 unchanged, and the ratio between width and length is changed as D = 1/4, 1/2, and 1/1. Second, keep the ratio between width and length D = 1/1 unchanged and the height to length ratio is changed as S = 1/10, 1/4, and 2/5. Several points are selected along the x-axis (Section A) and y-axis (Section B) to represent the principal direction displacement spectrum near the inclusion, five observation points are taken along the x-axis, including x/a = 0, 0.5, 1.0, 1.5, and 2.0, five observation points are taken along the y-axis, including y/a = 0, 0.5, 1.0, 1.5, and 2.0, respectively. As can be seen from Figure 6, the displacement amplitude increases and the displacement spectrum curve oscillates more intensely with the increase of inclusion thickness and width in the case of poroelastic soil, and the effect of scattering by inclusion thickness is greater than that by inclusion width. The oscillation degree of the spectrum curve in the range of 1ax is obviously greater than that the range of 1ax under P1 wave incidence. For example, when the ratio between width and length is D = 1/1, and the ratio between height and length is S = 1/10, 1/4, and 2/5, the displacement spectrum amplitudes are 3.8, 4.4, and 8, they all appear at the coordinate x = 0. When the inclusion is narrow and thin, the scattering effect of the inclusion by the P1 wave is very weak. For example, when S = 1/4 and D = 1/4, the displacement amplitude of each observation point is close to the free field, and the oscillation degree of the spectrum curve is also weak. When the frequency is low, the scattering effect of the inclusion on the P1 wave is not obvious. For example, when the incident frequency is η < 0.5, the displacement amplitude of each observation point is close to the free field. When the incident frequency is 1 < η < 2, the scattering effect of the inclusion body on the P1 wave is the most obvious, and the maximum displacement often appears here.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
5. Conclusions
Based on Biot’s theory, a high-accuracy indirect boundary element method is adopted to study the scattering of plane P1 wave by inclusion in a poroelastic half-space, the influence of incident frequency, incident angle, boundary permeability condition, porosity, inclusion depth, and geometric characteristics on elastic wave scattering is analyzed, some useful conclusions are obtained:(1)Due to the coupling of soil skeleton-pore water, the influence of inclusion on elastic wave scattering is significantly different between the poroelastic soil and the dry soil. When the porosity is n = 0.3, the surface displacement amplitude of the dry soil is larger than the poroelastic-drained and the poroelastic-undrained conditions. When the porosity is n = 0.36, the surface displacement amplitude of poroelastic soil is larger than that in the dry soil. The influence of the incident angle is also obvious; the difference between the displacement amplitude of the dry soil and the poroelastic soil decreases with the increase of the incident angle, and the spatial distribution will gradually become similar.(2)The depth of inclusion in poroelastic soil has a significant effect on the surface displacement amplitude and the spatial distribution characteristics, the scattering effect of inclusion on P1 wave increases gradually with the shallowing of the depth. The degree of increase can be up to three times under low frequency incident. For high frequency incidence, the fluctuation of surface displacement amplitude becomes more complex. When the depth is deep (H = 1.5) and the incident frequency is low (η = 0.5 and 1.0), the scattering effect of the inclusion on the plane P1 wave is very small, and the surface displacement amplitude is close to the free field.(3)When the inclusion is very narrow and thin, the scattering effect of the inclusion on P1 wave is very weak. For example, the spectrum curve of displacement amplitude at each observation point is close to the free field when S = 1/4 and D = 1/4, and the spectrum curve oscillation degree is also very weak. The surface displacement amplitude will produce a certain focusing effect under the influence of inclusion, the vibration degree of the spectrum curve at the points within the range of 1ax is significantly greater than that of points beyond the range of 1ax.
It should be noted that present method is more suitable for solving linear wave motion. While for nonlinear problems, the advantage of BEM is not so significant compared with the FEM or FDM.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was financially supported by the Innovation Foundation of Nanjing Institute of Technology (CKJB201909), Natural Science Foundation of Jiangsu Province (BK20150726), National Natural Science Foundation of China (51808283), and Natural Science Fund for Colleges and Universities in Jiangsu Province (18KJB580007).