Finite Difference Grid Research Articles

An Integrated Taylor Expansion and Least Squares Approach to Enhanced Acoustic Wave Staggered Grid Finite Difference Modeling

The staggered grid finite difference method has emerged as one of the most commonly used approaches in finite difference methodologies due to its high computational accuracy and stability. Inevitably, discretizing over time and space domains in finite difference methods leads to numerical artifacts. This paper introduces a novel approach that combines the widely used Taylor series expansion with the least squares method to effectively suppress numerical dispersion. We have derived the coefficients for the staggered grid finite difference method by integrating Taylor series expansions with the least squares method. To validate the effectiveness of our approach, we conducted analyses on accuracy, dispersion, and stability, alongside simple and complex numerical examples. The results indicate that our method not only inherits the capabilities of the original Taylor series and least squares methods in suppressing numerical dispersion across small and medium wavenumber ranges but also surpasses the original methods. Moreover, it demonstrates robust dispersion suppression capabilities at high wavenumber ranges.

Numerical Modeling of Scholte Wave in Acoustic-Elastic Coupled TTI Anisotropic Media

Numerical modeling of acoustic-elastic media is helpful for seismic exploration in the deepwater environment. We propose an algorithm based on the staggered grid finite difference to simulate wave propagation in the interface between fluid and transversely isotropic media, where the interface does not need to consider the boundary condition. We also derive the stability conditions of the proposed method. Scholte waves, which are generated at the seafloor, exhibit distinctly different propagation behaviors than body waves in ocean-bottom seismograms. Numerical examples are used to characterize the wavefield of Scholte waves and discuss the relationship between travel time and the Thomsen parameters. Thomsen parameters are assigned clear physical meanings, and the magnitude of their values directly indicates the strength of the anisotropy in the media. Numerical results show that the velocity of the Scholte wave is positively correlated with ε and negatively correlated with δ. And the curve of the arrival time of the Scholte wave as a whole is sinusoidal and has no symmetry in inclination. The velocity of the Scholte wave in azimuth is positively related to the polar angle. The energy of the Scholte wave is negatively correlated with the distance from the source to the fluid-solid interface. The above results provide a basis for studying oceanic Scholte waves and are beneficial for utilizing the information provided by Scholte waves.

Propagation characteristics of show sea seismic wave in different underwater media

The propagation characteristics of sea seismic wave (SSSW) in different underwater media are different. This study mainly focuses on the propagation characteristics of SSSW in different underwater media to improve the detection ability of SSSW mines. The propagation of SSSW in different underwater media, including clay, sandy clay, sandstone, mudstone, hard stone, and dolomite, is simulated using the staggered grid finite difference method. The dispersive energy maps of different underwater media are extracted using the f-k dispersive energy extraction method. Compared with the results of Scholte wave dispersion characteristic equations, the simulation results are consistent with the theoretical results. At the same time, the attenuation law of the SSSW train is obtained. The results provide a theoretical basis for developing SSSW mine fuses.

Simulate the elastic wavefields in media with an irregular surface topography based on staggered grid finite difference

Abstract Wave equation forward modeling is a useful method to study the propagation regulation of seismic wavefields. Finite difference (FD) is one of the most extensively employed numerical approaches for computing wavefields in earthquake and exploration seismology. However, the FD approach relying on regular grids often struggles to calculate wavefields in regions featuring surface topographies. The elastic wave equation can more accurately describe the propagation of seismic wavefields in elastic media compared to the acoustic wave equation. We introduce a new FD scheme to calculate the elastic wavefields in an isotropic model with a surface topography. The novel approach can use a conventional staggered grid FD (SGFD) approach based on regular grids. A new elastic model with a horizontal surface is first obtained from the nearby surface's elastic properties and the undulating terrain elevation. We subsequently employ a topography-related strategy to eliminate the effects of surface topographies on the seismic wavefields in models with irregular surface topographies. The merits of our proposed scheme lie in its ability to stable numerically compute wavefields in models with irregular surface topographies without altering the conventional SGFD relying on regular grids. To validate the effectiveness and practicality of our method, we utilize elastic models featuring complex surface topographies. Numerical experiments demonstrate that our approach efficiently calculates elastic wavefields in isotropic media with irregular topographies based on conventional SGFD.

Enhancing computational efficiency in 3-D seismic modelling with half-precision floating-point numbers based on the curvilinear grid finite-difference method

SUMMARY Large-scale and high-resolution seismic modelling are very significant to simulating seismic waves, evaluating earthquake hazards and advancing exploration seismology. However, achieving high-resolution seismic modelling requires substantial computing and storage resources, resulting in a considerable computational cost. To enhance computational efficiency and performance, recent heterogeneous computing platforms, such as Nvidia Graphics Processing Units (GPUs), natively support half-precision floating-point numbers (FP16). FP16 operations can provide faster calculation speed, lower storage requirements and greater performance enhancement over single-precision floating-point numbers (FP32), thus providing significant benefits for seismic modelling. Nevertheless, the inherent limitation of fewer 16-bit representations in FP16 may lead to severe numerical overflow, underflow and floating-point errors during computation. In this study, to ensure stable wave equation solutions and minimize the floating-point errors, we use a scaling strategy to adjust the computation of FP16 arithmetic operations. For optimal GPU floating-point performance, we implement a 2-way single instruction multiple data (SIMD) within the floating-point units (FPUs) of CUDA cores. Moreover, we implement an earthquake simulation solver for FP16 operations based on curvilinear grid finite-difference method (CGFDM) and perform several earthquake simulations. Comparing the results of wavefield data with the standard CGFDM using FP32, the errors introduced by FP16 are minimal, demonstrating excellent consistency with the FP32 results. Performance analysis indicates that FP16 seismic modelling exhibits a remarkable improvement in computational efficiency, achieving a speedup of approximately 1.75 and reducing memory usage by half compared to the FP32 version.

Numerical simulation of borehole compressional wave and shear wave in 3D vug formation

The staggered grid finite difference method (SGFDM) of a monopole source is used to simulate a three-dimensional vug reservoir model and study the effect of acoustic logging responses of different vug models on radial probing depth. The results show that the first arrival of the head wave peak of the compressional wave (P-wave) and the shear wave (S-wave) is unrelated to the radius of the vug, and the amplitude of the head wave peak of the P-wave and S-wave decreases as the vug volume increases. Compared with the volume change of the vug, the radial distance from the vug wall has little influence, while the vertical source distance has large influence on the P-wave and S-wave. When there are multiple vugs in the model, the amplitudes of the P-wave and S-wave head wave peaks change sinusoidally with the angle between the vugs. The ellipsoidal vug model with the same volume has a greater influence on the P-wave and S-wave than the spherical vug model. In the ellipsoidal vug model, the axial vug size has a greater impact on the first arrival of the head wave peak, while the radial vug size significantly influences the amplitude of the head wave peak. Finally, we validate the numerical simulation conclusions by comparing them with actual logging data responses in complex formations, demonstrating the practical value of the elastic wave response simulations for vugs.

Dynamic stress response in reservoirs stimulated by fluctuating fluid pressure during pulsating hydraulic fracturing

Determining the most effective scheme for fracturing operations depends on accurate calculations of reservoir stress. Traditional fracturing theories usually focus on stable fluid injection modes, which employ static assumptions in modelling and analysis. However, pulsating hydraulic fracturing (PHF) introduces fluctuating fluid pressure to stimulate reservoirs, which requires considering dynamic factors and exploring the dynamic stress response within the reservoir. Furthermore, reservoirs are often classified as homogeneous continuous or porous media materials, but the boundaries between these classifications are sometimes apparent. This paper investigated the development of both an elastodynamics model and a poroelastodynamics model for calculating dynamic stress responses during PHF. The staggered grid finite difference method addressed these models, with the outer boundary employing the perfectly matched layer (PML) method to simulate an infinite reservoir. The numerical simulation results were extensively verified and compared with analytical solutions, and the differences between various models and their relevant conditions were thoroughly analyzed. Analyzing fluid pressure fluctuation parameters revealed noteworthy insights into reservoir stress response characteristics. Results indicate that as the Biot coefficient and permeability increase, the disparities between the elastodynamics and poroelastodynamics models become prominent. In practical terms, high Biot coefficient or high permeability reservoirs necessitate using the poroelastodynamics model for accurate reservoir stress analysis. As the frequency of fluid pressure fluctuation increases, the amplitude of response of reservoir circumferential stress also increases, while the response amplitude of radial stress and fluid pressure decreases. These findings offer valuable theoretical guidance for parameter design in engineering applications during PHF.

ELASTOHYDRODYNAMIC EFFECTS IN LUBRICATION OF JOURNAL BEARINGS UNDER TURBULENT CONDITIONS

The current paper is aimed to investigate the deformability effects of the bearing surface on the lubrication performance of cylindrical journal bearing with finite length under turbulent conditions. The bush is coated with a thin resilient layer with a smooth surface, the radial distortions of which are of the same order of magnitude as the film thickness. The elastohydrodynamic problem is studied under elasticity conditions in accordance with the plain strain hypothesis (Winkler model). The modified Reynolds equation is carried out on the base of the turbulent lubrication model of Hirs (bulk flow model). The partial differential equation is solved numerically by successive over-relaxation technique on a finite difference grid, considering different values of coating materials, lubricant viscosity, and journal speed. Results for the performance characteristics of the bearing are given for several selected values of Reynolds number and elasticity parameters of the soft liner on the bush.

A time-domain SGFD-FK hybrid method for 2D teleseismic elastic wave modeling and inversion

AbstractFull waveform inversion (FWI) has proved to be a reliable tool for high-resolution imaging of lithospheric structures at various depths down to the upper mantle. However, when the size of the model is large, the computational burden is significative and applications are restricted to low frequencies. To tackle this issue, we developed a new 2D time-domain hybrid method to simulate high-frequency teleseismic body waves propagating through a local heterogeneous elastic Earth model: The frequency-wavenumber (FK) integration method is coupled with the staggered grid finite difference method (SGFD). The FK method is used to compute the wavefield due to obliquely incident plane P and SV waves in a 1D multilayered half space that excites a local heterogeneous region. Inside this region, the velocity-stress staggered grid FD method (SGFD) is used to accurately deal with wave propagation in heterogeneous media. Spurious waves that might be generated at the boundaries of the local region are avoided using convolutional perfectly matched layers (CPML). This new hybrid method inherits the low-memory requirements of the FK method and the accuracy, efficiency and easy implementation of the SGFD. The new hybrid method is benchmarked against the analytical FK method for some canonical models and shows good agreement with analytical solutions. Subsequently, our modeling tool is incorporated into a full waveform inversion algorithm adapted for teleseismic configurations to invert the incident P wave and its coda. The inversion is carried out using a gradient approach that is efficiently implemented via the adjoint-state method. The results suggest that our hybrid method and FWI algorithm represent a valuable tool for 2D forward and inverse regional applications using teleseismic data sets.

A linear, second-order accurate, positivity-preserving and unconditionally energy stable scheme for the Navier–Stokes–Poisson–Nernst–Planck system

Electrokinetic phenomena involving the transport of ions in fluids generally appears in biological system, hydrodynamics and geodynamics. The underlining mathematical models involve the nonlinear coupling between the fluid motion and ion transportation. Performing efficient, accurate and structure preserving computation for these electrohydrodynamic models is challenging due to their multiphysics and complex nature. The Navier–Stokes–Poisson–Nernst–Planck (NSPNP) system is a widely used electrokinetic model in the existing literature. In this paper, a linear, second-order accurate, positivity-preserving and unconditionally energy stable scheme is developed for the NSPNP system. First, a square root transformation is employed to ensure the positivity of ion concentrations. Next, a specific ordinary differential equation is introduced to deal with the nonlinear coupling terms satisfying the “zero-energy-contribution” feature, which is combined with the scalar auxiliary variable method for nonlinear potentials to achieve a decoupled energy stable scheme. Therefore, only several linearly independent equations need to be solved at each time step. The staggered grid finite difference method is used for the spatial discretizations, and the unconditional energy stability of the fully discrete scheme is established. Finally, the unique solvability of the scheme is proved. Numerical experiments are carried out to demonstrate the accuracy and performance of the proposed scheme, and the effect of electric Coulomb force on the flow and motion of ions are also illustrated.

Simulation and analysis of elastic waves in partially saturated double-porosity media based on finite difference method

Double-porosity poroelastic model takes into account the effect of mesoscopic flow induced by rock heterogeneity on dispersion and attenuation of elastic waves, and has obtained good application results in the quantitative explanation of seismic data in heterogeneous reservoirs. Wavefield simulation based on double-porosity model not only helps visualize the propagation characteristics of the elastic waves but also lays the foundation for seismic imaging. In this work, we perform wavefield simulation and analysis based on the Santos-Rayleigh model which incorporates mesoscopic and global flow in a partially-saturated double-porosity medium. Specifically, the mesoscopic flow mechanism is represented with a Zener viscoelastic model. The comparison shows that the Zener model can accurately capture the propagation characteristics of fast P-wave, but fails to describe the attenuation characteristics of slow P3 wave in the low-frequency band. It implies that Zener viscoelastic model and slow wave modes follow different mechanisms. Then the staggered grid finite-difference method is used to simulate wave propagation in a double-porosity medium, and the stiff problem is solved with a time-splitting algorithm, which can significantly improve computational efficiency. Based on the above methods, the correctness of our algorithm is verified with derived analytical solution for a P-wave source in a uniform partially saturated poroelastic medium. Analytical and numerical solutions are in good agreement and mean error is 0.33%. We provide some examples of wavefield snapshots and seismograms in homogeneous and layered heterogeneous media at seismic and ultrasonic frequencies. The simulation results demonstrate the strong attenuation of fast P-wave and no change of S-wave in the seismic band due to mesoscopic flow mechanism, which is consistent with the theoretical prediction of double-porosity model. Moreover, the energy of fast P-wave is concentrated in solid phase while slow waves are stronger in fluid phase. This work contributes to the understanding of broadband elastic wave propagation in a heterogeneous partially saturated porous medium and can be applied to the reservoir imaging with broadband geophysical data.

Setochno-kharakteristicheskiy metod povyshennogo poryadka dlya sistem giperbolicheskikh uravneniy s kusochno-postoyannymi koeffitsientami

A new approach is considered for increasing the order of accuracy of the grid-characteristic method in the region of coefficient jumps. The approach is based on piecewise polynomial interpolation for schemes of the second and third orders of accuracy for the case where the interface between the media is consistent with a finite-difference grid. The method is intended for numerical simulation of the propagation of dynamic wave disturbances in heterogeneous media. Systems of hyperbolic equations with variable coefficients are used to describe the considered physical processes. The description of the numerical method and the results of its testing are given.

Fast 3D Magnetotelluric Forward Modeling Based on a Regularization Approach in Spherical Coordinate

The earth’s curvature shouldn’t be ignored for the large-scale magnetotelluric (MT) modeling. A forward code of spherical coordinate MT method based on staggered grid finite difference method is developed. In addition, the forward problem is usually solved numerically, mostly based on iterative solution. Krylov subspace method is widely used to solve electromagnetic forward modeling problems in frequency domain. However, as the period increases, the convergence of these methods becomes significantly slower (especially for long-period MT data). We added a scaled grad-div operator to the double curl equation of the electric field to perform the nondivergence condition (CCGD). The results show that CCGD method is more efficient than other methods, which is an efficient method for our code.

3D Forward and Inversion Research on Frequency Domain Electromagnetic Tipper Data

The tipper is one of the main parameters of the electromagnetic exploration data in the frequency domain and has an important role and significant advantage in the exploration under complex geological conditions. In this paper, the 3D forward modeling is carried out based on the staggered grid finite difference method, and the analytical solution of the layered medium is used to verify its correctness; the L-BFGS algorithm is used to perform the 3D inversion of the tipper data. Through the three-dimensional complex model, the forward response characteristics and spatial distribution of the tipper are analyzed. The inversion results show that the inversion using the inversion data has less influence on static displacement and better inversion results.

An improved acoustic/elastic interface approach for 2D staggered grid finite-difference modeling of Rayleigh waves in the presence of surface topography

Rayleigh waves are generated along the free-surface and their propagation can be strongly influenced by surface topography. Therefore, a critical aspect of Rayleigh-wave study is to understand the wave-propagation behavior in the presence of surface topography. The acoustic/elastic boundary approach was incorporated into a ‘stair-case’ mesh for Rayleigh-wave modeling in the presence of surface topography. The conventional acoustic/elastic boundary approach, however, fails to completely fulfill the free-surface boundary condition at some grid nodes, leading to unphysical simulation results. To solve this problem, we develop a stable acoustic/elastic boundary approach that fully satisfies the free-surface boundary condition by adjusting the free boundary conditions at problematic nodes in the mesh. The accuracy is demonstrated by modeling tests in two-dimensional isotropic models in the presence of surface topography. And compared to the improved vacuum formulation, the proposed method takes less time. Our approach optimizes the free-surface boundary condition and enables the high-precision numerical simulation of Rayleigh waves.

Higher-order thin layer method (HTLM) based wavefield modeling approach

Traditional free vibration-based forward modeling approaches provide only the theoretical modal solution that lacks information regarding modal amplitudes. Moreover, they ignore the influence of source-receiver layout on the wavefield. Numerical wavefield modeling approaches can generate a dispersion spectrum mimicking the actual field test; however, numerical methods are not efficient in terms of computational time and produce numerical artifacts. Furthermore, they can be contaminated with leaky waves for soil strata with a high Poisson's ratio. This study presents a new wavefield modeling approach for developing full dispersion spectra of stratified media incorporating source-receiver configuration. The proposed method calculates the theoretical modal amplitude using the higher-order thin layer method and determines the vertical displacement responses at specified virtual receiver positions in the frequency-space domain. Afterward, it searches for spectral values at each frequency over a broad range of trial phase velocities to create a full dispersion spectrum. The accuracy and computational efficacy of the present method are demonstrated on different soil profiles adopted from literature, including normally dispersive media, irregularly dispersive media with high-velocity contrast, and a crustal level model. Results indicate the proposed approach can clearly distinguish the relative spectral amplitude of each mode and dominant dispersion trend. Calculated seismograms and dispersion images are compared with the existing 2D staggered grid finite difference method. Although both methods yield similar results, the current method is at least two orders (x102) faster than the staggered grid finite difference method. Unlike numerical approaches, the proposed method does not require additional computational resources to model materials with a high Poisson's ratio. Furthermore, it minimizes the contamination of leaky waves in the dispersion spectrum when modeling a high Poisson's ratio. The proposed method does not produce any numerical artifacts and can be used to generate accurate synthetic seismograms as well as dispersion spectra. In addition, the proposed wavefield modeling approach can be effectively used to develop an accurate and faster inversion scheme by including the entire dispersion spectrum.

Simulation of Heat Storage and Heat Regeneration in Phase Change Material

The present study explores numerically the energy storage and energy regeneration during Melting and Solidification processes in Phase Change Materials (PCM) used in Latent Heat Thermal Energy Storage (LHTES) systems. Transient two-dimensional (2-D) conduction heat transfer equations with phase change have been solved utilizing the Explicit Finite Difference Method (FDM) and Grid Generation technique. A Fortran computer program was built to solve the problem. The study included four different Paraffin's. The effects of container geometrical shape, which included cylindrical and square sections of the same volume and heat transfer area, the container volume or mass of PCM, variation of mass flow rate of heat transfer fluid (HTF), and temperatures difference between PCM and HTF were all investigated. Results showed that the PCMs in a cylindrical container melt and solidify quicker than the square container. The increase in mass flow rate and/or temperature difference decreases the time required for complete phase change. Paraffin's solidify quicker than they melt and store more energy than they release

Finite-difference modeling with topography using 3D viscoelastic parameter-modified free-surface condition

The free-surface boundary condition is a crucial aspect in the numerical modeling of (visco)elastic wave equations, especially when using a finite-difference (FD) method in the presence of surface topography. The parameter-modified method is a widely used approach to solve this problem. In this regard, the adaptive Poisson’s ratio parameter-modified method has proven to be effective in accurately simulating seismic surface waves within the FD discretization framework. Based on the equivalent medium theory, vacuum approximation, and mathematical limit, we develop a viscoelastic parameter-modified (VPM) method for the implementation of the free-surface boundary condition in the 3D viscoelastic wave equation. Our approach modifies the viscoelastic constitutive relation and density at the free surface and provides a formulation in terms of displacement and stress. We determine that our VPM method is more general than the viscoelastic stress-image method as it includes the latter as a limit case when the Poisson’s ratio equals zero. The presented free-surface method is represented implicitly within the FD grid, and we provide implementation details when simulating surface waves with topography in the standard staggered-grid FD scheme. We support our theory’s feasibility and accuracy through numerical examples.

NUMERICAL ESTIMATION OF SOIL DENSITY BY THE METHOD OF FINITE DIFFERENCE GRIDS IN MATHEMATICAL MODELING OF COMPACTION OF SUBSIBLE SOILS BY DEEP EXPLOSIONS

Link for citation: Tarasenko E.O. Numerical estimation of soil density by the method of finite difference grids in mathematical modeling of compaction of subsible soils by deep explosions. Bulletin of the Tomsk Polytechnic University. Geo Аssets Engineering, 2023, vol. 334, no. 5, рр.103-108. In Rus. Relevance. In order to ensure the duration of operation and strength of buildings and their expected design and construction on subsiding loess, it is necessary to carry out a numerical grade of the density characteristics of soils. Subsiding loesses are widespread in the world. They can be observed in Russia, China, the countries of Central Asia, etc. On the territory of Russia, their occurrence is noted at least at 17 % of the entire territory of the country. Almost ubiquitous distribution of loess soils occurs on the territory of the North Caucasus region (80–85 % of area). There is a large capacity of them (up to 50 m). This work is devoted to numerical modeling of the grade of density characteristics of subsiding soils using the finite difference grid method in the framework of a mathematical model of loess compaction by deep explosions. The use of the instrument of computational mathematics made it possible to solve the problem of estimating the compacted subsidence soil density depending on the depth of the explosive charge. The main aim of the research is to conduct a numerical assessment of soil density within the framework of mathematical modeling of compaction by the method of deep explosions of subsiding loess soils. Methods: application of the instrument of computational mathematics, differential equations, methods and results of solid state physics, engineering geology, soil mechanics. Experimental technique. Numerical modeling of the assessment of the increase in the compacted soil by deep explosions based on the grid method was carried out. Results. Based on the grid method using an explicit finite-difference scheme, discrete linear dynamic systems were constructed for numerical calculation of compacted subsidence soil density by deep explosions. The paper describs the layer-by-layer approximation to the solution of a differential equation in partial derivatives with given initial and boundary conditions, which mathematically describes the technological process of compaction of subsiding soil for the cases of surface ejection and camouflage compaction. The stability condition for the obtained finite-difference schemes is indicated. A computational experiment was carried out. It demonstrated the adequacy of the proposed method for estimating the compacted soil density to experimental data. Findings. The constructed numerical solutions of problems within the framework of mathematical modeling of the compaction of subsiding soils by deep explosions indicate the possibility of their practical application as calculation methods to reduce the cost of carrying out production work on the compaction of loess. The methods proposed in the work for determining soil density can be implemented in a number of areas of applied research in soil mechanics.

P- and S-wave energy current density vectors dot product imaging condition of source time-reversal imaging

SUMMARY Source time-reversal imaging (TRI) based on decoupled elastic wave equation can utilize vector P- and S-wave time differences and achieve high-precision source location in complex geological models. The imaging condition is critical for TRI. However, because of the orthogonally polarized properties of P and S waves, traditional vector dot product imaging condition directly applied to TRI will decrease the effective imaging values. In contrast, the energy current density vectors of P and S waves represent the propagation directions of the wavefields and are almost parallel. Their dot product can result in the maximum imaging energy. Based on this principle, we propose a P- and S-wave energy current density vectors dot product imaging condition (PSEDPIC), which uses the propagation direction information of P and S waves at the source point to suppress imaging artefacts generated by waves with inconsistent propagation directions. Numerical tests reveal that PSEDPIC can (1) reduce the image artefacts, (2) improve the imaging spatial resolution and (3) allow a shallower imaging region. In addition, if the numerical simulation algorithm used in TRI can reconstruct the seismic wavefield accurately in the presence of surface topography, the impact of an observation system with elevation differences on imaging can be eliminated automatically. For this reason, we use the curvilinear grid finite-difference method to directly reconstruct the wavefield in TRI to solve the problem of data elevation correction. The test results of 3-D synthetic and field data for microseismic monitoring demonstrate the effectiveness of the proposed method.