Parallel Simulation of 3D Wave Propagation by Domain Decomposition

Elastic wave propagation in a segmented X-ray computed tomography model of a concrete specimen

The problem of tsunami warning and the use of the grid-characteristic method developed for studying wave processes in heterogeneous media is addressed. A multilayer geological model with curved boundaries and contrasting elastic parameters is considered. Various seismic survey methods, as well as different geometries of the boundaries, different values of rock densities, and various longitudinal and shear wave velocities, are examined. The boundary value problem is numerically solved jointly for the elastic and acoustic wave equations (seismic wave propagation in rocks and water layer, respectively). We use the numerical grid-characteristic method with combined structured curvilinear and regular computational grids. The wave velocity fields (velocity as derivative of displacement) and stress fields are calculated and analyzed. Synthetic seismograms are calculated for the coastal zone. It is found that some specific features of seismograms and some wave modes can be used for tsunami warning. The optimal parameters of seismic surveying are determined. Numerical simulation of wave propagation effects is used for these purposes. The proposed numerical grid-characteristic method can be used in forward numerical modeling in order to interpret seismograms recorded in the coastal zone.

Physical mechanisms implied in the propagation of ultrasound waves in trabecular bones are not yet clearly understood. The relative role of scattering and absorption in the attenuation process is unknown. To dissociate these two effects, a simulation of wave propagation in 3D volumes of trabecular bone is performed, and compared to experimental results. The simulation algorithm accounts for scattering and volumetric absorption into both the saturating fluid and bone. The simulation software uses a finite difference approach based on the Virieux numerical scheme. An incident plane wave is propagated on a volume of bone of approximately 5×5×8 mm3. These volumes were reconstructed from high resolution micro-CT experiments. Several types of simulations are computed: bone described as a fluid matrix saturated by fluid (i.e., neglecting the shear wave propagation into the matrix) and bone described as a solid matrix saturated by fluid, with or without absorption. From the transmitted signals, attenuation and speed of sound through the specimens are calculated and compared to experimental values obtained on the same specimens. Finally, the contribution of scattering to the total attenuation is estimated, and the importance of mode conversions as well as absorption.

Accurate simulations of seismic wave propagation in complex geological structures with great and rapid variations of topography are of primary interest for environmental and industrial applications. Unfortunately, difficulties arise for such complex environments, due essentially to the existence of shadow zones, head waves, diffractions and edge effects. Usually, methods and codes are tested against "validated" ones, but one might wonder which method/code ultimately approaches the "real" solution. An original approach for seismics is to compare synthetic seismic data to controlled laboratory data for a well-described configuration, in order to analyze the respective limitations of each method/code. This is one of the objectives of the BENCHIE project, which brings together laboratories in France, Norway and Russia. In this presentation we will present some preliminary results provided by both laboratory experiments conducted in a tank and numerical simulations of wave propagation. The laboratory data have been obtained by zero-offset acquisitions at different ultrasonic frequencies on the Marseille model which is made up of anticlines, fault and truncated pyramid. The numerical results have been obtained by two methods: the Spectral-Element Method and the Tip-Wave Superposition Method.

In numerical simulation of soil media, the problem of wave propagation is of a great importance. Proposed in this paper is a new formulation of infinite elements which can be used in numerical simulation of wave propagation problems in infinite domains. In numerical simulation of wave propagation, the finite elements are used to model the near field, whereas the infinite elements are used to represent the behaviour of the far field. Formulations and various implementation aspects of the proposed infinite elements are illustrated. The accuracy and efficiency of the proposed approach is considered by comparing the obtained results with analytical and other numerical results. For better explanation, a couple of examples were analyzed such as one dimensional wave propagation problems arising from the Heaviside step function and impulse functions. In order to get a more complete insight, two dimensional wave propagations in soil medium were considered and the results are presented accordingly. Finally, a soil layer subjected to seismic excitation was analysed. The new infinite element was developed using the User Programmable Features of the ANSYS software, which enables creating new elements within the robust ANSYS core algorithm. The main advantage of the proposed infinite elements is that they can be used directly within a finite framework with minor modifications such as the Jacobian matrix and added absorbing properties. The performed analysis using the infinite elements shows very promising results and provides a good tool for simulation of boundary conditions.

Although the earth is 3-dimensional (3-D), numerical simulations of wave propagation through laterally heterogeneous media are easier to formulate and more practical to use in 2-D. In this thesis, schemes to model seismic wave propagation through laterally varying structures with 2-D numerical algorithms are developed and applied to earthquake and explosion problems. In Chapter 1, 2-D source expressions that have the same radiation patterns as their 3-D counterparts are derived which can rectify the following three problems: the use of 2-D simulations generally results in line source tails on what would be impulsive arrivals in 3-D, 1/√R rather than 1/R amplitude decay for body waves, and no decay rather than 1/√R amplitude decay for surface waves. Because this technique approximately transforms waves from a cartesian 2-D grid to a cylindrically symmetric 3-D world, slightly anisotropic geometrical spreading in 2-D better approximates isotropic spreading in 3-D than simple isotropic spreading in 2-D does. In Section 1.7, a correction to the explosive source expression reduces energy traveling vertically out of the source region, but leaves unchanged the energy traveling laterally out of the source region. In some cases, this correction will significantly improve the results of using a 2-D grid to simulate elastic wave propagation from an explosive point source. In Chapter 2, synthetic seismograms are constructed for the strong motions of the 1968 Borrego Mountain earthquake recorded at EI Centro. A good fit to the data results from using the laterally varying model determined by a detailed refraction survey and the source parameters determined by teleseismic waveform modeling. Shallow faulting is no longer necessary to explain the long-period surface-wave development. Synthetic seismograms calculated for the 1971 San Fernando earthquake show strong effects due to lateral variation in sediment thickness in the San Fernando valley and the Los Angeles basin. Using previously determined basin structure and teleseismically determined source parameters, two-dimensional SH and P-SV finite difference calculations can reproduce the amplitude and duration of the strong motion velocities recorded across the basins in Los Angeles in the period range from 1 to 10 seconds. The edges of basins nearest the seismic source show ground motion amplification up to a factor of three, and tend to convert direct shear waves into Love and Rayleigh waves that travel within the basins. The computed motions are sensitive to the mechanism and location of earthquakes. A strike-slip earthquake on the Newport-Inglewood fault zone, for example, would produce different patterns of peak velocity and duration of shaking across the San Fernando and Los Angeles basins. In Chapter 3, the effect of shallow station structure and lateral velocity variation are investigated for records of the Amchitka explosion Milrow. The differences between the Meuller-Murphy, Heimberger-Hadley, and von Seggern-Blandford reduced displacement potential (RDP) source representations are small compared to the differences between using various possible velocity structures. Creager and Jordan (1986) propose that penetration of subducting slabs under the Kurile Islands and other subduction zones to depths of at least 1000 km is necessary to explain the t ravel time anomalies of deep earthquakes. Such penetration would also affect the amplitudes and waveforms of the body waves from these earthquakes. In Chapter 4, synthetic seismograms appropriate for a record section in a plane perpendicular to the strike of the slab are presented using a coupled finite-difference and Kirchhoff method. An inferred shear-wave version of the compressional-wave velocity structure of Creager and Jordan (1986) produces an amplitude decrease up to a factor of four and waveform broadening up to 20 seconds for SH arrivals with a take-off angle pointing straight down t he slab. Slabs that extend only 300 km below the earthquake but are half as thick and twice as anomalously fast as Creager and Jordan 's (1986) velocity model will roughly preserve the travel time variation pattern, and show less waveform broadening, but produce first arrivals that are emergent. Slabs that become thicker with depth show less waveform broadening. Reconciliation of the amplitude, waveform distort ion, and timing of body waves from deep events is necessary to understand the geometry of slabs near and below the 6.50 km discontinuity.

Numerical methods are widely used in seismic exploration to simulate wave propagation; however, the algorithms are based on various assumptions. The accuracy of numerical simulations is of particular interest in the case of realistic geologic setups. The direct comparison of numerical results can have limitations, and an alternative approach can be the comparison of synthetic results with experimental data, obtained for a small-scale physical model in laboratory conditions. Laboratory experiments are repeatable and provide high-quality data for a known configuration. We have developed a possible workflow to adapt the numerical simulations and the laboratory experiments to each other, such that the two can be easily compared with high accuracy. The model is immersed in a water tank, and a conventional pulse-echo technique is used to collect the reflection data in zero-offset and offset configurations. We use a spectral-element method for the numerical modeling. The model geometry is implemented using a nonstructured mesh, and the computational cost can be optimized using larger elements and higher-order basis functions. The real source transducer characteristics are implemented based on a new approach: laboratory characterization of the impulse response, followed by an inversion step to obtain a numerically equivalent source. The comparison of the zero-offset synthetic and laboratory results reveals an excellent fit in terms of arrival time, phase, and amplitude. Minor amplitude mismatches may be attributed to the noise recorded in the laboratory data and to the possible inaccuracy of the proposed source implementation. Comparison of the simulated and laboratory offset traces also exhibits a good fit in general, but with significantly less accuracy for some arrivals than in the zero-offset case. This can be mainly attributed to the inaccuracies of the transducer positions during the laboratory measurements combined with the strong topography of the model.

Context. In stratified atmospheres, acoustic waves can only propagate if their frequency is higher than the cutoff value. The determination of the cutoff frequency is fundamental for several topics in solar physics, such as evaluating the contribution of the acoustic waves to the chromospheric heating or the application of seismic techniques. However, different theories provide different cutoff values. Aims. We developed an alternative method to derive the cutoff frequency in several standard solar models, including various quiet-Sun and umbral atmospheres. The effects of magnetic field and radiative losses on the cutoff are examined. Methods. We performed numerical simulations of wave propagation in the solar atmosphere using the code MANCHA. The cutoff frequency is determined from the inspection of phase-difference spectra computed between the velocity signal at two atmospheric heights. The process is performed by choosing pairs of heights across all the layers between the photosphere and the chromosphere to derive the vertical stratification of the cutoff in the solar models. Result. The cutoff frequency predicted by the theoretical calculations departs significantly from the measurements obtained from the numerical simulations. In quiet-Sun atmospheres, the cutoff shows a strong dependence on the magnetic field for adiabatic wave propagation. When radiative losses are taken into account, the cutoff frequency is greatly reduced and the variation of the cutoff with the strength of the magnetic field is lower. The effect of the radiative losses in the cutoff is necessary to understand recent quiet-Sun and sunspot observations. In the presence of inclined magnetic fields, our numerical calculations confirm that the cutoff frequency is reduced as a result of the reduced gravity experienced by waves that propagate along field lines. An additional reduction is also found in regions with significant changes in the temperature, which is due to the lower temperature gradient along the path of field-guided waves. Conclusions. Our results show solid evidence that the cutoff frequency in the solar atmosphere is stratified. The cutoff values are not correctly captured by theoretical estimates. In addition, most of the widely used analytical cutoff formulae neglect the effect of magnetic fields and radiative losses, whose role is critical for determining the evanescent or propagating nature of the waves.

Numerical simulation of wave propagation is performed through 31 volumes of trabecular bone. These volumes were reconstructed from high synchrotron microtomography experiments and were used as the input geometry in a three-dimensional (3D) finite-difference simulation tool developed in our laboratory. The simulation algorithm accounts for propagation in both the saturating fluid and bone but absorption is not taken into account. Numerical predictions are consistent with experimental observations in trabecular bones : linear frequency dependence of attenuation, quasi-linear increase of attenuation and speed of sound with the bone volume fraction, negative phase velocity dispersion in most of the specimens, propagation of fast and slow waves depending on the orientation of the trabecular network with respect to the direction of propagation of the ultrasound. Moreover, the predicted attenuation is in close agreement with the experimental data obtained for the same specimens. Coupling numerical simulation with real 3D bone microarchitecture provides a powerful tool to investigate the physics of ultrasound propagation in trabecular structures.

SummaryBiot theory predicts wave velocities in a saturated granular medium using the pore geometry, viscosity, densities, and elastic moduli of the solid skeleton and pore fluid, neglecting the interaction between constituent particles and local flow, which becomes essential as the wavelength decreases. Here, a hydro‐micromechanical model, for direct numerical simulations of wave propagation in saturated granular media, is implemented by two‐way coupling the lattice Boltzmann method (LBM) and the discrete element method (DEM), which resolve the pore‐scale hydrodynamics and intergranular behavior, respectively. The coupling scheme is benchmarked with the terminal velocity of a single sphere settling in a fluid. In order to mimic a small amplitude pressure wave entering a saturated granular medium, an oscillating pressure boundary on the fluid is implemented and benchmarked with the one‐dimensional wave equation. The effects of input waveforms and frequencies on the dispersion relations in 3D saturated poroelastic media are investigated with granular face‐centered‐cubic crystals. Finally, the pressure and shear wave velocities predicted by the numerical model at various effective confining pressures are found to be in excellent agreement with Biot analytical solutions, including his prediction for slow compressional waves.

We propose a domain decomposition method based on a combination of the finite element and finite difference methods for numerical simulation of elastic wave propagation. It is not limited by model geometries and easy to be applied to seismic wave modelling with surface topography. It is also convenient to be implemented by parallel computational programs. The results of numerical experiments validate the accuracy of the algorithm.

Ultrasound measurements are used for monitoring the structural health of concrete. The elastic waves propagating in concrete structures are attenuated by scattering and intrinsic attenuation. Both quantities have to be determined for a detailed material characterization and for realistic numerical simulations of wave propagation in such structures. We present in this paper a numerical approach to determine the intrinsic attenuation of concrete. For this purpose, a realistic digital three-dimensional concrete specimen is created that accounts for correct volume fractions of aggregates of different sizes embedded into a mortar matrix. Using a viscoelastic finite difference scheme, we perform simulations of the ultrasound wave field in these digital models with different attenuation levels. The numerical observations are compared with results from ultrasonic measurements on a sample of saturated concrete to estimate the intrinsic attenuation. Embedded sensors are used for those long-term laboratory experiments. For this saturated concrete, we obtain Q=130 in a frequency band around 60 kHz.

This paper deals with the development of a computational code for the numerical simulation of wave propagation through domains with a complex geometry consisting in both solids and moving fluids. The emphasis is on the numerical simulation of ultrasonic flowmeters (UFMs) by modeling the wave propagation in solids with the equations of linear elasticity (ELE) and in fluids with the linearized Euler equations (LEEs). This approach requires high performance computing because of the high number of degrees of freedom and the long propagation distances. Therefore, the numerical method should be chosen with care. In order to minimize the numerical dissipation which may occur in this kind of configuration, the numerical method employed here is the nodal discontinuous Galerkin (DG) method. Also, this method is well suited for parallel computing. To speed up the code, almost all the computational stages have been implemented to run on graphical processing unit (GPU) by using the compute unified device architecture (CUDA) programming model from NVIDIA. This approach has been validated and then used for the two-dimensional simulation of gas UFMs. The large contrast of acoustic impedance characteristic to gas UFMs makes their simulation a real challenge.

In seismic numerical simulations of wave propagation, it is very important for us to consider surface topography and attenuation, which both have large effects (e.g., wave diffractions, conversion, amplitude/phase change) on seismic imaging and inversion. An irregular free surface provides significant information for interpreting the characteristics of seismic wave propagation in areas with rugged or rapidly varying topography, and viscoelastic media are a better representation of the earth’s properties than acoustic/elastic media. In this study, we develop an approach for seismic wavefield simulation in 2D viscoelastic isotropic media with an irregular free surface. Based on the boundary-conforming grid method, the 2D time-domain second-order viscoelastic isotropic equations and irregular free surface boundary conditions are transferred from a Cartesian coordinate system to a curvilinear coordinate system. Finite difference operators with second-order accuracy are applied to discretize the viscoelastic wave equations and the irregular free surface in the curvilinear coordinate system. In addition, we select the convolutional perfectly matched layer boundary condition in order to effectively suppress artificial reflections from the edges of the model. The snapshot and seismogram results from numerical tests show that our algorithm successfully simulates seismic wavefields (e.g., P-wave, Rayleigh wave and converted waves) in viscoelastic isotropic media with an irregular free surface.

Numerical simulation of wave propagation by using a hybrid method with an arbitrary order accuracy in both spatial and temporal approximations