Wave Propagation and Scattering around a Radially Inhomogeneous Cylindrical Inclusion in a Full Space

: The geological structure, such as inclusions, may strongly affect the wave propagation and underground motions during earthquakes. Previous studies mainly focus on geological inclusion with the homogeneous medium. In this paper, the propagation and scattering of incident plane SH waves in and around an inhomogeneous cylindrical inclusion with a radially-varying modulus is studied. In terms of a radial wave function expansion, a rigorous analytical approach is formulated for general computation for the elastodynamic problem. A comprehensive set of numerical examples are presented to illustrate the sensitivity of the underground motion to the rigidity proﬁle of the geological inclusion.


Introduction
The influence of local irregular geological structure on ground and underground motions is an important research topic in seismology and earthquake engineering [1,2].The underground inclusion, as one kind of local irregular geological structures, has different properties, such as rigidity from the surrounding medium.Therefore, it may affect the earthquake wave propagation and cause nearby engineering structures, such as tunnels to suffer stronger seismic action [3][4][5][6].However, there are still no relevant provisions to follow in the engineering seismic code.Thus, the propagation and scattering of seismic waves in sites containing inclusions will affect the sustainability of the safety of engineering structures.
To reveal the propagation mechanism of seismic waves around the inclusion, a lot of numerical and analytical studies have been carried out.The numerical research based on the boundary type method is the most fruitful due to its advantage to accurately satisfy the boundary radiation condition of infinite space.For example, Dravinski [7] used the boundary integral method to solve the problem of elastic inclusion in a half space.Hadley [8] used the dynamic full plane boundary element method to study the scattering of seismic waves by inclusions in a non-uniform elastic half space.Heymsfield [9] solved the two-dimensional scattering of SH waves by rock inclusions in the soil layer by using the direct boundary element method.Dravinski and Sheikhhasani [10] used the same method to study the scattering of plane harmonic SH waves by arbitrarily shaped multilayer inclusions embedded in a half space.Dong et al. [11] also solved the inclusion problem in an elastic half space using different integral equation methods.Rus and Gallego [12] proposed a direct differentiation-based boundary integral equation method for sensitivity analysis of inclusions in elastodynamic media.Parvanova et al. [13] used a direct displacement-based boundary element method to study the dynamic response of a medium with multiple inclusions under anti plane strain conditions.Liu et al. [14] studied the scattering of plane SH waves by a group of inclusions in an elastic half space based on the indirect boundary element method.Panji et al. [15] proposed a direct time-domain numerical method based on the half space Green's function to study the scattering of SH waves by arbitrary shaped inclusions.Mojtabazadeh-Hasanlouei et al. [16] analyzed the underground multiple inclusions model under the propagation of SH wave based on the time-domain boundary element method established by the half space Green's function.Qi et al. [17] calculated the displacement field and stress field of a half space with elliptical inclusion under a line source load using the Green's function method, and analyzed the factors affecting the dynamic stress concentration factor (DSCF) around the elliptical inclusion.
Numerical methods can deal with more real and complex cases, but their accuracy often needs to be verified by analytical solutions.Pao and Mow [18] first established the model of scattering of SH waves by a cylindrical inclusion in a full space, and derived the wave function series solution to it.Later, Smerzini et al. [19] used the wave function expansion method to analyze the scattering of SH waves by a cylindrical inclusion in a half space.Jiang et al. [20] studied the dynamic response of a shallow circular inclusion under incident SH waves in a radially inhomogeneous half space by applying complex function theory and multipolar coordinate system.Recently, Tokmechi et al. [21] proposed the wave function series solution for SH wave scattering by a semicircular inclusion in a half space.
It is worth noting that the above studies all assume that the inclusions are of a homogeneous and an elastic medium.However, the actual inclusions are often inhomogeneous.Previous studies have shown that the propagation law of waves in the inhomogeneous medium is different from that in the homogeneous medium [22][23][24].Accordingly, the wave field in a full space containing inhomogeneous inclusions may be different as well.To reveal the mechanism of an inhomogeneous inclusion modifying the wave field, the wave function series solution to SH wave scattering by a cylindrical inclusion whose shear modulus varies continuously in the form of power function in the radial direction is derived in this paper.The contribution of this paper is twofold.One is to fill the gap of analytical solution to SH wave scattering by an inhomogeneous inclusion; the other is to reveal the influence of the inhomogeneous rigidity distribution characteristics of the inclusion on the seismic field around it through frequency-domain and time-domain parameter analysis.

Model and Excitation
The two-dimensional model studied in this paper is shown in Figure 1.It represents a cylindrical inclusion with radius b in a full space.The inclusion is composed of two parts defined as region 1 for r < a and region 2 for a ≤ r ≤ b, respectively.It is assumed that the materials of the core region 1 of the inclusion and the full space denoted as region 3  where the inclusion is located are both elastic, homogeneous and isotropic with constant density ρ j , shear modulus µ j and shear wave velocity c j (j = 1 or 3, indicating region 1 or region 3 ).In contrast, the material in the region 2 of the inclusion is elastic, isotropic but inhomogeneous with a continuously varying rigidity in the radial direction.To facilitate the analytical solution, the shear modulus µ 2 and shear wave velocity c 2 are assumed as power functions of radial distance r and the mass density ρ 2 remains constant.The continuously varying rigidity profile of the inhomogeneous region of the inclusion can be expressed as follows: where β is the inhomogeneity exponent, and µ 0 and c 0 are the shear modulus and shear wave velocity at the inner surface of region 2 , respectively.It should be noted that β can be any positive or negative number depending on the actual inhomogeneity profile and there is no rigorous limit of its value.
where β is the inhomogeneity exponent, and μ0 and c0 are the shear modulus and shear wave velocity at the inner surface of region ②, respectively.It should be noted that β can be any positive or negative number depending on the actual inhomogeneity profile and there is no rigorous limit of its value.To focus on the effect of continuous change of wave velocity across the interfaces among the three regions, it is assumed that c2(a) = c1 and c2(b) = c3.The variation of the wave velocity along the x-axis in the cross section at y = 0 of the model is shown in Figure 2. One can see that the wave velocity profile of the inhomogeneous inclusion is dependent on the key parameter of c1/c3.The excitation of the model is a plane SH wave of unit amplitude with harmonic vibration in the z direction.To facilitate the solution, a Cartesian coordinate system (x, y) and a polar coordinate system (r, θ) are established with the origin O at the center of the inclusion.The definition of x-axis, y-axis, and the angle θ can be seen in Figure 1.To focus on the effect of continuous change of wave velocity across the interfaces among the three regions, it is assumed that c 2 (a) = c 1 and c 2 (b) = c 3 .The variation of the wave velocity along the x-axis in the cross section at y = 0 of the model is shown in Figure 2. One can see that the wave velocity profile of the inhomogeneous inclusion is dependent on the key parameter of c 1 /c 3 .
where β is the inhomogeneity exponent, and μ0 and c0 are the shear modulus and shear wave velocity at the inner surface of region ②, respectively.It should be noted that β can be any positive or negative number depending on the actual inhomogeneity profile and there is no rigorous limit of its value.To focus on the effect of continuous change of wave velocity across the interfaces among the three regions, it is assumed that c2(a) = c1 and c2(b) = c3.The variation of the wave velocity along the x-axis in the cross section at y = 0 of the model is shown in Figure 2. One can see that the wave velocity profile of the inhomogeneous inclusion is dependent on the key parameter of c1/c3.The excitation of the model is a plane SH wave of unit amplitude with harmonic vibration in the z direction.To facilitate the solution, a Cartesian coordinate system (x, y) and a polar coordinate system (r, θ) are established with the origin O at the center of the inclusion.The definition of x-axis, y-axis, and the angle θ can be seen in Figure 1.The excitation of the model is a plane SH wave of unit amplitude with harmonic vibration in the z direction.To facilitate the solution, a Cartesian coordinate system (x, y) and a polar coordinate system (r, θ) are established with the origin O at the center of the inclusion.The definition of x-axis, y-axis, and the angle θ can be seen in Figure 1.

Construction of Wavefields
The wave field w 1 in region 1 should satisfy the Helmholtz equation in the polar coordinate system (r, θ): where k 1 = ω 1 /c 1 is the shear wave number.
By solving Equation (3) with the method of separating variables, the standing wave field in region 1 can be constructed as follows: where A n is the wave field coefficient to be determined, and J n (•) is the first kind Bessel function of order n.
According to Zhang et al. [25], the wave fields in regions 2 and 3 can be directly expressed as: where B n , C n and D n are coefficients to be determined.The expressions of U n (r) and V n (r) are consistent with Equations ( 18) and ( 19) in Zhang et al. [25], respectively.H n (1) (•) is the first kind Hankel function of order n.

Boundary Conditions
The displacement and stress continuity conditions on the two circular interfaces at r = a and b among the three regions need to be satisfied, i.e.,

Determination of Unknown Coefficients of Wave Fields
Substituting Equations ( 4) and ( 5) into ( 7), one has: Using the orthogonality of complex exponential function, i.e., 2π 0 e inθ e iqθ dθ = 2π, n + q = 0 0, n + q = 0 (12) Equation ( 11) can be simplified as: Substituting Equations ( 4) and ( 5) into (8) and also taking advantage of Equation ( 12), the following relational equation can be obtained as: Similarly, substituting Equations ( 4) and ( 5) into ( 9) and ( 10), one has: Combining Equations ( 13) and ( 14) to eliminate A n , the following relational equation can be obtained as: where Similarly, combining Equations ( 15) and ( 16) to eliminate D n , and taking advantage of the Wronskian formula of Bessel function, the following relational equation can be obtained as: where Combining Equations ( 17) and ( 20), the coefficients B n and C n can be derived as: Taking Equations ( 23) and ( 24) into ( 13) and (15), respectively, the coefficients A n and D n can be obtained as: Thus, all coefficients of wave fields have been obtained in an explicit form through rigorous mathematical derivation.After a truncation of the infinite term series in Equations ( 4)-( 6) into finite terms, the wave field in each region can be calculated.In order to ensure the accuracy of the calculation results, convergence tests have been conducted and a frequency dependent truncation of the series term at n ≧ 3 + 4η is found to be enough for high accuracy of the results.The convergence test examples are shown in Figure 3 for four representative positions (O, P 1 , P 2 and P 3 in Figure 1).For the calculation cases involved in this paper, a high enough calculation accuracy can be achieved by truncating the series term to n = 15 for brevity.
Equations ( 4)-( 6) into finite terms, the wave field in each region can be calculated.In order to ensure the accuracy of the calculation results, convergence tests have been conducted and a frequency dependent truncation of the series term at n ≧ 3 + 4η is found to be enough for high accuracy of the results.The convergence test examples are shown in Figure 3 for four representative positions (O, P1, P2 and P3 in Figure 1).For the calculation cases involved in this paper, a high enough calculation accuracy can be achieved by truncating the series term to n = 15 for brevity.

Comparison with an Exact Solution by Pao and Mow
If the inhomogeneity component β is set to be 0 and the material parameters in regions ① and ② are set to be identical but different from that in region ③ (ρ1 = ρ2 ≠ ρ3， c1 = c2 ≠ c3), the model in this paper can theoretically be degenerated into the homogeneous inclusion model in Pao and Mow [18].Here is the process of the theoretical degeneration.

Comparison with an Exact Solution by Pao and Mow
If the inhomogeneity component β is set to be 0 and the material parameters in regions 1 and 2 are set to be identical but different from that in region 3 (ρ 1 = ρ 2 = ρ 3 , c 1 = c 2 = c 3 ), the model in this paper can theoretically be degenerated into the homogeneous inclusion model in Pao and Mow [18].Here is the process of the theoretical degeneration. When one can find: According to Zhang et al. [25], in Equation ( 5), one has: then Equations ( 18) and ( 21) can be reduced as: Substituting Equations ( 30) and (31) into Equations ( 23)-( 26), the coefficients A n , B n , C n and D n can be degenerated into: C n = 0 (33) Since the series terms in the Equations ( 16) and ( 20) in Pao and Mow [18] are n = 0, 1, 2, . . ., +∞ and those of the model in this paper are n = −∞, . . ., +∞, one can substitute Equation (32) into Equations ( 4) and ( 5) and unify them in the same series form as follows for the convenience of comparison.
Similarly, taking Equation (32) into Equation ( 6), one has: Note that Equations ( 35) and ( 36) are completely consistent with Equations ( 16) and ( 20) in Pao and Mow [18], which verifies the correctness of the theory in this paper.

Comparison with the Exact Solution for Free Field
Setting the parameters of the radially inhomogeneous cylindrical inclusion model in this paper as β = 0.01, ρ 1 /ρ 3 = ρ 2 /ρ 3 = 1, c 1 /c 3 = 1, the model can be approximately degenerated into a homogeneous free field model in a full space.The displacement amplitudes |u| at the inner surface (r = a) and the outer surface (r = b) of the inclusion are calculated to be approximately 1, as shown in Figure 4.This agreement with that predicted by the exact solution for a free field in a full space verifies the accuracy of the method in this paper as well.
  Since the series terms in the Equations (1.16) and (1.20) in Pao and Mow [18] are n = 0,1,2, …, +∞ and those of the model in this paper are n = −∞, …, +∞, one can substitute Equation (32) into Equations ( 4) and ( 5) and unify them in the same series form as follows for the convenience of comparison.

Comparison with the Exact Solution for Free Field
Setting the parameters of the radially inhomogeneous cylindrical inclusion model in this paper as β = 0.01, ρ1/ρ3 = ρ2/ρ3 = 1, c1/c3 = 1, the model can be approximately degenerated into a homogeneous free field model in a full space.The displacement amplitudes |u| at the inner surface (r = a) and the outer surface (r = b) of the inclusion are calculated to be approximately 1, as shown in Figure 4.This agreement with that predicted by the exact solution for a free field in a full space verifies the accuracy of the method in this paper as well.For the soft inclusion, the wave velocity in the region 1 is set as c 1 /c 3 = c 2 (a)/c 3 = 0.5, while c 1 /c 3 = c 2 (a)/c 3 =2 is taken for the hard inclusion.

Effect of the Rigidity Profile of the Inclusion on the Frequency Domain Displacement and DSCF
To reveal the influence of the inhomogeneous rigidity or wave velocity profile on the seismic response of the inclusion, the surface displacement amplitudes as well as the dynamic stress concentration factors of both soft and hard inclusions at different dimensionless frequencies are calculated, as shown in Figures 5-7.The dimensionless parameters in the calculation model are set as a/b = 0.1, ρ1/ρ3 = ρ2/ρ3 = 1, and c2 (b)/c3 = 1.For the soft inclusion, the wave velocity in the region ① is set as c1/c3 = c2 (a)/c3 = 0.5, while c1/c3 = c2 (a)/c3 =2 is taken for the hard inclusion.

Effect of the Rigidity Profile of the Inclusion on the Frequency Domain Displacement and DSCF
To reveal the influence of the inhomogeneous rigidity or wave velocity profile on the seismic response of the inclusion, the surface displacement amplitudes as well as the dynamic stress concentration factors of both soft and hard inclusions at different dimensionless frequencies are calculated, as shown in Figures 5-7.The dimensionless parameters in the calculation model are set as a/b = 0.1, ρ1/ρ3 = ρ2/ρ3 = 1, and c2 (b)/c3 = 1.For the soft inclusion, the wave velocity in the region ① is set as c1/c3 = c2 (a)/c3 = 0.5, while c1/c3 = c2 (a)/c3 =2 is taken for the hard inclusion.It is obvious from Figures 5-7 that the right side of soft inclusion experiences a more intense response than that of the hard one in terms of both |u| and DSCF.In contrast, the left sides of both inclusions have a similar dynamic response.To explain this phenomenon, the wave propagation process in the time domain is illustrated as follows.It is obvious from Figures 5-7 that the right side of soft inclusion experiences a more intense response than that of the hard one in terms of both |u| and DSCF.In contrast, the left sides of both inclusions have a similar dynamic response.To explain this phenomenon, the wave propagation process in the time domain is illustrated as follows.

Propagation Process of the Ricker Wavelet in the Full Space Containing Soft and Hard Inclusions
The propagation process of the Ricker wavelet in a full space containing a soft or hard inclusion is shown in Figures 8 and 9, respectively.The dimensionless parameters in the calculation model are consistent with those in the frequency domain analysis.Figure 8 shows the propagation process of incident Ricker wavelet in a full space with a soft inclusion.It is set that the Ricker wavelet reaches the soft inclusion at t = 7.5 s (Figure 8a).When the Ricker wavelet enters the soft inclusion, the wave in the soft inclusion gradually lags behind the wave in the full space due to the difference in wave velocity (c2 < c3).Therefore, the overall wavefront in the soft inclusion becomes concave and curved backward as shown in Figure 8b-e.When the Ricker wavelet passes through the soft core of the inclusion, the core acts as a secondary wave source and generates a scattered wave that begins to propagate outward as shown in Figure 8f.A superposition between the scattered wave and the direct wave can be observed in Figure 8g-i.
The wave progress in Figure 9 for the hard inclusion is different from that in Figure 8 after a similar start at t = 7.5 s.When the Ricker wavelet enters the hard inclusion, the wave in the inclusion propagates faster than that in the full space due to the larger wave velocity c2 of the hard inclusion.The overall waveform is convex and curved forward as shown in Figure 9b-e.Comparing Figure 9f-i to Figure 8f-i, it can be found that the scattered wave generated by the hard inclusion is less obvious than that generated by the soft inclusion.The superposition of energy between the scattered wave by the hard inclusion and the direct wave is not so concentrated comparing to that in the case of soft inclusion.This explains the more intense response of the soft inclusion relative to the hard one.Figure 8 shows the propagation process of incident Ricker wavelet in a full space with a soft inclusion.It is set that the Ricker wavelet reaches the soft inclusion at t = 7.5 s (Figure 8a).When the Ricker wavelet enters the soft inclusion, the wave in the soft inclusion gradually lags behind the wave in the full space due to the difference in wave velocity (c 2 < c 3 ).Therefore, the overall wavefront in the soft inclusion becomes concave and curved backward as shown in Figure 8b-e.When the Ricker wavelet passes through the soft core of the inclusion, the core acts as a secondary wave source and generates a scattered wave that begins to propagate outward as shown in Figure 8f.A superposition between the scattered wave and the direct wave can be observed in Figure 8g-i.
The wave progress in Figure 9 for the hard inclusion is different from that in Figure 8 after a similar start at t = 7.5 s.When the Ricker wavelet enters the hard inclusion, the wave in the inclusion propagates faster than that in the full space due to the larger wave velocity c 2 of the hard inclusion.The overall waveform is convex and curved forward as shown in Figure 9b-e.Comparing Figure 9f-i to Figure 8f-i, it can be found that the scattered wave generated by the hard inclusion is less obvious than that generated by the soft inclusion.The superposition of energy between the scattered wave by the hard inclusion and the direct wave is not so concentrated comparing to that in the case of soft inclusion.This explains the more intense response of the soft inclusion relative to the hard one.

Acceleration Time Histories at Three Representative Positions of Soft and Hard Inclusions
To be more useful in engineering, a procedure based on the transfer function will be adopted in this section to obtain the acceleration response in the time domain.There are three major steps to obtain the results from the frequency domain to the time domain.Firstly, the input displacement time history can be obtained after twice numerical integration on the incident acceleration time history.Secondly, the input displacement Fourier spectrum obtained with the fast Fourier transform technique (FFT) should be multiplied by the frequency domain displacement transfer functions to obtain the response Fourier spectrum at the targeted locations.Finally, through the inverse fast Fourier transform (iFFT) and twice numerical differentiation, the earthquake acceleration response in the time domain can be obtained.In this way, the amplification or de-amplification of the wave at each frequency is taken into account in the site responses.The transfer functions and acceleration time histories of three representative positions (P 1 , O and P 2 in Figure 1 For a specific seismic wave time history, such as Deyang wave and El Centro wave widely used in the literature, it is more than sufficient to compute the transfer functions in its major frequency band within 20 Hz.Assuming the incident Deyang wave and El Centro wave both have a PGA of 0.1 g, their acceleration time histories and Fourier spectrums are both plotted in Figures 10 and 11, respectively.

Acceleration Time Histories at Three Representative Positions of Soft and Hard Inclusions
To be more useful in engineering, a procedure based on the transfer function will be adopted in this section to obtain the acceleration response in the time domain.There are three major steps to obtain the results from the frequency domain to the time domain.Firstly, the input displacement time history can be obtained after twice numerical integration on the incident acceleration time history.Secondly, the input displacement Fourier spectrum obtained with the fast Fourier transform technique (FFT) should be multiplied by the frequency domain displacement transfer functions to obtain the response Fourier spectrum at the targeted locations.Finally, through the inverse fast Fourier transform (iFFT) and twice numerical differentiation, the earthquake acceleration response in the time domain can be obtained.In this way, the amplification or de-amplification of the wave at each frequency is taken into account in the site responses.The transfer functions and acceleration time histories of three representative positions (P1, O and P2 in Figure 1 For a specific seismic wave time history, such as Deyang wave and El Centro wave widely used in the literature, it is more than sufficient to compute the transfer functions in its major frequency band within 20 Hz.Assuming the incident Deyang wave and El Centro wave both have a PGA of 0.1 g, their acceleration time histories and Fourier spectrums are both plotted in Figures 10 and 11, respectively.

Acceleration Time Histories at Three Representative Positions of Soft and Hard Inclusions
To be more useful in engineering, a procedure based on the transfer function will be adopted in this section to obtain the acceleration response in the time domain.There are three major steps to obtain the results from the frequency domain to the time domain.Firstly, the input displacement time history can be obtained after twice numerical integration on the incident acceleration time history.Secondly, the input displacement Fourier spectrum obtained with the fast Fourier transform technique (FFT) should be multiplied by the frequency domain displacement transfer functions to obtain the response Fourier spectrum at the targeted locations.Finally, through the inverse fast Fourier transform (iFFT) and twice numerical differentiation, the earthquake acceleration response in the time domain can be obtained.In this way, the amplification or de-amplification of the wave at each frequency is taken into account in the site responses.The transfer functions and acceleration time histories of three representative positions (P1, O and P2 in Figure 1       To take into account the frequency dependent modification of motions, an acceleration time history synthesizing procedure in Dai et al. [26] is adopted.Figure 13 shows the acceleration time histories at three positions of soft (Figure 13a-c) and hard (Figure 13df) inclusions under the Deyang (NS) wave recorded during the 2008 Wenchuan earthquake in China.It can be found that the PGA at positions P1, O and P2 of the soft inclusion is 1.0, 2.0 and 1.5 times that of the original Deyang (NS) wave (Figure 10a), while the PGA at positions P1, O and P2 of the hard inclusion is 1.0, 0.5 and 0.4 times that of the original wave.Similarly, Figure 14 shows the acceleration time histories of three positions of soft (Figure 14a-c) and hard (Figure 14d-f) inclusions under El Centro wave.It can be found that the PGA at positions P1, O and P2 of the soft inclusion is 1.3, 2.8 and 1.7 times that of To take into account the frequency dependent modification of motions, an acceleration time history synthesizing procedure in Dai et al. [26] is adopted.Figure 13 shows the acceleration time histories at three positions of soft (Figure 13a-c) and hard (Figure 13d-f) inclusions under the Deyang (NS) wave recorded during the 2008 Wenchuan earthquake in China.It can be found that the PGA at positions P 1 , O and P 2 of the soft inclusion is 1.0, 2.0 and 1.5 times that of the original Deyang (NS) wave (Figure 10a), while the PGA at positions P 1 , O and P 2 of the hard inclusion is 1.0, 0.5 and 0.4 times that of the original wave.To take into account the frequency dependent modification of motions, an acceleration time history synthesizing procedure in Dai et al. [26] is adopted.Figure 13 shows the acceleration time histories at three positions of soft (Figure 13a-c) and hard (Figure 13df) inclusions under the Deyang (NS) wave recorded during the 2008 Wenchuan earthquake in China.It can be found that the PGA at positions P1, O and P2 of the soft inclusion is 1.0, 2.0 and 1.5 times that of the original Deyang (NS) wave (Figure 10a), while the PGA at positions P1, O and P2 of the hard inclusion is 1.0, 0.5 and 0.4 times that of the original wave.Similarly, Figure 14 shows the acceleration time histories of three positions of soft (Figure 14a-c) and hard (Figure 14d-f) inclusions under El Centro wave.It can be found that the PGA at positions P1, O and P2 of the soft inclusion is 1.3, 2.8 and 1.7 times that of Similarly, Figure 14 shows the acceleration time histories of three positions of soft (Figure 14a-c) and hard (Figure 14d-f) inclusions under El Centro wave.It can be found that the PGA at positions P 1 , O and P 2 of the soft inclusion is 1.3, 2.8 and 1.7 times that of the original El Centro wave (Figure 10a), while the PGA at positions P 1 , O and P 2 of the hard inclusion is 1.3, 0.7 and 0.3 times that of the original wave.Overall, the amplification or reduction effect on the incident El Centro wave is stronger than that on Deyang (NS) wave due to the difference in the frequency content between the two incident waves.
Sustainability 2022, 14, x FOR PEER REVIEW 14 of 16 the original El Centro wave (Figure 10a), while the PGA at positions P1, O and P2 of the hard inclusion is 1.3, 0.7 and 0.3 times that of the original wave.Overall, the amplification or reduction effect on the incident El Centro wave is stronger than that on Deyang (NS) wave due to the difference in the frequency content between the two incident waves.

Conclusions
An analytical solution to SH wave scattering by a cylindrical inclusion with inhomogeneous elastic modulus in a full space is derived based on the wave function expansion method in this paper.The frequency domain results show that the rigidity profile of the inclusion has a great influence on both the displacement and stress around the inclusion.Generally, a softer inclusion will experience larger displacement and stress due to the superposition between scattered and direct waves, especially in its core and shaded side.Time domain results also illustrate that the soft inclusion may amplify the PGA of acceleration by more than 150% in the shaded side of the inclusion.This warrant special attention paid the modification of earthquake input for underground structures built near inhomogeneous inclusions.Extension of the single layer inclusion model herein to a multilayer one is helpful for studying the inclusion problem with a more generous inhomogeneous profile.

Conclusions
An analytical solution to SH wave scattering by a cylindrical inclusion with inhomogeneous elastic modulus in a full space is derived based on the wave function expansion method in this paper.The frequency domain results show that the rigidity profile of the inclusion has a great influence on both the displacement and stress around the inclusion.Generally, a softer inclusion will experience larger displacement and stress due to the superposition between scattered and direct waves, especially in its core and shaded side.Time domain results also illustrate that the soft inclusion may amplify the PGA of acceleration by more than 150% in the shaded side of the inclusion.This warrant special attention paid on the modification of earthquake input for underground structures built near inhomogeneous inclusions.Extension of the single layer inclusion model herein to a multilayer one is helpful for studying the inclusion problem with a more generous inhomogeneous profile.

Figure 1 .
Figure 1.2D model for plane SH wave scattering by an inhomogeneous cylindrical inclusion in a full space.

Figure 2 .
Figure 2. Variation of wave velocity as x/b for the cross section at y = 0 of the inhomogeneous inclusion.

Figure 1 .
Figure 1.2D model for plane SH wave scattering by an inhomogeneous cylindrical inclusion in a full space.

Figure 1 .
Figure 1.2D model for plane SH wave scattering by an inhomogeneous cylindrical inclusion in a full space.

Figure 2 .
Figure 2. Variation of wave velocity as x/b for the cross section at y = 0 of the inhomogeneous inclusion.

Figure 2 .
Figure 2. Variation of wave velocity as x/b for the cross section at y = 0 of the inhomogeneous inclusion.

4 .
Results and Discussions in Both Frequency and Time Domains 4.1.Effect of the Rigidity Profile of the Inclusion on the Frequency Domain Displacement and DSCF To reveal the influence of the inhomogeneous rigidity or wave velocity profile on the seismic response of the inclusion, the surface displacement amplitudes as well as the dynamic stress concentration factors of both soft and hard inclusions at different dimensionless frequencies are calculated, as shown in Figures 5-7.The dimensionless parameters in the calculation model are set as a/b = 0.1, ρ 1 /ρ 3 = ρ 2 /ρ 3 = 1, and c 2 (b)/c 3 = 1.

Figure 5 .
Figure 5.Comparison of displacement amplitudes |u| and dynamic stress concentration factors (DSCF) between soft (red solid line) and hard (blue dashed line) inclusions at η = 1.0.

Figure 6 .
Figure 6.Comparison of displacement amplitudes |u| and dynamic stress concentration factors (DSCF) between soft (red solid line) and hard (blue dashed line) inclusions at η = 2.0.

Figure 5 .
Figure 5.Comparison of displacement amplitudes |u| and dynamic stress concentration factors (DSCF) between soft (red solid line) and hard (blue dashed line) inclusions at η = 1.0.

Figure 5 .
Figure 5.Comparison of displacement amplitudes |u| and dynamic stress concentration factors (DSCF) between soft (red solid line) and hard (blue dashed line) inclusions at η = 1.0.

Figure 6 .
Figure 6.Comparison of displacement amplitudes |u| and dynamic stress concentration factors (DSCF) between soft (red solid line) and hard (blue dashed line) inclusions at η = 2.0.Figure 6.Comparison of displacement amplitudes |u| and dynamic stress concentration factors (DSCF) between soft (red solid line) and hard (blue dashed line) inclusions at η = 2.0.

Figure 7 .
Figure 7.Comparison of displacement amplitudes |u| and dynamic stress concentration factors (DSCF) between soft (red solid line) and hard (blue dashed line) inclusions at η = 3.0.

Figure 8 .
Figure 8.Time domain snapshots of the full space with a soft inclusion (c 1 /c 3 = 1/2) under the Ricker wavelet.

Figure 9 .
Figure 9.Time domain snapshots of the full space with a hard inclusion (c1/c3 = 2) under the Ricker wavelet.

Figure 9 .
Figure 9.Time domain snapshots of the full space with a hard inclusion (c 1 /c 3 = 2) under the Ricker wavelet.

Figure 10 .
Figure 10.Acceleration time history (a) and spectrum (b) of the Deyang (NS) wave.

Figure 11 .
Figure 11.Acceleration time history (a) and spectrum (b) of the El Centro wave.

Figure 12
Figure 12 shows the displacement transfer functions at three representative positions O, P1 and P2 for both soft (c1/c3 = 1/2) and hard (c1/c3 = 2) inclusions under incident SH wave at varying frequencies.It can be seen from Figure 12a for the soft inclusion case that the displacement transfer function amplitude |u| at position P1 is around 1 and |u| at position O reaches 2 for most frequencies, while |u| at position P2 oscillates between 1 and 2. These values means that there is no amplification of displacement amplitudes at position P1, there is an amplification of 200% at position O and below 200% at position P2 for most frequencies.When it comes to the case of a hard inclusion in Figure 12b, there is no obvious amplification in transfer function at position P1 (|u| ≈ 1) but an obvious reduction at position O (|u| ≈ 0.5) and position P2 (|u| ≈ 0.25) for most frequencies.

Figure 10 .
Figure 10.Acceleration time history (a) and spectrum (b) of the Deyang (NS) wave.
) around both the soft and hard inclusions are presented.The calculation parameters are set as a = 100 m, b = 1000 m, c3 = 1000 m/s, ρ1 = ρ2 = ρ3 = 2700 kg/m 3 .For a specific seismic wave time history, such as Deyang wave and El Centro wave widely used in the literature, it is more than sufficient to compute the transfer functions in its major frequency band within 20 Hz.Assuming the incident Deyang wave and El Centro wave both have a PGA of 0.1 g, their acceleration time histories and Fourier spectrums are both plotted in Figures10 and 11, respectively.

Figure 10 .
Figure 10.Acceleration time history (a) and spectrum (b) of the Deyang (NS) wave.

Figure 11 .
Figure 11.Acceleration time history (a) and spectrum (b) of the El Centro wave.

Figure 12
Figure 12 shows the displacement transfer functions at three representative positions O, P1 and P2 for both soft (c1/c3 = 1/2) and hard (c1/c3 = 2) inclusions under incident SH wave at varying frequencies.It can be seen from Figure 12a for the soft inclusion case that the displacement transfer function amplitude |u| at position P1 is around 1 and |u| at position O reaches 2 for most frequencies, while |u| at position P2 oscillates between 1 and 2. These values means that there is no amplification of displacement amplitudes at position P1, there is an amplification of 200% at position O and below 200% at position P2 for most frequencies.When it comes to the case of a hard inclusion in Figure 12b, there is no obvious amplification in transfer function at position P1 (|u| ≈ 1) but an obvious reduction at position O (|u| ≈ 0.5) and position P2 (|u| ≈ 0.25) for most frequencies.

Figure 11 .
Figure 11.Acceleration time history (a) and spectrum (b) of the El Centro wave.

Figure 12
Figure12shows the displacement transfer functions at three representative positions O, P 1 and P 2 for both soft (c 1 /c 3 = 1/2) and hard (c 1 /c 3 = 2) inclusions under incident SH wave at varying frequencies.It can be seen from Figure12afor the soft inclusion case that the displacement transfer function amplitude |u| at position P 1 is around 1 and |u| at position O reaches 2 for most frequencies, while |u| at position P 2 oscillates between 1 and 2. These values means that there is no amplification of displacement amplitudes at position P 1 , there is an amplification of 200% at position O and below 200% at position P 2 for most frequencies.When it comes to the case of a hard inclusion in Figure12b, there is no obvious amplification in transfer function at position P 1 (|u| ≈ 1) but an obvious reduction at position O (|u| ≈ 0.5) and position P 2 (|u| ≈ 0.25) for most frequencies.

Figure 12 .
Figure 12.Transfer functions of motions at three representative positions of soft (a) and hard (b) inclusions under SH wave.

Figure 13 .
Figure 13.Acceleration time histories of three representative positions for soft (a-c) and hard (d-f) inclusions under Deyang (NS) wave.

Figure 12 .
Figure 12.Transfer functions of motions at three representative positions of soft (a) and hard (b) inclusions under SH wave.

Sustainability 2022 , 16 Figure 12 .
Figure 12.Transfer functions of motions at three representative positions of soft (a) and hard (b) inclusions under SH wave.

Figure 13 .
Figure 13.Acceleration time histories of three representative positions for soft (a-c) and hard (d-f) inclusions under Deyang (NS) wave.

Figure 13 .
Figure 13.Acceleration time histories of three representative positions for soft (a-c) and hard (d-f) inclusions under Deyang (NS) wave.

Figure 14 .
Figure 14.Acceleration time histories of three representative positions for soft (a-c) and hard (d-f) inclusions under El Centro wave.

Figure 14 .
Figure 14.Acceleration time histories of three representative positions for soft (a-c) and hard (d-f) inclusions under El Centro wave.