# Use of Finite Point Method for Wave Propagation in Nonhomogeneous Unbounded Domains

Authors: S. Moazam; Z. Nalbantoglu; M. Celikag
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Wave propagation in an unbounded domain surrounding the stimulation resource is one of the important issues for engineers. Past literature is mainly concentrated on the modelling and estimation of the wave propagation in partially layered, homogeneous, and unbounded domains with harmonic properties. In this study, a new approach based on the Finite Point Method (FPM) has been introduced to analyze and solve the problems of wave propagation in any nonhomogeneous unbounded domain. The proposed method has the ability to use the domain properties by coordinate as an input. Therefore, there is no restriction in the form of the domain properties, such as being periodical as in the case of existing similar numerical methods. The proposed method can model the boundary points between phases with trace of errors and the results of this method satisfy both conditions of decay and radiation.


---

## Body

## 1. Introduction

Wave propagation in the unbounded domains is one of the important engineering issues. To solve this problem, many researches have been carried out and the ones that are most relevant to the proposed method are referred to in this study. Wang et al. [1] tried to simulate wave propagation in domains with nonhomogeneous cross-anisotropic properties. In another research, Wang et al. [2] simulated the wave parameters, such as stress and displacement, in a nonhomogeneous transversely isotropic half-space subjected to a uniform vertical circular load. Daros [3] presented a solution for SH-waves in a nonhomogeneous anisotropic media. Ke et al. [4] worked on simulations of Love waves in a nonhomogeneous, saturated, porous layered half-space with linearly varying properties. In 2006, Bazyar and Song [5] tried to simulate time-harmonic response of nonhomogeneous elastic unbounded domain using Scaled Boundary Finite Element Method. However, the aforementioned studies could not lead to a general method that can simulate unbounded domain with any kind of properties. Furthermore, each one of the mentioned methods has its own limitations in the shape and form of stimulating function problems, which are avoided in the proposed method.In recent years, due to the rapid development in computer power, researchers from various disciplines have developed a special interest in numerical solution to many problems. Numerical modelling of problems has been the major focus of the researchers and undoubtedly wave propagation has been an important part of this numerical research approach. Over the years, physicists also acknowledged the nature of masses not just as particles but also as waves [6] and this emphasizes the importance of the wave propagation modelling.The research by Boroomand and Mossaiby [7] introduced a new method, which is based on Finite Element Method, to estimate the wave propagation in a homogeneous unbounded domain. Lately, Moazam [8] used Finite Point Method (FPM) and Finite Difference Method (FDM) to develop the method of Boroomand and Mossaiby [7] to a meshless version to avoid problems caused by element mesh in simulating the wave propagation in unbounded domain. Then, Moazam et al. [9] further developed the method by Moazam [8] so that it could be used for homogeneous domains due to arbitrary stimulation and also it can reduce the estimation errors.The aim of this study is to investigate the effect of wave propagation due to the vibration of heavy machinery on the surrounding region of the nonhomogeneous arbitrary domain using a new FPM-based numerical method.The assumptions and concepts of them will be explained in the following parts of this study.


## 2. Methodology

### 2.1. Elastic Wave Propagation in Unbounded Domain

The Finite Point Method was used with specific form of generalized coordinate that was developed by Moazam et al. [9] and will be explained as follows. The assumptions and concepts which are used in this study will be explained in the following parts of this study.The wave equations are given below:(1)STDSU-ρU¨=Fx,y,t,x,y∈R2,in which U and U¨ are the magnitude of displacement of wave and its second derivative in respect to time, respectively; S is a differential equation that signifies the relative deformation; D is a matrix of material properties; ρ is the unit weight of the domain; and finally, F is the stimulation function of the domain (a Dirac-Delta function in the specified direction and time with sinusoidal form).One of the uses of the above formula is the elastic wave propagation in which all functions and operators are written in vector format. To solve this equation, a Cartesian coordinate system is adopted, where the center of this coordinate system is used as the stimulation point. IfU is considered as(2)U=ueiωt,then u is a Fourier transformation of U and ω is the value of the stimulation frequency. Consequently, these values are substituted in (1) to obtain (3)STDSu+ρω2u=f,where f is a Fourier transformation of stimulation function of F.According to the stimulation function shape, to solve this problem, symmetric and antisymmetric displacement condition can be used in the domain. Then, the equation is given as follows:(4)STDSu+ρω2u=0,x∈0,∞×0,∞.This study considers the importance of the reliability of the domain properties which can solve the problem. Therefore, the stimulation f can be applied as a boundary condition and thereby (4) can be classified as part of homogeneous equations group with constant coefficients. As a result, one of the significant properties of the differential equations with constant coefficients, such as proportionality, is given in (5)ux,y=Aeαx+βy,where A is a constant vector and α and β are two constant undefined scalars. Equation (6) is derived from the exponential function properties in the x and y direction:(6)ux+nLx,y+mLy=Aeα(x+nLx)+β(y+mLy)=μ1nμ2mux,y,where Lx and Ly are arbitrary specified values in x and y direction and m and n are positive numbers.By substituting (5) into (4), (7) can be obtained:(7)LAeαx+βy=0orL·A=0,where L is a matrix including values based on α and β. The nullspace of a matrix is equivalent to the matrix when it reaches to zero:(8)L=0.According to the characteristic of (7), α and β are the main factors relating to the issue which will be discussed in Results and Discussions of this paper.One of these variables can be calculated in terms of the other one:α=f(β) or β=g(α). It must be noted that, depending on the degree of characteristic of equation, there may be more than one answer to each of these equations.The homogenous solution of this equation may be obtained by using the superposition of spectral solutions. For example,β=g(α) makes [7](9)u=∫α∑iAieαx+βiydα=∫α∑iAieαx+fi(α)ydα.The inner sigma in the overall nullspace of L matrix and the overall integration gives the possible values of α.


### 2.2. Decay and Radiation Condition

Decay condition of amplitude means decreasing the amplitude while increasing the distance from the stimulation point ((n→∞,m→∞)⇒u→0 (6)); that is,(10)μ1<1,μ2<1.One should note that μ1 and μ2 can have complex values; thus, (10) is a circle with the radius of one in Gaussian coordinate as shown in Figure 1.Figure 1
According to decay condition, the domain ofμ1 and μ2 in Gaussian coordinate is a circle with the radius of 1.In wave propagation, problems like radiation condition should be considered. Therefore, given the physical nature of such a problem, the energy emitted towards infinity represents the energy returned from infinity and this changed the shape of the wave and the prerequisite as follows:(11)U=Ae(a+ib)x+(c+id)y+iωt=Aeax+cyei(bx+dy+ωt)a<0,c<0,b<0,d<0.The first two inequalities, a<0 and c<0, are satisfying the amplitude reduction condition. The last part of (11), ei(bx+dy+ωt), is showing the direction of wave movement. The proportional coefficients, μ1 and μ2, can be found in(12)μ1=e-aLx+ibLx=e-aLxe-ibLx,μ2=e-cLy+idLy=e-cLye-idLy,where Lx and Ly are both arbitrary lengths parallel to x and y directions. Therefore, they should be positive values.Each complex number can be shown asreiθ=Re+iIm [10], where r is called the modulus of this complex number, θ is the argument of the complex number, Re is the real part of the complex number, and Im is the value of imaginary part of the complex number.Therefore-bLx and -dLy are in the role of argument in the definition of proportion coefficients. Thus the acceptable range of proportion coefficients of μ1 and μ2 can be seen in Figure 2 or as below set:(13){∀(μ1,μ2)∣Reμ1·Imμ1<0,Reμ2·Imμ2<0}.Figure 2
According to decay and radiation condition, the domain ofμ1 and μ2 in Gaussian coordinate is a circle with the radius of 1.


### 2.3. Finite Point Method

In recent years, Finite Point Method has been developed as one of the numerical methods to solve differential equation problems. Since it is a meshless method, there is no need to carry out mesh generation [11–13] and it is known as a method for avoiding the errors which occur as a result of element networks [7, 9]. Using FPM solution in (4), where a series of regular and equal intervals are connected to each other in both horizontal and vertical directions (unit size), the equation can be given as(14)ux,y≈u^x,y=∑i=1mαifix,y=fTα,where u^ is a set of point value estimations and f=f1f2f3⋯fmT is an appropriate set of basic functions. This set of functions has to be selected from the complete expansion of binomial to the power of natural numbers which is shown as follows:(15)0throw:x+y0=1,1strow:x+y1=1x+1y,2ndrow:x+y2=1x2+2xy+1y2,3rdrow:x+y3=1x3+3x2y+3xy2+1y3,4throw:x+y4=1x4+4x3y+6x2y2+4xy3+1y4,5throw:x+y5=1x5+5x4y+10x3y2+10x2y35throw:x+y5+5xy4+1y5.Since there are nine points in each cloud of computations selected in this method, the selected set should have at least nine terms [14]. Different sets of nine and more terms were tried to be checked for their suitability to be used for the present research. Eventually, nine terms as shown in Figure 3, known as the set with the least error, were used.Figure 3
Selected set of terms to be used for the FPM basic functions.Therefore the set of basic functions can be shown as follows:(16)f=1,x,y,x2,xy,y2,x2y,xy2,x2y2T,m=9.In (14), the values of αi and the unknown values called “generalized coordinates” are the estimates obtained when their functions are determinate. Using this method, the value of αi is determined in the approximate location of the subscales. When the values are equal to the desired function, then the following equation is established:(17)u⌢j=uxj,yj=u-j,j=1,…,9.In (17), u-j is the point value of the function in the jth point. This equation can be developed based on (14) and it is given as (18)fT(xj,yj)α=u-j,j=1,…,9.Accordingly, with this system of equations where both sides are equal, a regular problem can be solved by using Finite Point Method without the need of using other methods, such as least square method [15, 16].Up to this point, the method is similar to FPM proposed by Moazam et al. [9] and Finite Element Method proposed by Boroomand and Mossaiby [7]. But as it was mentioned by these researchers, the methods were just applicable in the homogeneous media. Therefore, there should be some improvements in the numerical method so that it can be used in arbitrary media. This improvement has been discussed by Moazam et al. [9] and presented in this research using a variable Dx,y matrix instead of a constant D matrix similar to what has been shown in [17] by Boroomand and Mossaiby. In the present study, Finite Point Method was used instead of the Finite Element Method.Since the method used in this research, which is based on Finite Point Method, is a frequency domain method, the frequency of vibration in each point is supposed to be equal to the frequency of stimulation resource [8]. The maximum magnitude of wave is the value which can be measured in this method. If the stimulation resource is assumed to be in a sinusoidal form as given by(19)d=dmax⁡sinωt,then the velocity and acceleration can be shown as follows:(20)V=dmax⁡ωcos⁡ωt=Vmax⁡cos⁡ωt⟶Vmax⁡=dmax⁡ω,a=dmax⁡ω2sinωt=amax⁡sinωt⟶amax⁡=dmax⁡ω2.Therefore amax⁡ can be determined as dmax⁡ω2. Since this special form of FPM which was mentioned above can estimate the maximum value of wave magnitude on each point of the domain, the method can also estimate Vmax⁡ and amax⁡ according to (20).


## 2.1. Elastic Wave Propagation in Unbounded Domain

The Finite Point Method was used with specific form of generalized coordinate that was developed by Moazam et al. [9] and will be explained as follows. The assumptions and concepts which are used in this study will be explained in the following parts of this study.The wave equations are given below:(1)STDSU-ρU¨=Fx,y,t,x,y∈R2,in which U and U¨ are the magnitude of displacement of wave and its second derivative in respect to time, respectively; S is a differential equation that signifies the relative deformation; D is a matrix of material properties; ρ is the unit weight of the domain; and finally, F is the stimulation function of the domain (a Dirac-Delta function in the specified direction and time with sinusoidal form).One of the uses of the above formula is the elastic wave propagation in which all functions and operators are written in vector format. To solve this equation, a Cartesian coordinate system is adopted, where the center of this coordinate system is used as the stimulation point. IfU is considered as(2)U=ueiωt,then u is a Fourier transformation of U and ω is the value of the stimulation frequency. Consequently, these values are substituted in (1) to obtain (3)STDSu+ρω2u=f,where f is a Fourier transformation of stimulation function of F.According to the stimulation function shape, to solve this problem, symmetric and antisymmetric displacement condition can be used in the domain. Then, the equation is given as follows:(4)STDSu+ρω2u=0,x∈0,∞×0,∞.This study considers the importance of the reliability of the domain properties which can solve the problem. Therefore, the stimulation f can be applied as a boundary condition and thereby (4) can be classified as part of homogeneous equations group with constant coefficients. As a result, one of the significant properties of the differential equations with constant coefficients, such as proportionality, is given in (5)ux,y=Aeαx+βy,where A is a constant vector and α and β are two constant undefined scalars. Equation (6) is derived from the exponential function properties in the x and y direction:(6)ux+nLx,y+mLy=Aeα(x+nLx)+β(y+mLy)=μ1nμ2mux,y,where Lx and Ly are arbitrary specified values in x and y direction and m and n are positive numbers.By substituting (5) into (4), (7) can be obtained:(7)LAeαx+βy=0orL·A=0,where L is a matrix including values based on α and β. The nullspace of a matrix is equivalent to the matrix when it reaches to zero:(8)L=0.According to the characteristic of (7), α and β are the main factors relating to the issue which will be discussed in Results and Discussions of this paper.One of these variables can be calculated in terms of the other one:α=f(β) or β=g(α). It must be noted that, depending on the degree of characteristic of equation, there may be more than one answer to each of these equations.The homogenous solution of this equation may be obtained by using the superposition of spectral solutions. For example,β=g(α) makes [7](9)u=∫α∑iAieαx+βiydα=∫α∑iAieαx+fi(α)ydα.The inner sigma in the overall nullspace of L matrix and the overall integration gives the possible values of α.


## 2.2. Decay and Radiation Condition

Decay condition of amplitude means decreasing the amplitude while increasing the distance from the stimulation point ((n→∞,m→∞)⇒u→0 (6)); that is,(10)μ1<1,μ2<1.One should note that μ1 and μ2 can have complex values; thus, (10) is a circle with the radius of one in Gaussian coordinate as shown in Figure 1.Figure 1
According to decay condition, the domain ofμ1 and μ2 in Gaussian coordinate is a circle with the radius of 1.In wave propagation, problems like radiation condition should be considered. Therefore, given the physical nature of such a problem, the energy emitted towards infinity represents the energy returned from infinity and this changed the shape of the wave and the prerequisite as follows:(11)U=Ae(a+ib)x+(c+id)y+iωt=Aeax+cyei(bx+dy+ωt)a<0,c<0,b<0,d<0.The first two inequalities, a<0 and c<0, are satisfying the amplitude reduction condition. The last part of (11), ei(bx+dy+ωt), is showing the direction of wave movement. The proportional coefficients, μ1 and μ2, can be found in(12)μ1=e-aLx+ibLx=e-aLxe-ibLx,μ2=e-cLy+idLy=e-cLye-idLy,where Lx and Ly are both arbitrary lengths parallel to x and y directions. Therefore, they should be positive values.Each complex number can be shown asreiθ=Re+iIm [10], where r is called the modulus of this complex number, θ is the argument of the complex number, Re is the real part of the complex number, and Im is the value of imaginary part of the complex number.Therefore-bLx and -dLy are in the role of argument in the definition of proportion coefficients. Thus the acceptable range of proportion coefficients of μ1 and μ2 can be seen in Figure 2 or as below set:(13){∀(μ1,μ2)∣Reμ1·Imμ1<0,Reμ2·Imμ2<0}.Figure 2
According to decay and radiation condition, the domain ofμ1 and μ2 in Gaussian coordinate is a circle with the radius of 1.


## 2.3. Finite Point Method

In recent years, Finite Point Method has been developed as one of the numerical methods to solve differential equation problems. Since it is a meshless method, there is no need to carry out mesh generation [11–13] and it is known as a method for avoiding the errors which occur as a result of element networks [7, 9]. Using FPM solution in (4), where a series of regular and equal intervals are connected to each other in both horizontal and vertical directions (unit size), the equation can be given as(14)ux,y≈u^x,y=∑i=1mαifix,y=fTα,where u^ is a set of point value estimations and f=f1f2f3⋯fmT is an appropriate set of basic functions. This set of functions has to be selected from the complete expansion of binomial to the power of natural numbers which is shown as follows:(15)0throw:x+y0=1,1strow:x+y1=1x+1y,2ndrow:x+y2=1x2+2xy+1y2,3rdrow:x+y3=1x3+3x2y+3xy2+1y3,4throw:x+y4=1x4+4x3y+6x2y2+4xy3+1y4,5throw:x+y5=1x5+5x4y+10x3y2+10x2y35throw:x+y5+5xy4+1y5.Since there are nine points in each cloud of computations selected in this method, the selected set should have at least nine terms [14]. Different sets of nine and more terms were tried to be checked for their suitability to be used for the present research. Eventually, nine terms as shown in Figure 3, known as the set with the least error, were used.Figure 3
Selected set of terms to be used for the FPM basic functions.Therefore the set of basic functions can be shown as follows:(16)f=1,x,y,x2,xy,y2,x2y,xy2,x2y2T,m=9.In (14), the values of αi and the unknown values called “generalized coordinates” are the estimates obtained when their functions are determinate. Using this method, the value of αi is determined in the approximate location of the subscales. When the values are equal to the desired function, then the following equation is established:(17)u⌢j=uxj,yj=u-j,j=1,…,9.In (17), u-j is the point value of the function in the jth point. This equation can be developed based on (14) and it is given as (18)fT(xj,yj)α=u-j,j=1,…,9.Accordingly, with this system of equations where both sides are equal, a regular problem can be solved by using Finite Point Method without the need of using other methods, such as least square method [15, 16].Up to this point, the method is similar to FPM proposed by Moazam et al. [9] and Finite Element Method proposed by Boroomand and Mossaiby [7]. But as it was mentioned by these researchers, the methods were just applicable in the homogeneous media. Therefore, there should be some improvements in the numerical method so that it can be used in arbitrary media. This improvement has been discussed by Moazam et al. [9] and presented in this research using a variable Dx,y matrix instead of a constant D matrix similar to what has been shown in [17] by Boroomand and Mossaiby. In the present study, Finite Point Method was used instead of the Finite Element Method.Since the method used in this research, which is based on Finite Point Method, is a frequency domain method, the frequency of vibration in each point is supposed to be equal to the frequency of stimulation resource [8]. The maximum magnitude of wave is the value which can be measured in this method. If the stimulation resource is assumed to be in a sinusoidal form as given by(19)d=dmax⁡sinωt,then the velocity and acceleration can be shown as follows:(20)V=dmax⁡ωcos⁡ωt=Vmax⁡cos⁡ωt⟶Vmax⁡=dmax⁡ω,a=dmax⁡ω2sinωt=amax⁡sinωt⟶amax⁡=dmax⁡ω2.Therefore amax⁡ can be determined as dmax⁡ω2. Since this special form of FPM which was mentioned above can estimate the maximum value of wave magnitude on each point of the domain, the method can also estimate Vmax⁡ and amax⁡ according to (20).


## 3. Results and Discussions

### 3.1. Developing Formulas for Plain Strain Wave with FPM

The real nodal value vectoru is given below:(21)u≈Nu-,where u- is the nodal wave value vector and N is a matrix containing used shape functions in the FPM.As discussed in Section2.3, for shape functions in the N matrix, the set of functions given in (16) is used.Substituting (21) in the main equation (4) and using least square method which can be used in the FPM, (22) can be obtained:(22)K-ω2Mu-=ftK-ω2M=Kd⟹Kdu-=ft.This resultant equation (22) is separated according to the Dirichlet and Neumann boundary conditions, same as what has been done in [7–9]. Hence, the formulation will be specialized for the FPM method.In Section3.2, Case Studies, the specialized formulations for the plain strain waves in unbounded domain are represented in the form of (23) to (26).


### 3.2. Case Studies

The developed method explained in the present study is capable of estimating the wave propagation in unbounded domains that have nonhomogeneous properties. To evaluate the ability of this method in modelling domains with nonhomogeneous characteristics, discrete Green’s functions were used to model some case studies of elastic shear wave problems in this section.One of the most popular problems occurs while using numerical modelling of the wave propagation in unbounded domains with more than one material scattering at the boundary line between each pair of the materials [18]. Therefore, to evaluate the capability of this method to avoid such problems, the case study should have at least two different materials. The case study which is going to be solved in this part is a problem of elastic shear wave propagation in an unbounded domain with two different materials. The condition of these two distinct materials is as shown in Figure 4. According to the assumed shear modulus and the geometry of the two materials, the central part of the material is stiffer than the surrounding material. The stiffer central part was assumed to be circular to have two axes of symmetry to reduce the time consumption by modelling just a quarter of the domain.Figure 4
Circular part with stiff material.To show the ability of the method to model the domains with two materials having different ratios of shear modulus, it was decided to solve four different problems. Since the materials will have different ratios of shear modulus with respect to each other, the shear modulus of material 1 (the central part) was assumed to be constant at 100 MPa and the shear modulus of the surrounding part, which is the complier in comparison with the central part, was assumed to have different values in each problem. The list of the shear modulus for surrounding part that was used in this method is presented in Table1.Table 1
Shear modulus of surrounding part and the ratio of the shear modulus of two different parts.
Cases
Shear modulus of complier part (MPa)
Shear modulus of complier material
Shear modulus of stiffer material
Case  1
1
0.01
Case  2
10
0.10
Case  3
20
0.20
Case  4
50
0.50Moreover, to be able to compare the results of these four different problems, the radius of stiffer central part (Figure4) was assumed to be constant with the value of 0.5 m.The problem is a simplified version of more general three-dimensional cases represented by the following equation:(23)σij,j+ρω2ui=0,u1=u=0,u2=v=0,u3=w,which can also be shown as the following differential equation:(24)∂σ31∂u1+∂σ32∂u2+ρω2u3=0,σ31=Gu1,u2∂u3∂u1,σ32=Gu1,u2∂u3∂u2.Shear modulus of the material is shown with the notation of G(u1,u2) which is considered variable with respect to the coordinate. Therefore, the following equation is obtained:(25)S=∂∂u1∂∂u2,D(u1,u2)=G(u1,u2)1001.With such a definition of the problem, the essential boundary conditions are those in which w is prescribed and the Neumann boundary conditions are those in which GnTSw, with n, a unit normal to the boundary, is prescribed. For evaluation of discrete Green’s functions, the Neumann boundary conditions are considered as (26)G(x,y)∂w∂u2u2=0,u1≠0=0,G(x,y)∂w∂u1u1=0,u2≠0=0.The harmonic motivation term is applied in u3 direction and the motivation at u1=0 and u2=0. In this example, it is considered that the geometry of material 1 and material 2 is as shown in Figure 4. As previously explained, the shear modulus of the domain is changing based on the material; hence, there are two different shear moduli in the present example. Also the domain density is assumed to be constant in the whole domain. The frequency is selected as ω=1. The numerical solution is performed over an area containing 30 × 30 points. The number of integration points is selected to be 30 for integration on r and 30 for integration on θ. The convergence studies for this method are done and presented in [9]. As a summary of this convergence study, it could be said that the minimum number of points selected for integration on r and for integration on θ to have acceptable convergence is equal to the number of points along x and y. Therefore, since the network of points along each axis of x and y is 30 for the example given in this study, the number of points for integration on r and 30 for integration on θ decided to be 30.Figures5 and 6 represent the results of different problems which were solved with the presented method in this research. In Figure 5, the real part of the response of the domain against a sinusoidal Dirac-Delta motivation at the coordinate center is represented. As it can be seen, there are four different results in Figure 5 depicted which are related to the four different problems mentioned in Table 1. The imaginary parts of all results are presented in Figure 6 in the same order. For better understanding, all of the results are presented in two forms of 3D-graph and contour. The real and imaginary parts of the solution are plotted over an area of 30 × 30 points.Figure 5
Real part of a30×30 points domain response to the sinusoidal Dirac-Delta motivation in the form of 3D-graph and contour with a circular stiffer part with the shear modulus of 100 MPa and radius of 0.5 m (the horizontal axis shows the distance from motivation point). (a) G2=1 (MPa), (b) G2=10 (MPa), (c) G2=20 (MPa), and (d) G2=50 (MPa).
(a)
(b)
(c)
(d)Figure 6
Imaginary part of a30×30 points domain response to the sinusoidal Dirac-Delta motivation in the form of 3D-graph and contour with a circular stiffer part with the shear modulus of 100 MPa and radius of 0.5 m (the horizontal axis shows the distance from motivation point). (a) G2=1 (MPa), (b) G2=10 (MPa), (c) G2=20 (MPa), and (d) G2=50 (MPa).
(a)
(b)
(c)
(d)As it is clear, the results shown in the real part of the answers, Figure5, are the amplitude of the domain response on each point to the point sinusoidal motivation on the center point of the stiffer part. Also the imaginary part of the results shows the value of phase difference between stimulation point’s wave and the wave value on each point of the domain.According to Figures5 and 6, there are no scattering effects on the boundary points of the stiffer central part and this shows the ability of the method to model the domain with more than one material.Moreover, the results of this method have the same change rate due to the change of shear modulus of the domain material. It means the more the value of shear modulus of the domain material, the more the wave length resulting in the domain caused by the same motivation.Also the satisfaction of both decay and radiation conditions clearly can be seen in Figures5 and 6, because in none of the graphs the reflection from boundary solution parts and nonsense amplitude change is seen.


### 3.3. Comparison of Proposed Method with Generalized Haskell Matrix Method

The numerical simulation of elastic wave propagation in nonhomogeneous media based on generalized Haskell matrix method by Jiangfeng and Youming [19] is selected and compared with the method presented in this study. The selected method is able to model the problems which were supposed to be solved with the method presented in this study. By comparing the results obtained from both methods, the efficiency of the present method will be analysed and discussed.In order to have a better view of the comparison, the section graphs of wave, instead of whole three-dimensional graphs of results, are given in Figures7 to 9. Therefore, a straight section through the stimulation point in the domain is shown. Since the problems are all symmetric, a straight section through center point, which represents the whole domain response, is selected.Figure 7
Real part of a domain response to the sinusoidal Dirac-Delta Dirichlet motivation in the form of point graph for the first problem.In Figures7 to 9, the comparison between results of generalized Haskell matrix method [19] and the method which is presented in this paper is shown.According to Figures7, 8, and 9, it can be clearly seen that the method which is presented in this study gives smoother curves than the ones generated by generalized Haskell matrix method despite the fact that wave value does not change suddenly in one material.Figure 8
Real part of a domain response to the sinusoidal Dirac-Delta Dirichlet motivation in the form of point graph for the second problem.Figure 9
Real part of a domain response to the sinusoidal Dirac-Delta Dirichlet motivation in the form of point graph for the third problem.


## 3.1. Developing Formulas for Plain Strain Wave with FPM

The real nodal value vectoru is given below:(21)u≈Nu-,where u- is the nodal wave value vector and N is a matrix containing used shape functions in the FPM.As discussed in Section2.3, for shape functions in the N matrix, the set of functions given in (16) is used.Substituting (21) in the main equation (4) and using least square method which can be used in the FPM, (22) can be obtained:(22)K-ω2Mu-=ftK-ω2M=Kd⟹Kdu-=ft.This resultant equation (22) is separated according to the Dirichlet and Neumann boundary conditions, same as what has been done in [7–9]. Hence, the formulation will be specialized for the FPM method.In Section3.2, Case Studies, the specialized formulations for the plain strain waves in unbounded domain are represented in the form of (23) to (26).


## 3.2. Case Studies

The developed method explained in the present study is capable of estimating the wave propagation in unbounded domains that have nonhomogeneous properties. To evaluate the ability of this method in modelling domains with nonhomogeneous characteristics, discrete Green’s functions were used to model some case studies of elastic shear wave problems in this section.One of the most popular problems occurs while using numerical modelling of the wave propagation in unbounded domains with more than one material scattering at the boundary line between each pair of the materials [18]. Therefore, to evaluate the capability of this method to avoid such problems, the case study should have at least two different materials. The case study which is going to be solved in this part is a problem of elastic shear wave propagation in an unbounded domain with two different materials. The condition of these two distinct materials is as shown in Figure 4. According to the assumed shear modulus and the geometry of the two materials, the central part of the material is stiffer than the surrounding material. The stiffer central part was assumed to be circular to have two axes of symmetry to reduce the time consumption by modelling just a quarter of the domain.Figure 4
Circular part with stiff material.To show the ability of the method to model the domains with two materials having different ratios of shear modulus, it was decided to solve four different problems. Since the materials will have different ratios of shear modulus with respect to each other, the shear modulus of material 1 (the central part) was assumed to be constant at 100 MPa and the shear modulus of the surrounding part, which is the complier in comparison with the central part, was assumed to have different values in each problem. The list of the shear modulus for surrounding part that was used in this method is presented in Table1.Table 1
Shear modulus of surrounding part and the ratio of the shear modulus of two different parts.
Cases
Shear modulus of complier part (MPa)
Shear modulus of complier material
Shear modulus of stiffer material
Case  1
1
0.01
Case  2
10
0.10
Case  3
20
0.20
Case  4
50
0.50Moreover, to be able to compare the results of these four different problems, the radius of stiffer central part (Figure4) was assumed to be constant with the value of 0.5 m.The problem is a simplified version of more general three-dimensional cases represented by the following equation:(23)σij,j+ρω2ui=0,u1=u=0,u2=v=0,u3=w,which can also be shown as the following differential equation:(24)∂σ31∂u1+∂σ32∂u2+ρω2u3=0,σ31=Gu1,u2∂u3∂u1,σ32=Gu1,u2∂u3∂u2.Shear modulus of the material is shown with the notation of G(u1,u2) which is considered variable with respect to the coordinate. Therefore, the following equation is obtained:(25)S=∂∂u1∂∂u2,D(u1,u2)=G(u1,u2)1001.With such a definition of the problem, the essential boundary conditions are those in which w is prescribed and the Neumann boundary conditions are those in which GnTSw, with n, a unit normal to the boundary, is prescribed. For evaluation of discrete Green’s functions, the Neumann boundary conditions are considered as (26)G(x,y)∂w∂u2u2=0,u1≠0=0,G(x,y)∂w∂u1u1=0,u2≠0=0.The harmonic motivation term is applied in u3 direction and the motivation at u1=0 and u2=0. In this example, it is considered that the geometry of material 1 and material 2 is as shown in Figure 4. As previously explained, the shear modulus of the domain is changing based on the material; hence, there are two different shear moduli in the present example. Also the domain density is assumed to be constant in the whole domain. The frequency is selected as ω=1. The numerical solution is performed over an area containing 30 × 30 points. The number of integration points is selected to be 30 for integration on r and 30 for integration on θ. The convergence studies for this method are done and presented in [9]. As a summary of this convergence study, it could be said that the minimum number of points selected for integration on r and for integration on θ to have acceptable convergence is equal to the number of points along x and y. Therefore, since the network of points along each axis of x and y is 30 for the example given in this study, the number of points for integration on r and 30 for integration on θ decided to be 30.Figures5 and 6 represent the results of different problems which were solved with the presented method in this research. In Figure 5, the real part of the response of the domain against a sinusoidal Dirac-Delta motivation at the coordinate center is represented. As it can be seen, there are four different results in Figure 5 depicted which are related to the four different problems mentioned in Table 1. The imaginary parts of all results are presented in Figure 6 in the same order. For better understanding, all of the results are presented in two forms of 3D-graph and contour. The real and imaginary parts of the solution are plotted over an area of 30 × 30 points.Figure 5
Real part of a30×30 points domain response to the sinusoidal Dirac-Delta motivation in the form of 3D-graph and contour with a circular stiffer part with the shear modulus of 100 MPa and radius of 0.5 m (the horizontal axis shows the distance from motivation point). (a) G2=1 (MPa), (b) G2=10 (MPa), (c) G2=20 (MPa), and (d) G2=50 (MPa).
(a)
(b)
(c)
(d)Figure 6
Imaginary part of a30×30 points domain response to the sinusoidal Dirac-Delta motivation in the form of 3D-graph and contour with a circular stiffer part with the shear modulus of 100 MPa and radius of 0.5 m (the horizontal axis shows the distance from motivation point). (a) G2=1 (MPa), (b) G2=10 (MPa), (c) G2=20 (MPa), and (d) G2=50 (MPa).
(a)
(b)
(c)
(d)As it is clear, the results shown in the real part of the answers, Figure5, are the amplitude of the domain response on each point to the point sinusoidal motivation on the center point of the stiffer part. Also the imaginary part of the results shows the value of phase difference between stimulation point’s wave and the wave value on each point of the domain.According to Figures5 and 6, there are no scattering effects on the boundary points of the stiffer central part and this shows the ability of the method to model the domain with more than one material.Moreover, the results of this method have the same change rate due to the change of shear modulus of the domain material. It means the more the value of shear modulus of the domain material, the more the wave length resulting in the domain caused by the same motivation.Also the satisfaction of both decay and radiation conditions clearly can be seen in Figures5 and 6, because in none of the graphs the reflection from boundary solution parts and nonsense amplitude change is seen.


## 3.3. Comparison of Proposed Method with Generalized Haskell Matrix Method

The numerical simulation of elastic wave propagation in nonhomogeneous media based on generalized Haskell matrix method by Jiangfeng and Youming [19] is selected and compared with the method presented in this study. The selected method is able to model the problems which were supposed to be solved with the method presented in this study. By comparing the results obtained from both methods, the efficiency of the present method will be analysed and discussed.In order to have a better view of the comparison, the section graphs of wave, instead of whole three-dimensional graphs of results, are given in Figures7 to 9. Therefore, a straight section through the stimulation point in the domain is shown. Since the problems are all symmetric, a straight section through center point, which represents the whole domain response, is selected.Figure 7
Real part of a domain response to the sinusoidal Dirac-Delta Dirichlet motivation in the form of point graph for the first problem.In Figures7 to 9, the comparison between results of generalized Haskell matrix method [19] and the method which is presented in this paper is shown.According to Figures7, 8, and 9, it can be clearly seen that the method which is presented in this study gives smoother curves than the ones generated by generalized Haskell matrix method despite the fact that wave value does not change suddenly in one material.Figure 8
Real part of a domain response to the sinusoidal Dirac-Delta Dirichlet motivation in the form of point graph for the second problem.Figure 9
Real part of a domain response to the sinusoidal Dirac-Delta Dirichlet motivation in the form of point graph for the third problem.


## 4. Conclusions

In this study, an improved version of the method previously presented by Moazam et al. [9] is proven to have the capability of modelling the wave propagation in an unbounded, nonhomogeneous domain. The ability of the method was checked with using more than one material with different properties. According to the results obtained and the comparison of the findings, the following are the conclusions.(i)
Since there is no scattering in the results, it can be said that the method developed in this study can clearly model a multiphase domain as well as a homogeneous domain.(ii)
Due to the results depicted in Figures5 and 6 which were obtained for different materials with different shear modulus, the method is capable of showing the change of wave length occurring with the change of shear modulus of the domain.(iii)
The results obtained in Figures5 to 6 indicate that the method is able to model the unbounded domain without any effect of reflecting waves from outgoing boundaries.(iv)
According to the results obtained, the more the distance from the stimulation point, the less the wave displacement value indicating that the developed method has capability of showing decay in the wave propagation modelling.(v)
According to the results obtained in the case study, the method is able to model the multiphase domain with trace of error in the phase boundary. Also the results indicate that decay and radiation conditions are being met.(vi)
When present method is compared with the generalized Haskell matrix method, the results showed approximately the same wave lengths for different problems. However, the curves obtained by using the method presented in this study are more smooth and closer to the real behaviour than those obtained from the generalized Haskell matrix method.



---

*Source: 914207-2015-03-10.xml*