Characterization of reservoir properties like porosity and permeability in reservoir models typically relies on history matching of production data, well pressure data, and possibly other fluid-dynamical data. Calibrated (history-matched) reservoir models are then used for forecasting production and designing effective strategies for improved oil and gas recovery. Here, we perform assimilation of both flow and deformation data for joint inversion of reservoir properties. Given the coupled nature of subsurface flow and deformation processes, joint inversion requires efficient simulation tools of coupled reservoir flow and mechanical deformation. We apply our coupled simulation tool to a real underground gas storage field in Italy. We simulate the initial gas production period and several decades of seasonal natural gas storage and production. We perform a probabilistic estimation of rock properties by joint inversion of ground deformation data from geodetic measurements and fluid flow data from wells. Using an efficient implementation of the ensemble smoother as the estimator and our coupled multiphase flow and geomechanics simulator as the forward model, we show that incorporating deformation data leads to a significant reduction of uncertainty in the prior distributions of rock properties such as porosity, permeability, and pore compressibility. Copyright © 2015 John Wiley & Sons, Ltd.

Damage induced by microcracking affects not only the mechanical behaviour of geomaterials but also their hydraulic properties. Evaluating these impacts is important for many engineering applications, such as the safety assessment of radioactive waste disposal facilities. This paper presents a new constitutive model accounting simultaneously for the impact of damage on hydraulic and mechanical properties of unsaturated poroplastic geomaterials. The hydro-mechanical coupling is formulated by means of the thermodynamic framework for partially saturated media, extended by taking into account isotropic damage and plasticity. State and complementary laws are governed by the so-called plastic effective stress and equivalent pore pressure. Assuming a bimodal pore size distribution for cracked porous media, the hydraulic part (water retention curve and hydraulic conductivity) is modelled using phenomenological functions of damage variable. The participation of damage on both mechanical and hydraulic part enables this model to describe bilateral couplings between them. This coupled model is then validated against a number of experimental data obtained from Callovo-Oxfordian argillite, which is the possible host rock for a radioactive waste disposal in France. Parametric studies are also carried out to check the consistency and to better demonstrate the bilateral couplings in the model. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, an implementation of fractional plastic flow rule in the framework of implicit and explicit procedures is under consideration. The fractional plastic flow rule is obtained from a generalisation of the classical plastic flow rule utilising fractional calculus. The key feature of this new concept is that in general, the non-associative flow is obtained without necessity of additional potential assumption. If needed, the model can cover the anisotropy induced by plastic deformation. Illustrative examples showing the unusual flexibility of this model are also presented. Copyright © 2015 John Wiley & Sons, Ltd.

In engineering practices, different numerical methods for fluid flow simulation and solid deformation/stress simulation are adopted to model fluid–structure interaction problems in porous media. Cell-centered finite volume method is widely used in fluid flow simulation, while the solid deformation/stress simulation is usually accomplished by using the Galerkin vertex-centered finite element method, which leads to the incompatibility between cell variables with nodal variables. Therefore, the data transfer between cell variables and nodal variables is inevitable. Consequently, this kind of transfer will lead to extra artificial error. Hence, the major concern is how to minimize the error due to cell to node projections. In this paper, a problem of pore pressure diffusion within a one-dimensional heterogeneous porous medium is investigated. We present a new projection scheme and corresponding error formula, where the error control factor is introduced. The new projection scheme is based on piecewise linear interpolations. Results demonstrate that if the error control factor is chosen properly, the error due to the projection from cell to node can be controlled effectively, and the most desired zero error can be achieved. Finally, we analyze some practical cases in consideration of permeability contrast and mesh uniformity. Copyright © 2015 John Wiley & Sons, Ltd.

A new constitutive model for soft structured clays is developed based on an existing model called S-CLAY1S, which is a Cam clay type model that accounts for anisotropy and destructuration. The new model (E-SCLAY1S) uses the framework of logarithmic contractancy to introduce a new parameter that controls the shape of the yield surface as well as the plastic potential (as an assumed associated flow rule is applied). This new parameter can be used to fit the coefficient of earth pressure at rest, the undrained shear strength or the stiffness under shearing stress paths predicted by the model. The improvement to previous constitutive models that account for soil fabric and bonding is formulated within the contractancy framework such that the model predicts the uniqueness of the critical state line and its slope is independent of the contractancy parameter. Good agreement has been found between the model predictions and published laboratory results for triaxial compression tests. An important finding is that the contractancy parameter, and consequently the shape of the yield surface, seems to change with the degree of anisotropy; however, further study is required to investigate this response. From published data, the yield surface for isotropically consolidated clays seems ‘bullet’ or ‘almond’ shaped, similar to that of the Cam clay model; while for anisotropically consolidated clays, the yield surface is more elliptical, like a rotated and distorted modified Cam clay yield surface. Copyright © 2015 John Wiley & Sons, Ltd.

The failure of a discrete elastic-damage axial system is investigated using both a discrete and an equivalent continuum approach. The Discrete Damage Mechanics approach is based on a microstructured model composed of a series of periodic elastic-damage springs (axial Discrete Damage Mechanics lattice system). Such a discrete damage system can be associated with the finite difference formulation of a Continuum Damage Mechanics evolution problem. Several analytical and numerical results are presented for the tensile failure of this axial damage chain under its own weight.

The nonlocal Continuum Damage Mechanics models examined in this paper are mainly built from a continualization procedure applied to centered or uncentered finite difference schemes. The asymptotic expansion of the first-order upward difference equations leads to a first-order nonlocal model, whereas the asymptotic expansion of the centered finite difference equations leads to a second-order nonlocal Eringen's approach. To complete this study, a phenomenological nonlocal gradient approach is also examined and compared with the first continualization methods.

A comparison of the discrete and the continuous problems for the chains shows the effectiveness of the new micromechanics-based nonlocal Continuum Damage modeling, especially for capturing scale effects. For both continualized approaches, the length scale of the nonlocal models depends only on the cell size, while for the so-called phenomenological approach, the length scale may depend on the loading parameter. This apparent load-dependent length scale, already discussed in the literature with numerical arguments, is found to be sensitive to the postulated structure of the nonlocal model calibrated according to a lattice approach. Copyright © 2015 John Wiley & Sons, Ltd.

Asymptotic behaviour of soil deserves particular attention: If soil is deformed with a proportional strain path, the resulting stress path approaches asymptotically a proportional stress path. In this arcticle, we review existing experimental evidence on this phenomenon and discuss it in the frame of barodesy. Here, the presented relation is a modification of a barodetic expression and includes Jáky's relation, inhibits tensile stress and is able to predict asymptotic stress ratios based on experimental findings. The proposed relation is compared with experimental data as well as with the so-called stress-dilatancy relations and other constitutive relations proposed so far. Copyright © 2015 John Wiley & Sons, Ltd.

The cohesive-frictional nature of cementitious geomaterials raises great interest in the discrete element method (DEM) simulation of their mechanical behavior, where a proper bond failure criterion is usually required. In this paper, the failure of bond material between two spheres was investigated numerically using DEM that can easily reproduce the failure process of brittle material. In the DEM simulations, a bonded-grain system (composed of two particles and bond material in between) was discretized as a cylindrical assembly of very fine particles connecting two large end spheres. Then, the bonded-grain system was subjected to compression/tension, shear, rolling and torsion loadings and their combinations until overall failure (peak state) was reached. Bonded-grain systems with various sizes were employed to investigate bond geometry effects. The numerical results show that the compression strength is highly affected by bond geometry, with the tensile strength being dependent to a lesser degree. The shear, rolling and torsion strengths are all normal force dependent; i.e., with an increase in the normal force, these strengths first increase at a declining rate and then start to decrease upon the normal force exceeding a critical value. The combined actions of shear force, rolling moment and torque lead to a spherical failure envelope in a normalized loading space. The fitted bond geometry factors and bond failure envelopes obtained numerically in this three-dimensional study are qualitatively consistent with those in previous two-dimensional experiments. The obtained bond failure criterion can be incorporated into a future bond contact model. Copyright © 2015 John Wiley & Sons, Ltd.

Water pipe cooling has been widely used for the temperature control and crack prevention of massive concrete structures such as high dams. Because both under-cooling and over-cooling may reduce the efficiency of crack prevention, or even lead to great harm to structures, we need an accurate and robust numerical tool for the prediction of cooling effect. Here, a 3D *discrete FEM Iterative Algorithm* is introduced, which can simulate the concrete temperature gradient near the pipes, as well as the water temperature rising along the pipes. On the basis of the heat balance between water and concrete, the whole temperature field of the problem can be computed exactly within a few iteration steps. Providing the pipe meshing tool for building the FE model, this algorithm can take account of the water pipe distribution, the variation of water flow, water temperature, and other factors, while the *traditional equivalent algorithm* based on semi-theoretical solutions can only solve problems with constant water flow and water temperature. The validation and convergence are proved by comparing the simulated results and analytical solutions of two standard second-stage cooling problems. Then, a practical concrete block with different cooling schemes is analyzed and the influences of cooling factors are investigated. In the end, detailed guidance for pipe system optimization is provided. Copyright © 2015 John Wiley & Sons, Ltd.

We propose a discrete element model for brittle rupture. The material consists of a bidimensional set of closed-packed particles in contact. We explore the isotropic elastic behavior of this regular structure to derive a rupture criterion compatible to continuum mechanics. We introduce a classical criterion of mixed mode crack propagation based on the value of the stress intensity factors, obtained by the analysis of two adjacent contacts near a crack tip. Hence, the toughness becomes a direct parameter of the model, without any calibration procedure. We verify the consistency of the formulation as well as its convergence by comparison with theoretical solutions of tensile cracks, a pre-cracked beam, and an inclined crack under biaxial stress. Copyright © 2015 John Wiley & Sons, Ltd.

This paper presents a numerical formulation of a three dimensional embedded beam element for the modeling of piles, which incorporates an explicit interaction surface between soil and pile. 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 model assumes perfect adherence between beam and soil along the interaction surface. The paper presents a comparison of the results obtained by means of the present formulation and by means of a previously formulated embedded pile element without interaction surface, as well as reference semi-analytical solutions and a fully 3D finite element (FE) model. It is seen that the proposed embedded element provides a better convergence behavior than a previously formulated embedded element and is able to reproduce key features of a full 3D FE model. Copyright © 2015 John Wiley & Sons, Ltd.

The fundamental solution of a continuous line source, injecting fluid at a constant rate over the thickness of a poroelastic reservoir bounded by infinite impermeable elastic layers, is derived in this paper. This idealized problem has applications in hydrogeology and in petroleum engineering, as it can be used to assess the mechanical perturbations caused by injection or withdrawal of fluid in the subsurface through a vertical well. Construction of the solution takes advantage of the uncoupling of the pore pressure field, which, in this particular case, is given by the classical singular solution of the diffusion equation for an infinite line source. The mechanical fields then are determined by solving an elasticity-like problem with a body force field proportional to the time-dependent pore pressure gradient. On account of the problem symmetries, the Navier equations of elasticity reduce to two uncoupled partial differential equations for the radial and vertical (axial) displacement components, which are solved by a twofold application of Fourier and Hankel transforms. The solution exhibits different regimes at small, intermediate, and large times. When the diffusion radius, proportional to the square root of time, is smaller than or comparable to the thickness of the permeable layer, the induced deformation is confined to a region with a characteristic dimension of the same order as the diffusion radius. At large time, when the diffusion radius is large compared with the permeable layer thickness, the deformation rate in the reservoir is essentially oedometric (uniaxial). The different regimes of solutions are justified with a conceptual model based on identifying the evolving characteristics of complementary interior and exterior domain problems. The derived solution can serve as a valuable benchmark for coupled reservoir simulators. It also provides insights in to such problems as waterflooding, shearing at reservoir/cap rock interfaces, and stress redistribution around producing wells. Copyright © 2015 John Wiley & Sons, Ltd.

This paper addresses the open question of the effect of the third invariant in the definition of the strength criterion of porous media. This is here achieved in the context of a porous medium made up of a solid phase, in which strength criterion appears as an explicit function of Lode's angle. Application is provided from the Gurson hollow sphere model for which the macroscopic strength exhibits a lack of symmetry for axisymmetric macroscopic strain rate loadings. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, rate-dependent plasticity is employed to regularize, for the scope of mesh independency, the numerical solution in strain localization process of multiphase geomaterials. Towards this goal, the viscoplastic Duvaut–Lions and the viscoplastic Perzyna models are implemented in an existing finite element code for multiphase porous media. Perzyna's model is then extended according to the non-local approach, allowing for handling weakly rate-sensitive materials for which viscoplasticity is not sufficient to meet numerical requirements. Assuming quasi-static and isothermal conditions, the models are validated and applied in the numerical simulation of an undrained plane strain biaxial compression test on initially water-saturated dense sand taken from the existing literature. An overview of the significant factors such as loading velocity and soil permeability in conjunction with viscosity on the formation of localization pattern is presented. The interaction of the two internal lengths introduced by non-locality and viscosity is also discussed. Finally, the obtained results are analysed and compared emphasizing the capabilities and shortcomings of each regularization technique. Copyright © 2015 John Wiley & Sons, Ltd.

Hydraulic fracturing involves the initiation and propagation of fractures in rock formations by the injection of pressurized fluid. The largest use of hydraulic fracturing is in enhancing oil and gas production. Tiltmeters are sometimes used in the process to monitor the generated fracture geometry by measuring the fracture-induced deformations. Fracture growth parameters obtained from tiltmeter mapping can be used to study the effectiveness of such stimulations. In this work, we present a novel scheme that uses the ensemble Kalman Filter (EnKF) to assimilate tiltmeter data using a simple process model to describe the evolution of fracture growth parameters, and an observation model that maps the fracture geometry with the observed tilt. The forward observation model is based on the analytical solution for computing the displacements and tilts due to a point source displacement discontinuity in an elastic half-space developed by Okada . The displacement and tilts for any given fracture geometry are then obtained by numerical integration of this solution, by considering multiple point sources to be located at the quadrature points. The proposed method is validated using synthetic data sets generated from polygon and elliptical shaped fracture geometries. Finally, real data from a field site, where asymmetry was measured from the intersections of the hydraulic fracture with offset boreholes, have been analyzed. Preliminary results show that, in addition to extracting the fracture dip, orientation, and volume, the procedure is able to satisfactorily predict fracture growth parameters when the fracture is relatively close to the tiltmeter array and provides some insight into the development of asymmetry when the measurements are relatively far from the fracture plane. Copyright © 2015 John Wiley & Sons, Ltd.

The paper presents detailed FE simulation results of concrete elements under mixed-mode failure conditions according to the so-called shear-tension test by Nooru-Mohamed, characterized by curved cracks. A continuous and discontinuous numerical two-dimensional approach was used. In order to describe the concrete's behaviour within continuum mechanics, two different constitutive models were used. First, an elasto-plastic model with isotropic hardening and softening was assumed. In a compression regime, a Drucker–Prager criterion with a non-associated flow rule was used. In turn, in a tensile regime, a Rankine criterion with an associated flow rule was adopted. Second, an isotropic damage constitutive model was applied with a single scalar damage parameter and different definitions of the equivalent strain. Both constitutive laws were enriched by a characteristic length of micro-structure to capture properly strain localization. As an alternative approach, the extended finite element method was used. Our results were compared with the experimental ones and with results of other FE simulations reported in the literature. Copyright © 2015 John Wiley & Sons, Ltd.

Large-scale engineering computing using the discontinuous deformation analysis (DDA) method is time-consuming, which hinders the application of the DDA method. The simulation result of a typical numerical example indicates that the linear equation solver is a key factor that affects the efficiency of the DDA method. In this paper, highly efficient algorithms for solving linear equations are investigated, and two modifications of the DDA programme are presented. The first modification is a linear equation solver with high efficiency. The block Jacobi (BJ) iterative method and the block conjugate gradient with Jacobi pre-processing (Jacobi-PCG) iterative method are introduced, and the key operations are detailed, including the matrix-vector product and the diagonal matrix inversion. Another modification consists of a parallel linear equation solver, which is separately constructed based on the multi-thread and CPU-GPU heterogeneous platforms with OpenMP and CUDA, respectively. The simulation results from several numerical examples using the modified DDA programme demonstrate that the Jacobi-PCG is a better iterative method for large-scale engineering computing and that adoptive parallel strategies can greatly enhance computational efficiency. Copyright © 2015 John Wiley & Sons, Ltd.

We pay a revisit to some classical geomechanics problems using a novel computational multiscale modelling approach. The multiscale approach employs a hierarchical coupling of the finite element method (FEM) and the discrete element method. It solves a boundary value problem at the continuum scale by FEM and derives the material point response from the discrete element method simulation attached to each Gauss point of the FEM mesh. The multiscale modelling framework not only helps successfully bypass phenomenological constitutive assumptions as required in conventional modelling approaches but also facilitates effective cross-scale interpretation and understanding of soil behaviour. We examine the classical retaining wall and footing problems by this method and demonstrate that the simulating results can be well validated and verified by their analytical solutions. Furthermore, the study sheds novel multiscale insights into these classical problems and offers a new tool for geotechnical engineers to design and analyse geotechnical applications based directly upon particle-level information of soils. Copyright © 2015 John Wiley & Sons, Ltd.

The paper presents the finite volume formulation and numerical solution of finite strain one-dimensional consolidation equation. The equation used in this study utilises a nonlinear continuum representation of consolidation with varying compressibility and hydraulic conductivity and thus inherits the material and geometric nonlinearity. Time-marching explicit scheme has been used to achieve transient solutions. The nonlinear terms have been evaluated with the known previous time step value of the independent variable, that is, void ratio. Three-point quadratic interpolation function of Lagrangian family has been used to evaluate the face values at discrete control volumes. It has been shown that the numerical solution is stable and convergent for the general practical cases of consolidation. Performance of the numerical scheme has been evaluated by comparing the results with an analytical solution and with the piecewise piecewise-linear finite difference numerical model. The approach seems to work well and offers excellent potential for simulating finite strain consolidation. Further, the parametric study has been performed on soft organic clays, and the influence of various parameters on the time ate consolidation characteristics of the soil is shown. Copyright © 2015 John Wiley & Sons, Ltd.

An analytical solution of the plane strain problem of the deformation of a homogeneous, isotropic, poroelastic layer of uniform thickness overlying a homogeneous, isotropic, elastic half-space due to two-dimensional seismic sources buried in the elastic half-space has been obtained. The integral expressions for the displacements, stresses and pore pressure have been obtained using the stress function approach by applying suitable boundary conditions at the free surface and the interface. The solution obtained is in the Laplace–Fourier transform domain. The case of a vertical dip-slip line dislocation for the oceanic crust model of Earth is studied in detail. Schapery's formula is used for the Laplace inversion and the extended Simpson's formula for the Fourier inversion. Diffusion of pore pressure in the layer is studied numerically. Contour maps showing the pore pressure in the poroelastic layer have been plotted. The effect of the compressibility of the solid and fluid constituents on pore pressure has also been studied. Copyright © 2015 John Wiley & Sons, Ltd.

The mechanical properties of calcarenites are known to be significantly affected by water saturation: both stiffness and strength decrease for wetting in the short term and for chemical dissolution in the long term. Both processes mainly affect bonds among grains: immediately after inundation depositional bonds fall in suspension, whereas diagenetic bonds dissolve more slowly. In this paper, the authors started from the micro-structural analysis of the weathering processes to conceive a strain hardening hydro-chemo-mechanical coupled elastoplastic constitutive model. The concept of extended hardening rules is here enriched: weathering functions have been determined by employing a micro to macro simplified upscaling procedure. Chemical damage is incorporated into the formulation by means of a scalar damage function. Its evolution is also described by using a multiscale approach. A new term is added to the strain rate tensor in order to incorporate the dissolution induced chemical deformations developing once the soft rock is turned into a granular material. A calibration procedure for the constitutive parameters is suggested, and the model is validated by using both coupled and uncoupled chemo-mechanical experimental test results. Copyright © 2015 John Wiley & Sons, Ltd.

This paper presents a numerical solution for the analysis of the axisymmetric thermo-elastic problem in transversely isotropic material due to a buried heat source by means of extended precise integral method. By virtue of the Laplace–Hankel transform applied into the basic governing equations, an ordinary differential matrix equation is achieved, which describes the relationship between the generalized stresses and displacements in transformed domain. An extended precise integration method is introduced to solve the aforementioned matrix equation, and the actual solution in the physical domain is acquired by inverting the Laplace–Hankel transform. Numerical examples are carried out to demonstrate the accuracy of the proposed method and elucidate the influence of the character of transverse isotropy, the anisotropy of linear expansion coefficient, the anisotropy of thermal diffusivity, and medium's stratification on the thermo-elastic response. Copyright © 2015 John Wiley & Sons, Ltd.

A time-domain viscous-spring transmitting boundary is presented for transient dynamic analysis of saturated poroelastic media with linear elastic and isotropic properties. The *u–U* formulation of Biot equation in cylindrical coordinate is adopted in the derivation. By this general viscous-spring boundary, the effective stress and pore fluid pressure on the truncated boundary of the computational area are replaced by a set of continuously distributed spring and dashpot elements, of which the parameters are defined assuming an infinite permeability and considering the two dilatational waves. Numerical examples demonstrate good absorption of both the two cylindrical dilatational waves by the proposed ‘drained’ boundary. For general two-dimensional wave propagation problems, acceptable accuracy can still be achieved by setting the proposed boundary relatively far away from the scatter. Numerical comparison shows that the results obtained by using this boundary are more accurate for all permeability values than those by the traditional viscous-spring or viscous boundaries established for *u–U* formulation. Copyright © 2015 John Wiley & Sons, Ltd.

Behavior of unsaturated soils is influenced by many factors, and the influences of these factors are usually coupled together. Suction-controlled triaxial (SCTX) tests are considered to allow researchers to investigate influences of individual variables on unsaturated soils under specified stress path with controls of stresses, pore water, and air pressures. In the past 50 years, SCTX testing method has been established as a standard approach to characterize constitutive behavior of unsaturated soils. Most important concepts for modern unsaturated soil mechanics were developed upon results from the SCTX tests. Among these, one of the most important contributions in the constitutive modeling of elasto-plastic behavior for unsaturated soils is the Barcelona basic model (BBM) proposed by Alonso *et al.* in 1990. The BBM successfully explained many features of unsaturated soils and received extensive acceptance.

However, the SCTX tests are designed based upon the divide-and-conquer approach in which an implicit assumption is used: soil behavior is stress-path independent. However, it is well-established that unsaturated soil behavior is elasto-plastic and stress-path dependent. It is found that the SCTX tests in fact cannot control the stress path of an unsaturated soil during loading. This incapability, in combination with complicated loading/collapse behavior of unsaturated soils, makes the SCTX tests for characterizing unsaturated soil questionable. This paper discusses the limitations of the SCTX tests in the characterization of unsaturated soils. A possible solution to the problem was proposed based on a newly developed modified state surface approach. The discussions are limited for isotropic conditions. Copyright © 2015 John Wiley & Sons, Ltd.

Micromechanical characterization of multiple-porosity and multiple-permeability fluid-saturated porous materials from the properties of their single-porosity constituents is, to date, an open problem in our poromechanics society. This paper offers an in-depth view to this problem by considering the thermodynamic potential energy density, consistent with Biot's original definition, together with the general thought experiment, which allows for independent control of the sample's confining stress and distinct fluid pore pressures within its individual porosity networks. The complete set of well-known poroelastic constants, namely, Biot–Willis effective stress, Skempton's pore pressure, and specific storage coefficients, as well as drained, undrained, and Biot moduli for a fluid-saturated porous material, is herein identified with the reformulated theory. In particular, Gassmann relation for the bulk compressibility of the fluid-saturated material is accordingly upgraded to the case being addressed in this study.

The practical implications of the theory are showcased through a class of analytical solutions to the time-dependent poroelastic responses of shale to compression, when the hierarchical structure of its porous networks are accounted for at different levels of complexity and inter-porosity exchange effects. For this purpose, the laboratory setup of a quasi-2D compression test is considered, in which disk-shaped fluid-saturated samples of shale are allowed to drain laterally, while being sealed and confined from the top and bottom. A general closed-form solution to this problem is derived in the Laplace space, and the inverse numerical results for the cases of single-porosity, double-porosity, triple-porosity, and quadruple-porosity shale are discussed in the time domain. Copyright © 2015 John Wiley & Sons, Ltd.

The paper deals with the numerical solution of Biot's equations of coupled consolidation obtained by a mixed formulation combining continuous Galerkin finite-element and multipoint flux approximation finite-volume methods. The solution algorithm relies on the recently developed fixed-stress solution scheme, in which first the flow problem and then the mechanical one are addressed iteratively. We show that the algorithm can be interpreted as a particular block triangular preconditioning strategy applied within a Richardson iteration. The key component to the success of the preconditioner is the sparse approximation to the Schur complement based on a pressure space mass matrix scaled by a weighting factor that depends element-wise on the inverse of a suitable bulk modulus. The accuracy of the method is assessed, making use of well-known analytical solutions from the literature. Numerical results demonstrate robustness and low computational cost of the fixed-stress scheme in accurately capturing the two-way coupling between deformation and pressure. Copyright © 2015 John Wiley & Sons, Ltd.

A constitutive model for simulation of the behavior of unsaturated interfaces is presented here. The model is an extension of an existing critical state compatible interface model for dry and saturated interfaces that was already proposed by one of the authors [Lashkari, A. 2013. Int. J. Numer. Anal. Meth. Geomech. **37**(8): 904–931]. For a proper simulation of the behavior of partially saturated interfaces, the extended model is formulated in terms of two pairs of work conjugate stress–strain-like variables. The modified model simulations are compared with the existing data of dry, unsaturated, and saturated interfaces. For each interface type, it is shown that the proposed model can capture the essential elements of the behavior using a unique set of parameters. Copyright © 2015 John Wiley & Sons, Ltd.

The numerical simulation of rapid landslides is quite complex mainly because constitutive models capable of simulating the mechanical behaviour of granular materials in the pre-collapse and post-collapse regimes are still missing. The goal of this paper is to introduce a constitutive model capable of capturing the response of dry granular flows from quasi-static to dynamic conditions, in particular when the material experiences a sort of solid-to-fluid phase transition. An ideal assembly of identical spheres under simple shear conditions is considered. In the constitutive model, void ratio and granular temperature have been chosen as state variables, and both shear and normal stresses are computed as the sum of two contributions: the quasi-static one and the collisional one. The former is determined by using a perfect elasto-plastic model including the critical state concept, while the latter is derived from the kinetic theory of granular gases. The evolution of the granular temperature, fundamentally governing the material phase transition, is obtained by imposing the kinetic fluctuating energy balance. The constitutive relationship has been integrated, under both constant pressure and constant volume conditions, and the influence of shear strain rate, initial void ratio and normal pressure on the mechanical response has been investigated. Copyright © 2015 John Wiley & Sons, Ltd.

An analytical approach using the three-dimensional displacement of a soil is investigated to provide analytical solutions of the horizontal response of a circular pile subjected to lateral soil movements in nonhomogeneous soil. The lateral stiffness coefficient of the pile shaft in nonhomogeneous soil is derived from the rocking stiffness coefficient that is obtained from the analytical solution, taking into account the three-dimensional displacement represented in terms of scalar potentials in the elastic three-dimensional analysis. The relationship between horizontal displacement, rotation, moment, and shear force of a pile subjected to lateral soil movements in nonhomogeneous soil is obtainable in the form of the recurrence equation. For the relationship between the lateral pressure and the horizontal displacement, it is assumed that the behavior is linear elastic up to lateral soil yield, and the lateral pressure is constant under the lateral soil yield. The interaction factors between piles subjected to both lateral load and moment are calculated, taking into account the lateral soil movement. The formulation of the lateral displacement and rotation of the pile base subjected to lateral loads in nonhomogeneous soils is presented by taking into account the Mindlin equation and the equivalent thickness for soil layers in the equivalent elastic method. For lateral movement, lateral pressure, bending moment, and interaction factors, there are small differences between results obtained from the 1-D and the 3-D displacement methods except a very flexible pile. Copyright © 2015 John Wiley & Sons, Ltd.

This paper presents a superposition method expanded for computing impedance functions (IFs) of inclined-pile groups. Closed-form solutions for obtaining horizontal, vertical, and rocking IFs, estimated by using pile-to-pile interaction factors, are proposed. IFs of solitary inclined piles, crossed IFs, and explicit incorporation of compatibility conditions for pile-head movements are also appropriately taken into consideration. All of these factors should be known in advance and will be computed and shown for the most relevant cases. The accuracy of the proposed closed-form solutions is verified for 2 × 2 and 3 × 3 square inclined-pile groups embedded in an isotropic viscoelastic homogeneous half-space soil medium, with hysteretic damping. The pile-to-pile interaction factors are computed by means of a three-dimensional time-harmonic boundary elements–finite elements coupling formulation. The results indicate that the IFs obtained from the proposed method are in good agreement with those obtained from the coupling formulation. Furthermore, crossed vertical-rocking IFs of solitary piles need to be appropriately considered for obtaining rocking IFs when the number of piles is small. Copyright © 2015 John Wiley & Sons, Ltd.

Compressive loading of granular materials causes inter-particle forces to develop and evolve into force chains that propagate through the granular body. At high-applied compressive stresses, inter-particle forces will be large enough to cause particle fracture, affecting the constitutive behavior of granular materials. The first step to modeling particle fracture within force chains in granular mass is to understand and model the fracture of a single particle using actual three-dimensional (3D) particle shape. In this paper, the fracture mode of individual silica sand particles was captured using 3D x-ray radiography and Synchrotron Micro-computed Tomography (SMT) during in situ compression experiments. The SMT images were used to reconstruct particle surfaces through image processing techniques. Particle surface was then imported into Abaqus finite element (FE) software where the experimental loading setup was modeled using the extended finite element method (XFEM) where particle fracture was compared to experimental fracture mode viewed in radiograph images that were acquired during experimental loading. Load-displacement relationships of the FE analysis were also compared with experimental measurements. 3D FE modeling of particle fracture offers an excellent tool to map stress distribution and monitors crack initiation and propagation within individual sand particles. Copyright © 2015 John Wiley & Sons, Ltd.

In the numerical modeling of fluid flow in heterogeneous geological media, large material contrasts associated with complexly intersected material interfaces are challenging, not only related to mesh discretization but also for the accurate realization of the corresponding boundary constraints. To address these challenges, we developed a discontinuous approach for modeling fluid flow in heterogeneous media using the numerical manifold method (NMM) and the Lagrange multiplier method (LMM) for modeling boundary constraints. The advantages of NMM include meshing efficiency with fixed mathematical grids (covers), the convenience of increasing the approximation precision, and the high integration precision provided by simplex integration. In this discontinuous approach, the elements intersected by material interfaces are divided into different elements and linked together using the LMM. We derive and compare different forms of LMMs and arrive at a new LMM that is efficient in terms of not requiring additional Lagrange multiplier topology, yet stringently derived by physical principles, and accurate in numerical performance. To demonstrate the accuracy and efficiency of the NMM with the developed LMM for boundary constraints, we simulate a number of verification and demonstration examples, involving a Dirichlet boundary condition and dense and intersected material interfaces. Last, we applied the developed model for modeling fluid flow in heterogeneous media with several material zones containing a fault and an opening. We show that the developed discontinuous approach is very suitable for modeling fluid flow in strongly heterogeneous media with good accuracy for large material contrasts, complex Dirichlet boundary conditions, or complexly intersected material interfaces. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, the onset of mechanical instability in time-sensitive elasto-viscoplastic solids is theoretically analyzed at the constitutive level and associated with the occurrence of ‘spontaneous accelerations’ under stationary external perturbations. For this purpose, a second-order form of Perzyna's constitutive equations is first derived by time differentiation, and a sufficient stability condition is identified for general mixed loading programs. These loading conditions are in fact the most general in both laboratory tests and real boundary value problems, where a combination of certain stress and strain components is known/prescribed.

The theoretical analysis leads to find precise stability limits in terms of material hardening modulus. In the case of constitutive relationships with isotropic strain-hardening, no instabilities are possible while the hardening modulus is larger than the so-called ‘controllability modulus’ defined for (inviscid) elasto-plastic materials. It is also shown that the current stress/strain rate may also directly influence the occurrence of elasto-viscoplastic instability, which is at variance with elasto-plastic inviscid media. Copyright © 2015 John Wiley & Sons, Ltd.

This paper presents an analytical solution for the lateral dynamic response of a pipe pile in a saturated soil layer. The wave propagations in the saturated soil and the pipe pile are simulated by Biot's three-dimensional poroelastic theory and one-dimensional elastic theory, respectively. The governing equations of soil are solved directly without introducing potential functions. The displacement response and dynamic impedances of the pipe pile are obtained based on the continuous conditions between the pipe pile and both the outer and inner soil. A comparison with an existing solution is performed to verify the proposed solution. Selected numerical results for the lateral dynamic responses and impedances of the pipe pile are presented to reveal the lateral vibration characteristics of the pile-soil system. Copyright © 2015 John Wiley & Sons, Ltd.

Computational fluid dynamics and discrete element method (CFD–DEM) is extended with the volume of fluid (VOF) method to model free-surface flows. The fluid is described on coarse CFD grids by solving locally averaged Navier–Stokes equations, and particles are modelled individually in DEM. Fluid–particle interactions are achieved by exchanging information between DEM and CFD. An advection equation is applied to solve the phase fraction of liquid, in the spirit of VOF, to capture the dynamics of free fluid surface. It also allows inter-phase volume replacements between the fluid and solid particles. Further, as the size ratio (SR) of fluid cell to particle diameter is limited (i.e. no less than 4) in coarse-grid CFD–DEM, a *porous sphere* method is adopted to permit a wider range of particle size without sacrificing the resolution of fluid grids. It makes use of more fluid cells to calculate local porosities. The developed solver (*cfdemSolverVOF*) is validated in different cases. A dam break case validates the CFD-component and VOF-component. Particle sedimentation tests validate the CFD–DEM interaction at various Reynolds numbers. Water-level rising tests validate the volume exchange among phases. The *porous sphere* model is validated in both static and dynamic situations. Sensitivity analyses show that the SR can be reduced to 1 using the *porous sphere* approach, with the accuracy of analyses maintained. This allows more details of the fluid phase to be revealed in the analyses and enhances the applicability of the proposed model to geotechnical problems, where a highly dynamic fluid velocity and a wide range of particle sizes are encountered. Copyright © 2015 John Wiley & Sons, Ltd.

In the present paper, a constitutive model for the description of the dissipation in the concrete is provided. The theoretical description is based on a micromorphic model in which the microstructure is constituted by a kinematical scalar descriptor *φ* whose time derivative is linked to a dissipative potential. The scalar *φ* can be interpreted as the relative displacement between two opposite faces of the microcracks, and our physical interpretation of dissipation is indeed linked to the friction force (in a mixed Coulomb-type and viscous-type behavior) between them. To evaluate the effects of bending on the dissipation, the 3D model is then reduced by means of standard Saint-Venant's procedure in case of combined compression and bending over a cylindrical domain. A qualitative analysis of the reduced ODEs model is then provided. Numerical results showing comparison between different types of dissipative force and between pure compression and combined compression and bending are included in a dedicated section. Finally, the proposed model and our physical interpretation of the dissipation are supported by some experimental data concerning standard concrete and a concrete enriched by adding to the mixture a filler constituted by micro-particles capable of improving the dissipative behavior of the material. Measured data show very good fit with our theoretical previsions and provide a sufficiently sound basis for further deepening of the theoretical description of the considered phenomena. Copyright © 2015 John Wiley & Sons, Ltd.

A finite element algorithm for frictionless contact problems in a two-phase saturated porous medium, considering finite deformation and inertia effects, has been formulated and implemented in a finite element programme. The mechanical behaviour of the saturated porous medium is predicted using mixture theory, which models the dynamic advection of fluids through a fully saturated porous solid matrix. The resulting mixed formulation predicts all field variables including the solid displacement, pore fluid pressure and Darcy velocity of the pore fluid. The contact constraints arising from the requirement for continuity of the contact traction, as well as the fluid flow across the contact interface, are enforced using a penalty approach that is regularised with an augmented Lagrangian method. The contact formulation is based on a mortar segment-to-segment scheme that allows the interpolation functions of the contact elements to be of order *N*. The main thrust of this paper is therefore how to deal with contact interfaces in problems that involve both dynamics and consolidation and possibly large deformations of porous media. The numerical algorithm is first verified using several illustrative examples. This algorithm is then employed to solve a pipe-seabed interaction problem, involving large deformations and dynamic effects, and the results of the analysis are also compared with those obtained using a node-to-segment contact algorithm. The results of this study indicate that the proposed method is able to solve the highly nonlinear problem of dynamic soil–structure interaction when coupled with pore water pressures and Darcy velocity. Copyright © 2015 John Wiley & Sons, Ltd.

Integrating ground heat exchanger elements into concrete piles is now considered as an efficient energy solution for heating/cooling of buildings. In addition to the static load of buildings, the concrete piles also undergo a cycle of thermal deformation. In the case of single energy pile, calculation methods already exist and permit to perform a proper geotechnical design. In the case of energy pile group, the thermo-mechanical interactions within the group are more complex. Very few experimental results on the energy pile group are available so that numerical analysis can be an interesting way to provide complementary results about their behavior. This paper deals with a numerical analysis including a comparison between a single energy pile and an energy pile group with different boundary conditions at the pile head. In order to take into account the stress reversal induced by the thermal expansions and contractions, a cyclic elastoplastic constitutive model is introduced at the soil–pile interface. The analysis aims to give some insights about the long-term cyclic interaction mechanisms in the energy pile group. Based on this qualitative study, some guidance can be brought for the design of energy piles in the case where group effects should be considered. Copyright © 2015 John Wiley & Sons, Ltd.

The displacement discontinuity method (DDM) is frequently used in geothermal and petroleum applications for modeling the behavior of fractures in linear-elastic rocks. The DDM requires O(*N*^{2}) memory and O(*N*^{3}) floating point operations (where *N* is the number of unknowns) to construct the coefficient matrix and solve the linear system of equations by direct methods. Therefore, the conventional implementation of the DDM is not computationally efficient for very large systems of cracks, often limiting its application to small-scale problems. This work presents an approach for solving large-scale fracture problems using the fast multipole method (FMM). The approach uses both the DDM and a kernel-independent version of the FMM along with a preconditioned generalized minimal residual algorithm to accelerate the solution of linear systems of equations using desktop computers. Using the fundamental solutions for constant displacement discontinuity in a two-dimensional elastic medium, several numerical examples involving fracture networks representing fractured reservoirs are treated. Numerical results show good agreement with analytical solutions and demonstrate the efficiency of the FMM implementation of the DDM for large-scale simulations. Copyright © 2015 John Wiley & Sons, Ltd.

Sand production is a complex physical process that depends on the external stress and flow rate conditions as well as on the state of the material. Models developed for the prediction of sand production are usually solved numerically because of the complexity of the governing equations. Testing of new sand production models can very well be performed through calibration with laboratory experiments, which by construction possess geometric symmetry facilitating explicit mathematical analysis. We introduce an erosion model that is built upon the physics (poro-mechanical coupling of the fluid-solid system) usually incorporated in erosion models for the prediction of sand production. Around this model, we set up a mathematical framework in which sand production models because of erosion can be tested and calibrated without having to resort to complex numerical work or specialised software. The model is validated by data of volumetric sand production from a hollow cylinder test on synthetic sandstone. Generalisations of the model, which are naturally incorporated in the same framework and have useful phenomenological features, are discussed. Copyright © 2015 John Wiley & Sons, Ltd.

Modeling the flow in highly fractured porous media by finite element method (FEM) has met two difficulties: mesh generation for fractured domains and a rigorous formulation of the flow problem accounting for fracture/matrix, fracture/fracture, and fracture/boundary fluid mass exchanges. Based on the recent theoretical progress for mass balance conditions in multifractured porous bodies, the governing equations for coupled flow and deformation in these bodies are first established in this paper. A weak formulation for this problem is then established allowing to build a FEM. Taking benefit from recent development of mesh-generating tools for fractured media, this weak formulation has been implemented in a numerical code and applied to some typical problems of hydromechanical coupling in fractured porous media. It is shown that in this way, the FEM that has proved its efficiency to model hydromechanical phenomena in porous media is extended with all its performances (calculation time, couplings, and nonlinearities) to fractured porous media. Copyright © 2015 John Wiley & Sons, Ltd.

This paper is devoted to multi-scale modeling of elastic–plastic deformation of a class of geomaterials with a polycrystalline microstructure. We have extended and improved the simplified polycrystalline model presented in [Zeng T. *et al.*, 2014. *Mech. Mater.***69**(1):132–145]. A rigorous and fully consistent self-consistent (SC) scheme is proposed to describe the interaction among plastic mineral grains. We have also deeply discussed the numerical issues related to the numerical implementation of the proposed micromechanical model. The efficiency of the proposed model and the related numerical procedure is evaluated in several representative cases. We have compared the numerical results respectively obtained from the fully SC model and two simplified ones. It is found that the SC model produces a softer stress–strain response than that of the simplified models. The comparisons between the estimation of overall behavior of a granite in different loading conditions and experimental data are also conducted. Copyright © 2015 John Wiley & Sons, Ltd.

Cavity expansion theory assists in the interpretation of *in situ* tests including the cone penetration test and pressuremeter test. In this paper, a cavity expansion analysis is presented for unsaturated silty sand exhibiting hydraulic hysteresis. The similarity technique is used in the analysis. The soil stress–strain behaviour is described by a bounding surface plasticity model. Results of oedometric compression tests, isotropic compression tests and triaxial shear tests for both saturated and unsaturated states are used to calibrate the model. The void ratio, suction, degree of saturation and effective stress are fully coupled in the analysis. The influence of where the initial hydraulic state is located on the soil–water characteristic curve on the cavity wall pressure is investigated and found to be significant. Also, the effects of three different drainage conditions (constant suction, constant moisture content and constant contribution of suction to the effective stress) on cavity wall pressure are studied. It is found that the drainage condition in which the contribution of suction to the effective stress is constant offers a good approximation to the other two. This may simplify interpretation of *in situ* tests. When testing occurs quickly, meaning a constant moisture content condition prevails, a constant contribution of suction condition can be assumed without loss of significant accuracy. The contribution of suction assumed in the interpretation can be taken as being equal to the *in situ* value, although this discovery may not be applicable to all soil types, constitutive models and soil–water characteristic curves. Copyright © 2015 John Wiley & Sons, Ltd.

This paper presents an analytical layer element solution to axisymmetric thermal consolidation of multilayered porous thermoelastic media containing a deep buried heat source. By applying the Laplace–Hankel transform to the state variables involved in the basic governing equations of porous thermoelasticity, the analytical layer elements that describe the relationship between the transformed generalized stresses and displacements of a finite layer and a half-space are derived. The global stiffness matrix equation is obtained by assembling the interrelated layer elements, and the real solutions in the physical domain are achieved by numerical inversion of the Laplace–Hankel transform after obtaining the solutions in the transformed domain. Finally, numerical calculations are performed to demonstrate the accuracy of this method and to investigate the influence of heat source's types, layering, and the porous thermoelastic material parameters on thermal consolidation behavior. Copyright © 2015 John Wiley & Sons, Ltd.

The Sandia GeoModel is a continuum elastoplastic constitutive model that captures many features of the mechanical response for geological materials over a wide range of porosities and strain rates. Among the specific features incorporated into the formulation are a smooth compression cap, isotropic/kinematic hardening, nonlinear pressure dependence, strength differential effect, and rate sensitivity. This study attempts to provide enhancements regarding computational tractability, domain of applicability, and robustness of the model. A new functional form is presented for the yield and plastic potential functions. This reformulation renders a more accurate, robust, and efficient model as it eliminates spurious solutions attributed to the original form. In addition, we achieve a high-performance implementation, because the local iterative method is allowed to recast residual vectors with a uniform dimensionality. The model is also furnished with a smooth, elliptical tension cap to account for the tensile yielding. Moreover, an efficient algorithm is introduced, which decreases the computational cost by differentiating the updated shear yield surface from the cap surface based on the trial relative stress state. Finally, various numerical examples including a large-scale boundary value problem are presented to demonstrate the fidelity of the modified model and to analyze its numerical performance. Copyright © 2015 John Wiley & Sons, Ltd.

This paper presents a simple three-dimensional (3D) Distinct Element Method (DEM) for numerical simulation of the mechanical behavior of bonded sands. First, a series of micro-mechanical tests on a pair of aluminum rods glued together by cement with different bond sizes were performed to obtain the contact mechanical responses of ideally bonded granular material. Second, a 3D bond contact model, which takes into account the influences of bond sizes, was established by extending the obtained 2D experimental results to 3D case. Then, a DEM incorporating the new contact model was employed to perform a set of drained triaxial compression tests on the DEM bonded specimens with different cement contents under different confining pressures. Finally, the mechanical behavior of the bonded specimens was compared with the available experimental results. The results show that the DEM incorporating the simple 3D bond contact model is able to capture the main mechanical behavior of bonded sands. The bonded specimen with higher cement content under lower confining pressure exhibits more pronounced strain softening and shear dilatancy. The peak and residual strengths, the apparent cohesion and peak/residual friction angles, and the position and slope of the critical state line increase with increase in cement content. Microscopically, bond breakage starts when the system starts to dilate and the maximum rate of bond breakage coincides with the maximum rate of dilation. Bond breakage is primarily due to tension-shear failure and the percentage of such failures is independent of both confining pressure and cement content. Copyright © 2015 John Wiley & Sons, Ltd.

In contrast to the traditional approach that computes the reliability index in the uncorrelated standard normal space (u-space), the reliability analysis that is simply realized in the original space (x-space, non-Gaussian type) would be more efficient for practical use, for example, with the Low and Tang's constrained optimization approach. On the other hand, a variant of Hasofer, Lind, Rackwits and Fiessler algorithm for first-order reliability method is derived in this paper. Also, the new algorithm is simply formulated in x-space and requires neither transformation of the random variables nor optimization tools. The algorithm is particularly useful for reliability analysis involving correlated non-Gaussian random variables subjected to implicit limit state function. The algorithm is first verified using a simple example with closed-form solution. With the aid of numerical differentiation analysis in x-space, it is then illustrated for a strut with complex support and for an earth slope with multiple failure modes, both cases involving implicit limit state surfaces. Copyright © 2015 John Wiley & Sons, Ltd.

Methane hydrate (MH) is a new energy resource in the 21^{st} century. But the dissociation of MH from sediments during the MH exploration or oil/gas exploration under a hydrate layer accompanied by the softening of soils and formation of excess pore gas pressure may lead to ground failures and environmental disasters. In this study, experiments on modeling the weakening and failure of the sediment by heat-induced dissociation of tetrahydrofuran (THF) hydrate were presented. The failure mode of gas outburst was observed. Gas outbursts is a process where gas and soils in hydrate-dissociation zone burst out after the continuous skeleton of over-layer is fractured during the expansion of the dissociation zone and the formation of gas zone and excess pore gas pressure. An analytical method is presented by decoupling heat transfer and soil deformation. The geometrical and mechanical similarities for gas outburst are obtained. An empirical criterion for the occurrence of outburst is proposed using the theory of thermal conduction, rigid plastic mechanics, and the experimental data. Copyright © 2015 John Wiley & Sons, Ltd.

When an underwater tunnel is excavated, the groundwater may flow into the tunnel. The seepage forces consequently induced can have important effects on the effective stresses around the tunnel. Moreover, the influences of the free surface of a shallow underwater tunnel should also be considered. In this research, an analytical solution is presented to calculate the seepage-induced effective stresses around a shallow underwater tunnel in an elastic half plane. The solution uses the complex variable method and consists of conformally mapping the half plane with a hole onto a transformed circular ring. The coefficients of the various terms in the Laurent series expansions of the stress functions in the transformed region can be obtained from the boundary conditions. The total stress distribution around a shallow underwater tunnel can be calculated by the potentials in the half plane. The effective stress can be obtained by subtracting the pore pressure from the total stress. The analytical solution is validated by numerical simulations and can be used to perform both the short-term and long-term analyses. By using the proposed solution, it is found that the circumferential effective stresses around the tunnel increase greatly because of seepage, and they increase with the increase of water depth in both the undrained and drained conditions. Copyright © 2015 John Wiley & Sons, Ltd.

The purpose of this paper is to analyse the development and the evolution of the fracture process zone during fracture and damage in quasi-brittle materials. A model taking into account the material details at the mesoscale is used to describe the failure process at the scale of the heterogeneities. This model is used to compute histograms of the relative distances between damaged points. These numerical results are compared with experimental data, where the damage evolution is monitored using acoustic emissions. Histograms of the relative distances between damage events in the numerical calculations and acoustic events in the experiments exhibit good agreement. It is shown that the mesoscale model provides relevant information from the point of view of both global responses and the local failure process. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, the problem of propagation of localized deformation associated with formation of macrocracks/shear bands is studied in both tensile and compressive regimes. The main focus here is on enhancement of the constitutive law with embedded discontinuity to provide a discrete representation of the localization phenomenon. This has been accomplished by revising the formulation and coupling it with the level-set method for tracing the propagation path. Extensive numerical studies are conducted involving various fracture modes, ranging from brittle to frictional, and the results are compared with the experimental data as well as those obtained using XFEM methodology. Copyright © 2015 John Wiley & Sons, Ltd.

This paper presents a simplified finite element analysis technique, the ‘Press-Replace’ technique, to model pile penetration problems in geotechnical engineering, particularly, pile jacking. The method is employed in standard finite element analysis software. The method involves a straining and a consequent geometry update phase. First, a cone penetration test in (undrained) clay is modelled and compared with the results of analytical, semi-analytical and more advanced finite element techniques. The model sensitivity for the step size and mesh is investigated using a hypoplastic constitutive model. An optimum way of modelling based on the numerical performance is shown. Copyright © 2015 John Wiley & Sons, Ltd.

This paper discusses the excess pore-air and pore-water pressure dissipations and the average degree of consolidation in the 2D plane strain consolidation of an unsaturated soil stratum using eigenfunction expansion and Laplace transformation techniques. In this study, the application of a constant external loading on a soil surface is assumed to immediately generate uniformly or linearly distributed initial excess pore pressures. The general solutions consisting of eigenfunctions and eigenvalues are first proposed. The Laplace transform is then applied to convert the time variable t in partial differential equations into the Laplace complex argument s. Once the domain is obtained, a simplified set of equations with variable s can be achieved. The final analytical solutions can be computed by taking a Laplace inverse. The proposed equations predict the two-dimensional consolidation behaviour of an unsaturated soil stratum capturing the uniformly and linearly distributed initial excess pore pressures. This study investigates the effects of isotropic and anisotropic permeability conditions on variations of excess pore pressures and the average degree of consolidation. Additionally, isochrones of excess pore pressures along vertical and horizontal directions are presented. It is found that the initial distribution of pore pressures, varying with depth, results in considerable effects on the pore-water pressure dissipation rate whilst it has insignificant effects on the excess pore-air pressure dissipation rate. Copyright © 2015 John Wiley & Sons, Ltd.

We present hereafter the formulation of a Timoshenko finite element straight beam with internal degrees of freedom, suitable for nonlinear material problems in geomechanics (e.g., beam type structures and deep pile foundations). Cubic shape functions are used for the transverse displacements and quadratic for the rotations. The element is free of shear locking, and we prove that one element is able to predict the exact tip displacements for any complex distributed loadings and any suitable boundary conditions. After the presentation of the virtual power and the weak form formulations, the construction of the elementary stiffness matrix is detailed. The analytical results of the static condensation method are provided. It is also proven that the element introduced by Friedman and Kosmatka in , with shape functions depending on material properties, is derived from the new beam element. Validation is provided using linear and material nonlinear applications (reinforced concrete column under cyclic loading) in the context of a multifiber beam formulation. Copyright © 2015 John Wiley & Sons, Ltd.

Large deformations and discontinuous problems can be calculated using the discontinuous deformation analysis (DDA) method by solving time steps, and this method is suitable for simulating the seismic dynamic response of engineering rock mass structures. However, the boundary setting must be carefully analyzed. In this paper, four boundary settings for the DDA method are investigated. First, the contributions to the DDA equations for nonreflecting boundaries (including the viscous boundary and the viscoelastic boundary) are deduced based on the Newmark method. Second, a free-field boundary is introduced in the DDA method with boundary grid generation and coupling calculation algorithms to accurately simulate external source wave motion, such as earthquakes. Third, seismic input boundary treatments are intensively examined, and the force input method is introduced based on nonreflecting boundaries. Finally, the static-dynamic unified boundary is implemented to ensure consistent boundary transformation. The boundary setting method in the DDA method is discussed, and the suggested treatments are used to analyze the seismic dynamic response of underground caverns. Copyright © 2015 John Wiley & Sons, Ltd.

Geomechanical analysis needed to treat anisotropic (specifically trigonal symmetry) grains of carbonates jumbled together to form an overall isotropic polycrystalline poroelastic material is summarized. Poroelastic effects enter the problem via fractured anisotropic solids saturated by fluids. Copyright © 2014 John Wiley & Sons, Ltd.

Shrink–swell soils can cause distresses in buildings, and every year, the economic loss associated with this problem is huge. This paper presents a comprehensive system for simulating the soil–foundation–building system and its response to daily weather conditions. Weather data include rainfall, solar radiation, air temperature, relative humidity, and wind speed, all of which are readily available from a local weather station or the Internet. These data are used to determine simulation flux boundary conditions. Different methods are proposed to simulate different boundary conditions: bare soil, trees, and vegetation. A coupled hydro-mechanical stress analysis is used to simulate the volume change of shrink–swell soils due to both mechanical stress and water content variations. Coupled hydro-mechanical stress-jointed elements are used to simulate the interaction between the soil and the slab, and general shell elements are used to simulate structural behavior. All the models are combined into one finite element program to predict the entire system's behavior. This paper first described the theory for the simulations. A site in Arlington, Texas, is then selected to demonstrate the application of the proposed system. Simulation results are shown, and a comparison between measured and predicted movements for four footings in Arlington, Texas, over a 2-year period is presented. Finally, a three-dimensional simulation is made for a virtual residential building on shrink–swell soils to identify the influence of various factors. Copyright © 2015 John Wiley & Sons, Ltd.

Accurate prediction of the interactions between the nonlinear soil skeleton and the pore fluid under loading plays a vital role in many geotechnical applications. It is therefore important to develop a numerical method that can effectively capture this nonlinear soil-pore fluid coupling effect. This paper presents the implementation of a new finite volume method code of poro-elasto-plasticity soil model. The model is formulated on the basis of Biot's consolidation theory and combined with a perfect plasticity Mohr-Coulomb constitutive relation. The governing equation system is discretized in a segregated manner, namely, those conventional linear and uncoupled terms are treated implicitly, while those nonlinear and coupled terms are treated explicitly by using any available values from previous time or iteration step. The *implicit–explicit* discretization leads to a linearized and decoupled algebraic system, which is solved using the fixed-point iteration method. Upon the convergence of the iterative method, fully nonlinear coupled solutions are obtained. Also explored in this paper is the special way of treating traction boundary in finite volume method compared with FEM. Finally, three numerical test cases are simulated to verify the implementation procedure. It is shown in the simulation results that the implemented solver is capable of and efficient at predicting reasonable soil responses with pore pressure coupling under different loading situations. Copyright © 2015 John Wiley & Sons, Ltd.

Fast closure of rock fractures has been commonly observed in the initial stage of fluid flow experiments at environmental temperatures under low or moderate normal stresses. To fully understand the mechanisms that drive this fast closure, the evolution of local stresses acting on contacting asperities on the fracture surfaces prior to fluid flow tests needs to be evaluated. In this study, we modeled numerically the asperity deformation and failure processes during initial normal loading, by adopting both elastic and elastic–plastic deformation models for the asperities on a real rock fracture with measured surface topography data, and estimated their impact on initial conditions for fluid flow test performed under laboratory conditions. Compared with the previous models that simulate the normal contact of a fracture as the approach of two rigid surfaces without deformations, our models of deformable asperities yielded smaller contact areas and higher stresses on contacting asperities at a given normal stress or normal displacement. The results show that the calculated local stresses were concentrated on the contacts of a few major asperities, resulting in crushing of asperity tips. With these higher contact stresses, however, the predicted closure rates by pressure solution are still several orders of magnitude lower than that of the experimental measurements at the initial stage of fluid flow test. This indicates that single pressure solution may not likely to be the principal compaction mechanism for this fast closure, and that the damages on contacting asperities that occur during the initial normal loading stage may play an important role. Copyright © 2015 John Wiley & Sons, Ltd.

The backfilling materials of borehole heat exchangers (BHE), particularly the grout material, must provide a suitable thermal contact and ensure durability to the induced thermal stresses because of the heat loading. In this paper, the thermal stresses that occurred in BHEs because of heat injection or extraction is investigated with an analytical solution of a hollow cylinder model that is adapted for time-dependent heat loading, the geometry of a BHE, and the thermo-mechanical properties of surrounding ground conditions. Firstly, the hollow cylinder model is solved with the considered boundary conditions in 2D plane stress. Secondly, the temperature differences at the inner and outer circles of the cylinder are evaluated with the heat line source models for continuous and discontinuous loading to observe the impact of the heat loading schedule. The developed analytical solution for thermal stress investigation is validated with numerical models. It is demonstrated that the analytical solutions agree well with numerical results for two types of BHE configurations (co-axial and single U-shaped pipes). Furthermore, the calculated maximum stresses are compared with the tensile strength of grout materials obtained from Brazilian tests. It is predicted that the thermal contraction of the grout, partially constrained by the surrounding rock, generates tensile stresses that may lead to cracking in the BHE. According to the results, the stiffness of rock has a primary role on the developed tensile stresses, and the relationship between the thermal conductivity of the ground and of the grout induces a proportional impact on the magnitude of thermal stresses. Copyright © 2015 John Wiley & Sons, Ltd.

This study presents two three-parameter failure criteria for cohesive-frictional materials based on the Mohr–Coulomb failure function. One proposed failure criterion can be linked to Mogi's empirical formula and incorporates the well-known Von-Mises, Drucker–Prager, and Linear Mogi criteria as special cases. Another one with smooth and convex cross sections contains a general Lode dependence in the deviatoric plane and includes the Matsuoka–Nakai and Lade–Duncan Lode dependences as special cases. The effect of the intermediate principal stress on the strength of the material can be taken into account in both criteria. The proposed criteria are numerically calibrated against polyaxial data sets of many different types of rocks and concrete_{.} The comparison results show that the performance of the proposed criteria is excellent, and the failure criterion with a general Lode dependence performs better than the other one for concrete. Copyright © 2015 John Wiley & Sons, Ltd.