NUMERICAL SOLUTION OF A BOUNDARY VALUE PROBLEM WITH EFFECTIVE BOUNDARY CONDITIONS FOR CALCULATION OF GRAVITY Текст научной статьи по специальности «Математика»

Математические заметки СВФУ Январь—март, 2021. Том 28, № 1

UDC 519.63+531.26

NUMERICAL SOLUTION OF A BOUNDARY VALUE PROBLEM WITH EFFECTIVE BOUNDARY CONDITIONS FOR CALCULATION OF GRAVITY D. Kh. Ivanov and P. N. Vabishchevich

Abstract. Forward modeling of gravity field on the base of boundary-value problem solution is a promising technique against traditional summation methods. To calculate gravity of a body with known physical and geometrical properties, one can firstly solve a boundary-value problem for gravitational potential and then calculate its gradient. This approach is more common, but it requires two operations. Another approach is to solve a boundary value problem formulated directly for gravity itself. The main difficulties for methods based on boundary-value problem solution are proper boundary condition, domain size and domain discretization. For the first approach there are plenty of works dealing with these cases in contrast with the second approach. In this paper, authors discuss and compare boundary conditions of two types: the Dirichlet and Robin, in terms of approximation accuracy for the second approach using the finite element method. Calculation results are presented for a test problem, when the gravitational field is produced by a homogeneous body in the shape of a right rectangular prism. A more effective boundary condition is a Robin type condition derived from a simple asymptotic approximation of the gravitational field by the equivalent point mass.

DOI: 10.25587/SVFU.2021.74.56.008

Keywords: gravitational field, Poisson equation, approximate boundary condition, finite element method.

Introduction

The study of gravitational potential, its analysis and interpolation is an integral part of applied problem of gravity exploration [1]. The gravity exploration is to determine the distribution of heterogeneous masses in the upper layers of the Earth's crust from observed gravitational anomalies on or near the Earth's surface. The forward problem of gravimetry is the calculation of gravitational potential and its derivatives produced by a body with known geometric properties and a given density contrast [2].

The gravitational potential of a three dimensional body inside a bounded domain D c R3 with a given density function p(x, y, z) can be calculated using the

This work was supported by the Mega Grant of the Russian Federation Government (No. 14. Y26.31. 0013), the grant of Russian Foundation for Basic Research (No. 19-31-50044 and No. 2001-00207) and the Grant of the President of the Russian Federation for State Support of Young Scientists (No. MK-1131.2020.1).

© 2021 D. Kh. Ivanov and P. N. Vabishchevich

following integral expression

U(P)=7J^^dv, P = (x, y, z) £ R3, (1)

D

where y = 6.674 ■ 10-11 m3 kg-1 s-2 is the gravitational constant, Q = (x',y',z') denotes an integration point, dv = dx'dy'dz' and

r(Q, P) = ((x - x')2 + (y - y')2 + (z - z')2)1/2 .

The gradient of the potential expressed as follows

9:=vu = 1jp{Q)v7(chr)dv (2)

D

defines gravitational field, alternative name is just gravity. We note that g = (gx,gy,gz), where gx and gy constitute the horizontal components of the gravitational field, and gz is the vertical component of the gravitational field. The gravitational potential satisfies Poisson's equation

AU = -4n7p(P), P e R3, (3)

where zero density contrast is implied outside D. An additional condition for the equation (3) is usually as follows

lim U(P) = 0, Q e D, (4)

which originally comes from physics and provides a bounded solution. Applying differentiation procedure on the solutions, one can derive desired component of the gravity.

Over the past two decades more and more attention is paid to development of effective numerical methods for solving forward problems of gravimetry. The existing methods can be divided into two main groups: the first group includes methods based on the definition of a gravity quantity according to (1), (2), the other includes methods based on solution of the boundary value problem for partial differential equation (3), (4).

The first group is traditional and uses closed-form, analytical solutions. Analytical representation of the volume integrals (1), (2) for bodies in simple geometrical forms with specific density functions has a vast history. The earliest work in this direction is [3], where a numerical algorithm was given for exact computation of the gravitational field produced by a homogeneous polygonal body. Works [4-7] consider more complex case with polyhedral bodies with depth-dependent density contrasts. Thorough list of literature can be found therein. A more universal tool for approximate numerical integration is based on the use of cubature formulae [8-11]. The approach based on the use of cubature formulae is easier to numerically implement since only simple algebraic operations, such as addition and multiplication, are used avoiding singularity of the integrand. For bodies with complex geometry and

arbitrary density contrast approximate integration technique is based on dividing the domain into small elements in a simple geometric shape with a given density function, for which there is the analytical solution.

Calculation of the gravitational data produced by arbitrary three dimensional bodies requires large computing power and resources. A great attention is paid to numerical algorithms for calculation of gravitational data from the solution of boundary value problems for the Poisson equation (3), (4). Prerequisite for that is a rapid development of numerical methods for approximate solution of multidimensional boundary value problems for partial differential equations and intense development of computational technologies. The common problem of this kind of approach refers to finiteness of a computational domain and approximation of the boundary condition. Computational costs are directily connected with the size of a computational domain and a mesh. Gravitational potential is computed on the entire computational mesh at the same time when extracting required data at observation points. For a large number of observation points the use of such computational procedures provides a more efficient calculation of gravitational anomalies than the forward integration method, as it is shown in [10,12].

Here we give a brief overview of different techniques based on solving boundary value problems in gravity forward modelling. The most popular methods are finite difference method [13,14], finite volume [14-17] and finite element methods [10,12, 17-19]. Also we need to mention works [20,21] with application of spectral element method. In [22] a forward modelling of gravitational data is considered using modified finite difference method with radial basis functions in an arbitrary cloud of nodes. In traditional methods, computational mesh is regular usually with elements in form of the right rectangular prism. In finite volume and finite element methods it is possible to use unstructured meshes with elements as tetrahedrons or hexahedrons for more accurate modelling geometrical properties of the material body.

In gravitational exploration, the basic observed gravitational data is the vertical component of gravitational field, gz (vertical gravity). To compute gz one can firstly solve the boundary value problem for the gravitational potential, then the value of gz at observation points can be obtained on the basis of numerical differentiation procedures (see for example [17]). An alternative approach is to state the problem directly for the vertical gravity itself. The proper equation is derived from taking derivative of the gravitational potential equation (3) along vertical axis. It was used in works [14,15]. Need to mention that for such approach challenges occur with respect to derivatives of discontinuous density function. The possible treatment of that issue is a weak formulation of the boundary value problem for the Poisson's equation as it was implemented for the finite difference method in [14] and for the finite volume method in [14,15].

To solve boundary value problems of gravimetry in a finite domain and to get more accurate results, we need set a boundary condition that approximates the behavior of the gravitational field at far distance. According to asymptotic attenuations of the gravitational potential and its gradient as r-1 and r-2, respec-

tively, a homogeneous Dirichlet boundary condition can be applied for both terms. In [17] zero values are set for the gravitational potential on the boundary, in [14,19] for the vertical gravity. Another effective condition is Robin type condition, it was applied in [12] for the potential and compared in [10]. For two dimensional case these conditions were compared in [23]. Need to mention that the use of the Robin condition presented in [12] for two dimensional case is incorrect since the gravitational potential should be considered as logarithmic function which does not decays in radial direction. In [14] a Robin condition was considered for the vertical gravity and compared with homogeneous Dirichlet condition. These conditions was claimed to give same results for the aspect ratio of the domain size to the body size equal to 10. According to a far-field approximation of the vertical gravity, in [15] a more accurate Robin boundary condition was realized separately on each facet of a prismatic computation domain. In that work the aspect ration was set to 6. These two works lack investigation of the domain size effect on accuracy.

It is worth mentioning that a different approach for taking into account the boundary condition at infinity was implemented in works [20,21]. The idea is the use of spectral element method adding infinite elements in an outer single layer around the mesh. Such approach applied to the potential was comparised with pure finite element method in [23]. Consideration of this method is out of the scope of the present paper.

In this paper, modelling of the vertical gravity calculation produced by a material body is considered in a sense of numerical solution of the boundary value problem stated directly for the vertical gravity in a bounded domain. Instead of a strong statement of the problem, a weak formulation is considered to alleviate the vertical derivative of the density function for the use of finite element method. Boundary conditions are considered in two types: Dirichlet and Robin, which are homogeneous Dirichlet condition, asymptotic Dirichlet condition, Robin condition with a constant variable and asymptotic Robin condition. Here, the asymptotic conditions are obtained from simple asymptotic approximation of the gravitational potential by a point source with equivalent mass located at the barycenter of the body. An analysis is focused on investigation of the effect of the domain size on the approximation accuracy for each boundary condition. For numerical implementation, standard piecewise linear finite elements are applied. Computational domain is discretized by tetrahedron elements. After analyzing the boundary conditions on uniform meshes, we test the use of a non-uniform regular mesh. For that we select an inner uniform mesh containing the body, outer of which a non-uniform mesh is constructed such that the element size increases from the body in geometric progression.

1. Boundary value problem and approximate boundary conditions

Let O = [—L, L]3 c R3 denote a large enough bounding domain, dO is its boundary. Let D c O denote a volume of a material body with a given density

distribution p(x, y, z). The vertical axis z is pointed down, the origin coincides with the center of see Fig. 1. Outside D the density function p is set to zero.

Fig. 1. Computational domain Q, material body D and observation surface T.

By taking derivative of equation (2) with respect to z, we obtain the following Poisson's equation for the vertical gravity gz

dp(P)

A gz = -47T7-

dz

P e fi,

(5)

which should be supplemented by some approximate boundary condition. Assume that we need to determine values of gravity in the observation manifold r, a horizontal surface as shown in Fig. 1 .

In faraway distance gravitational potential of any compact body can be asymptotically approximated by an effect of the gravitating point source with equivalent mass located at the barycenter. Hence the gravitational potential at far distance is approximated as follows

U(P)*> 7-,

r

where m = J pdv is the total mass of the body, D

r = ((x - xo)2 + (y - yo)2 + (z - zo)2)1/2

is a distance between the barycenter (x0,y0, z0) and the observation point (x,y,z). From that we can derive approximation for the vertical gravity:

gz (P)

(p ),

Qznp(p) _7m

z - zo

Thus, at large distances it vanishes quickly as 1/r2, therefore one can simply apply the homogeneous Dirichlet condition for equation (5). If the total mass m and the barycenter (x0 ,y0 ,z0) are given, the verical gravity on the boundary dfi can be approximated by the function gZmp (P), named as the asymptotic Dirichlet boundary condition. The first type of Dirichlet condition is more rough and its accuracy strongly depends on remoteness of the domain boundary. Obviously, in contrast to that, the second type gives more accurate results, but it depends on the mass value and the barycenter coordinates.

3

For a function 1/r , k > 1, the following expression takes place

d ( 1 \ k 1

For potential (k = 1) we have As in [12] one can apply

Or r

^- + -U = 0 on <90, On r

where n is the unit outward normal vector at the boundary. Although the transition from radial partial derivative to normal partial derivetive on the boundary of an arbitrary prismatic domain should be done more accurately. In this case, taking derivative of U = Ym/r over any direction n we have

On

which becomes

dU ^r ■ n jj _ q dn r2

Here r denotes the vector-radius from (x0,y0,z0) to (x,y,z). Thus, such boundary condition is more recommended for the boundary-value problem for the potential in a prismatic domain. We do not give an investigation in this area since we are focusing on calculation of the vertical gravity by the equation (5). If k = 2 for the gravity amplitude we have

M + Hto| = „.

Or r

This leads to consideration of a simple Robin boundary condition for the vertical gravity:

dgz , „

--1- agz = 0,

On

where the parameter a = O(1/r). In [14] a =1, but it should depends on the domain size L. According to the aforementioned procedure for the potential, a higher accurate condition can be provided after taking a proper derivative of function gZmp. In this case the parameter a(x,y, z) depends on spatial variables. As given in [15], this parameter has the next expression

3(x—xq) _ _j_

*{x,y,z)= { n = ±ey,

n = ±ez,

r2 z — zo ' z'

where ev is the unit vector of the axis v = x,y, z. Introducing a vector

/ 3(x - xo) 3(y - yo) 3(z - zo) 1 \

a(x,y,z)=[---,---,-----, (6)

r2 r2 r2 z - zo

the asymptotic Robin boundary condition has the following expression

dgz , „

—- + a • n gz = 0. on

Features of the stated boundary value problem for the equation (5) using the formulated boundary condition are associated with the sign of the coefficient a ■ n at the boundary of the prismatic domain fi. Standard restrictions [24] on non-negativity of this coefficient at z = ±L may not be fulfilled despite natural assumption that z0 e fi.

Combining above mentioned cases, we highlight next boundary conditions:

1. Exact Dirichlet condition with analytical solution gf to test a method in the prismatic domain:

gz(P)= gf (P), P e dfi;

2. Homogeneous Dirichlet condition:

gz(P) = 0, P e dfi;

3. Asymptotic Dirichlet condition:

gz (P ) = gZmp(P ), P e dfi;

4. Robin condition with a constant parameter a = const = O(1/L):

^ + agz = 0, P G Oïl;

5. Asymptotic Robin condition with the variable parameter a(x,y,z):

^ +a-n gz = 0, P€dil. on

The exact boundary condition (condition 1) will help to compare the other boundary conditions 2-5. In further, we will refer the condition 4 as the constant Robin condition.

2. Numerical implementation

Preliminaries for finite element method [25] start with introduction of Hilbert space Vh C L2(H) of first order finite elements. For any two elements from that space scalar product and norm are defined in the standard way:

(u, v) = J uv dv, ||u|| = (u, u)1/2 Vu, v e Vh. n

Basis functions are piecewise linear functions defined on cells of a mesh built over the computational domain H. For simplicity a nonuniform regular mesh can be constructed with tetrahedron shaped cells. Heterogeneity of the mesh consists from refinement of element sizes near the material body and coarsening element sizes near the computational boundary. In result, we can count on a higher accuracy in calculation of the gravitational data.

Instead of the strong statement of the boundary value problem with chosen boundary condition we consider a weak formulation for equation (5). Firstly, integration of the right hand side of (5) with a variable function v G Vh gives

fzvdv = j ±{pv)dv-jP^-dv.

dp \ f dp id — ,v J = I —Dr' ~~

n n n

Applying the divergence theorem to the first component in the right-hand side, it turns to the integral over the boundary dO:

d

J ~Q~(fyv)dv = J pv cos(nz)ds = 0,

dv

n an

which vanishes due to zero density at the boundary.

For the Dirichlet boundary conditions the variational formulation is as follows: find gz € Vh such that

dv\ w _

as

('Vgz ,Vv) = -4^7 - Vv€Vc where V0h = {v G Vh : v = 0 on Oiî}. By defining new notation

(u, v)an = J uvds, an

the first boundary conditions at the boundary nodes (for P G dO) give

(gz, v)dn = (gE,v)dn, (gz,v)dn = 0, (gz,v)dn = (grmp,v)an Vv g Vh, respectively.

Taking into account the Robin boundary conditions (conditions 4 and 5) the variational formulation is slightly modified:

( dv\

(Vgz, Vv) +a(gz,v)dn = -47r7 Vv G Vh,

for boundary condition 4 and

( dv\

(Vgz,Vv) + (ct-n gz,v)dn = -47r7 ( J Vv G Vh,

for boundary condition 5.

Let Th denote the mesh and Pi, i = 1,... ,N, be nodes of the mesh, and the mesh elements represent tetrahedrons. Let dTh be a part of the mesh nodes corresponding to the set of boundary nodes.

In our numerical examples we use uniform and nonuniform regular meshes. According to the rapid attenuation of the gravitational potential and its derivatives the nonuniform regular mesh can be constructed such that element sizes get larger from the barycenter. We select an area D', D c D' C O, where more refine uniform mesh is constructed. Outside this area a nonuniform regular mesh is built such

that the size of the elements increases in geometrical progression along each spatial direction.

In the space Vh with piecewise linear basis functions i = 1,... , N, an approximate solution u has the following form:

N

u(P )=£ ulVl (P), (7)

i= 1

where ui = u(Pi). The basis functions satisfy next condition

v(Pi =j. 10, i = j.

Variational formulation for equation (5) with Dirichlet boundary conditions leads to solution of the system of linear algebraic equations

Au = b,

A = (aij), aij = (V^i, Vvj), j = 1,...,N,

u=(ui), b= -4tt7(p, ^P-), i = l,...,N.

dz

For the Robin boundary conditions, a(vi)dTh or (a ■ n )dTh are added to the stiffness matrix A.

Governing stiffness matrix A is symmetric, positively defined and sparse. Since direct solution of a big system of linear algebraic equations consumes large computational resources, we use iterative method such as conjugate gradient method (CG) with incomplete Cholesky factorization (ICC) as the preconditioner [26]. The vertical gravity at the observation points is calculated by interpolation using mesh nodes values of u by (7).

In this paper, the software for numerical solution of stated variational problems is implemented on the base of the open source computational platform FEniCS [27, https://fenicsproject.org/]. Computational meshes are constructed by gmsh [28, http://gmsh.info/].

3. Numerical experiments

We consider the material body inside the right rectangular prism of dimensions 1 x 1 x 0.5 (km) with a homogeneous density distribution p = 2000kg m-3. The origin of the coordinate system coincides with the barycenter of the body and the prism's center, the coordinate axes coincide with the main axes of the prism coincide, see Fig. 1. Observation data are collected in the square surface T of 2 km side at height H above the origin either in the whole horizontal section S.

First, computational tetrahedral mesh Th is uniform regular with the element size 1/12 km along each space direction. In this case, the maximum length between two nodes of the mesh is h = y/Z/12. For the observation surface an uniform mesh is constructed also with the element size 1/12 km. Depending on the height H, its nodes can coincide with the nodes of the main mesh.

For the homogeneous rectangular prism the analytical expression of the vertical gravity at the point P is given as [29,30]

2 2 2 , \ gf(P) =7 Vijk ( zk tan-1 X%yj -xt In (yj + rijk)-yj ln(xi+rijk)), (8)

i=i j=i k=i V ZkTljk ' where xi, yj, zk are coordinates of the prism ends in a reference system corresponding to the observation point, herein

rijk = (x2 + yj2 + z2)1/2, j = (-1)i(-1)j (-1)k.

The analytical solution on the observation surface at the height H =1 km is illustrated in the Fig. 2 (left), 1 mGal = 0.00001 ms-2.

Fig. 2. Analytical vertical gravity at the height H = 1km (left), L2 errors in the observation section for different parameters a in the boundary condition 4 (right).

3.1. Optimal value for the constant parameter in the Robin boundary condition. Here, we numerically analyze an effect of the constant parameter a in the boundary condition 4 on the approximation accuracy. Due to estimation a ~ 2/r, values for this parameter are considered in [0.5/L, 2.5/L]. For every value a, we calculate L2-norm error of an approximate solution in the horizontal section S of the computational domain ÎÎ at H =1 km:

1/2

£2 = (/(ffz - £)2 ds))

S

For the observation surface T corresponding error is denoted by e2(r).

For series of half-length of £2: L = 2,3,4 and 5 km, the dependence of the error e2 on the parameter a is given in the Fig. 2 (right). It is seen that for each L there is an optimal value a located in [1.5/L, 1.7/L]. It is reasonable to compare this result with the average value of the variable coefficient in the asymptotic Robin boundary condition 5, according to (6), over the boundary of This average is calculated via taking an integral over facets of the prism £:

i i

1 ra.nds = ±. /7 1.585

24L2 J 3L \J J 1 + x2 + y2 I L

on

The optimal value for the constant a is close to a. This boundary condition will be further used with a = 1.6/L.

3.2. Effect of the boundary conditions and the computational domain size. For every parameter L = 2,3,4, 5 km and boundary conditions 1-5 the boundary value problem is solved. Table 1 lists a runtime corresponding to solution of the system of linear algebraic equations with iterative method CG with preconditioner ICC, a number of iterations needed for the convergence with a stopping criteria as relative error 10-6. Besides that there are L2-norm error e2 and the maximum absolute error over the section S at z = 1 km, and also the error e2(r) in the selected surface r. Fig. 3 shows dependences of errors e2 and on the domain half-size L.

Absolute error of computed values of the vertical gravity Sa = gz — gf and relative error S = |Sa|/gf at discrete nodes on the surface at H =1 km are illustrated in Fig. 4-7. As expected, the approximation accuracy is getting better with increase in the computational domain size for each boundary condition. Moreover, the result of the exact boundary condition 1 determines the error of finite element approximation of the stated problem for equation (5), and it is seen to be stable.

The worst result belongs to the boundary condition 2 (the homogeneous Dirich-let condition). For this condition, when L = 5km (the most distant boundary) the L2-norm error in the observation surface r is 3.3 times bigger than with the other conditions. The next accurate condition is the boundary condition 4 (the constant Robin condition). From Table 1 and Fig. 3 it is seen that the accuracy of this condition decreases and almost matches with the exact boundary condition 1 for the largest distant boundary.

The boundary condition 3 (the asymptotic Dirichlet condition) and 5 (the asymptotic Robin condition) show approximately same results. Even for short distant boundary the approximation accuracy of both conditions almost matches with the result of the exact boundary condition 1. Need to mention, from the Table 1 it is seen that the asymptotic Dirichlet condition gives less L2-norm error on the observation area r eighter maximum norm error on the section S than the others even with the exact boundary condition. Although this condition differs from condition 5 in approximation near the boundary, the latter has a better result, see relative errors, S, in Fig. 4-7.

To analyze the effect of mass value in condition 3, Fig. 8 presents errors e2, e2(r) and for different mass values m. It is seen that L2-norm error on the whole section S does not exceed the accuracy of the exact boundary condition. Nevertheless, there is some area depending on the domain size, where maximum error and L2-norm error on the observation area r can be less than with the exact condition. For small domain sizes these conditions are more sensible to the mass value, smoothening out for large sizes.

Hence, when the mass m of the material body is unknown, for given barycen-ter (x0,y0,z0) the asymptotic Robin boundary condition with variable parameter a(x, y, z) defined as (6) is the most preferable condition to approximate the vertical gravity. If the barycenter is also unknown, one can apply the Robin boundary con-

Table 1. Characteristic results of vertical gravity calculation: BC — boundary condition, t — computational time in seconds, K — iterations

mesh BC t K t2 e2(r) £-00

1 0.19 33 0.030 0.028 0.024

2 0.19 33 1.967 1.106 0.618

L = 2 ,N = 117649 3 0.27 33 0.036 0.019 0.021

4 0.26 44 0.104 0.070 0.064

5 0.25 44 0.034 0.028 0.024

1 0.91 46 0.030 0.028 0.024

2 0.91 46 0.918 0.348 0.216

L = S,N = 389017 3 0.92 46 0.030 0.027 0.023

4 1.17 60 0.060 0.038 0.033

5 1.20 61 0.031 0.028 0.024

1 2.63 56 0.030 0.028 0.024

2 2.64 56 0.526 0.160 0.097

L = 4,^ = 912673 3 2.64 56 0.030 0.028 0.023

4 3.63 78 0.044 0.032 0.026

5 3.66 79 0.030 0.028 0.024

1 6.20 68 0.030 0.028 0.024

2 6.21 68 0.341 0.093 0.060

L = 5 ,N = 1771561 3 6.21 68 0.030 0.028 0.024

4 8.21 91 0.037 0.030 0.025

5 8.68 96 0.030 0.028 0.024

L, km

Fig. 3. L2 error (left) and error (right) of approximations over the horizontal section at H = 1 km for different L in each boundary conditions 1—5.

dition with the constant parameter a = 1.6/L, if the boundary of the computational domain is sufficiently far from the center at least five times.

-0.45 -0.90 -1.35 -1.80 -2.25

il

i-2 mGal 3.3

k. 4

II

-4.2 -5.1 -6.0 -6.9 7.8 -8.7 9.6

0

x, km

-1.0 -0.5 0.0 0.5 1.0 x, km

0

x, km

-1.0 -0.5 0.0 0.5 1.0 x, km

-0.45 -0.90 -1.35 -1.80 -2.25

024 x, km

0

x, km

1.0 -0.5 0.0 0.5 1.0 x, km

-1.0 -1.5 -2.0 -2.5

024 x, km

-1.0 -0.5 0.0 0.5 1.0 x, km

-0.45 -0.90

-1.80 -2.25

-2 0 2 4 x, km

0

x, km

4

-2

4

-2

4

x, km

-1

-2

4

x, km

-2

-3

-2 mGal

4

4

-2

4

4

4

4

-2

-2

-2

-2

4

x

-3

-3

-2 mGal

4

4

4

x, km

-4 -5.0

-2.5 0.0 2.5 x, km

x 10-2 mGal 0.45 00 -0.45 -0.90 -1.35 -1.80 -2.25

0-2 mGal i m-1.8

-2.4 -3.0 -3.6 -4.2 -4.8 -5.4 -6.0

-2.5 0.0 2.5 5.0 x, km

4 2

m,k 0 -2 -4

4 2

m,k 0 -2 -4

-0.45 -0.90 -1.35 -1.80 -2.25

-1.2 -1.6 -2.0 -2.4

-0.45 -0.90 -1.35 -1.80 -2.25

-2.5 0.0 x, km

-2.5 0.0 2.5 x, km

-1.0 -0.5 0.0 0.5 1.0 x, km

-2

x, km

-2 mGal

4

-4

4

x, km

-2 mGal

m/m m/m

L, km

\

—\—

\ / y

v \ / y

\ \ / >

/

........^

fet

VCX- ___^

0.8 1.0 ^ 1.2 1.4

Fig. 8. £2 (left), £2(r) (center) and (right) for different mass values m in the

boundary conditions 3.

3.3. The use of nonuniform regular meshes. By the aforementioned facts, we can obtain a good approximation of the vertical gravity using the asymptotic Robin boundary condition. Accuracy improvement can be achieved by refining the mesh. However to avoid huge computational costs it is better to perform an adaptive strategy of mesh construction applying nonuniform regular or unstructured meshes.

We consider the computational domain with L = 4 km. Nonuniform regular mesh is constructed as follows. Each edge of £ is divided into three intervals: [-4, -1], [-1,1] and [1,4]. The middle part is uniformly discretized by step 1/24km (two times finer than before). Each outer part is divided into M subintervals with increasing sizes in geometric progression from the center. The size of the first subin-terval near the center is set as q/24, where the progression factor q > 1 is found from the relation

g qM ~ 1 24 q - 1

We consider three cases: 1) q = 1.133, M = 18; 2) q = 1.08, M = 24; 3) q = 1.052, M = 30. Mesh elements for one side of £ are presented in Fig. 9. Total number of nodes is N = 614125, 912673 and 1295029, respectively. In the second mesh, the number of nodes is the same as in the case with the uniform mesh.

Fig. 9. Nonuniform regular mesh: left — q = 1.133, M = 18; middle — q = 1.08, M = 24; right — q = 1.052, M = 30.

Table 2 and Fig. 10 show computation results by the asymptotic Robin boundary condition. There is a noticable increase in the approximation accuracy for all

Table 2. Computation on nonuniform mesh: t — solution time, K — iterations

q M N t, s K t2 e2(r) £oo

1.133 18 614125 3.47 104 0.028 0.018 0.015

1.08 24 912673 5.53 106 0.014 0.012 0.010

1.052 30 1295029 7.46 104 0.011 0.010 0.008

Fig. 10. Nonuniform mesh for L = 4 km: upper row — q = 1.133, M = 18; middle row — q = 1.08, M = 24; lower row — q = 1.052, M = 30.

cases.

4. Conclusion

Calculation of gravity data on the basis of solution the Poisson's equation is a promising alternative to numerical integration techniques. The transition to a finite computational domain requires a boundary condition. The main attention is paid to the expression of that condition and the size of the computational domain to reduce

its effects on the accuracy.

In this paper, two classes of approximate boundary conditions (Dirichlet and Robin) where distinguished to solve the boundary value problem stated directly to the vertical gravity using the finite element method. Clarification of the asymptotic Dirichlet and Robin boundary conditions is based on approximation by the field of the equivalent point mass located at the barycenter. In addition to that conditions, the homogeneous Dirichlet boundary condition and Robin boundary conditions with the constant variable were considered. The use of the finite element method is motivated to handle the problem associated with the derivative of the density function by considering thr weak formulation of the boundary value problem.

The results of numerical calculations for the prismatic body showed that the asymptotic Dirichlet boundary condition gives a good accuracy, depending on the observation area and the error type, even for approximate mass values. Sensibility of that parameter decreases with increase in the domain size. The asymptotic Robin boundary condition with the variable coefficient gives almost the same accuracy as the exact boundary condition. This condition seems to be more preferable since it is less parameterized.

The results on the nonuniform regular meshes with elements, whose sizes grow in geometric progression from the body, have demonstrated the possibility of improvement the approximate accuracy without significantly increasing computational costs.

REFERENCES

1. Barnett C. T., "Theoretical modeling of the magnetic and gravitational fields of an arbitrarily shaped three-dimensional body," Geophys., 41, No. 6, 1353—1364 (1976).

2. Blakely R. J., Potential Theory in Gravity and Magnetic Applications, Camb. Univ. Press, Cambridge (1996).

3. Brenner S. C. and Scott L. R., The Mathematical Theory of Finite Element Methods, Springer, New York (2008).

4. Butler S. L. and Sinha G., "Forward modeling of applied geophysics methods using Comsol and comparison with analytical and laboratory analog models," Comput. Geosci., 42, 168—176 (2012).

5. Cai Y. and Wang C.-Y., "Fast finite-element calculation of gravity anomaly in complex geological regions," Geophys. J. Int., 162, No. 3, 696-708 (2005).

6. Casenave F., Metivier L., Pajot-Metivier G., and Panet I., "Fast computation of general forward gravitation problems," J. Geodesy, 90, No. 7, 655-675, (2016).

7. Chakravarthi V., Ramamma B., and Reddy T. V., "Gravity anomaly modeling of sedimentary basins by means of multiple structures and exponential density contrast-depth variations: A space domain approach," J. Geol. Soc. India, 82, No. 5, 561-569 (2013).

8. Cordell L., "Gravity analysis using an exponential density-depth function; San Jacinto Graben, California," Geophys., 38, No. 4, 684-690 (1973).

9. D'Urso M. G., "Analytical computation of gravity effects for polyhedral bodies," J. Geodesy, 88, No. 1, 13-29 (2014).

10. D'Urso M. G. and Trotta S., "Gravity anomaly of polyhedral bodies having a polynomial density contrast," Surv. Geophys., 38, No. 4, 781-832 (2017).

11. D'Urso M. G., "Gravity effects of polyhedral bodies with linearly varying density," Celest. Mech. Dyn. Astronomy, 120, No. 4, 349-372 (2014).

12. Evans L. C., Partial Differential Equations, 2nd ed., Amer. Math. Soc. (2010) (Grad. Stud. Math.; vol. 19).

13. Farquharson C. G. and Mosher C. R. W., "Three-dimensional modelling of gravity data using finite differences," J. Appl. Geophys., 68, No. 3, 417-422 (2009).

14. Garcia-Abdeslem J., "The gravitational attraction of a right rectangular prism with density varying with depth following a cubic polynomial," Geophys., 70, No. 6, J39—J42 (2005).

15. Geuzaine C. and Remacle J.-F., "Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities," Int. J. Numer. Methods Eng., 79, No. 11, 1309—1331 (2009).

16. Gharti H. N., Tromp J., and Zampini S., "Spectral-infinite-element simulations of gravity anomalies," Geophys. J. Int., 215, No. 2, 1098-1117 (2018).

17. Gupta H., Encyclopedia of Solid Earth Geophysics, Springer, Dordrecht (2011).

18. Guzman S., Forward Modeling and Inversion of Potential Field Data Using Partial Differential Equations, MSc Thes., Colorado School of Mines (2015).

19. Haber E., Holtham E., and Davis K., "Large-scale inversion of gravity gradiometry with differential equation," in: SEG Technical Program Expanded Abstracts 2014, pp. 1302-1307, Soc. Explor. Geophys. (2014).

20. Haji T. K., Faramarzi A., Metje N., Chapman D., and Rahimzadeh F., "Challenges associated with finite element methods for forward modelling of unbounded gravity fields," UK Assoc. Comput. Mech. Conf. Pap., 16, 1-4 (2019).

21. Haji T. K., Faramarzi A., Metje N., Chapman D., and Rahimzadeh F., "Development of an infinite element boundary to model gravity for subsurface civil engineering applications," Int. J. Numer. Anal. Methods Geomech., 44, No. 3, 418-431 (2020).

22. Hoschl V. and Burda M., "Computation of gravity anomalies caused by three-dimensional bodies of arbitrary shape and arbitrary varying density," Stud. Geophys. Geodaet., 25, No. 4, 315-320 (1981).

23. Howell L. E., Forward Modeling the Gravitational Field Using a Direct Solution of Poisson's Equation, MSc Thes., Colorado School of Mines (2013).

24. Jahandari H. and Farquharson C. G., "Forward modeling of gravity data using finite-volume and finite-element methods on unstructured grids," Geophys., 78, No. 3, G69-G80 (2013).

25. Jiang L., Liu J., Zhang J., and Feng Z., "Analytic expressions for the gravity gradient tensor of 3D prisms with depth-dependent density," Surv. Geophys., 39, No. 3, 337-363 (2018).

26. Kwok Y.-K., "Gravity gradient tensors due to a polyhedron with polygonal facets," Geophys. Prospect., 39, No. 3, 435-443 (1991).

27. Logg A., Mardal K. A., and Wells G., Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Springer-Verl., Berlin; Heidelberg (2012) (Lect. Notes Comput. Sci. Eng.; vol. 84).

28. Long J. and Farquharson C. G., "Three-dimensional forward modelling of gravity data using mesh-free methods with radial basis functions and unstructured nodes," Geophys. J. Int., 217, No. 3, 1577-1601 (2019).

29. Maag E., Capriotti J., and Li Y., "3D gravity inversion using the finite element method," in: SEG Technical Program Expanded Abstracts 2017, pp. 1713-1717, Soc. Explor. Geophys. (2017).

30. Martin R., Chevrot S., Komatitsch D., Seoane L., Spangenberg H., Wang Y., Dufrechou G., Bonvalot S., and Bruinsma S., "A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data-application to the western Pyrenees," Geophys. J. Int., 209, No. 1, 406-424 (2017).

31. May D. A. and Knepley M. G., "Optimal, scalable forward models for computing gravity anomalies," Geophys. J. Int., 187, No. 1, 161-177 (2011).

32. Nagy D., "The gravitational attraction of a right rectangular prism," Geophys., 31, No. 2, 362-371 (1966).

33. Nagy D., Papp G., and Benedek J., "The gravitational potential and its derivatives for the prism," J. Geodesy, 74, No. 7-8, 552-560 (2000).

34. Okabe M., "Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies," Geophys., 44, No. 4, 730-741 (1979).

35. Paul M. K., "The gravity effect of a homogeneous polyhedron for three-dimensional interpretation," Pure Appl. Geophys., 112, No. 3, 553-561 (1974).

36. Plouff D., "Gravity and magnetic fields of polygonal prisms and application to magnetic terrain corrections," Geophys., 41, No. 4, 727-741 (1976).

37. Pohanka V., "Optimum expression for computation of the gravity field of a homogeneous polyhedral body," Geophys. Prospect., 36, No. 7, 733-751 (1988).

38. Hamayun and Prutkin, I. and Tenzer, R., "The optimum expression for the gravitational potential of polyhedral bodies having a linearly varying density distribution, J. Geodesy, 83, No. 12, 1163 (2009).

39. Rao, C. V. and Raju, M. L. and Chakravarthi, V., "Gravity modelling of an interface above which the density contrast decreases hyperbolically with depth," J. Appl. Geophys., 34, No. 1, 63-67 (1995).

40. Ren Z., Chen C, Zhong Y., Chen H., Kalscheuer T, Maurer H., Tang J., and Hu X., Exact gravity field for polyhedrons with polynomial density contrasts of arbitrary orders, arXiv: 1810.11768 (2018).

41. Ren Z., Chen C., Zhong Y., Chen H., Kalscheuer T., Maurer H., Tang J., and Hu X., "Recursive analytical formulae of gravitational fields and gradient tensors for polyhedral bodies with polynomial density contrasts of arbitrary non-negative integer orders," Surv. Geophys., 41, 695-722 (2020).

42. Ren Z., Tang J., Kalscheuer T., and Maurer H., "Fast 3-D large-scale gravity and magnetic modeling using unstructured grids and an adaptive multilevel fast multipole method," J. Geophys. Res., Solid Earth, 122, No. 1, 79-109 (2017).

43. Saad Y., Iterative Methods for Sparse Linear Systems, SIAM (2003).

44. Talwani M., Worzel J. L., and Landisman M., "Rapid gravity computations for two-dimensional bodies with application to the Mendocino submarine fracture zone," J. Geophys. Res., 64, No. 1, 49-59 (1959).

45. Tsoulis D., "Analytical computation of the full gravity tensor of a homogeneous arbitrarily shaped polyhedral source using line integrals," Geophys., 77, No. 2, F1-F11 (2012).

46. Zhang J. and Jiang L., "Analytical expressions for the gravitational vector field of a 3-D rectangular prism with density varying as an arbitrary-order polynomial function," Geophys. J. Int., 210, No. 2, 1176-1190 (2017).

Received November 25, 2020 Revised January 22, 2021 Accepted February 26, 2021

Dulus Kh. Ivanov

M. K. Ammosov North-Eastern Federal University, Institute of Mathematics and Informatics, 48 Kulakovsky Street, Yakutsk 677000, Russia; Yakutsk Branch of the Regional Scientific and Educational Mathematical Center "Far Eastern Center of Mathematical Research", 48 Kulakovsky Street, Yakutsk 677000, Russia i,am.djoos@gmail.com

Petr N. Vabishchevich

Nuclear Safety Institute of RAS,

52 B. Tulskaya Street, Moscow 115191, Russia;

Academy of Science of the Republic of Sakha (Yakutia),

33 Lenin Ave., Yakutsk 677007, Russia

vabishchevich@gmail.com