Transient wave propagation in a one-dimensional poroelastic column

Summary Biot's theory of porous media governs the wave propagation in a porous, elastic solid infiltrated with fluid. In this theory, a second compressional wave, known as the slow wave, has been identified. In this paper, Biot's theory is applied to a one-dimensional continuum. Despite the simplicity of the geometry, an exact solution of the full model, and a detailed analysis of the phenomenon, so far have not been achieved. In the present approach, an analytical solution in the Laplace transform domain is obtained showing clearly two compressional waves. For the special case of an inviscid fluid, a closed form exact solution in time domain is obtained using an analytical inverse Laplace transform. For the general case of a viscous fluid, solution in time domain is evaluated using the Convolution Quadrature Method of Lubich. Of all the inverse methods previously investigated, it seems that only the method of Lubich is efficies and stable enough to handle the highly transient cases such as impact and step loadings. Using properties of three widely different real materials, the wave propagating behavior, in terms of stress, pore pressure, displacement, and flux, are examined. Of most interest is the identification of second compressional wave and its sensitivity of material parameters.     

A new visco- and elastodynamic time domain Boundary Element formulation

The usual time domain Boundary Element Method (BEM) contains fundamental solutions which are convoluted with time-dependent boundary data and integrated over the boundary surface. If the fundamental solution is known, e.g., in Elastodynamics, the temporal convolution can be performed analytically when the boundary data are approximated by polynomial shape functions in time and in the boundary elements. This formulation is well known, but the resulting time-stepping BEM procedure produces instabilities and high numerical damping, when the time step size is chosen too small and too large, respectively. Moreover, in case of viscoelastic or poroelastic domains, the fundamental solution is known only in the frequency domain such that the time history of a response can only be obtained by an inverse transformation of the frequency domain results. Here, a new approach for the evaluation of the convolution integrals, the so-called “Operational Quadrature Methods” developed by LUBICH, is presented. In this formulation, the convolution integral is numerically approximated by a quadrature formula whose weights are determined by the Laplace transform of the fundamental solution and a linear multistep method. Hence, the frequency domain fundamental solution can be used without the need of an inverse transformation. Therefore, the extension to viscoelastic problems succeeds using the elastic-viscoelastic correspondence principle.     

Existence of longitudinal waves in anisotropic poroelastic solids

In anisotropic fluid-saturated porous solids, four waves can propagate along a general phase direction. However, solid particles in different waves may not vibrate in mutually orthogonal directions. In the propagation of each of these waves, the displacement of pore–fluid particles may not be parallel to that of solid particles. The polarization for a wave is the direction of aggregate displacement of the particles of the two constituents of a porous aggregate. These polarizations, for different waves, are not mutually orthogonal. Out of the four waves in anisotropic poroelastic medium, two are termed as quasi-longitudinal waves. The prefix ‘quasi’ refers to their polarization being nearly, but not exactly, parallel to the direction of propagation. The existence of purely longitudinal waves in an anisotropic poroelastic medium is ensured by the stationary characters of two expressions. These expressions involve the elastic (stiffness and coupling) coefficients of a porous aggregate and the components of phase direction. Necessary and sufficient conditions for the existence of longitudinal waves are discussed for different anisotropic symmetries. Conditions are also discussed for the existence of the apparent longitudinal waves, i.e., the propagation of wave motion with the particle displacement parallel to the ray direction instead of the phase direction. A graphical solution of a numerical example is shown to check the existence of these apparent longitudinal waves for general directions of phase propagation.     

Transient dynamic analysis of a fluid-saturated porous gradient elastic column

The dynamic response of a fluid-saturated porous gradient elastic column to a transient disturbance is determined analytically and numerically. The basic dynamic theory of a fluid-saturated poroelastic medium due to Biot is modified by replacing the classical linear elastic model of the solid skeleton by the simple gradient elastic model of Mindlin with just one elastic constant (internal length scale) in addition to the classical ones. Thus, the new theory, which is presently restricted to the one-dimensional case, can take into account the microstructural effects of the solid skeleton. After the establishment of appropriate boundary and initial conditions, the one-dimensional dynamic column problem is solved analytically with the aid of the Laplace transform with respect to time. The time domain response is finally obtained by a numerical inversion of the transformed solution. The effect of the solid microstructure on the response is assessed and discussed.     

An accelerated symmetric time-domain boundary element formulation for elasticity

Wave propagation phenomena occur often in semi-infinite regions. It is well known that such problems can be handled well with the boundary element method (BEM). However, it is also known that the BEM, with its dense matrices, becomes prohibitive with respect to storage and computing time. Focusing on wave propagation problems, where a formulation in time domain is preferable, the mentioned limit of the method becomes evident. Several approaches, amongst them the adaptive cross approximation (ACA), have been developed in order to overcome these drawbacks mainly for elliptic problems.The present work focuses on time dependent elastic problems, which are indeed not elliptic. The application of the presented fast boundary element formulation on such problems is enabled by introducing the well known Convolution Quadrature Method (CQM) as time stepping scheme. Thus, the solution of the time dependent problem ends up in the solution of a system of decoupled Laplace domain problems. This detour is worth since the resulting problems are again elliptic and, therefore, the ACA can be used in its standard fashion.The main advantage of this approach of accelerating a time dependent BEM is that it can be easily applied to other fundamental solutions as, e.g., visco- or poroelasticity.     

Nonlinear surface waves on an elastic solid

The effect of quadratic elastic nonlinearity on the propagation of surface Rayleigh waves on an isotropic elastic solid is examined. Using the method of multiple scales an approximate solution is obtained which is uniformly valid in both spatial directions as well as in time. An arbitrary wave profile is considered and an integro-differential equation is derived for the Fourier transform of the displacement on the boundary. In the case of a quasi-monochromatic wave explicit expressions are derived for the variations of the amplitudes of the fundamental and second and third harmonics along the boundary.     

An explanation for the anomalous ultrasonic slow wave in underwater sand

An explanation is given for the propagation time of the well-known anomalous ultrasonic slow wave observed in water-saturated sand using a three-layer elastic model. The rapid increase of elastic properties of sand with depth causes conversion of near-grazing underwater acoustic waves into multiple coupled shear and compressional waves.     

Initial conditions contribution in a BEM formulation based on the convolution quadrature method

This work presents a two‐dimensional boundary element method (BEM) formulation for the analysis of scalar wave propagation problems. The formulation is based on the so‐called convolution quadrature method (CQM) by means of which the convolution integral, presented in time‐domain BEM formulations, is numerically substituted by a quadrature formula, whose weights are computed using the Laplace transform of the fundamental solution and a linear multistep method. This BEM formulation was initially developed for scalar wave propagation problems with null initial conditions. In order to overcome this limitation, this work presents a general procedure that enables one to take into account non‐homogeneous initial conditions, after replacing the initial conditions by equivalent pseudo‐forces. The numerical results included in this work show the accuracy of the proposed BEM formulation and its applicability to such kind of analysis. Copyright © 2006 John Wiley & Sons, Ltd.     

The interaction of two inclined cracks with dynamic stress wave loadings

To gain insight into the phenomenon of the interaction of stress waves with material defects and the linkage of two cracks, the transient response of two semi-infinite inclined cracks subjected to dynamic loading is examined. The solutions are obtained by the linear superposition of fundamental solutions in the Laplace transform domain. The fundamental solution is the exponentially distributed traction on crack faces proposed by Tsai and Ma [1]. The exact closed form solutions of stress intensity factor histories for these two inclined cracks subjected to incident plane waves and diffracted waves are obtained explicitly. These solutions are valid for the time interval from initial loading until the first wave scattered at one crack tip returns to the same crack tip after being diffracted by another crack tip. The result shows that the contribution of diffracted waves to stress intensity factors is much less than the incident waves. The probable crack propagation direction is predicted from the fracture criterion of maximum circumferential tensile stress. The linkage of these two cracks is also investigated in detail.     

An inverse-source problem for maximization of pore-fluid oscillation within poroelastic formations

This paper discusses a mathematical and numerical modeling approach for identification of an unknown optimal loading time signal of a wave source, atop the ground surface, that can maximize the relative wave motion of a single-phase pore fluid within fluid-saturated porous permeable (poroelastic) rock formations, surrounded by non-permeable semi-infinite elastic solid rock formations, in a one-dimensional setting. The motivation stems from a set of field observations, following seismic events and vibrational tests, suggesting that shaking an oil reservoir is likely to improve oil production rates. This maximization problem is cast into an inverse-source problem, seeking an optimal loading signal that minimizes an objective functional – the reciprocal of kinetic energy in terms of relative pore-fluid wave motion within target poroelastic layers. We use the finite element method to obtain the solution of the governing wave physics of a multi-layered system, where the wave equations for the target poroelastic layers and the elastic wave equation for the surrounding non-permeable layers are coupled with each other. We use a partial-differential-equation-constrained-optimization framework (a state-adjoint-control problem approach) to tackle the minimization problem. The numerical results show that the numerical optimizer recovers optimal loading signals, whose dominant frequencies correspond to amplification frequencies, which can also be obtained by a frequency sweep, leading to larger amplitudes of relative pore-fluid wave motion within the target hydrocarbon formation than other signals.     

Potential discontinuity simulations with the numerical Green’s function in steady and transient problems

This paper applies the numerical Green’s function (NGF) boundary element formulation (BEM) first in standard form to solve the Laplace equation and then, coupled to the operational quadrature method (OQM), to solve time domain problems (TD-BEM). Both involve the analysis of potential discontinuities in the respective scalar model simulation. The implementation of the associated Green’s function acting as the fundamental solution is advantageous since element discretization of actual discontinuity surfaces are no longer required. In the OQM the convolution integral is substituted by a quadrature formula, whose weights are computed using the fundamental solution in the Laplace domain, producing the direct solution to the problem in the time domain. Applications of the NGF to problems involving the Laplace equation and its transient counterpart are presented for two-dimensional potential flow examples, confirming that the formulation is stable and accurate.     

Scalar wave propagation in 2D: a BEM formulation based on the operational quadrature method

This work presents a boundary element method formulation for the analysis of scalar wave propagation problems. The formulation presented here employs the so-called operational quadrature method, by means of which the convolution integral, presented in time-domain BEM formulations, is substituted by a quadrature formula, whose weights are computed by using the Laplace transform of the fundamental solution and a linear multistep method. Two examples are presented at the end of the article with the aim of validating the formulation.     

Dynamic stress intensity factor (Mode I) of a penny-shaped crack in an infinite poroelastic solid

This paper considers the transient stress intensity factor (Mode I) of a penny-shaped crack in an infinite poroelastic solid. The crack surfaces are impermeable. By virtue of the integral transform methods, the poroelastodynamic mixed boundary value problems is formulated as a set of dual integral equations, which, in turn, are reduced to a Fredholm integral equation of the second kind in the Laplace transform domain. Time domain solutions are obtained by inverting Laplace domain solutions using a numerical scheme. A parametric study is presented to illustrate the influence of poroelastic material parameters on the transient stress intensity. The results obtained reveal that the dynamic stress intensity factor of poroelastic medium is smaller than that of elastic medium and the poroelastic medium with a small value of the potential of diffusivity shows higher value of the dynamic stress intensity factor.     

Soil-structure interaction and wave propagation problems in 2D by a Duhamel integral based approach and the convolution quadrature method

In this paper, a new methodology for analyzing wave propagation problems, originally presented and checked by the authors for one-dimensional problems [18], is extended to plane strain elastodynamics. It is based on a Laplace domain boundary element formulation and Duhamel integrals in combination with the convolution quadrature method (CQM) [13], [14]. The CQM is a technique which approximates convolution integrals, in this case the Duhamel integrals, by a quadrature rule whose weights are determined by Laplace transformed fundamental solutions and a multi-step method. In order to investigate the accuracy and the stability of the proposed algorithm, some plane wave propagation and interaction problems are solved and the results are compared to analytical solutions and results from finite element calculations. Very good agreement is obtained. The results are very stable with respect to time step size. In the present work only multi-region boundary element analysis is discussed, but the presented technique can easily be extended to boundary element – finite element coupling as will be shown in subsequent publications.     

Dynamic stress concentration studies by boundary integrals and Laplace transform

The dynamic stress field and its concentrations around holes of arbitrary shape in infinitely extended bodies under plane stress or plane strain conditions are numerically determined. The material may be linear elastic or viscoelastic, while the dynamic load consists of plane compressional waves of harmonic or general transient nature. The method consists of applying the Laplace transform with respect to time to the governing equations of motion and formulating and solving the problem numerically in the transfomed domain by the boundary integral equation method. The stress field can then be obtaind by a numerical inversion of the trasformed solution. The correspondence principle is invoked for the case of viscoelastic material behavious. The method is simplified for the case of harmonic waves where no numerical inversion is involved.     

A regularized collocation boundary element method for linear poroelasticity

This article presents a collocation boundary element method for linear poroelasticity, based on the first boundary integral equation with only weakly singular kernels. This is possible due to a regularization of the strongly singular double layer operator, based on integration by parts, which has been applied to poroelastodynamics for the first time. For the time discretization the convolution quadrature method (CQM) is used, which only requires the Laplace transform of the fundamental solution. Furthermore, since linear poroelasticity couples a linear elastic with an acoustic material, the spatial regularization procedure applied here is adopted from linear elasticity and is performed in Laplace domain due to the before mentioned CQM. Finally, the spatial discretization is done via a collocation scheme. At the end, some numerical results are shown to validate the presented method with respect to different temporal and spatial discretizations.     

A finite element Laplace transform solution technique for the wave equation

A technique is described for the solution of the wave equation with time dependent boundary conditions. The finite element solution accompanied by the numerical Laplace inversion process seems to be an efficient procedure to treat such problems. The programming involved is straightforward in the sense that numerical Laplace inversion routines can be directly used as a time integration procedure after obtaining standard finite element differential equation solutions in the transformed domain. Some results are presented for one- and two- dimensional applications, such as wave propagation in longitudinal bars and wave propagation in harbours.     

基于快速多极子基本解方法(FMM-MFS)的弹性波二维散射模拟研究

针对弹性波二维散射问题,发展一种新的快速多极子基本解方法（FMM-MFS）。方法基于单层位势理论,通过在虚边界上设置膨胀波线源和剪切波线源以构造散射波场,从而避免了奇异性的处理和边界单元离散;结合快速多极子展开技术（FMM）,大幅度降低了计算量和存储量,突破了传统方法难以处理大规模散射问题的瓶颈。以全空间孔洞对P、SV波的二维散射为例,给出了具体求解步骤,并在个人计算机上实现了上百万自由度问题的快速精确计算。在方法效率和精度检验基础上,分别以单孔洞和随机孔洞群对平面波（P、SV波）的散射为例进行计算模拟,揭示了孔洞（群）周围弹性波散射的若干重要规律。     

Effects of fiber ellipticity and orientation on dynamic stress concentrations in porous fiber-reinforced composites

Interaction of time harmonic fast longitudinal and shear incident plane waves with an elliptical fiber embedded in a porous elastic matrix is studied. The novel features of Biot dynamic theory of poroelasticity along with the classical method of eigen-function expansion and the pertinent boundary conditions are employed to develop a closed form series solution involving Mathieu and modified Mathieu functions of complex arguments. The complications arising due to the non-orthogonality of angular Mathieu functions corresponding to distinct wave numbers in addition to the problems associated with appearance of additional angular dependent terms in the boundary conditions are all avoided by expansion of the angular Mathieu functions in terms of transcendental functions and subsequent integration, leading to a linear set of independent equations in terms of the unknown scattering coefficients. A MATHEMATICA code is developed for computing the Mathieu functions in terms of complex Fourier coefficients which are themselves calculated by numerically solving appropriate sets of eigen-systems. The analytical results are illustrated with numerical examples in which an elastic fiber of elliptic cross section is insonified by a plane fast compressional or shear wave at normal incidence. The effects of fiber cross sectional ellipticity, angle of incidence (fiber two-dimensional orientation), and incident wave polarization (P, SV, SH) on dynamic stress concentrations are studied in a relatively wide frequency range. Limiting cases are considered and fair agreements with well-known solutions are established.