第1页 / 共17页
第2页 / 共17页
第3页 / 共17页
第4页 / 共17页
第5页 / 共17页
第6页 / 共17页
第7页 / 共17页
第8页 / 共17页
论文研究 - 广义有限差分法（GFDM）在地球物理测试建模中的应用.pdf
The Application of the Generalized Finite Difference Method (GFDM) for Modelling Geophysical Test
Abstract
Keywords
1. Introduction
2. Fundamentals of the GFDM
2.1. Explicit Formulae
2.2. Equation of Motion, Boundary Conditions and PML
2.3. Explicit Generalized Differences
2.4. Stability of the GFD Scheme and Stability of PML Regions
2.5. Heterogeneous Formulation
3. Geophysical Investigation Modelling
3.1. Model of a Cross-Hole Test
3.2. Model of Seismic Refraction
4. Discussion
4.1. Model of a Cross-Hole Test
4.2. Model of Seismic Refraction
5. Conclusion
Acknowledgements
Conflicts of Interest
References
Journal of Geoscience and Environment Protection, 2019, 7, 1-17
http://www.scirp.org/journal/gep
ISSN Online: 2327-4344
ISSN Print: 2327-4336
The Application of the Generalized Finite
Difference Method (GFDM) for Modelling
Geophysical Test
Angel Muelas1, Eduardo Salete1, Juan José Benito1, Francisco Ureña2, Luis Gavete3, Miguel Ureña1
1UNED, ETS Ingenieros Industriales, C/Juan del Rosal, Madrid, Spain
2UCLM, IMACI, Ciudad Real, Spain
3UPM, ETSIM, Madrid, Spain
How to cite this paper: Muelas, A., Salete,
E., Benito, J. J., Ureña, F., Gavete, L., &
Ureña, M. (2019). The Application of the
Generalized Finite Difference Method
(GFDM) for Modelling Geophysical Test.
Journal of Geoscience and Environment
Protection, 7, 1-17.
https://doi.org/10.4236/gep.2019.74001
Received: February 20, 2019
Accepted: April 7, 2019
Published: April 10, 2019
Copyright © 2019 by author(s) and
Scientific Research Publishing Inc.
This work is licensed under the Creative
Commons Attribution International
License (CC BY 4.0).
http://creativecommons.org/licenses/by/4.0/
Open Access
Abstract
The possibility of using a nodal method allowing irregular distribution of
nodes in a natural way is one of the main advantages of the generalized finite
difference method (GFDM) with regard to the classical finite difference me-
thod. Moreover, this feature has made it one of the most-promising meshless
methods because it also allows us to reduce the time-consuming task of mesh
generation and the numerical solution of integrals. This characteristic allows
us to shape geological features easily whilst maintaining accuracy in the re-
sults, which can be a source of great interest when dealing with this kind of
problems. Two widespread geophysical investigation methods in civil engi-
neering are the cross-hole method and the seismic refraction method. This
paper shows the use of the GFDM to model the aforementioned geophysical
investigation tests showing precision in the obtained results when comparing
them with experimental data.
Keywords
Meshless Methods, Generalized Finite Difference Method, Geophysics
1. Introduction
Geotechnical geophysical tools are routinely used to image the subsurface of the
Earth and determine site geology, stratigraphy, and rock quality. Commonly
employed geophysical investigation methods include seismic refraction, seismic
reflection, MASW, cross-hole seismic tomography, electrical resistivity, GPR
(Ground Penetrating Radar), electromagnetics, gravity, etc.
DOI: 10.4236/gep.2019.74001 Apr. 10, 2019
1
Journal of Geoscience and Environment Protection
A. Muelas et al.
DOI: 10.4236/gep.2019.74001
The cross-hole seismic technique determines the compressional (P-) and/or
shear (S-) wave velocity of materials at different depths. Cross-hole testing takes
advantage of generating and recording (seismic) body waves, both the P-waves
and S-waves, at selected depth intervals where the source and receiver(s) are
maintained at equal elevations for each measurement.
P-waves are generated with a sparker or small explosive device such that along
the assumed straight-ray propagation path the seismic impulse compresses and
rarefies the materials radially toward the receiver borehole(s). S-waves generated
in cross-hole testing may be split into two wave types, each with different par-
ticle motions—SV- and SH-waves, respectively.
The seismic refraction method is used to estimate depths and seismic veloci-
ties of geological layers. Seismic velocities can assist in the interpretation of the
geological profile as well as the evaluation of the rippability of bedrock.
In the seismic refraction method, a seismic source (a hammer hitting on a
plate, an explosive, etc.) is used to generate body waves, most commonly,
P-waves.
The wave arrival time is measured by a series of evenly spaced geophones lo-
cated on the ground surface.
The seismic waves propagate downwards through the ground until they are
reflected or refracted off subsurface layers. Refracted waves are detected by ar-
rays of a number of geophones spaced at regular intervals of 1 - 10 meters, de-
pending on the desired depth penetration of the survey.
Some head waves enter a high velocity medium near the critical angle and
travel in the high velocity medium nearly parallel to the interface between layers.
Then, the wave refracted along that interface will overtake the direct wave at
some distance from the source. This point at which the refracted wave overtakes
the direct wave arrival is known as the critical distance and is used to estimate
the depth to the refracting surface.
The GFDM has been proven to be a successful method to model seismic wave
propagation using regular or irregular meshes. Researchers as Jensen (1972),
Liszka and Orkisz (1980), Orkisz (1998) and Perrone and Kao (1975) have con-
tributed to developing the GFDM in different aspects of its applications. Benito,
Gavete, & Ureña (2001) have developed the explicit formulae and have also ap-
plied this meshless method to solve the advection diffusion equation (Ureña,
Benito, & Gavete, 2011) to solve parabolic and hyperbolic equations (Benito,
Ureña, & Gavete, 2007) and to solve second order non-linear elliptic partial dif-
ferential equations (Gavete et al., 2017).
The application of GFDM to the solution of the problem of seismic wave
propagation was introduced by Benito et al. (2017) where a GFD scheme for
elastic, homogeneous isotropic medium was proposed and the stability and dis-
persion studied. A perfectly matched layer (PML) absorbing boundary was in-
cluded in the numerical model by Benito et al. (2013), its stability was analyzed
in (Salete et al., 2017) and the influence of the topography and domain irregular-
2
Journal of Geoscience and Environment Protection
A. Muelas et al.
ities was shown in (Benito et al., 2015). An adaptive process to minimize the in-
fluence of the irregularity of the star nodes distribution on the dispersion was
proposed in (Salete et al., 2011). New schemes for seismic wave propagation in
Kelvin-Voight viscoelastic media and some improvements in the heterogeneous
schemes to deal with irregular interfaces were proposed in (Benito et al., 2013),
(Benito et al., 2018).
Moreover, the application of the finite-difference method to the seismic wave
propagation and earthquake problems has been also tackled by authors as Boore,
Levander, Carcione and Moczo (1998) among others.
In this paper, the GFDM has been applied to model two widespread geophys-
ical tests, namely the cross-hole method and the Seismic refraction method try-
ing to analyze the possibility of using it in other geophysical problems.
The paper is organized as follows. In Section 2 the formulation of GFDM to
solve seismic wave propagation problem containing the stability condition and
the scheme in PML has been included. The geophysical investigation modelling
is shown in Section 3 and the obtained results are analyzed and evaluated in Sec-
tion 4. Finally, in Section 5, some conclusions are obtained.
2. Fundamentals of the GFDM
2.1. Explicit Formulae
Let
2RΩ ⊂
be a domain and
{
M x
1,
=
x
, N
}
⊂ Ω
a discretization of the domain Ω with N points. Each point of the discretization
M is denoted as a node.
For each one of the nodes of the domain, where the value of U is unknown,
with the
Es-star is defined as a set of selected points
E
s
)
central node x0 ∈ M and
is a set of points located in the
s M
neighborhood of x0. In order to select the points, different criteria as the four
quadrants or distance can be used (Benito et al., 2001).
⊂
x x
1,
0
M
x
s
x
i
∈
1,
{
(
i
}
,
,
=
,
=
We consider an Es-star with the central node x0, where (x0, y0) are the coordi-
nates of the central node, (xi, yi) are the coordinates of the ith node in the Es-star,
and hi = xi − x0 and ki = yi − y0. If U0 = U(x0) is the value of the function at the
central node of the star and Ui = U(xi) are the function values at the rest of the
= , then, according to the Taylor series expansion it is
nodes, with
known that:
1,
s
i
,
U U h
i
=
+
0
i
+
1
2
2
h
i
U
∂
0
x
∂
2
U
∂
2
x
∂
+
k
i
0
+
k
U
∂
0
y
∂
2
U
∂
2
y
∂
2
i
0
+
2
h k
i
i
2
U
∂
0
x y
∂ ∂
+
,
i
=
1,
,
s
(1)
If in Equation (1) the higher than second order terms are ignored, a second
order approximation for the Ui function is obtained. This is indicated as ui. It is
then possible to define the function B(u) as:
3
Journal of Geoscience and Environment Protection
DOI: 10.4236/gep.2019.74001
A. Muelas et al.
( )
B u
=
∑
s
i
1
=
(
u
0
−
u
i
)
+
(2)
+
k
i
u
∂
0
y
∂
+
2
h k
i
i
2
u
∂
0
x y
∂ ∂
2
2
ω
i
h
i
u
∂
0
x
∂
2
u
∂
0
2
y
∂
+
1
2
2
h
i
2
u
∂
0
2
x
∂
+
k
2
i
(
)
i
i
=
w w h k
,
i
are positive symmetrical weighting functions, having the
where
property as defined in Lancaster and Salkauskas (1986) that are functions de-
creasing monotonically in magnitude as the distance to the center of the weight-
ing function increases. See also Levin (1998) for more detail.
1
dist
Some weighting functions as potential
can be
used (Benito et al., 2001), where n ∈ . If the norm given by Equation (2) is
minimized with respect to the partial derivatives, the following system of linear
equations is obtained (for each node of the discretization)
or exponential
e n dist
− ⋅
n
2
=AD b (3)
The matrix A is a positive definite matrix and the approximation is of second
order,
(
2
h kθ
,i
2
i
)
(Gavete et al., 2017).
The solution of system Equation (3) is unique and it is a linear combination of
the function values obtained at the nodes. Then
+ ∑
s
i
1
=
1
−
D A b
0 0
m u
= −
=
i
i
m u
(4)
with
(5)
The explicit expressions of matrices A, D, b and coefficients m0, mi may be
m
= ∑
i
s
m=
i
1
0
seen in (Benito et al., 2001; Benito et al., 2015; Benito et al., 2007).
2.2. Equation of Motion, Boundary Conditions and PML
The equations of motion for a perfectly elastic, homogeneous and isotropic me-
dium in the domain
(6)
The Ui are the components of the displacements, ρ is the density and λ and µ
U
µ
U
k ik
,
i kk
,
+
+
i tt
,
2RΩ ⊂
U
ρ
are
)
(
λ µ
=
are Lame’s elastic coefficients.
In this paper two types of boundary conditions are considered: Dirichlet
boundary conditions and free surface. On the free surface, the following condi-
tions we impose
σ
ij
⋅
n
j
=
)
(
λ δ µ ⋅
kk
ij
+
2
ij
n
j
=
( )
g t
i
(7)
where gi(t) is equal to zero when there are not forces on the surface.
PML x-dir can be conceptually assumed up by a single transformation of the
original equation. Then, wherever and x derivative appears in the wave equa-
tions, it is replaced in the form
∂
x
∂
(8)
∂
x
∂
→
1
( )
x x
δ
ω
1
+
i
DOI: 10.4236/gep.2019.74001
4
Journal of Geoscience and Environment Protection
where
i = − in Equation (8).
1
Substituting the expression Equation (8) in Equation (6) the equations for
A. Muelas et al.
PML part (Benito et al., 2013), are:
2
U
δ
U
2
δ
+
+
U
ρ
i tt
,
i t
,
=
(
)
λ µ
+
U
i
+
U
µ
i kk
,
(9)
j
,
ji
2.3. Explicit Generalized Differences
The explicit scheme formulae for the second spatial derivatives with second or-
der approximation (see Equation (4)), are
∑
s
l
1
=
(10)
n
0
m u
i
n
l
m u
i
2
h k
,
i
+ Θ
= −
0
u
i
,
l
jk
0
jk
+
2
i
jk
t n t
= ∆
where hl = xl − x0 and kl = yl − y0 the superscript n denotes time step, the super-
scripts 0 and l refer to the central node and the rest of nodes of the star respec-
tively and s is the number of nodes in the star (in this work s = 8 and the star
nodes are selected by using the distance criteria) (Benito et al., 2001).
0
jkm are the coefficients that multiply the approximate values of the functions
0n
iu ) for the time n∆t in the generalized finite difference
Ui at the central node (
l
jkm are the coefficients that
explicit expressions for the space derivatives, and
multiply the approximate values of the functions Ui at the other nodes of the star
( n l
iu ). In all these expressions the cross-terms are equal.
The replacement in Equation (1) of the explicit expressions obtained for the
spatial derivatives (Equation (10)) and the central-difference formula for the
time derivative lead to the explicit differences scheme:
n
1
0
+
u
i
=
2
n
0
u
i
−
n
1
0
−
u
i
+
(
(
µ
+
−
n
0
m u
i
0
kk
+
2
)
t
∆
ρ
∑
s
l
1
=
+
)(
(
λ µ
)
n
m u
l
kk
l
j
−
n
0
m u
?
i
0
ji
+
∑
s
l
1
=
n
m u
l
ji
l
j
)
(
t
+ Θ ∆
2
)
,
2
h k
,
i
2
i
(11)
and the scheme in GFDM for the PML part can be obtained in a similar way:
)
)(
(
λ µ
n
0
m u
i
n
m u
1
0
+
u
i
∑
s
l
1
=
2
δ
0
u
i
t
Δ
(
)
=
−
+
+
−
+
2
l
ji
0
ji
(
)
(
l
j
n
n
2
2
)
t
∆
ρ
+
1
t
Δ
δ
(
µ
−
1
+
n
0
m u
i
0
kk
+
∑
s
l
1
=
n
m u
l
kk
l
j
)
(
1
− −
t
Δ
δ
)
n
1 0
−
u
i
+ Θ
(
t
∆
2
)
,
2
h k
,
i
2
i
(12)
where Θ represents the truncation error of the scheme.
In order to implement the free surface conditions, Equation (7), we need to
add the same number of nodes outside of our domain as there exists on the free
boundary and place them near it. Substituting the first order derivatives that ap-
pear in the free surface conditions Equation (7) by the explicit expressions (Equ-
ation (10)) the system of 2p equations is obtained:
p
n
m u
i
l
n
m u
i
n
m u
n
m u
−
+
+
−
+
p
k
p
k
p
j
l
k
l
k
l
j
∑
N
l
1
=
)
(13)
{
(
λ
(
+ −
s
l
1
=
∑
∑
N
l
1
=
)
) }
(
=
µ
⋅
n
j
n
m u
p
i
p
j
+
n
m u
l
i
l
j
( )
g t
i
DOI: 10.4236/gep.2019.74001
where the 2p unknowns are the displacements in the p added nodes. These un-
5
Journal of Geoscience and Environment Protection
A. Muelas et al.
knowns appear in the summation.
2.4. Stability of the GFD Scheme and Stability of PML Regions
The stability of the schemes has already been studied in (Ureña et al., 2011) and
its formulation leads to
t
∆ <
(
2
α β
+
2
4
+
)
) (
0
m
11
+
0
m
22
(
0
m
11
+
0
m
22
2
)
+
(
0
m
12
2
)
(14)
where α and β are the velocities of the P-waves and S-waves respectively.
The stability analysis in the continuum of the PML has already been studied in
(Salete et al., 2011).
The dispersion for the phase and group velocities have also been studied in
depth in (Gavete et al., 2017).
2.5. Heterogeneous Formulation
If the density and Lamé coefficients λ and µ are functions of spatial variables, in
agreement with Moczo (1998), in seismological problems the heterogeneous ap-
proach can be complicated and, then, so-called heterogeneous formulation is
preferred. Therefore, the same formulae are used for all points in the domain
and material discontinuities are accounted for spatial variation of the material
parameters (Benito et al., 2013).
3. Geophysical Investigation Modelling
The following is the application of the GFDM to model of a cross-hole and
model of seismic refraction studies of the aforementioned seismic tests.
3.1. Model of a Cross-Hole Test
A geotechnical investigation has been carried out to elaborate a geotechnical
profile on a construction site. Three boreholes with core recovery and an addi-
tional cross-hole test were performed in order to obtain a seismic velocity vs.
depth profile. Based on the geotechnical investigation aforementioned, the sub-
soil comprises an upper layer of medium dense to dense sand up to 6-meter
depth, overlying a very dense silty sand layer. The λ parameter, the dynamic
Shear Modulus and the compression and shear wave velocity vs. depth are
shown in Figure 1.
The dynamic parameters shown in Figure 1 have been assigned to the subsoil
as a function of soil depth.
The cross-hole test has been analyzed by in a GFDM model, applying a vertic-
al pulse (to generate S-waves) and a horizontal pulse (to generate P-waves) at
different depths. The step pulse is applied for a time increment of t = 1e−5 sec.
The model is 16 m (vertical) × 6 m (horizontal), with a regular node spacing of
0.05 m. The number of nodes is 38,400. The star is made up of eight nodes plus
the central node. The selection of nodes associated to each central node is based
6
Journal of Geoscience and Environment Protection
DOI: 10.4236/gep.2019.74001
A. Muelas et al.
Compression and shear wave velocity vs. Depth
100
300
600
200
400
500
700
800
900 1000
vp (m/s)
vs (m/s)
λ (MPa) - Gdyn (MPa) vs. Depth
0
100
200
300
400
500
600
700
800
900 1000
0
)
m
(
h
t
p
e
D
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
λ (MPa)
Gdyn (MPa)
)
m
(
h
t
p
e
D
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Figure 1. Dynamic soil properties (Borehole CH03).
on the distance criterion. The weighting function adopted is the potential func-
tion 1/(dist)3. A free surface boundary condition has been imposed to the upper
and right limits of the model. An absorbing boundary condition (Perfectly
Matched Layer) has been imposed to the lower and left limits of the model.
Figure 2 shows the moment at which the front S-wave reaches the geophone
located at 9-meter depth in the GFDM model. The shear, head or Von-Schmidt
and surface or Rayleigh waves can clearly be seen in Figure 2.
The vertical displacement over time for the geophone located at 9-meter depth
is shown in Figure 3. The shear wave velocity is obtained from the time of wave
arrival, taking into account the distance between the source and the geophone (3
meters).
By applying a horizontal pulse at the source location, the P-waves propagation
has also been analyzed. The pressure or dilatational, shear, head and surface
waves can be seen in Figure 4.
The curve that appears in Figure 5 shows the horizontal displacement ux over
time for the geophone located at 9-meter depth. The P-waves velocity is ob-
tained from the time of wave arrival, similarly to the S-waves case.
3.2. Model of Seismic Refraction
The seismic refraction method is usually performed in quarries or open-pit
mines to make an estimate of the rock volume that may be extracted from the
site and the most appropriate excavation method. For example, the seismic re-
fraction method carried out in a quarry in Valparaiso (Chile) to quantify the
weathered and fresh rock volume. The layout of the seismic refraction profiles is
shown in Figure 6. As a result of the geophysical study, the longitudinal profile
shown in (Figure 7) was obtained.
Basically, the subsoil is composed of three layers (sandy oil, weathered rock
and fresh rock).
The purpose of the model is to simulate the real geophysical investigation,
measuring the travel time for the wave to reach the geophones located on the
7
Journal of Geoscience and Environment Protection
DOI: 10.4236/gep.2019.74001
A. Muelas et al.
Figure 2. Cross-hole test model-front shear waves reaching the 9-meter depth geophone.
Figure 3. Vertical displacement vs. time for geophone located at 9 m depth.
ground surface along the profile. The real data obtained from the geophysical
investigation are shown in Figure 8. The different dromochrones (distance vs.
travel time curves) have been drawn for the different pulse locations. Three
slopes may be distinguished in most of these curves, corresponding to the three
soil layers aforementioned: sandy soil, weathered rock and fresh rock. Three
representative dromochrones have been selected (Figure 9) to enhance these
three slopes.
The slope of the first section of the three curves is very similar, leading to a
seismic velocity of 350 - 500 m/s for the sandy soil layer. Nevertheless, the point
where the curve changes its slope depends on the thickness of the soil layer. On
curve 2, the time when the refracted wave reaches the direct wave is shorter than
the one for curves 1 and 3. This difference leads to a smaller thickness of layer 1
in the neighborhood of x = 30 m. Regarding the second section of the curves, the
8
Journal of Geoscience and Environment Protection
DOI: 10.4236/gep.2019.74001
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
