# Wave Propagation in Unbounded Domains under a Dirac Delta Function with FPM

Authors: S. Moazam; B. Boroomand; S. Naimi; M. Celikag
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Wave propagation in unbounded domains is one of the important engineering problems. There have been many attempts by researchers to solve this problem. This paper intends to shed a light on the finite point method, which is considered as one of the best methods to be used for solving problems of wave propagation in unbounded domains. To ensure the reliability of finite point method, wave propagation in unbounded domain is compared with the sinusoidal unit point stimulation. Results indicate that, in the case of applying stimulation along one direction of a Cartesian coordinate, the results of finite point method parallel to the stimulation have less error in comparison with the results of finite element method along the same direction with the same stimulation.


---

## Body

## 1. Introduction

The rapid development of computers and computation power within the last decade encouraged researchers from different disciplines to show more interest in the usage of numerical methods. Wave propagation is one of those numerical modeling problems which have been a major focus of some of the researchers. However, it is only in recent years that physicists had acknowledged the nature of masses not just as particles but also as waves [1] and they emphasized the importance of the wave propagation modeling. The methods used to solve the differential equations have been categorized into two groups: with or without mesh network [2]. Previous studies showed that using mesh in modeling the wave propagation may cause wave to emanate lead [3]. According to Fatahpour [3], this lead is caused by the shape of the elements and their positioning with respect to each other. In addition, Gerdes and Ihlenburg [4] and Harari and Nogueira [5] highlighted the effects of the shape function problem of the elements used for modeling wave propagation in unbounded domains in their studies. Furthermore, the finite element modeling of wave propagation resulted in the phase difference problems of response, numerical approximation, and pollution error [4].Accordingly, based on the problems of network in wave propagation, there are two methods that can be used for solving the wave propagation problems. Meshless method offers solutions despite the problems associated with its use, such as singularity of stiffness matrices, nonstability, and difficulties in ensuring the accuracy of the number of points in the domain. On the other hand, finite difference method, which is one of the oldest numerical methods, can also be used to solve the problems caused by meshes in modeling wave propagations. This method is limited due to the need for a regular grid of points in an infinite domain. However, these problems can be resolved by using a special storage combination and replication in other parts of the environment.The following methods are often preferred to the ones mentioned earlier since they are very successful for large quantity of numerical modeling of unbounded domains.(a)
Methods which are based on boundary integral equations: according to Kirsch [6], this method has some limitations associated with the properties of domain, such as homogeneous, isotropic, and linear. This method can further be classified into two subgroups: direct and indirect integral equations that are dealing with the physical [7] and mathematical aspects [8], respectively. Therefore, boundary element method has been used successfully to solve the problems with unbounded equations [9–12]. However, there are disadvantages of this method which include the inaccessibility to basic functions of different problems, such as nonhomogeneous domains and complicated calculations that sometimes trigger the singularity of integrals.
(b)
As demonstrated by the first monograph in the world [13], dynamic and transient infinite elements have been developed to solve wave propagation and a broad range of scientific and engineering problems [14, 15]. Zhao et al. established the coupled method of finite and dynamic infinite elements [16, 17] for solving wave scattering problems associated with many real scientific and engineering problems involving semi-infinite and infinite domains, for example,
(i)
dynamic concrete gravity dam-foundation interaction and dynamic embankment dam-foundation interaction problems during earthquakes [18, 19],
(ii)
seismic free field distributions along the surfaces of natural canyons [20, 21],
(iii)
dynamic interactions between three-dimensional framed structures and their foundations [22],
(iv)
dynamic interactions between concrete retaining walls and their foundations [23]. In addition, Zhao and Valliappan also developed the coupled method of finite and transient infinite elements for solving transient seepage flow, heat transfer, and mass transport problems involving semi-infinite and infinite domains [24–26].
(c)
Nonreflecting boundary conditions are shape based by placing B on virtual boundary around the stimulation reservoirΩ (Figure 1) in such a way to allow for the waves to go outward without any reflection inside. Therefore, it is costly to simulate the full infinite domain.
At a glance, this kind of simulation seems easy and simple to perform. But research conducted for the past thirty years has shown that such boundary simulation is hard to perform. In addition, the limited numerical solutions available so far also indicate existence of possible problems with such boundary simulations [27–29] and researchers do not have a consensus on this matter [30]. Therefore, recent studies are aimed at achieving better developed stimulations [31–34].
Absorbing layer or perfectly matched layer method was first introduced by Berenger in 1994 [35] upon completion of the nonreflecting boundaries. Recently, extensive studies have been conducted on how to develop this method for 2- and 3-dimensional domains [36].
(d)
Dynamic solution of unbounded domains using finite element method was first introduced by Boroomand and Mossaiby [37]. In this research the method is further developed to solve the wave leading problem caused by element arrangement and shape functions.Figure 1
Stimulation and nonreflecting boundary.


## 2. Materials and Methods

### 2.1. Elastic Wave Propagation in Unbounded Domain

In this research work, finite point and finite element methods were used to study the wave propagation in unbounded domain [37]. The wave equations are given below:
(1)
S
T
D
S
U
-
ρ
U
¨
=
F
(
x
,
y
,
t
)
,
(
x
,
y
)
∈
R
2
,
where U and U
¨ are the value of wave function and the second derivative of wave function of 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, while the center of this coordinate system is used as the stimulation point. IfU is considered as
(2)
U
=
u
e
i
ω
t
,
then u is a Fourier transformation of U, i
=
-
1, and ω is the value of the stimulation frequency. Consequently, these values are substituted in (2) to obtain
(3)
S
T
D
S
u
+
ρ
ω
2
u
=
f
.
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)
S
T
D
S
u
+
ρ
ω
2
u
=
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)
u
(
x
,
y
)
=
A
e
α
x
+
β
y
.
A, α, and β are constant vector and two undefined scalars. The following equation is derived from the exponential function properties in the x- and y-direction:
(6)
u
(
x
+
n
L
x
,
y
+
m
L
y
)
=
A
e
α
(
x
+
n
L
x
)
+
β
(
y
+
m
L
y
)
=
(
μ
1
)
n
(
μ
2
)
m
u
(
x
,
y
)
,
where L
x and L
y are arbitrary specified values in x- and y-direction and m and n are positive numbers.By substituting (5) into (4), the following is obtained:
(7)
L
A
e
α
x
+
β
y
=
0
or
L
·
A
=
0
.
L is a matrix including values based on α and β. The nullspace of a matrix is equivalent to the matrix when it reaches zero:
(8)
|
L
|
=
0
.
According to the characteristic of (7), α and β are the main factors relating to the issue discussed in the Results and Discussions part 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
(
α
) like the following [37]:
(9)
u
=
∫
α
∑
i
A
i
e
α
x
+
β
i
y
d
α
=
∫
α
∑
i
A
i
e
α
x
+
f
i
(
α
)
y
d
α
.
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 amplitude with 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 a Gaussian coordinate. 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 as well as the prerequisite [37] as follows:
(11)
U
=
A
e
(
a
+
i
b
)
x
+
(
c
+
i
d
)
y
+
i
ω
t
=
A
e
a
x
+
c
y
e
i
(
b
x
+
d
y
+
ω
t
)
=
A
e
(
a
+
)
x
+
(
c
+
i
d
)
ω
t
a
<
0
,
b
<
0
,
c
<
0
,
d
<
0
.


### 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 [38–40] and it is known to be the best method for avoiding the errors which occur as a result of element networks [7, 37]. Using finite point solution in (4), where 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
(12)
u
(
x
,
y
)
≈
u
^
(
x
,
y
)
=
∑
i
=
1
m
α
i
f
i
(
x
,
y
)
=
f
T
α
.
u
^ is a set of point value estimation and f
=
[
f
1
f
2
f
3
⋯
f
m
]
T is an appropriate set of basic functions. In this paper, functions are selected as follows:
(13)
f
=
[
1
,
x
,
y
,
x
2
,
x
y
,
y
2
,
x
2
y
,
x
y
2
,
x
2
y
2
]
T
,
m
=
9
.In (12), 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, the following equation is established:
(14)
u
⌢
j
=
u
(
x
j
,
y
j
)
=
u
-
j
,
j
=
1
,
…
,
9
.
In the above equation u
-
j is the point value of the function in the jth point. This equation can be developed based on (12) and it is given as
(15)
f
T
(
x
j
,
y
j
)
α
=
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 [41, 42].


## 2.1. Elastic Wave Propagation in Unbounded Domain

In this research work, finite point and finite element methods were used to study the wave propagation in unbounded domain [37]. The wave equations are given below:
(1)
S
T
D
S
U
-
ρ
U
¨
=
F
(
x
,
y
,
t
)
,
(
x
,
y
)
∈
R
2
,
where U and U
¨ are the value of wave function and the second derivative of wave function of 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, while the center of this coordinate system is used as the stimulation point. IfU is considered as
(2)
U
=
u
e
i
ω
t
,
then u is a Fourier transformation of U, i
=
-
1, and ω is the value of the stimulation frequency. Consequently, these values are substituted in (2) to obtain
(3)
S
T
D
S
u
+
ρ
ω
2
u
=
f
.
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)
S
T
D
S
u
+
ρ
ω
2
u
=
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)
u
(
x
,
y
)
=
A
e
α
x
+
β
y
.
A, α, and β are constant vector and two undefined scalars. The following equation is derived from the exponential function properties in the x- and y-direction:
(6)
u
(
x
+
n
L
x
,
y
+
m
L
y
)
=
A
e
α
(
x
+
n
L
x
)
+
β
(
y
+
m
L
y
)
=
(
μ
1
)
n
(
μ
2
)
m
u
(
x
,
y
)
,
where L
x and L
y are arbitrary specified values in x- and y-direction and m and n are positive numbers.By substituting (5) into (4), the following is obtained:
(7)
L
A
e
α
x
+
β
y
=
0
or
L
·
A
=
0
.
L is a matrix including values based on α and β. The nullspace of a matrix is equivalent to the matrix when it reaches zero:
(8)
|
L
|
=
0
.
According to the characteristic of (7), α and β are the main factors relating to the issue discussed in the Results and Discussions part 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
(
α
) like the following [37]:
(9)
u
=
∫
α
∑
i
A
i
e
α
x
+
β
i
y
d
α
=
∫
α
∑
i
A
i
e
α
x
+
f
i
(
α
)
y
d
α
.
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 amplitude with 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 a Gaussian coordinate. 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 as well as the prerequisite [37] as follows:
(11)
U
=
A
e
(
a
+
i
b
)
x
+
(
c
+
i
d
)
y
+
i
ω
t
=
A
e
a
x
+
c
y
e
i
(
b
x
+
d
y
+
ω
t
)
=
A
e
(
a
+
)
x
+
(
c
+
i
d
)
ω
t
a
<
0
,
b
<
0
,
c
<
0
,
d
<
0
.


## 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 [38–40] and it is known to be the best method for avoiding the errors which occur as a result of element networks [7, 37]. Using finite point solution in (4), where 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
(12)
u
(
x
,
y
)
≈
u
^
(
x
,
y
)
=
∑
i
=
1
m
α
i
f
i
(
x
,
y
)
=
f
T
α
.
u
^ is a set of point value estimation and f
=
[
f
1
f
2
f
3
⋯
f
m
]
T is an appropriate set of basic functions. In this paper, functions are selected as follows:
(13)
f
=
[
1
,
x
,
y
,
x
2
,
x
y
,
y
2
,
x
2
y
,
x
y
2
,
x
2
y
2
]
T
,
m
=
9
.In (12), 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, the following equation is established:
(14)
u
⌢
j
=
u
(
x
j
,
y
j
)
=
u
-
j
,
j
=
1
,
…
,
9
.
In the above equation u
-
j is the point value of the function in the jth point. This equation can be developed based on (12) and it is given as
(15)
f
T
(
x
j
,
y
j
)
α
=
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 [41, 42].


## 3. Results and Discussions

The points consist of operator identified in (4) and point 5 in Figure 2 and then using numerical results equation (16a) and (16b) will be given as(16a)
u
5
=
E
v
1
8
(
1
-
ν
)
+
E
u
4
1
-
ν
2
-
E
v
7
8
(
1
-
ν
)
+
E
u
2
2
(
1
+
ν
)
+
[
E
(
ν
-
3
)
1
-
ν
2
+
α
]
u
5
+
E
u
8
2
(
1
+
ν
)
-
E
v
3
8
(
1
-
ν
)
+
E
u
6
1
-
ν
2
+
E
v
9
8
(
1
-
ν
)
,
(16b)
v
5
=
E
u
1
8
(
1
-
ν
)
+
E
v
4
1
-
ν
2
-
E
u
7
8
(
1
-
ν
)
+
E
v
2
2
(
1
+
ν
)
+
[
E
(
ν
-
3
)
1
-
ν
2
+
α
]
v
5
+
E
v
8
2
(
1
+
ν
)
-
E
u
3
8
(
1
-
ν
)
+
E
v
6
1
-
ν
2
+
E
u
9
8
(
1
-
ν
)
.Figure 2
Arrangement of points in the operation related to each point.E and ν are modulus of elasticity and Poisson’s ratio of the domain that will allow the wave propagation. From these two equations, the method discussed in Section 2 with respect to (6) is given as(17)
[
α
+
E
(
ν
-
3
)
1
-
ν
2
+
E
(
1
-
ν
2
)
μ
1
+
E
μ
1
1
-
ν
2
+
E
2
(
1
+
ν
)
μ
2
+
E
μ
2
2
(
1
+
ν
)
E
8
(
1
-
ν
)
(
1
μ
1
μ
2
-
μ
1
μ
2
-
μ
2
μ
1
+
μ
1
μ
2
)
E
8
(
1
-
ν
)
(
1
μ
1
μ
2
-
μ
1
μ
2
-
μ
2
μ
1
+
μ
1
μ
2
)
α
+
E
(
ν
-
3
)
1
-
ν
2
+
E
(
1
+
ν
)
μ
1
+
E
μ
1
1
+
ν
+
E
2
(
1
-
ν
2
)
μ
2
+
E
μ
2
2
(
1
-
ν
2
)
]
.Hence, the interrelation with the previous studies is determined by using characteristic of (17), where α is the same as ρ
ω
2 in (4) and the values of μ
1 and μ
2 are based on the infinite element methods.The solutions developed [37] in Section (a) are the exact solution in Section (b). The solutions in Section (c) are used in Figures 3, 4, 5, and 6 by using ρ
ω
2
=
100.Real part ofx-direction response. (a) Using method in [37]. (b) Exact response. (c) Present method.
(a)
(b)
(c)Imaginary part ofx-direction response. (a) Using method in [37]. (b) Exact response. (c) Present method.
(a)
(b)
(c)Real part ofy-direction response. (a) Using method in [37]. (b) Exact response. (c) Present method.
(a)
(b)
(c)Imaginary part ofy-direction response. (a) Using method in [37]. (b) Exact response. (c) Present method.
(a)
(b)
(c)


## 4. Conclusion

The method developed by Boroomand and Mossaiby [37] is described together with the finite point method which is an alternative method to the finite element method. This paper shows the ability of the new method in solving problems for the infinite domains with homogeneous properties. Semifinite media can also be observed with this method using appropriate boundary conditions.Based on the numerical results, the following are the points concluded.(i)
Discreet Green’s functions [37] can easily be estimated with finite point method.
(ii)
Results obtained through this new method have pollution error like the basic finite element method. However, in the finite element usage, wave lengths in the pressure direction are increased when compared with exact solutions. On the other hand, the finite point method increases the wave lengths in the shear direction when compared with exact solutions.
(iii)
There are two kinds of wave propagation problems commonly encountered in the engineering practice [2, 3]. One of them is the wave radiation problem (the machine foundation vibration is an example of this kind) [2] and the other one is the wave scattering problem (the seismic response of a structure is an example of this kind) [3]. Since the source of vibration should be obtained as a boundary condition, only the first kind of wave propagation, which is wave radiation, can be observed in this method.



---

*Source: 470346-2014-02-06.xml*