Absorbing boundary conditions for dynamic analysis of fluid-saturated porous media

Based on Biot's two-phase mixture theory and the paraxial approximation, the absorbing boundary conditions in the time domain for u-w, u-U and u-p formulations are presented for the dynamic analysis of fluid-saturated porous media. These absorbing boundaries are equivalent to the viscous boundaries in the fundamental mode. The expressions for the energy ratio and reflection coefficient between the reflected and incident waves along the absorbing boundary are given. The numerical results show that the proposed absorbing boundaries can greatly suppress the spuriously reflected wave and model the far field of the system. These results also dynamic analysis of infinite fluid-saturated porous media. for the transient dynamic analysis of infinite fluid-saturated porous media.     

General absorbing boundary conditions for dynamic analysis of fluid-saturated porous media

General absorbing boundary conditions based on Biot's two-phase mixture theory and paraxial approximation is presented for the dynamic analysis of fluid-saturated porous media with isotropic, transverse isotropic, and anisotropic properties. For the last two cases, the equivalent Lame's constants, under conditions of uniqueness, are introduced to facilitate the analytical solutions. The numerical results show that the proposed absorbing boundary can greatly suppress spuriously reflected waves and efficiently model the far field of the system with sufficient accuracy.     

Comparison of Galerkin-type discretization techniques for two-phase flow in heterogeneous porous media

A comparison of Standard Galerkin, Petrov-Galerkin, and Fully-Upwind Galerkin methods for the simulation of two-phase flow in heterogeneous porous media is presented. On the basis of the coupled pressure-saturation equations, a generalized formulation for all three finite element methods is derived and analysed. For flow in homogeneous media, the Petrov-Galerkin method gives excellent results. But this method fails miserably for problems with heterogeneous media. This is because it is not able to capture correctly processes that take place at interfaces when, for instance, the capillary pressure-saturation relationship after Brooks and Corey is assumed. The Fully-Upwind Galerkin method is superior to the Petrov-Galerkin approach because it is able to give correct results for flow in homogeneous and heterogeneous media for the two models of van Genuchten and Brooks-Corey. The widely used formulation which is correct for the homogeneous case cannot be used for heterogeneous media. Instead the straightforward approach of gradp_c in combination with a chord-slope technique must be utilized.     

Application of a second-order paraxial boundary condition to problems of dynamics of circular foundations on a porous layered half-space

A half-space finite element and a consistent transmitting boundary in a cylindrical coordinate system are developed for analysis of rigid circular (or cylindrical) foundations in a water-saturated porous layered half-space. By means of second-order paraxial approximations of the exact dynamic stiffness for a half-space in plane-strain and antiplane-shear conditions, the corresponding approximation for general three-dimensional wave motion in a Cartesian coordinate system is obtained and transformed in terms of cylindrical coordinates. Using the paraxial approximations, the half-space finite element and consistent transmitting boundary are formulated in a cylindrical coordinate system. The development is verified by comparison of dynamic compliances of rigid circular foundations with available published results. Examination of the advantage of the paraxial condition vis-á-vis the fixed condition shows that the former achieves substantial gain in computational effort. The developed half-space finite element and transmitting boundary can be employed for accurate and effective analysis of foundation dynamics and soil–structure interaction in a porous layered half-space.     

A high-order time-domain transmitting boundary for cylindrical wave propagation problems in unbounded saturated poroelastic media

Based on the u–p formulation of Biot equation with an assumption of zero permeability coefficient, a high-order transmitting boundary is derived for cylindrical elastic wave propagation in infinite saturated porous media. By this transmitting boundary the total stresses on the truncated boundaries of a numerical model, such as a finite element model, are replaced by a set of spring, dashpot and mass elements, with some additionally introduced auxiliary degrees of freedom. The transmitting boundaries are incorporated into the DIANA SWANDYNE II program and an unconditionally stable implicit time integration algorithm is adopted. Despite the assumption made in the derivation of the transmitting boundary, numerical examples show that it can provide highly accurate results for cylindrical elastic wave propagation problems in infinite saturated porous medium in case the u–p formulation is applicable. Although the direct applications of the proposed transmitting boundary to general two dimensional wave problems in infinite saturated porous media are not highly accurate, acceptable accuracy can still be achieved by placing the transmitting boundary at relatively large distance from the wave source.     

基于Comsol的饱和多孔介质动力方程的数值模拟及应用

两相饱和多孔介质的动力响应问题在地震工程领域具有重要的研究意义,由于涉及到固相和液相的动力耦合,使得该问题的求解尤为复杂。本文利用Comsol在求解多场耦合问题上的优点,针对Biot饱和多孔介质u-U耦合形式下的波动方程特征,经过一系列微分算子运算和矩阵变换得到导数形式下的波动方程,基于Comsol Multiphysics提供的广义偏微分方程模式对变形后的波动方程进行求解,并把改进后的无限元边界应用到无限域动力问题的模拟中。通过与饱和多孔介质动力响应的解析解进行对比,验证模型求解技术的可行性和正确性,并在此基础上讨论饱和土地基中空沟隔振效果与饱和土体参数孔隙率、泊松比的关系。通过研究分析,可以为饱和土地基中空沟隔振设计提供一些有价值的参考。     

Gaussian beams in inhomogeneous media: A review

The paper outlines the most important results of the paraxial complex geometrical optics (CGO) in respect to Gaussian beams diffraction in the smooth inhomogeneous media and discusses interrelations between CGO and other asymptotic methods, which reduce the problem of Gaussian beam diffraction to the solution of ordinary differential equations, namely: (i) Babich’s method, which deals with the abridged parabolic equation and describes diffraction of the Gaussian beams; (ii) complex form of the dynamic ray tracing method, which generalizes paraxial ray approximation on Gaussian beams and (iii) paraxial WKB approximation by Pereverzev, which gives the results, quite close to those of Babich’s method. For Gaussian beams all the methods under consideration lead to the similar ordinary differential equations, which are complex-valued nonlinear Riccati equation and related system of complex-valued linear equations of paraxial ray approximation. It is pointed out that Babich’s method provides diffraction substantiation both for the paraxial CGO and for complex-valued dynamic ray tracing method. It is emphasized also that the latter two methods are conceptually equivalent to each other, operate with the equivalent equations and in fact are twins, though they differ by names. The paper illustrates abilities of the paraxial CGO method by two available analytical solutions: Gaussian beam diffraction in the homogeneous and in the lens-like media, and by the numerical example: Gaussian beam reflection from a plane-layered medium.     

Transmitting boundary for water-saturated transversely isotropic strata based on the u–U formulation

A new transmitting boundary in a cylindrical coordinate system has been developed for modeling the elastic waves radiating out to an infinite boundary in water-saturated transversely isotropic soil strata over a rigid bedrock. The saturated soil strata are assumed to consist of a porous material and modeled as a transversely isotropic two-phase medium, based on the u−U formulation. The newly developed transmitting boundary is combined with the finite-elements model of the near-field region, using the same u−U formulation, and applied to the study of the dynamics of a rigid circular foundation in porous isotropic or transversely isotropic layered strata, either fully or partly saturated with water. The verification and application examples give valuable insights into new and interesting aspects of the dynamic behavior of rigid circular foundations in fully or partly saturated two-phase ground in terms of permeability, transverse anisotropy, and ground-water table level.     

流体饱和多孔隙介质弹性波方程边界元解法研究

基于流体饱和多孔隙各向同性介质模型,本文首先推导了流体饱和多孔隙介质中弹性波传播的频率域系统动力方程及边界积分方程,然后给出了流体饱和多孔隙介质弹性波方程的基本解,最后,利用本文给出的边界元方法对流体饱和多孔隙各向同性介质中的弹性波传播进行了数值模拟.结果表明：不论是从固相位移,还是液相位移的地震合成记录都能看到明显的慢速P波,本文提出的流体饱和多孔隙介质弹性波边界元法是有效可行的.     

Propagation of seismic waves through liquefied soils

To predict the earthquake response of saturated porous media it is essential to correctly simulate the generation, redistribution, and dissipation of excess pore water pressure during and after earthquake shaking. To this end, a reliable numerical tool requires a dynamic, fully coupled formulation for solid–fluid interaction and a versatile constitutive model. Presented in this paper is a 3D finite element framework that has been developed and utilized for this purpose. The framework employs fully coupled dynamic field equations with a u–p–U formulation for simulation of pore fluid and solid skeleton interaction and a SANISAND constitutive model for response of solid skeleton. After a detailed verification and validation of the formulation and implementation of the developed numerical tool, it is employed in the seismic response of saturated porous media. The study includes examination of the mechanism of propagation of the earthquake-induced shear waves and liquefaction phenomenon in uniform and layered profiles of saturated sand deposits.     

饱和度对波在土层交界面的反射、透射系数的影响

把均匀流体概念引入Biot两相多孔介质动力理论中，用Biot的两相多孔介质模型模拟不完全饱和土层，给出SV波，P波从不完全饱和土层入射到弹性土层时，在土层交界面上反射，透射系数的表达式，结果表明与完全饱和相比，饱和度发生很小的变化就会对交界面上反射，透射系数产生很大的影响。今后应该重视饱和度变化对地震动力响应的影响。     

Dimensional analysis of two-phase flow including a rate-dependent capillary pressure–saturation relationship

The macroscopic modelling of two-phase flow processes in subsurface hydrosystems or industrial applications on the Darcy scale usually requires a constitutive relationship between capillary pressure and saturation, the P_c(S_w) relationship. Traditionally, it is assumed that a unique relation between P_c and S_w exists independently of the flow conditions as long as hysteretic effects can be neglected. Recently, this assumption has been questioned and alternative formulations have been suggested. For example, the extended P_c(S_w) relationship by Hassanizadeh and Gray [Hassanizadeh SM, Gray WG. Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries. Adv Water Resources 1990;13(4):169–86] proposes that the difference between the phase pressures to the equilibrium capillary pressure is a linear function of the rate of change of saturation, thereby introducing a constant of proportionality, the coefficient τ. It is desirable to identify cases where the extended relationship needs to be considered. Consequently, a dimensional analysis is performed on the basis of the two-phase balance equations. In addition to the well-known capillary and gravitational number, the dimensional analysis yields a new dimensionless number. The dynamic number Dy quantifies the ratio of dynamic capillary to viscous forces. Relating the dynamic to the capillary as well as the gravitational number gives the new numbers DyC and DyG, respectively. For given sets of fluid and porous medium parameters, the dimensionless numbers Dy and DyC are interpreted as functions of the characteristic length and flow velocity. The simulation of an imbibition process provides insight into the interpretation of the characteristic length scale. The most promising choice for this length scale seems to be the front width. We conclude that consideration of the extended P_c(S_w) relationship may be important for porous media with high permeability, small entry pressure and high coefficient τ when systems with a small characteristic length (e.g. steep front) and small characteristic time scale are under investigation.     

Modelling of viscoacoustic wave propagation in transversely isotropic media using decoupled fractional Laplacians

A new wave equation is derived for modelling viscoacoustic wave propagation in transversely isotropic media under acoustic transverse isotropy approximation. The formulas expressed by fractional Laplacian operators can well model the constant-Q (i.e. frequency-independent quality factor) attenuation, anisotropic attenuation, decoupled amplitude loss and velocity dispersion behaviours. The proposed viscoacoustic anisotropic equation can keep consistent velocity and attenuation anisotropy effects with that of qP-wave in the constant-Q viscoelastic anisotropic theory. For numerical simulations, the staggered-grid pseudo-spectral method is implemented to solve the velocity–stress formulation of wave equation in the time domain. The constant fractional-order Laplacian approximation method is used to cope with spatial variable-order fractional Laplacians for efficient modelling in heterogeneous velocity and Q media. Simulation results for a homogeneous model show the decoupling of velocity dispersion and amplitude loss effects of the constant-Q equation, and illustrate the influence of anisotropic attenuation on seismic wavefields. The modelling example of a layered model illustrates the accuracy of the constant fractional-order Laplacian approximation method. Finally, the Hess vertical transversely isotropic model is used to validate the applicability of the formulation and algorithm for heterogeneous media.     

Boundary integral formulation and two-dimensional fundamental solutions for dynamic behavior analysis of unsaturated soils

In this paper the coupled equations governing the dynamic behavior of unsaturated soils are derived based on the poromechanics theory within the framework of the suction-based mathematical model presented by Gatmiri (1997) [Gatmiri B. Analysis of fully coupled behavior of unsaturated porous medium under stress, suction and temperature gradient. Final report of CERMES-EDF, 1997] and Gatmiri et al. (1998) [Gatmiri B, Delage P, Cerrolaza M, UDAM: a powerful finite element software for the analysis of unsaturated porous media. Adv Eng Software 1998; 29(1): 29–43]. In this formulation, the solid skeleton displacements, water pressure and air pressure are presumed to be independent variables. The Boundary Integral formulations as well as fundamental solutions for such a dynamic u–p_w–p_a theory are presented in this paper for the first time. The boundary integral equations are derived via the use of the weighted residuals method in a way that permits an easy discretization and implementation in a Boundary Element code. Also, the associated two dimensional (2D) fundamental solutions for such deformable porous medium with linear elastic behavior are derived in Laplace transform domain using the method of Hörmander. Finally, some numerical results are presented to show the accuracy of the proposed solutions. The derived results are verified analytically by comparison with the previously introduced corresponding fundamental solutions in elastodynamic limiting case.     

Field scale impacts of spatially correlated relative permeability in heterogeneous multiphase systems

Two-dimensional numerical simulations of two-phase (DNAPL-water) flow in spatially correlated random fields demonstrate the influence of nonwetting phase (NWP) relative permeability–saturation (k_r,N–S_W) relationships correlated to porous media intrinsic permeability (k). Both the volume of porous media invaded by the NWP and the length of time during which the NWP is migrating are under predicted if k_r,N–k correlation is not accounted for in the model formulation. Not accounting for the k_r,N–k correlation resulted in under predicting the volume of porous media invaded by up to approximately 10%, which is likely not significant for many practical applications. However, not accounting for the k_r,N–k correlation resulted in under predicting field scale migration times by up to a factor of 4, which is likely significant in that the migration times are on the order of years to several decades for the DNAPL (1,2-DCE) considered in this study. The under prediction of migration times was greater for lower permeability aquifers.     

On the transport of moisture in polythermal glaciers

Abstract

We reconsider the problem of formulation of a model for polythermal glaciers, focussing attention in particular on the temperate zone where ice and water can coexist at the melting temperature. The energy equation for the ice-water mixture in this zone introduces a moisture flux, and a constitutive law for this flux is required. By analogy with the flow through a porous medium, we use Darcy's law (i.e. the second momentum equation of a two-phase flow model with “porous” geometry), and then require a mechanical constitutive relation relating the water pressure p _w to the average ice pressure p _i . Experience in two phase flows suggests that p _w =p _i may be problematical, and experience in soil mechanics suggests it is inaccurate. A constitutive relation is therefore presented based on work of Nye (1976), and its effect on the well-posedness of the model is examined. Considerations of the sort presented here have clear relevance in the formulation of similar problems in other geophysical situations, notably mantle convection.     

Paraxial ray methods for anisotropic inhomogeneous media

A new formalism of surface-to-surface paraxial matrices allows a very general and flexible formulation of the paraxial ray theory, equally valid in anisotropic and isotropic inhomogeneous layered media. The formalism is based on conventional dynamic ray tracing in Cartesian coordinates along a reference ray. At any user-selected pair of points of the reference ray, a pair of surfaces may be defined. These surfaces may be arbitrarily curved and oriented, and may represent structural interfaces, data recording surfaces, or merely formal surfaces. A newly obtained factorization of the interface propagator matrix allows to transform the conventional 6 × 6 propagator matrix in Cartesian coordinates into a 6 × 6 surface-to-surface paraxial matrix. This matrix defines the transformation of paraxial ray quantities from one surface to another. The redundant non-eikonal and ray-tangent solutions of the dynamic ray-tracing system in Cartesian coordinates can be easily eliminated from the 6 × 6 surface-to-surface paraxial matrix, and it can be reduced to 4 × 4 form. Both the 6 × 6 and 4 × 4 surface-to-surface paraxial matrices satisfy useful properties, particularly the symplecticity. In their 4 × 4 reduced form, they can be used to solve important boundary-value problems of a four-parametric system of paraxial rays, connecting the two surfaces, similarly as the well-known surface-to-surface matrices in isotropic media in ray-centred coordinates. Applications of such boundary-value problems include the two-point eikonal, relative geometrical spreading, Fresnel zones, the design of migration operators, and more.     

H-adaptivity applied to liquefiable soil in nonlinear analysis of soil–pile interaction

This paper is concerned with application of the h-adaptive finite element method to dynamic analysis of a pile in liquefiable soil considering large deformation. In finite element analysis of pile behavior in liquefiable soil during an earthquake, especially considering large deformation of liquefied ground, error due to discretization in the zone near the pile becomes very large. Our purpose was to refine the approximation of the finite element method. The updated Lagrangian formulation and a cyclic elasto-plastic model based on the kinematic hardening rule were adopted to deal with the nonlinearity of the soil. The mixed finite element and finite difference methods together with the u-p formulation and Biot's two-phase mixture theory were used. To improve the accuracy and increase the efficiency of finite element analysis, an h-adaptive scheme that included a posteriori error estimation and h-version mesh refinement was applied to the analysis. The calculated results of effective stress were smoothed locally by the extrapolation method and smoothed stress was used to calculate the L₂ norm of the effective stress error in the last step of the calculation of each time increment. The mesh was refined by a fission procedure based on the indication of the error estimate As a numerical example, a soil–pile interaction system loaded cyclically was analyzed by our method.     

Eigenrays in 3D heterogeneous anisotropic media,Part II: Dynamics

This paper is the second in a sequel of two papers and dedicated to the computation of paraxial rays and dynamic characteristics along the stationary rays obtained in the first paper. We start by formulating the linear, second‐order, Jacobi dynamic ray tracing equation. We then apply a similar finite‐element solver, as used for the kinematic ray tracing, to compute the dynamic characteristics between the source and any point along the ray. The dynamic characteristics in our study include the relative geometric spreading and the phase correction due to caustics (i.e. the amplitude and the phase of the asymptotic form of the Green's function for waves propagating in 3D heterogeneous general anisotropic elastic media). The basic solution of the Jacobi equation is a shift vector of a paraxial ray in the plane normal to the ray direction at each point along the central ray. A general paraxial ray is defined by a linear combination of up to four basic vector solutions, each corresponds to specific initial conditions related to the ray coordinates at the source. We define the four basic solutions with two pairs of initial condition sets: point–source and plane‐wave. For the proposed point–source ray coordinates and initial conditions, we derive the ray Jacobian and relate it to the relative geometric spreading for general anisotropy. Finally, we introduce a new dynamic parameter, similar to the endpoint complexity factor, presented in the first paper, used to define the measure of complexity of the propagated wave/ray phenomena. The new weighted propagation complexity accounts for the normalized relative geometric spreading not only at the receiver point, but along the whole stationary ray path. We propose a criterion based on this parameter as a qualifying factor associated with the given ray solution. To demonstrate the implementation of the proposed method, we use several isotropic and anisotropic benchmark models. For all the examples, we first compute the stationary ray paths, and then compute the geometric spreading and analyse these trajectories for possible caustics. Our primary aim is to emphasize the advantages, transparency and simplicity of the proposed approach.