Simulation of elastic wave propagation in geological media: Intercomparison of three numerical methods

For wave propagation in heterogeneous media, we compare numerical results produced by grid-characteristic methods on structured rectangular and unstructured triangular meshes and by a discontinuous Galerkin method on unstructured triangular meshes as applied to the linear system of elasticity equations in the context of direct seismic exploration with an anticlinal trap model. It is shown that the resulting synthetic seismograms are in reasonable quantitative agreement. The grid-characteristic method on structured meshes requires more nodes for approximating curved boundaries, but it has a higher computation speed, which makes it preferable for the given class of problems.     

Numerical solution of seismic exploration problems in the Arctic region by applying the grid-characteristic method

The goal of this paper is the numerical solution of direct problems concerning hydrocarbon seismic exploration on the Arctic shelf. The task is addressed by solving a complete system of linear elasticity equations and a system of acoustic field equations. Both systems are solved by applying the grid-characteristic method, which takes into account all wave processes in a detailed and physically correct manner and produces a solution near the boundaries and interfaces of the integration domain, including the interface between the acoustic and linear elastic media involved. The seismograms and wave patterns obtained by numerically solving these systems are compared. The effect of ice structures on the resulting wave patterns is examined.     

Numerical simulation of seismic activity by the grid-characteristic method

Seismic activity in homogeneous and layered enclosing rock masses is studied. A numerical mechanical-mathematical model of a hypocenter is proposed that describes the whole range of elastic perturbations propagating from the hypocenter. Synthetic beachball plots computed for various fault plane orientations are compared with the analytical solution in the case of homogeneous rock. A detailed analysis of wave patterns and synthetic seismograms is performed to compare seismic activities in homogeneous and layered enclosing rock masses. The influence exerted by individual components of a seismic perturbation on the stability of quarry walls is analyzed. The grid-characteristic method is used on three-dimensional parallelepipedal and curvilinear structured grids with boundary conditions set on the boundaries of the integration domain and with well-defined contact conditions specified in explicit form.     

Combined grid-characteristic method for the numerical solution of three-dimensional dynamical elastoplastic problems

A combined method blending the advantages of smoothed particles hydrodynamics (SPH) and the grid-characteristic method (GCM) is proposed for simulating elastoplastic bodies. Various grid methods, including the GCM, have long been used for the numerical simulation of elastoplastic media. This method applies to the simulation of wave processes in elastic media, including elastic impacts, in which case an advantage is the use of moving tetrahedral meshes. Additionally, fracture processes can be simulated by applying various fracture criteria. However, this is a technically complicated task with the accuracy of the results degrading due to the continual updating of the grid. A more suitable approach to the simulation of processes involving substantial fractures and deformations is based on SPH, which is a meshless method. However, this method also has shortcomings: it produces spurious modes, and the simulation of oscillations requires particle refinement. Thus, two families of methods are available that are optimal as applied to two different groups of problems. However, a realworld problem can frequently be a mixed one, which requires a substantial tradeoff in the numerical methods applied. Aimed at solving such problems, a combined GCM-SPH method is developed that blends the advantages of two constituting techniques and partially eliminates their shortcomings.     

On the Calculation of Border and Contact Nodes by a Grid-Characteristic Method on Non-Periodic Tetrahedral Grids

A grid-characteristic method for the numerical simulation of wave processes in continuum mechanics was initially proposed, and has been successfully applied to periodic hexagonal computational grids. Later it was proposed to adapt this method to non-periodic triangle and tetrahedral grids, and wide computational experience has been gained. However, this approach encounters some difficulties in the calculation of border and contact points when applied to various grid configurations in areas with complex geometries. In this paper, limitations of the method which cause such problems are considered, and some improvements to overcome them are proposed.     

Detailed simulation of the pulsating detonation wave in the shock-attached frame

The paper is devoted to the numerical investigation of the stability of propagation of pulsating gas detonation waves. For various values of the mixture activation energy, detailed propagation patterns of the stable, weakly unstable, irregular, and strongly unstable detonation are obtained. The mathematical model is based on the Euler system of equations and the one-stage model of chemical reaction kinetics. The distinctive feature of the paper is the use of a specially developed computational algorithm of the second approximation order for simulating detonation wave in the shock-attached frame. In distinction from shock capturing schemes, the statement used in the paper is free of computational artifacts caused by the numerical smearing of the leading wave front. The key point of the computational algorithm is the solution of the equation for the evolution of the leading wave velocity using the second-order grid-characteristic method. The regimes of the pulsating detonation wave propagation thus obtained qualitatively match the computational data obtained in other studies and their numerical quality is superior when compared with known analytical solutions due to the use of a highly accurate computational algorithm.     

Numerical treatment of classic scattering from bounded objects

Classic scattering from objects of arbitrary shape must generally be treated by numerical methods. It has proven very difficult to describe scattering from general bounded objects without resorting to frequency-limiting approximations. The starting point of many numerical methods is the Helmholtz integral representation of a given wavefield. From that point several departures are possible for constructing computationally feasible approximate schemes. To date, attempts at direct solutions have been rare.One method (originated by P. Waterman) that attacks the exact numerical solution for a very broad class of problems begins with the Helmholtz integral representations for a point exterior and interior to the target in a partial wave expansion. After truncating the partial wave space, one arrives at a set of matrix equations useful in describing the field. This method is often referred to as the T-matrix method, null-field, or extended integral equation method. It leads to a unique solution of the exterior boundary integral equation by incorporating the interior solution (extinction theorem) as a constraint. In principle, there are no theoretical limitations on frequency, although numerical complications can arise and must be appropriately dealt with for the method to be computationally reliable.For submerged objects the formalism will be outlined for acoustical scattering from targets that are rigid; sound-soft and penetrable; elastic solids; elastic shells; and layered elastic objects. Finally, illustrations of several numerical examples for the above will be presented to emphasize specific response features peculiar to a variety of targets.     

3D dynamic responses of a multi-layered transversely isotropic saturated half-space under concentrated forces and pore pressure

Dynamic Green's function plays an important role in the study of various wave radiation, scattering and soil-structure interaction problems. However, little research has been done on the response of transversely isotropic saturated layered media. In this paper, the 3D dynamic responses of a multi-layered transversely isotropic saturated half-space subjected to concentrated forces and pore pressure are investigated. First, utilizing Fourier expansion in circumferential direction accompanied by Hankel integral transform in radial direction, the wave equations for transversely isotropic saturated medium in cylindrical coordinate system are solved. Next, with the aid of the exact dynamic stiffness matrix for in-plane and out-of-plane motions, the solutions for multi-layered transversely isotropic saturated half-space under concentrated forces and pore pressure are obtained by direct stiffness method. A FORTRAN computer code is developed to achieve numerical evaluation of the proposed method, and its accuracy is validated through comparison with existing solutions that are special cases of the more general problems addressed. In addition, selected numerical results for a homogeneous and a layered material model are performed to illustrate the effects of material anisotropy, load frequency, drainage condition and layering on the dynamic responses. The presented solutions form a complete set of Green's functions for concentrated forces (including horizontal load in x(y)-direction, vertical load in z-direction) as well as pore pressure, which lays the foundation for further exploring wave propagation of complex local site in a layered transversely isotropic saturated half-space by using the BEMs.     

Grid-characteristic method on unstructured tetrahedral meshes

The goal of this paper is to develop a grid-characteristic method intended for high-performance computer systems and implemented on unstructured tetrahedral hierarchical meshes with the use of a multiple time step and high-order interpolation, including interpolation with a limiter, piecewise parabolic interpolation, and monotone interpolation. The method is designed for simulating complex three-dimensional dynamical processes in heterogeneous media. It involves accurately stated contact conditions and produces physically correct solutions of problems in seismology and seismic exploration. Hierarchical meshes make it possible to take into account numerous inhomogeneous inclusions (cracks, cavities, etc.) and to solve problems in a real-life formulation. The grid-characteristic method enables the use of a multiple time step. As a result, the computation time is considerably reduced and the efficiency of the method is raised. The method is parallelized on a computer cluster with an optimal use of system resources.     

High-performance computer simulation of wave processes in geological media in seismic exploration

A class of problems arising in seismic exploration are investigated, namely, seismic signal propagation in multilayered geological rock and near-surface disturbance propagation in massive rock with heterogeneities, such as empty or filled fractures and cavities. Numerical solutions are obtained for wave propagation in such highly heterogeneous media, including those taking into account the plastic properties of the rock, which can be manifested near a seismic gap or a wellbore. All types of explosion-generated elastic and elastoplastic waves and waves reflected from fractures and the boundaries of the integration domain are analyzed. The identification of waves in seismograms recorded with near-surface receivers is addressed. The grid-characteristic method is used on triangular, parallelepipedal, and tetrahedral meshes with boundary conditions set on the rock-fracture interface and on free surfaces in explicit form. The numerical method proposed is suitable for the study of the interaction between seismic waves and heterogeneous inclusions, since it ensures the most correct design of computational algorithms on the boundaries of the integration domain and at media interfaces. A parallel software code implemented with the help of OpenMP and MPI was used to execute computations on parallelepipedal and tetrahedral grids.     

A simple and robust boundary treatment for the forced Korteweg–de Vries equation

In this paper, we propose a simple and robust numerical method for the forced Korteweg–de Vries (fKdV) equation which models free surface waves of an incompressible and inviscid fluid flow over a bump. The fKdV equation is defined in an infinite domain. However, to solve the equation numerically we must truncate the infinite domain to a bounded domain by introducing an artificial boundary and imposing boundary conditions there. Due to unsuitable artificial boundary conditions, most wave propagation problems have numerical difficulties (e.g., the truncated computational domain must be large enough or the numerical simulation must be terminated before the wave approaches the artificial boundary for the quality of the numerical solution). To solve this boundary problem, we develop an absorbing non-reflecting boundary treatment which uses outward wave velocity. The basic idea of the proposing algorithm is that we first calculate an outward wave velocity from the solutions at the previous and present time steps and then we obtain a solution at the next time step on the artificial boundary by moving the solution at the present time step with the velocity. And then we update solutions at the next time step inside the domain using the calculated solution on the artificial boundary. Numerical experiments with various initial conditions for the KdV and fKdV equations are presented to illustrate the accuracy and efficiency of our method.     

Numerical computation of sharp travelling waves of a degenerate diffusion–reaction equation arising in biofilm modelling

Many degenerate diffusion–reaction equations permit sharp travelling wave solutions that describe the propagation of an interface with finite speed. If the equation is at least double degenerate, the derivative of the travelling wave solution can blow up at the interface, which poses considerable challenges for the computation of the travelling wave speed. We propose a numerical method for this problem that is based on the idea to approximate the multiple degenerate problem by a family of simple degenerate problems. For the latter we propose an interval-bracketing algorithm based on the theory of Sanchez-Garduno and Maini. The travelling wave speed of the original problem is obtained as the limit of the travelling wave speeds of the auxiliary problems. The performance of the method is investigated in a numerical simulation experiment for a problem that arises in the mathematical modelling of biofilm processes.     

Energetic BEM-FEM coupling for the numerical solution of the damped wave equation

Time-dependent problems modeled by hyperbolic partial differential equations can be reformulated in terms of boundary integral equations and solved via the boundary element method. In this context, the analysis of damping phenomena that occur in many physics and engineering problems is a novelty. Starting from a recently developed energetic space-time weak formulation for the coupling of boundary integral equations and hyperbolic partial differential equations related to wave propagation problems, we consider here an extension for the damped wave equation in layered media. A coupling algorithm is presented, which allows a flexible use of finite element method and boundary element method as local discretization techniques. Stability and convergence, proved by energy arguments, are crucial in guaranteeing accurate solutions for simulations on large time intervals. Several numerical benchmarks, whose numerical results confirm theoretical ones, are illustrated and discussed.     

Explicit periodic travelling waves for a discrete lambda–omega reaction–diffusion system

In 1973, Kopell and Howard introduced a λ–ω reaction–diffusion system and found an explicit family of periodic travelling wave solutions lying on circles with radius less than 1. Since λ–ω systems represent universal models for studying chemical processes, and onset of turbulent behaviour, etc., explicit solutions of λ–ω systems with delays or discrete λ–ω systems can be of further help when the only method for obtaining other solutions is through numerical computation. There are now much investigations of various λ–ω systems. However, it is of interest to note that none attempts to find explicit travelling wave solutions. In this paper, we investigate the existence of explicit solutions for the simplest Euler scheme of a λ–ω system with delays or advancements which is described as a coupled pair of partial difference equations. We are able to provide necessary as well as sufficient conditions for the existence of numerical periodic travelling wave solutions. Additionally, we also provide some examples to show that our explicit solutions are qualitatively different from those found by Kopell and Howard and hence they may be of interests for specialists in the area of reaction–diffusion systems.     

Three-Dimensional Green's Function for a Layered Isotropic Half-space

A numerical approach to calculate the Green's function for a layered half space is presented. It is based on the precise integration method (PIM), which is an efficient and accurate numerical method for the solution of one order ordinary differential equations. In the numerical implementation, the layered half space is divided into numerous mini-layers; and the dual vector form of the wave motion equation is introduced to combine two adjacent mini-layers/layers. The advantages of the proposed algorithm are: (a) it overcomes the exponent overflow generally encountered with employing the transfer matrix method; (b) it avoids solving the intractable transcendental functions in the stiffness matrix method and the huge matrix calculation in the thin layer method; (c) it imposes no limit to the thickness of layered strata and ensures convergence at high-frequency range. (© 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Numerical study of wave propagation in porous media with the use of the grid-characteristic method

Elastic wave propagation in a porous medium is numerically studied by applying the grid-characteristic method. On the basis of direct measurements of reflected and transmitted wave amplitudes, the reflection and decay coefficients are investigated as depending on the degree of porosity (percentage of the pore volume) and on the type of the filling substance (solid, liquid, or nothing). The reflection and decay coefficients are shown to be closely related to the porosity of the medium, which can be used in geological applications (estimation of porosity) and engineering applications (acoustic response attenuation).     

Numerical solution of electromagnetic scattering from a large partly covered cavity

The paper focuses on the numerical study of electromagnetic scattering from two-dimensional (2D) large partly covered cavities, which is described by the Helmholtz equation with a nonlocal boundary condition on the aperture. The classical five-point finite difference method is applied for the discretization of the Helmholtz equation and a linear approximation is used for the nonlocal boundary condition. We prove the existence and uniqueness of the numerical solution when the medium in the cavity is y-direction layered or the number of the mesh points on the aperture is large enough. The fast algorithm proposed in Bao and Sun (2005) [2] for open cavity models is extended to solving the partly covered cavity problem with (vertically) layered media. A preconditioned Krylov subspace method is proposed to solve the partly covered cavity problem with a general medium, in which a layered medium model is used as a preconditioner of the general model. Numerical results for several types of partly covered cavities with different wave numbers are reported and compared with those by ILU-type preconditioning algorithms. Our numerical experiments show that the proposed preconditioning algorithm is more efficient for partly covered cavity problems, particularly with large wave numbers.     

Inverse problem for wave propagation in a perturbed layered half-space

This paper is concerned with the inverse medium scattering problem in a perturbed, layered, half-space, which is a problem related to the seismologial investigation of inclusions inside the earth’s crust. A wave penetrable object is located in a layer where the refraction index is different from the other part of the half-space. Wave propagation in such a layered half-space is different from that in a homogeneous half-space. In a layered half-space, a scattered wave consists of a free wave and a guided wave. In many cases, only the free-wave far-field or only the guided-wave far-field can be measured.We establish mathematical formulas for relations between the object, the incident wave and the scattered wave. In the ideal condition where exact data are given, we prove the uniqueness of the inverse problem. A numerical example is presented for the reconstruction of a penetrable object from simulated noise data.