Three‐dimensional nonlinear finite element analysis of pile groups in saturated porous media using a new transmitting boundary

The dynamic behaviour of pile groups subjected to an earthquake base shaking is analysed. An analysis is formulated in the time domain and the effects of material nonlinearity of soil, pile–soil–pile kinematic interaction and the superstructure–foundation inertial interaction on seismic response are investigated. Prediction of response of pile group–soil system during a large earthquake requires consideration of various aspects such as the nonlinear and elasto‐plastic behaviour of soil, pore water pressure generation in soil, radiation of energy away from the pile, etc. A fully explicit dynamic finite element scheme is developed for saturated porous media, based on the extension of the original formulation by Biot having solid displacement (u) and relative fluid displacement (w) as primary variables (u–w formulation). All linear relative fluid acceleration terms are included in this formulation. A new three‐dimensional transmitting boundary that was developed in cartesian co‐ordinate system for dynamic response analysis of fluid‐saturated porous media is implemented to avoid wave reflections towards the structure. In contrast to traditional methods, this boundary is able to absorb surface waves as well as body waves. The pile–soil interaction problem is analysed and it is shown that the results from the fully coupled procedure, using the advanced transmitting boundary, compare reasonably well with centrifuge data. Copyright © 2007 John Wiley & Sons, Ltd.     

Numerical simulation of fully saturated porous materials

Fully coupled, porous solid–fluid formulation, implementation and related modeling and simulation issues are presented in this work. To this end, coupled dynamic field equations with u?p?U formulation are used to simulate pore fluid and soil skeleton (elastic–plastic porous solid) responses. Present formulation allows, among other features, for water accelerations to be taken into account. This proves to be useful in modeling dynamic interaction of media of different stiffnesses (as in soil–foundation–structure interaction). Fluid compressibility is also explicitly taken into account, thus allowing excursions into modeling of limited cases of non‐saturated porous media. In addition to these features, present formulation and implementation models in a realistic way the physical damping, which dissipates energy. In particular, the velocity proportional damping is appropriately modeled and simulated by taking into account the interaction of pore fluid and solid skeleton. Similarly, the displacement proportional damping is physically modeled through elastic–plastic processes in soil skeleton. An advanced material model for sand is used in present work and is discussed at some length. Also explored in this paper are the verification and validation issues related to fully coupled modeling and simulations of porous media. Illustrative examples describing the dynamical behavior of porous media (saturated soils) are presented. The verified and validated methods and material models are used to predict the behavior of level and sloping grounds subjected to seismic shaking. Copyright © 2008 John Wiley & Sons, Ltd.     

Hybrid time integration and coupled solution methods for nonlinear finite element analysis of partially saturated deformable porous media at small strain

The goal of the paper is to determine the most efficient, yet accurate and stable, finite element nonlinear solution method for analysis of partially saturated deformable porous media at small strain. This involves a comparison between fully implicit, semi‐implicit, and explicit time integration schemes, with monolithically coupled and staggered‐coupled nonlinear solution methods and the hybrid combination thereof. The pore air pressure p_a is assumed atmospheric, that is, p_a=0 at reference pressure. The solid skeleton is assumed to be pressure‐sensitive nonlinear isotropic elastic. Coupled partially saturated ‘consolidation’ in the presence of surface infiltration and traction is simulated for a simple one‐dimensional uniaxial strain example and a more complicated plane strain slope example with gravity loading. Three mixed plane strain quadrilateral elements are considered: (i) Q4P4; (ii) stabilized Q4P4S; and (iii) Q9P4; “Q” refers to the number of solid skeleton displacement nodes, and “P” refers to the number of pore fluid pressure nodes. The verification of the implementation against an analytical solution for partially saturated pore water flow (no solid skeleton deformation) and comparison between the three time integration schemes (fully implicit, semi‐implicit, and explicit) are presented. It is observed that one of the staggered‐coupled semi‐implicit schemes (SIS(b)), combined with the fully implicit monolithically coupled scheme to resolve sharp transients, is the most efficient computationally. Copyright © 2015 John Wiley & Sons, Ltd.     

A robust numerical framework for simulating localized failure and fracture propagation in frictional materials

A computationally robust framework for simulating geomaterial failure patterns is presented in this paper. Finite element simulations which feature the use of embedded discontinuities to track material failure are known to suffer from convergence issues due to a lack of robustness. Oftentimes, complex time step-cutting schemes or arc-length methods are required in order to achieve convergence. This may invariably limit the complexity of constitutive models available for use in tracking nonlinear material behavior. To this end, we use an implicit–explicit integration scheme [Impl–Ex (Oliver et al. in Comput Methods Appl Mech Eng 195(52):7093–7114, 2006)] coupled with a novel constitutive model which allows for combined opening and shearing displacement in tension, as well as frictional sliding in compression. We show that this framework is suitable for capturing complex fracture patterns in geomaterial structures without the need for elaborate continuance schemes.     

饱和土地下源u-P形式解答动力响应计算

一直以来,由Biot孔隙弹性动力方程得到的饱和土地下源Green函数都是u-w形式（u为固相介质位移,w为流相相对于固相的平均位移）。应用两相介质纵波解耦理论,得到了饱和土半空间地下点源荷载的u-P形式（P为孔压）Green函数频域解答;克服了u-w形式Green函数在边界元（BEM）积分时的增根影响。再由Hankel反演,结合Somigliana表象积分,完成BEM计算。并以计算结果分析了地下集中力作用时,饱和土位移、孔压、排水量等动力特性,这对地铁等交通工程、地震工程、土-结构动力相互作用（SSI）的响应计算都具有较重要应用价值。     

Stabilized single-point 4-node quadrilateral element for dynamic analysis of fluid saturated porous media

An accurate and efficient low-order quadrilateral mixed u?Cp element suitable for dynamic analysis of fluid saturated porous media is presented. The element uses physical hourglass stabilization to facilitate single-point integration for the solid phase, and non-residual stabilization of the fluid phase to circumvent instability in the incompressible-impermeable limit due to the use of equal-order interpolation for the displacement and pressure fields. Element behavior is verified and demonstrated through several numerical examples.     

A stabilized extended finite element framework for hydraulic fracturing simulations

We present a stabilized extended finite element formulation to simulate the hydraulic fracturing process in an elasto‐plastic medium. The fracture propagation process is governed by a cohesive fracture model, where a trilinear traction‐separation law is used to describe normal contact, cohesion and strength softening on the fracture face. Fluid flow inside the fracture channel is governed by the lubrication equation, and the flow rate is related to the fluid pressure gradient by the ‘cubic’ law. Fluid leak off happens only in the normal direction and is assumed to be governed by the Carter's leak‐off model. We propose a ‘local’ U‐P (displacement‐pressure) formulation to discretize the fluid‐solid coupled system, where volume shape functions are used to interpolate the fluid pressure field on the fracture face. The ‘local’ U‐P approach is compatible with the extended finite element framework, and a separate mesh is not required to describe the fluid flow. The coupled system of equations is solved iteratively by the standard Newton‐Raphson method. We identify instability issues associated with the fluid flow inside the fracture channel, and use the polynomial pressure projection method to reduce the pressure oscillations resulting from the instability. Numerical examples demonstrate that the proposed framework is effective in modeling 3D hydraulic fracture propagation. Copyright © 2016 John Wiley & Sons, Ltd.     

Thermodynamic modeling of binary CH4–CO2 fluid inclusions

An equation of state (EOS) explicit in Helmholtz free energy has been improved to calculate the PVTx and vapor–liquid phase equilibrium properties of CH₄–CO₂ fluid mixture. This EOS, where four mixing parameters are used, is based on highly accurate EOSs recommended by NIST for pure components (CH₄ and CO₂) and contains a simple generalized departure function presented by Lemmon and Jacobsen (1999). Comparison with experimental data available indicates that the EOS can calculate both vapor–liquid phase equilibrium and volumetric properties of this binary fluid system with accuracy close to that of experimental data up to high temperature and pressure within full range of composition. The EOS of CH₄–CO₂ fluid, together with the updated Gibbs free energy model of solid CO₂ (dry ice), is applied to calculate the CH₄ content (x_CH4) and molar volume (V_m) of the CH₄–CO₂ fluid inclusion based on the assumption that the volume of an inclusion keeps constant during heating and cooling. $V_{\text{m}}-x_{{\text{CH}}_4}$ diagrams are presented, which describe phase transitions involving vapor, liquid and CO₂ solid phases of CH₄–CO₂ fluid inclusions. Isochores of CH₄–CO₂ inclusions at given $x_{{\text{CH}}_4}$ and V_m can be easily calculated from the improved EOS.     

Comparison of numerical formulations for the modeling of tensile loaded suction buckets

The tensile resistance of a suction bucket is investigated using three different numerical formulations. The first formulation utilizes the three-field u-p-U formulation accounting for solid and fluid displacements, u and U, as well as the pore-fluid pressure, p. The two other formulations comprise the simpler u-p formulation in its dynamic and quasi-static form, accounting only for solid displacement and pore-fluid pressure. As basis for comparison, the tensile resistance of a single suction bucket is investigated using a velocity-driven model for a wide range of velocities. It is found, that the quasi-static u-p formulation is sufficient for most relevant velocities.     

A FIC‐based stabilized mixed finite element method with equal order interpolation for solid–pore fluid interaction problems

A new mixed displacement‐pressure element for solving solid–pore fluid interaction problems is presented. In the resulting coupled system of equations, the balance of momentum equation remains unaltered, while the mass balance equation for the pore fluid is stabilized with the inclusion of higher‐order terms multiplied by arbitrary dimensions in space, following the finite calculus (FIC) procedure. The stabilized FIC‐FEM formulation can be applied to any kind of interpolation for the displacements and the pressure, but in this work, we have used linear elements of equal order interpolation for both set of unknowns. Examples in 2D and 3D are presented to illustrate the accuracy of the stabilized formulation for solid–pore fluid interaction problems. Copyright © 2016 John Wiley & Sons, Ltd.     

Plasticity with generalized hardening: constitutive modeling and computational aspects

In this work, an extended theory of plasticity with generalized hardening is proposed to describe the response of geomaterials under both mechanical and environmental processes, which include as special cases several elastoplastic constitutive equations proposed in the literature to model such processes as desaturation or suction hardening, thermal softening, chemo-mechanical coupling effects in fine-grained soils, as well as weathering of soft rocks. In the formulation of the theory, the coupling between mechanical and environmental processes takes place at two levels: first, as an additional direct contribution to the constitutive stress changes, taking place in both elastic and elastoplastic processes; and second, as a result of the evolution of the internal state variables induced by changes in the environmental process variables. This last effect is incorporated through a set of generalized hardening rules. As an example of application, the general formulation is specialized to the particular case of weak calcarenite rocks undergoing degradation processes due to plastic deformations, changes in degree of saturation (short-term debonding) and chemical dissolution of the bond material and the solid grains (long-term debonding). The resulting model is implemented in a FE code by means of an implicit generalized backward Euler algorithm, suitably modified to incorporate the full formalism of plasticity with generalized hardening. Results of numerical simulations carried out at the element level show the accuracy and efficiency properties of the proposed stress-point algorithm. The simulation of a representative initial-boundary value problem demonstrates the practical relevance of environmental degradation effects in practical applications, over periods of time comparable with the life cycle of most geotechnical structures.     

On nonequilibrium models of spontaneous countercurrent imbibition

Spontaneous displacement of the non-wetting phase by a wetting phase in a porous medium, known as spontaneous imbibition, is an important mechanism of oil recovery from fractured reservoirs. In this paper, we consider the nonequilibrium model, proposed by Aryana and Kovscek, where consitutive relationships for multiphase flow in porous media are functions of a locally moving time-average saturation, and allow relaxation time to be an explicit function of local saturation. We obtain asymptotic self-similar solutions for early and late times. At very early stages, the time-scale of the process characterizing the cumulative volume of displaced fluid is a power function with an exponent of \(\frac {1}{2}+\frac {1}{2r+1}\) where r is the inverse of pore size distribution index of the medium in question. Additionally, the cumulative volume of displaced fluid at late times is independent of relaxation time, and this volume approaches the square root of time asymptotically. Finally, the late-time solution for recovery is compared with experimental observations.     

Three dimensional elasto‐plastic interface for embedded beam elements with interaction surface for the analysis of lateral loading of piles

This paper presents a numerical formulation for a three dimensional elasto‐plastic interface, which can be coupled with an embedded beam element in order to model its non‐linear interaction with the surrounding solid medium. The formulation is herein implemented for lateral loading of piles but is able to represent soil‐pile interaction phenomena in a general manner for different types of loading conditions or ground movements. The interface is formulated in order to capture localized material plasticity in the soil surrounding the pile within the range of small to moderate lateral displacements. The interface is formulated following two different approaches: (i) in terms of beam degrees of freedoms; and (ii) considering the displacement field of the solid domain. Each of these alternatives has its own advantages and shortcomings, which are discussed in this paper. The paper presents a comparison of the results obtained by means of the present formulation and by other well‐established analysis methods and test results published in the literature. Copyright © 2016 John Wiley & Sons, Ltd.     

Higher resolution total velocity <Emphasis Type="Italic">Vt</Emphasis> and <Emphasis Type="Italic">Va</Emphasis> finite-volume formulations on cell-centred structured and unstructured grids

Novel cell-centred finite-volume formulations are presented for incompressible and immiscible two-phase flow with both gravity and capillary pressure effects on structured and unstructured grids. The Darcy-flux is approximated by a control-volume distributed multipoint flux approximation (CVD-MPFA) coupled with a higher resolution approximation for convective transport. The CVD-MPFA method is used for Darcy-flux approximation involving pressure, gravity, and capillary pressure flux operators. Two IMPES formulations for coupling the pressure equation with fluid transport are presented. The first is based on the classical total velocity Vt fractional flow (Buckley Leverett) formulation, and the second is based on a more recent Va formulation. The CVD-MPFA method is employed for both Vt and Va formulations. The advantages of both coupled formulations are contrasted. The methods are tested on a range of structured and unstructured quadrilateral and triangular grids. The tests show that the resulting methods are found to be comparable for a number of classical cases, including channel flow problems. However, when gravity is present, flow regimes are identified where the Va formulation becomes locally unstable, in contrast to the total velocity formulation. The test cases also show the advantages of the higher resolution method compared to standard first-order single-point upstream weighting.     

Compositional flow modeling using a multi-point flux mixed finite element method

We present a general compositional formulation using multi-point flux mixed finite element (MFMFE) method on general hexahedral grids. The mixed finite element framework allows for local mass conservation, accurate flux approximation, and a more general treatment of boundary conditions. The multi-point flux inherent in MFMFE scheme allows the usage of a full permeability tensor. The proposed formulation is an extension of single and two-phase flow formulations presented by Wheeler and Yotov, SIAM J. Numer. Anal. 44(5), 2082–2106 (35) with similar convergence properties. Furthermore, the formulation allows for black oil, single-phase and multi-phase incompressible, slightly and fully compressible flow models utilizing the same design for different fluid systems. An accurate treatment of diffusive/dispersive fluxes owing to additional velocity degrees of freedom is also presented. The applications areas of interest include gas flooding, CO ₂ sequestration, contaminant removal, and groundwater remediation.     

Applications of the SW96 formulation in the thermodynamic calculation of fluid inclusions and mineral-fluid equilibria

The SW96 formulation explicit in Helmholtz free energy proposed by Span and Wagner (1996) is the most accurate multifunction equation of state of CO₂ fluid, from which all thermodynamic properties can be obtained over a wide temperature-pressure range from 216.592 to 1100 K and from 0 to 8000 bar with or close to experimental accuracy. This paper reports the applications of the SW96 formulation in fluid inclusions and mineral-fluid equilibria. A reliable and highly efficient algorithm is presented for the saturated properties of CO₂ so that the formulation can be conveniently applied in the study of fluid inclusions, such as calculation of homogenization pressures, homogenization densities (or molar volumes), volume fractions of vapor phase and isochores. Meanwhile, the univariant curves of some typical decarbonation reactions of minerals are calculated with the SW96 formulation and relevant thermodynamic models of minerals. The computer code of the SW96 formulation can be obtained from the corresponding author.     

Experimental Na-K distribution between amphiboles and aqueous chloride solutions, and a mixing model along the richterite – K-richterite join

The cation exchange equilibrium has been investigated by hydrothermal experiments at 700 and 800°C at 200 MPa. To avoid equilibration problems of conventional exchange experiments, we synthesized amphiboles with an excess fluid allowing exchange between solid and fluid during the experiment. The exchangeable cations Na and K were provided as excess 1 to 2n chloridic solution. These exchange syntheses can be described by the reaction equation with ^(aq) for hydroxides and chlorides in aqueous solutions and _{( s )} and _{( p )}?=?start and product fluid. The amphiboles grew in presence of the exchange fluid and adjusted their stoichiometry in equilibrium with the fluid phase. The solid products consist of more than 99% amphibole (Na,K-richterite_ss) with traces of diopside and quartz. The amphiboles are up to 1?mm long and often ≈ 40 μm thick. Detailed EMP- and HRTEM-observations show that they are chemically homogeneous and structurally wellordered. The experimental results give consistent phase relations in the reciprocal ternary system Na-richterite–K-richterite–NaCl–KCl. We analysed the product fluid with AAS- and ICP-methods. The Na-K distribution coefficients between fluid and amphiboles of the richterite–K-richterite join are close to unity at 700°C and 800°C at 200 MPa. Small systematic deviations are explained by a symmetric solution model for the A-position of the amphiboles. Using ideal mixing for H₂O-NaCl-KCl fluids, a mixing model for the system richterite–K-richterite is presented. We suggest that the composition of richterite solid solutions can be used as a sensor for NaCl/KCl-ratios in metamorphic fluids.     

Oxygen and hydrogen isotope evidence for meteoric water infiltration during mylonitization and uplift in the Ruby Mountains-East Humboldt Range core complex,Nevada

Stable isotope analyses of rocks and minerals associated with the detachment fault and underlying mylonite zone exposed at Secret Creek gorge and other localities in the Ruby-East Humboldt Range metamorphic core complex in northeastern Nevada provide convincing evidence for meteoric water infiltration during mylonitization. Whole-rock ¹⁸O values of the lower plate quartzite mylonites (95% modal quartz) have been lowered by up to 10 per mil compared with structurally lower, compositionally similar, unmylonitized material. Biotite from these rocks has D values ranging from -125 to -175, compared to values of -55 to-70 in biotite from unmylonitized rocks. Mylonitized leucogranites have large disequilibrium oxygen isotope fractionations ( ^{quartz-feldspar} up to 8 per mil) relative to magmatic values ( ^{quartz-feldspar}1 to 2 per mil)). Meteoric water is the only major oxygen and hydrogen reservoir with an isotopic composition capable of generating the observed values. Fluid inclusion water from unstrained quartz in silicified breccia has a D value of-119 which provides a plausible estimate of the D of the infiltrating fluid, and is similar to the isotopic composition of present-day and Tertiary local meteoric water. The quartzite mylonite biotites would have been in equilibrium with such a fluid at temperatures of 480–620° C, similar to independent estimates of the temperature of mylonitization. The relatively high temperatures required for isotopic exchange between quartz and water, the occurrence of fluid inclusion trails and deformed veins in quartzite mylonites, and the spatial association of the low-¹⁸O, low-D rocks with the shear zone all constrain isotopic exchange to the mylonitic (plastic) deformation event. These observations suggest thata significant amount of meteoric water infiltrated the shear zone during mylonitization to depths of at least 5 to 10 km below the surface. The depth of penetration of meteoric fluids into the lower plate mylonites was at least 70 meters below the detachment fault. In contrast, the upper-plate unmylonitized fault slices are dominated by brittle fracture and are often intensely veined (carbonates) or silicified (volcanic rocks and breccias). The fluids associated with the veining and silicification were also meteoric as evidenced by low ¹⁸O values of the veins, which are often 10 per mil lower than the adjacent carbonate matrix, and the exceptionally low ¹⁸O values (down to-4.4) of the breccias. Several previous studies have documented the infiltration of meteoric fluids into the brittley deformed upper plate rocks of core complexes, but this study provides convincing evidence that surface fluids have penetrated lower plate rocks undergoing plastic deformation. It is proposed that infiltration took place as the shear zone began the transition from plastic flow to brittle fracture while the lower plate rocks were being uplifted. During this period, plastic flow and brittle fracture were operating simultaneously, perhaps allowing upper plate meteoric fluids to be seismically pumped down into the lower plate mylonites.     

Dynamics of unsaturated soils using various finite element formulations

Unsaturated soils are three‐phase porous media consisting of a solid skeleton, pore liquid, and pore gas. The coupled mathematical equations representing the dynamics of unsaturated soils can be derived based on the theory of mixtures. Solution of these fully coupled governing equations for unsaturated soils requires tremendous computational resources because three individual phases and interactions between them have to be taken into account. The fully coupled equations governing the dynamics of unsaturated soils are first presented and then two finite element formulations of the governing equations are presented and implemented within a finite element framework. The finite element implementation of all the terms in the governing equations results in the complete formulation and is solved for the first time in this paper. A computationally efficient reduced formulation is obtained by neglecting the relative accelerations and velocities of liquid and gas in the governing equations to investigate the effects of fluid flow in the overall behavior. These two formulations are used to simulate the behavior of an unsaturated silty soil embankment subjected to base shaking and compared with the results from another commonly used partially reduced formulation that neglects the relative accelerations, but takes into account the relative velocities. The stress–strain response of the solid skeleton is modeled as both elastic and elastoplastic in all three analyses. In the elastic analyses no permanent deformations are predicted and the displacements of the partially reduced formulation are in between those of the reduced and complete formulations. The frequency of vibration of the complete formulation in the elastic analysis is closer to the predominant frequency of the base motion and smaller than the frequencies of vibration of the other two analyses. Proper consideration of damping due to fluid flows in the complete formulation is the likely reason for this difference. Permanent deformations are predicted by all three formulations for the elastoplastic analyses. The complete formulation, however, predicts reductions in pore fluid pressures following strong shaking resulting in somewhat smaller displacements than the reduced formulation. The results from complete and reduced formulations are otherwise comparable for elastoplastic analyses. For the elastoplastic analysis, the partially reduced formulation leads to stiffer response than the other two formulations. The likely reason for this stiffer response in the elastoplastic analysis is the interpolation scheme (linear displacement and linear pore fluid pressures) used in the finite element implementation of the partially reduced formulation. Copyright © 2008 John Wiley & Sons, Ltd.