Discontinuous Galerkin method for wave propagation in elastic media with inhomogeneous inclusions

1 Setting of numerical experiments

In numerical experiments we considered three variants of a computational domain geometry. Figure 1 presents the description of the indicated domains. In these cases we considered a rectangular domain with the width of 7000 m and the depth of 3000 m. The dimensions of the collector were 2000 m by 100 m. The depth of the collector was 1500 m in all the cases (see Fig. 1).

Figure 1.

Geometrical parameters of the model.

The domain of the collector was specified in two different ways: the sequence of parallel macro-cracks with the media parameters in the domain of collector being

1. the same as the parameters of the host medium;

2. different from the parameters of the host medium.

For the sake of brevity, below we call them the first and the second models.

In order to simplify the interpretation of seismic responses in the proposed models, we also consider the response from the model of a homogeneous isotropic bounded stratum.

We used a linear elastic model of the medium with the following parameters: the velocity of longitudinal waves for the host medium was c_p = 3200 m/s, the velocity of transverse waves was c_s = 1780 m/s, the density of the medium was ρ = 2450 kg/m³. If the interface of the media exists (see Figs. 1b and 1c), the characteristics of the lower layer, i.e., the velocities and density, are greater by 25% and 4%, respectively. In the case of a homogeneous isotropic layer its parameters were c_p = 2400 m/s, c_s = 1330 m/s, ρ = 2450 kg/m³. The second model had the following parameters in the fractured domain of collector: c_p = 2880 m/s, c_s = 1600 m/s, ρ = 2450 kg/m³.

We consider each of the presented variants of representation of the fractured domain in more detail. One method of specification of the collector is an explicit description of cracks in the rock [4, 14, 21, 22]. In our numerical experiments with the discrete model we consider the case of parallel cracks. Similar studies of the wave response from a system of randomly distributed and intersecting cracks were presented in [23, 24]. Presumably, the second model is more correct.

We simulate cracks and fractures in collectors by the model of an infinitely thin fluid-saturated crack [14]. In this model we have no need in significant refinement of calculation grids and hence we can avoid the related computational expenses. The geometrical parameters of the discrete and combined models of the collector are presented in Fig. 2.

Figure 2.

Position of cracks in the cluster.

Constructing calculation grids, we used the Triangle and Metis libraries [9, 25]. The grid mesh size (see Fig. 3) was specified non-uniformly to increase the accuracy of calculations in the domain of crack cluster and to reduce computational cost.

Figure 3.

Example of construction of computation unstructured grid.

The initial condition was specified by a plane wave with the amplitude in the form of Berlage pulse [21] and the length 200 m.

These settings of numerical experiment used geometric parameters corresponding to calculations performed with the use of the grid-characteristic method of [22] considering a wave response from a system of subvertical macro-cracks. A detailed study of the response from a single subvertical macro-crack was presented in [14].

2 Glide contact condition for discontinuous Galerkin method

The discontinuous Galerkin method for unstructured grids [10] forms the base of the software computational complex used in this work [16, 19]. The method has a series of positive features [6, 17] and the authors consider the ability to implement a high order approximation in spatial coordinates and in time as the basic feature among those ones, this ability is much-needed due to the wave nature of the solved problem.

In this paper we use Lagrange polynomials of 4th degree as a system of basis polynomials. The integration in time was performed by the Dormand–Prince integrator with adaptive step in time. Software implementation was optimized for calculations on high-performance multiprocessor computer systems.

Let us obtain determining relations for the glide contact condition forming the base of the model of an infinitely-thin fluid-saturated crack. In discontinuous Galerkin method, relations for any contact conditions (full sticking, i.e., at an interface between two elastic media the velocity and the traction are continuous, sliding, Coulomb’s law of friction) are reduced at the end to solution of a one-dimensional Riemann problem [15] subject to some additional relations. In particular, in the case of sliding they are the continuity of the normal velocity component and the absence of shear stresses on the interface. Paper [32] presents a detailed description what is the origin of the necessity to solve Riemann’s problem in the solution of a linear system of elasticity equations with the use of discontinuous Galerkin’s method and how one can approximately reduce it to a one-dimensional problem. This paper presents the solution for the case of full sticking only. Omitting details, we give a brief description of the numerical method and obtain an expression for the numerical flux for the glide condition in the two-dimensional case, which the authors did not meet previously.

The two-dimensional system of elasticity equations for the case of isotropic space has the following matrix form [15] in stress and velocity variables:

where u is the vector of 5 unknown variables u = (σ_xx , σ_yy, σ_xy, ν_x, ν_y)^T. Here and below we assume summation over all repeating indices. The eigenvalues of the matrices A_pq and B_pq are

The integration domain is divided into triangles T^(m) and enumerated. We assume that the matrices A_pq and B_pq are constant inside of T^(m). The solution in the triangle T^(m) is approximated by a linear combination of (N + 1)(N + 2)/2 polynomial functions Φ_k(x, y) of degree not exceeding N, not dependent on time, and forming a basis with the support T^(m) and with time-dependent coefficients

Multiplying (2.1) by the basis function Φ_k and integrating over the triangle T^(m), we get

Further, applying the formula of integration by parts, we obtain

The second summand appears due to the discontinuity of the solution u_h and the matrices A_pq, B_pq on the boundary of the triangle T^(m) in the general case. Here Fph,j is the numerical flux over the j th side of the triangle in the global system of coordinates, by (∂T^(m))_j we have denoted the sides of the triangle T^(m), j = 1, 2, 3 (see Fig. 4).

Figure 4.

Element of calculation grid.

Let Tpqmj be the matrix of transformation to the coordinate system X′Y′ (orthogonal transformation) related to the j th side of the triangle T^(m) (see Fig. 4). In this case the flux Fph,j can be obtained approximately, i.e.,

where ulb is the solution to the one-dimensional Riemann problem (see Fig. 5). In Figure 5, the characteristics of the medium marked by the symbol ⁽⁻⁾ correspond to the triangle T^(m), the other marked by ⁽⁺⁾ correspond to the triangle adjacent with respect to the j th side, u^a,b,c,d are the unknown states in the Riemann problem in the coordinate system X′Y′. Writing down the left and right limits of the solution relative to the considered j th side (see Fig. 4), we have

Figure 5.

Solution to one-dimensional Riemann’s problem.

Applying the Rankine–Hugoniot jump condition to each characteristic and taking into account the continuity of the normal velocity component and the absence of shear stresses on the contact, we get the following system of equations:

where α_i are the unknown coefficients, r_i are the columns of the matrix R_pq composed from the eigenvectors of the matrices Apq(m) and Apq(mj), i.e.,

Solving system of equations (2.7), we obtain the following expression for u^b:

In the case of full sticking when the condition of continuity of all components of velocity and stress tensor holds true on the contact, the coefficient α₁ coincides with those for the case of sliding, i.e., the longitudinal waves spread in the same way in these cases. In the case of full sticking the coefficient α₂ has the form

It is seen that α₂ for sliding (see formula (2.9)) is obtained from formula (2.10) by setting the longitudinal velocity and the shear stress equal to zero in the adjacent cell, i.e., cs+ = σxy+ = 0, which is equivalent to a contact with an acoustic medium without transverse waves. This fact clarifies the physical sense of the model of an infinitely thin fluid-filled crack.

Combining expressions (2.4), (2.5), (2.6), (2.9), (2.10), we complete the construction of the numerical scheme, a more detailed derivation was presented in [10, 16, 32].

3 Results of numerical experiments

The results are presented as a series of wave patterns for fixed time moments and seismograms. This paper contains the results of only three of nine calculations (three variants of the calculation domain geometry by (model of homogeneous isotropic layer + two models of fractured collector)). Figure 6a shows calculated seismograms of the horizontal (left) and vertical (right) velocity components for the scheme without media interface (see Fig. 1a). For the sake of brevity, below we call them ν_x- and ν_y-seismograms, respectively.

Figure 6.

Calculation seismograms for the scheme without media interface.

Figure 7 presents the calculated ν_x-(left) and ν_y-seismograms (right) for schemes with the media interface below the cluster (see Figs. 1b and 1c).

Figure 7.

Calculation seismograms for specification of the cluster according to the second model.

Figure 8 presents the wave pattern for a fixed time moment in the case of a scheme without media interface (see Fig. 1a) corresponding to seismogram 6a.

Figure 8.

Example of wave patterns for the scheme without media interface in the model of a homogeneous isotropic medium.

Figure 9 presents the wave pattern corresponding to seismogram 7 for the scheme with the media interface below the cluster specified by the second method.

Figure 9.

Example of wave patterns for the scheme with the cluster specified by the second method with the media interface below the cluster.

All calculations described in the paper were performed on the cluster of the Information Computing Center of the Novosibirsk State University.

4 Analysis of calculated wave patterns and seismograms

Let us consider the seismogram presented in Fig. 6a in more detail. The ν_y-seismogram clearly shows longitudinal p-waves reflected from the upper (II) and lower (III) boundaries of the cluster. Sequences of double waves reflected from the day surface (VII) are positioned after them. Within the model considered here, we can calculate the thickness of the cluster based on the values of the velocities of elastic waves in the medium and the time interval between fixation of waves obtained from the seismogram and this thickness coincides with the original one with a high accuracy.

The ν_x-seismogram clearly shows two symmetric (due to the symmetry of geometric parameters of the problem) wave fronts of hyperbolic form (I and VI) related to longitudinal and exchange waves diffracted from the borders of the cluster, respectively. A similar effect can be seen in the ν_y-seismogram, however, in comparison with components II and III, its contribution is not great. In the seismogram presented in Figure 6b these effects (I, VI) are seen less, in this case the media interface between the cluster and the host medium is absent and hence the reflection occurs at border subvertical macro-cracks. The diffracted exchange wave coming from the borders of the cluster is also clearly seen in Fig. 7. It is worth noting not only a quantitative, but also a qualitative difference of the indicated fronts. Comparing (I, VI) in Figs. 6a and 6b, we can see redistribution of phase characteristics of signals. In both cases presented in Figs. 6a and 6b) the distance between the fronts of I provides an information on the position and horizontal size of the cluster.

Sequences of horizontal plane wave fronts (IV) formed as the result of multiple reflections between cracks of the cluster (exchange waves) and giving us an information on the fractured structure of the host medium are marked on the ν_x-seismogram shown in Fig. 6b. In practice, it is difficult to detect exchange waves due to dissipation existing in the medium. The dependence of parameters of the periodic structure of exchange waves on the length of the incident wave and geometrical characteristics of the fractured domain is a subject of separate study.

The ν_y-seismogram presented in Fig. 7 indicates a retardation of the wave reflected from the media interface (VIII) in the central domain (at the position of the cluster). This is explained by differences in the characteristics of the media filling the cluster and the environment.

All the conclusions presented in this section of the paper are confirmed by the corresponding pictures of wave fields (see Figs. 8 and 9), which, in comparison with seismograms, provide not only additional information, but also give a more clear picture of the situation. The authors intentionally do not rely on wave patterns in their arguments because in real life there is no possibility to get them in geological-prospecting work.

5 Results

1. The discontinuous Galerkin method on unstructured computation grids is adapted and implemented for simulation of the response from subvertical macro-cracks in carbonate rocks.

2. Formulas for numerical flux are obtained in the case of a glide contact condition for solving the contact break problem in the fluid-saturated crack model.

3. Wave patterns are obtained using this method at fixed time moments as well as seismograms from the system of subvertical macro-cracks.

4. Border waves caused by interaction of an incident seismic pulse with the borders of the fractured layer are obtained.

5. The feasibility of detailed study of wave processes initiated by an incident seismic pulse acting on the system of subvertical macro-cracks is shown as well as for solution of direct seismic exploration problems with the use of the discontinuous Galerkin method.

6. The models proposed here are compared and the feasibility to take into account the inter-crack wave interactions in the simulation of a fractured layer is shown.

It is worth to point out that the comparison of calculated seismograms with field ones can help in drawing conclusions concerning geometrical and mechanical properties of fractured collectors.