In this study, a fully coupled thermo-hydro-mechanical (THM) theoretical model, called the TSINGHUA-THERMOSOIL model, which is applicable in saturated soils, was developed based on a non-equilibrium thermodynamic approach (see Groot S, Mazur P. Non-Equilibrium thermodynamics. Amsterdam: North-Holland Pub. Co., 1962). This model simplifies the mathematical description of a THM coupling system to the modeling of a group of migration coefficients and energy functions by means of the non-equilibrium thermodynamic theory. The important concept of granular entropy was introduced for describing the granular fluctuation and the resulting transient elasticity. Therefore, non-elastic mechanical constitutive relations that do not require such concepts as yield surface and flow rule are obtained. The physical process of bound pore water transforming into free water when the temperature rises was taken into consideration for describing the phenomenon of the volumetric compression of normally consolidated and slightly over-consolidated clays under non-isothermal conditions. The behavior of thermal consolidation and drained shear tests at different temperatures for saturated Kaolin clay with different over-consolidation ratios (OCR) were simulated using this model. The effect of the OCR value, the heating rate and temperature cycling on the thermal consolidation and the effect of the temperature on the shear strength were analyzed and compared with existing experimental data. The results showed that the model proposed in this paper can accurately describe the thermal consolidation and the basic law of shearing for saturated soil, which preliminarily proves the validity of the model. Copyright © 2016 John Wiley & Sons, Ltd.

This paper presents a biconcave bond model to investigate the effect of the cementation between grains on the mechanical behavior of rock. The proposed model considers the shape of the bonds among particles that have a biconcave cement form, based on observations of microscopic rock images. The general equations of the proposed model are based on Dvorkin theory. The accuracy and efficiency of the bond model is improved in three ways. After the biconcave bond model is implemented in the discrete element method software Particle Flow Code in 2 Dimensions, a series of numerical uniaxial compression tests were performed to investigate the relationships between the micro- to macro-parameters. The simulations revealed that the biconcave bond model reflects the effect of micro-parameters, such as the elastic modulus and Poisson's ratio of the cement, on the macroscopic deformation of cemented granular material. Variations in the bond geometry caused extremely diverse macro-mechanical behaviors. Experimental results concerning rock demonstrate that the biconcave bond model accurately captures the mechanical behavior of intact rock and supports an innovative method for investigating the relationships between the micro- and macro-parameters of cemented granular material. Copyright © 2016 John Wiley & Sons, Ltd.

An extensive literature on the shear behavior of continuum–particulate interfaces has been developed during the last four decades. However, relatively limited work regarding the behavior of interfaces under different loading conditions has been published. This paper presents a discrete element modeling study, along with comparisons from experimental data, of interface behavior under axial and torsional drained loading conditions. Detailed studies allow for links between micro-scale particle behavior and observed global response to be developed and for the latter to be evaluated in light of particle–particle and particle–continuum interactions. The results of this study indicate that axial and torsional interface shear induce inherently different loading conditions, as shown by the different failure envelopes, stress paths, and induced soil volume changes and deformations. Furthermore, the results presented in this paper indicate that particle-level mechanisms, such as particle rotations and contact slippage, play different roles in axial and torsional shear. Coordination number, polar histograms, particle displacements, particle rotations, and local void ratio measurements provide further insights into the fabric evolution, loading conditions, and failure mechanisms induced by these two shear modes. This study expands the current understanding of interface behavior and discusses potential improvements to geotechnical systems that leverage the characteristics of different imposed loading conditions. Copyright © 2016 John Wiley & Sons, Ltd.

Based on the Biot's poroelastic theory and using scalar potential functions both the ring load and point load displacement Green's functions for a transversely isotropic saturated porous full-space composed of an upper half-space, a finite thickness middle layer and a lower half-space is analytically presented for the first time. It is assumed that each region consists of a different transversely isotropic material. The equations of poroelastodymanics in terms of the solid displacements and the pore fluid pressure are uncoupled with the help of two scalar potential functions, so that the governing equations for the potential functions are either a second order wave equation or a repeated wave-heat transfer equation of sixth order. With the aid of Fourier expansion with respect to circumferential direction and Hankel integral transforms with respect to the radial direction in cylindrical coordinate system, the response is determined in the form of line integrals in the real space, followed by theorem of inverse Hankel integral transforms. The solutions degenerate to a single phase elastic material, and the results are compared with previous studies, where an excellent agreement may be observed with the results provided in the literature. Some examples of displacement Green's functions are finally given to illustrate the solution. Copyright © 2016 John Wiley & Sons, Ltd.

This paper presents the theoretical development and methodological motivation of a single surface anisotropic hyperplasticity model. The model extends the isotropic family of models developed by Coombs and Crouch by: (i) introducing anisotropic shearing into the yield surface, (ii) relating two of the material constants to a single physical quantity and (iii) using a more physically realistic pressure sensitive elastic free energy function. This model overcomes the difficulty of determining the constants of the isotropic two-parameter surface by analytically relating them to a single experimentally measurable physical quantity, namely the normalised hydrostatic position of the Critical State. This provides a model with a Critical State surface that is constant throughout the loading process, invariant of the level of anisotropy inherent in the yield envelope. The model is compared with experimental data from triaxial tests on Lower Cromer Till, contrasted against the SANIclay model and the recent model of Yang *et al*. (2015) as well as being compared with rarely considered experimental data from hollow cylinder tests on London Clay. Copyright © 2016 John Wiley & Sons, Ltd.

The stability of integration is essential to numerical simulations especially when solving nonlinear problems. In this work, a continuum damage mechanics model proposed by the first author is implemented with an integration method named cutting plane algorithm (CPA) to improve the robustness of the simulation. This integration method is one type of return mapping algorithm that bypasses the need for computing the gradients. We compare the current integration method with the previous direct method, and the result shows that the cutting plane algorithm exhibits excellent performance under large loading rate conditions. To enhance accuracy of the new method, a control procedure is utilized in the implementation of the algorithm based on error analysis. Thereafter, the theory of poromechanics is utilized with the damage model to account for the effects of fluid diffusion. Laboratory tests simulated with finite element method illustrate distinct behaviors of shale with different loading rates and indicate the development of microcrack propagation under triaxial compression. Copyright © 2016 John Wiley & Sons, Ltd.

Discontinuity layout optimization (DLO) is a recently presented topology optimization method for determining the critical layout of discontinuities and the associated upper bound limit load for plane two-dimensional and three-dimensional (3D) problems. The modelling process (pre-processing) for DLO includes defining the discontinuities inside a specified domain and building the target function and the global constraint matrix for the optimization solver, which has great influence on the the efficiency of the computation processes and the reliability of the final results. This paper focuses on efficient and reliable pre-processing of the discontinuities within the 3D DLO and presents a multi-slicing strategy, which naturally avoids the overlapping and crossing of different discontinuities. Furthermore, the formulation of the 3D discontinuity considering a shape of an arbitrary convex polygon is introduced, permitting the efficient assembly of the global constraint matrix. The proposed method eliminates unnecessary discontinuities in 3D DLO, making it possible to apply 3D DLO for solving large-scale engineering problems such as those involving landslides. Numerical examples including a footing test, a 3D landslide and a punch indentation are considered, illustrating the effectiveness of the presented method. Copyright © 2016 John Wiley & Sons, Ltd.

The paper describes an approach to the simulation of arbitrary laboratory tests with homogeneous stress and strain fields applicable to most constitutive models common in geomechanics. The method by Bardet and Choucair [Int. J. Numer. Anal. Meth. Geomech. 15(1):1–19, 1991] is generalized for an arbitrary number of controlling and controlled variables. The approach is illustrated for time-dependent thermo–hydro–mechanical constitutive models. The purpose of the method is to define an endless variety of different laboratory tests declaratively by means of two matrices ** E** and

It is normally accepted that materials inside the shear band undergo severe rotation of the principal stress direction, which causes non-coaxiality between the principal stress and principal plastic strain rate. However, classical plasticity flow theory implicitly assumes that the principal stress and the principal plastic strain rate are coaxial; thus, it may not correctly predict the onset of the shear band. In addition, classical continuum does not contain any internal length scales; as a result, it cannot provide a reasonable shear band thickness. In this study, the original vertex non-coaxial plastic model based on the classical continuum is extended to the Cosserat continuum. The corresponding codes are implemented via the interface of the user defined element subroutine in ABAQUS. Through a simple shear test, the effectiveness of the user's codes is verified. Through a uniaxial compression test, the influence of non-coaxiality on the onset, the orientation, and the thickness of the shear band is investigated. Results show that the onset of the shear localization is delayed, and the thickness of the shear band is widened when the non-coaxial degree increases, while the orientation of the shear band is little affected by the non-coaxial degree. In addition, it is found that the non-coaxiality can weaken the micro-polar effect to some extent; nonetheless, the Cosserat non-coaxial model still has its advantage over the classical non-coaxial model in capturing the pre-bifurcation as well as the post-bifurcation behaviors of strain localization. Copyright © 2016 John Wiley & Sons, Ltd.

Although the potential contact force proposed by Munjiza overcomes the difficulties inherent in the traditional discrete element methods, the physical meaning of the potential is not clear and the contact force derived from the original potential function is strongly dependent on the mesh configuration. In this study, we redefine a potential function and propose a new contact force calculation method based on a unified standard. Moreover, the new potential function retains all the advantages of the original potential function but has less mesh dependency. Copyright © 2016 John Wiley & Sons, Ltd.

Contact descriptions for interfaces between structural elements (e.g. piles, anchors and tunnel-linings) and soils are widely used in geotechnical engineering. Most of these constitutive interface models, developed for sands, consider 2D conditions. For clays, only a limited number of 3D interface models exist. In this paper, the fine-grained continuum models based on hypoplastic theories are adapted for the constitutive modelling of fine-grained interfaces. To this end, we develop a general approach to convert the existing continuum soil models into an interface model adopting reduced stress and strain-rate vectors and redefining tensorial operations in such a way that the formulation of the existing continuum models can be used without much modification. The hypoplastic Cam-clay and the explicitly formulated hypoplastic models are adapted to model the interface behaviour. In addition to the reduced stress and strain tensor formulations, we introduce an additional variable reducing the strength and stiffness of the interfaces when compared with the strength and stiffness of the soil. For verification, experiments in constant volume and constant normal stress conditions have been simulated. A comparison of the available experimental data from the literature and the simulations is presented. It is shown that the new hypoplastic interface models can describe a number of important phenomena of clay–structure interfaces. Copyright © 2016 John Wiley & Sons, Ltd.

This paper deals with the numerical simulation of the quicksand phenomenon using a coupled Discrete Elements – Lattice Boltzmann hydromechanical model. After the presentation of the developed numerical model, simulations of ascending fluid flow through granular deposits are performed. The simulations show that the quicksand actually triggers for a hydraulic gradient very close to the critical hydraulic gradient calculated from the global analysis of classical soil mechanics, that is, when the resultant of the applied external pressure balances submerged weight of the deposit. Moreover, they point out that the quicksand phenomenon does not occur only for hydraulic gradients above the critical hydraulic gradient, but also in some cases with slightly lower gradients. In such cases, a more permeable zone is first gradually built at the bottom of the deposit through a grain rearrangement, which increases the hydraulic gradient in the upper zones and triggers the phenomenon. Copyright © 2016 John Wiley & Sons, Ltd.

We propose a numerical method that couples a cohesive zone model (CZM) and a finite element-based continuum damage mechanics (CDM) model. The CZM represents a mode II macro-fracture, and CDM finite elements (FE) represent the damage zone of the CZM. The coupled CZM/CDM model can capture the flow of energy that takes place between the bulk material that forms the matrix and the macroscopic fracture surfaces. The CDM model, which does not account for micro-crack interaction, is calibrated against triaxial compression tests performed on Bakken shale, so as to reproduce the stress/strain curve before the failure peak. Based on a comparison with Kachanov's micro-mechanical model, we confirm that the critical micro-crack density value equal to 0.3 reflects the point at which crack interaction cannot be neglected. The CZM is assigned a pure mode II cohesive law that accounts for the dependence of the shear strength and energy release rate on confining pressure. The cohesive shear strength of the CZM is calibrated by calculating the shear stress necessary to reach a CDM damage of 0.3 during a direct shear test. We find that the shear cohesive strength of the CZM depends linearly on the confining pressure. Triaxial compression tests are simulated, in which the shale sample is modeled as an FE CDM continuum that contains a predefined thin cohesive zone representing the idealized shear fracture plane. The shear energy release rate of the CZM is fitted in order to match to the post-peak stress/strain curves obtained during experimental tests performed on Bakken shale. We find that the energy release rate depends linearly on the shear cohesive strength. We then use the calibrated shale rheology to simulate the propagation of a meter-scale mode II fracture. Under low confining pressure, the macroscopic crack (CZM) and its damaged zone (CDM) propagate simultaneously (i.e., during the same loading increments). Under high confining pressure, the fracture propagates in slip-friction, that is, the debonding of the cohesive zone alternates with the propagation of continuum damage. The computational method is applicable to a range of geological injection problems including hydraulic fracturing and fluid storage and should be further enhanced by the addition of mode I and mixed mode (I+II+III) propagation. Copyright © 2016 John Wiley & Sons, Ltd.

The research presented in this paper deals with the numerical analysis of projectile impact on regular strength concrete (RSC), high-strength concrete (HSC), and engineered cementitious composites (ECC) using the Lattice Discrete Particle Model (LDPM). The LDPM is chosen in this study as it naturally captures the failure mechanisms at the length scale of coarse aggregate of concrete, and its capabilities include the accurate depiction of both intrinsic and apparent rate effects in concrete, as well as fiber reinforcement effects. The model is used to predict the experimental impact response performed by four independent testing laboratories, and for each data set the model parameters are calibrated and validated using a combination of uniaxial compression, triaxial compression, uniaxial strain compression, and dogbone tests. In the first study, perforation experiments on RSC and HSC for varied impact velocities are carried out, and the exit velocity is compared with the available experimental data. The second study focuses on ECC, where multiple impact of steel and plastic fiber reinforced concrete panels are explored. A third investigation is performed on four RSC panels with varied thicknesses and subjected to the same impact velocity. In this instance, the model is used to predict the penetration depths for the different cases. Finally, in the last study, the response of large-thickness infinite panels of sizes ranging from 300 mm to 700 mm under projectile impact is considered. Copyright © 2016 John Wiley & Sons, Ltd.

Bifurcation of unsaturated soils into a localized shear band is a ubiquitous failure mode of partially saturated soils. The density and degree of saturation have major impacts on the inception of localized deformations in unsaturated soils. Unsaturated fluid flow may dramatically change the density and degree of fluid saturation of unsaturated soils. Therefore, the unsaturated fluid flow is a potential trigger for shear banding in such materials. In this paper, we derive a simplified bifurcation condition of localized deformation in unsaturated soils under the local transient condition at finite strain. This transient bifurcation condition is implemented into a nonlinear finite element code to study the inception of localized deformation in unsaturated soil specimens. Numerical simulations are conducted to study the impact of soil fabrics of density, a ‘bonding’ variable, and intrinsic permeability on the inception of localized failures via the transient bifurcation criterion. Mesh sensitivity analysis is performed to demonstrate the viscosity effect of unsaturated fluid flow on the localized deformation. Numerical simulations demonstrate that the transient bifurcation condition can detect the localized deformation triggered by the internal unsaturated fluid flow process in unsaturated soils. Copyright © 2016 John Wiley & Sons, Ltd.

In the traditional numerical reservoir simulations, the internodal transmissibility is usually defined as the harmonic mean of the permeabilities of the adjacent grids. This definition underestimates the phase flux and the speed of the saturation front, especially for the strong heterogeneous case. In this article, the internodal transmissibility is recalculated according to the nodal analytic solution. The redefined internodal transmissibility can be used directly to calculate the multiphase flow in the numerical reservoir simulations. Numerical examples show that, compared to the traditional numerical methods, the proposed scheme makes the convergences much faster as the refinement parameter increases, and the accuracy is independent of the heterogeneity. Copyright © 2016 John Wiley & Sons, Ltd.

Slope stability optimization, in the presence of a band of a weak layer between two strong layers, is accounted for in complicated geotechnical problems. Classical optimization algorithms are not suitable for solving such problems as they need a proper preliminary solution to converge to a valid result. Therefore, it is necessary to find a proper algorithm which is capable of finding the best global solution. Recently a lot of metaheuristic algorithms have been proposed which are able to evade local minima effectively. In this study four evolutionary algorithms, including well-known and recent ones, such as genetic algorithm, differential evolution, evolutionary strategy and biogeography-based optimization (BBO), are applied in slope stability analysis and their efficiencies are explored by three benchmark case studies. Result show BBO is the most efficient among these evolutionary algorithms and other proposed algorithms applied to this problem. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, we consider numerical algorithms for modeling of the time-dependent coupling between the fluid flow and deformation in elastic porous media. Here, we employ a four-field formulation which uses the total stress, displacement, flux, and pressure as its primary variables and satisfies Darcy's law and linear elasticity in mixed weak form. We present four different iteratively coupled methods, known as drained, undrained, fixed-strain, and fixed-stress splits, in which the diffusion operator is separated from the elasticity operator and the two subproblems are solved in a staggered way while ensuring convergence of the solution at each time step. A-priori convergence results for each iterative coupling which differs from those found when using a traditional two-field or three-field formulation are presented. We also present some numerical results to support the convergence estimates and to show the accuracy and efficiency of the algorithms. Copyright © 2016 John Wiley & Sons, Ltd.

The dynamic behaviors of railway ballast under cyclic loading are simulated with discrete element method (DEM). Dilated polyhedra are constructed based on the Minkowski sum operator in order to resemble the irregular shapes of ballast particles. The polyhedral particle generation, contact detection between particles and contact laws are presented. Ballast box tests with periodic lateral boundaries are conducted to simulate the dynamics of the sleeper and ballast particles. The settlement and effective stiffness of ballast bed are investigated under cyclic loadings with five distinct frequencies. The settlement of ballast bed is significant in the first several cycles and increases with the number of cycles gradually. The higher frequency loading generates larger displacement in the same simulation time. The effective stiffness of ballast bed increases gradually. To study the effect of particles' sharpness, dilated polyhedra with different dilating radii and spherical particles are also developed. Simulation results show the sharper the ballast particles are, the smaller the produced settlement. Copyright © 2016 John Wiley & Sons, Ltd.

As a result of deposition process and particle characteristics, granular materials can be inherently anisotropic. Many researchers have strongly suggested that the inherent anisotropy is the main reason for the deformation non-coaxiality of granular materials. However, their relationships are not unanimous because of the limited understanding of the non-coaxial micro-mechanism. In this study, we investigated the influence of inherent anisotropy on the non-coaxial angle using the discrete element method. Firstly, we developed a new discrete element method approach using rough elliptic particles and proposed a novel method to produce anisotropic specimens. Secondly, the effects of initial specimen density and particle characteristics, such as particle aspect ratio *A*_{m}, rolling resistance coefficient *β*, and bedding plane orientation *δ*, were examined by a series of biaxial tests and rotational principal axes tests. Findings from the numerical simulations are summarized as follows: (1) the peak internal friction angle *ϕ*_{p} and the non-coaxial angle *i* both increase with the initial density, *A*_{m} and *β*, and they both increase initially and then decrease with *δ* in the range of 0–90°; (2) among the particle characteristics, the influence of *A*_{m} is the most significant; and (3) for anisotropic specimens, the non-coaxial angle can be calculated using the double slip and rotation rate model. Then, an empirical formula was proposed based on the simulation results to depict the relationship between the non-coaxial angle and the particle characteristics. Finally, the particle-scale mechanism of non-coaxiality for granular materials was discussed from the perspective of energy dissipation. Copyright © 2016 John Wiley & Sons, Ltd.

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.

This article is an attempt at providing an insight into the development of hypoplasticity (including barodesy, which is a recent development of hypoplasticity) as a theory elaborated since 1977, when the first version was published by the first author, until present. The multiplicity of the many versions published since then is hard to overlook. This article presents a review and insight into the evolution of a theory and the struggle to formulate a satisfactory constitutive law. Among the many proposed versions, we focus on those ones that can be seen as changes of paradigm. Copyright © 2016 John Wiley & Sons, Ltd.

A Lagrangian numerical approach for the simulation of rapid landslide runouts is presented and discussed. The simulation approach is based on the so-called Particle Finite Element Method. The moving soil mass is assumed to obey a rigid-viscoplastic, non-dilatant Drucker–Prager constitutive law, which is cast in the form of a regularized, pressure-sensitive Bingham model. Unlike in classical formulations of computational fluid mechanics, where no-slip boundary conditions are assumed, basal slip boundary conditions are introduced to account for the specific nature of the landslide-basal surface interface. The basal slip conditions are formulated in the form of modified Navier boundary conditions, with a pressure-sensitive threshold. A special mixed Eulerian–Lagrangian formulation is used for the elements on the basal interface to accommodate the new slip conditions into the Particle Finite Element Method framework. To avoid inconsistencies in the presence of complex shapes of the basal surface, the no-flux condition through the basal surface is relaxed using a penalty approach. The proposed model is validated by simulating both laboratory tests and a real large-scale problem, and the critical role of the basal slip is elucidated. Copyright © 2016 John Wiley & Sons, Ltd.

Three-dimensional particle morphology is a significant problem in the discrete element modeling of granular sand. The major technical challenge is generating a realistic 3D sand assembly that is composed of a large number of random-shaped particles containing essential morphological features of natural sands. Based on X-ray micro-computed tomography data collected from a series of image processing techniques, we used the spherical harmonics (SH) analysis to represent and reconstruct the multi-scale features of real 3D particle morphologies. The SH analysis was extended to some highly complex particles with sharp corners and surface cavities. We then proposed a statistical approach for the generation of realistic particle assembly of a given type of sand based on the principle component analysis (PCA). The PCA aims to identify the major pattern of the coefficient matrix, which is made up of the SH coefficients of all the particles involved in the analysis. This approach takes into account the particle size effect on the variation of particle morphology, which is observed from the available results of micro-computed tomography and QICPIC analyses of sand particle morphology. Using the aforementioned approach, two virtual sand samples were generated, whose statistics of morphological parameters were compared with those measured from real sand particles. The comparison shows that the proposed approach is capable of generating a realistic sand assembly that retains the major morphological features of the mother sand. Copyright © 2016 John Wiley & Sons, Ltd.

Thermo-hydro-mechanical responses around a cylindrical cavity drilled or excavated in a low-permeability formation are studied when the cavity is subjected to a time-dependent thermal loading. The cavity is considered backfilled after it is supported by casing or lining. Solutions of temperature, pore water pressure, stress, and displacement responses are analytically formulated based on Biot's consolidation theory with the assumption that the backfilling material, supporting material, and surrounding low-permeability formation are poroelastic media. The solution is expressed in Laplace space, and numerical inversion techniques are used to find field variables in the real-time domain. After the solution is verified with the numerical results, it is applied in a large-scale *in situ* heating test – PRACLAY heating test – for a predictive reference calculation and an extensive parametric study. Another medium-scale *in situ* heating test – ATLAS III heating test – is also analyzed using the solution, which provides reasonable agreement with measurements. The new analytical solution proves to be a convenient tool for a good understanding of the resulting coupled thermo-hydro-mechanical behavior and is therefore valuable for the interpretation of measured data in engineering practices and for a rational design of potential radioactive waste repositories. Copyright © 2016 John Wiley & Sons, Ltd.

The complexity of formulations for the hydromechanical coupled mechanics of porous media is typically minimised by simplifying assumptions such as neglecting the effect of inertia terms. For example, three formulations commonly employed to model practical problems are classified as fully dynamic, simplified dynamic and quasi-static. Thus, depending on the porous media conditions, each formulation will have advantages and limitations. This paper presents a comprehensive analysis of these limitations when solving one-dimensional fully saturated porous media problems in addition to a new solution that considers a more general loading situation. A phase diagram is developed to assist on the selection of which formulation is more appropriate and convenient regarding particular cases of porosity and hydraulic conductivity values. Non-dimensional formulations are proposed to achieve this goal. Results using the analytical solutions are compared against numerical values obtained with the finite element method, and the effect of porosity is investigated. Copyright © 2016 John Wiley & Sons, Ltd.

A simple method called anisotropic transformed stress (ATS) method is proposed to develop failure criteria and constitutive models for anisotropic soils. In this method, stress components in different directions are modified differently in order to reflect the effect of anisotropy. It includes two steps of mapping of stress. First, a modified stress tensor is introduced, which is a symmetric multiplication of stress tensor and fabric tensor. In the modified stress space, anisotropic soils can be treated to be isotropic. Second, a TS tensor is derived from the modified stress tensor for the convenience of developing anisotropic constitutive models to account for the effect of intermediate principal stress. By replacing the ordinary stress tensor with the TS tensor directly, the unified hardening model is extended to model the anisotropic deformation of soils. Anisotropic Lade's criterion is adopted for shear yield and shear failure in the model. The form of the original model formulations remains unchanged, and the model parameters are independent of the loading direction. Good agreement between the experimental results and predictions of the anisotropic unified hardening model is observed. Copyright © 2016 John Wiley & Sons, Ltd.

This paper revisits the variational limit equilibrium (LE) analysis of three-dimensional (3D) slope stability in the context of limit analysis (LA). It proves the kinematic admissibility of the 3D mechanism in LA, although it was derived from LE variational extremization. It also includes algorithms in the realm of LA that are associated with the variational mechanism. A comparison between the variational results and reported LA upper-bound or LE closed-form results is conducted. It demonstrates that the variationally derived mechanism consistently yields upper-bound solutions for 3D symmetrical slopes that are as accurate as those produced by postulated mechanisms in LA. However, the results are more critical than those derived from spherical failure mechanism in LE. The generalized log spiral 3D mechanism rigorously legitimizes the variational slope stability analysis in both frameworks of mechanics LE and LA. Stability charts were produced where the 3D factor of safety can be assessed for a constrained length of failure, while including factors like pore water pressure and seismic loading. The results presented within this study demonstrate the capabilities of the variational 3D solution and can be used to evaluate approximate methods, numerical or closed-form, developed in 3D slope stability analyses. Copyright © 2016 John Wiley & Sons, Ltd.

This paper presents a numerical evaluation of three non-coaxial kinematic models by performing Distinct Element Method (DEM) simple shear tests on specimens composed of elliptical particles with different aspect ratios of 1.4 and 1.7. The models evaluated are the double-shearing model, the double-sliding free-rotating model and the double slip and rotation rate model (DSR^{2} model). Two modes of monotonic and cyclic simple shear tests were simulated to evaluate the role played by the inherent anisotropy of the specimens. The main findings are supported by all the DEM simple shear tests, irrespective of particle shape, specimen density or shear mode. The evaluation demonstrates that the assumption in the double-shearing model is inconsistent with the DEM results and that the energy dissipation requirements in the double-sliding free-rotating model appear to be too restrictive to describe the kinematic flow of elliptical particle systems. In contrast, the predictions made by the DSR^{2} model agree reasonably well with the DEM data, which demonstrates that the DSR^{2} model can effectively predict the non-coaxial kinematic behavior of elliptical particle systems. Copyright © 2016 John Wiley & Sons, Ltd.

Thermal fracturing can play an important role in development of unconventional petroleum and geothermal resources. Thermal fractures can result from the nonlinear deformation of the rock in response to thermal stress related to cold water injection as well as heating. Before the rock reaches the final failure stage, material softening and bulk modulus degradation can cause changes in the thermo-mechanical properties of the solid. In order to capture this aspect of the rock fracture, a virtual multidimensional internal bond-based thermo-mechanical model is derived to track elastic, softening, and the failure stages of the rock in response to the temporal changes of its temperature field. The variations in thermo-mechanical properties of the rock are derived from a nonlinear constitutive model. To represent the thermo-mechanical behavior of pre-existing fractures, the element partition method is employed. Using the model, numerical simulation of 3D thermal fracture propagation in brittle rock is carried out. Results of numerical simulations provide evidence of model verification and illustrate nonlinear thermal response and fracture development in rock under uniform cooling. In addition, fracture coalescence in a cluster of fractures under thermal stress is illustrated, and the process of thermal fracturing from a wellbore is captured. Results underscore the importance of thermal stress in reservoir stimulation and show the effectiveness of the model to predict 3D thermal fracturing. Copyright © 2016 John Wiley & Sons, Ltd.

At present, several of the existing elastoplastic constitutive models are adapted for describing the stress–strain behavior of unsaturated soils. However, most of them present certain limitations in this field. These limitations can be related to the basic model and/or added unsaturated state variables and formulations. In this regard, inability to model the hydro-mechanical behavior in constant water (CW) conditions is an example of these limitations. In this paper, an advanced version of CJS model is selected for adaptation to the unsaturated states. Adaptation to unsaturated states is achieved in the framework of effective stress approach. Effective stress equation and unsaturated state variables are selected based on the recent research existing in the literature. The developed model is capable of describing the complex behavior of unsaturated soil in the CW condition in addition to predicting the behavior at failure and post–failure, nonlinear elastoplastic behavior at low levels of stress and strain (by selecting a very small elastic domain), as well as wetting and collapse behaviors. In order to validate the model, results of triaxial tests in CD and CW conditions are used. The validation results indicate the good capability of the proposed model. Behavior of the unsaturated soils during wetting is an important issue. For this reason, the model is also evaluated based on the results of wetting and collapse triaxial tests. A comparison between the tests and simulation results shows that the model is able to predict the soil behavior under the wetting path. Copyright © 2016 John Wiley & Sons, Ltd.

Geomaterials such as soils and rocks are inherently anisotropic and sensitive to temperature changes caused by various internal and external processes. They are also susceptible to strain localization in the form of shear bands when subjected to critical loads. We present a thermoplastic framework for modeling coupled thermomechanical response and for predicting the inception of a shear band in a transversely isotropic material using the general framework of critical state plasticity and the specific framework of an anisotropic modified Cam–Clay model. The formulation incorporates anisotropy in both elastic and plastic responses under the assumption of infinitesimal deformation. The model is first calibrated using experimental data from triaxial tests to demonstrate its capability in capturing anisotropy in the mechanical response. Subsequently, stress-point simulations of strain localization are carried out under two different conditions, namely, isothermal localization and adiabatic localization. The adiabatic formulation investigates the effect of temperature on localization via thermomechanical coupling. Numerical simulations are presented to demonstrate the important role of anisotropy, hardening, and thermal softening on strain localization inception and orientation. Copyright © 2016 John Wiley & Sons, Ltd.

An adaptive substepping explicit integration scheme is developed for a porosity-dependent hydro-mechanical model for unsaturated soils. The model is referred to as the modified **σ**–Θ model in this paper, which features the employment of the subloading surface plasticity and the stress–saturation approach. On numerical aspects, convex/nonconvex subloading surfaces in the **σ**–Θ space may result in incorrect loading–unloading decisions during the integration. A new loading–unloading decision method is developed here to solve the problem and then embedded into the explicit integration scheme for the modified **σ**–Θ model. In addition, to enhance the accuracy of the explicit integration, local errors from both hydraulic and mechanical components are included in the error control for each substep. A drift correction method is also developed to ensure the state point lies on the subloading surface in the **σ**–Θ space within a set error level. The performance of the loading–unloading decision method for the modified **σ**–Θ model is discussed through comparing it with the conventional loading–unloading decision method. The importance of involving the hydraulic component in the error control is also demonstrated. The accuracy and efficiency of the proposed adaptive substepping explicit integration scheme for the modified *p*–Θ model are also studied via several numerical examples.

Under the proportional strain loading path, particle assemblies may exhibit various failure modes. Besides the strain localization, the diffuse failure may also occur under certain conditions. The diffuse failure mode corresponds to a homogeneous occurrence of failure with stress states strictly included within the plastic limit condition. This paper emphasizes the influences of the density degree and the rolling resistance under the strain path. A contact model considering rolling friction is adopted in a discrete element method analysis as an approximate means to account for the effects of particle shape. Mechanical responses indicate that loose assemblies without the rolling resistance are more vulnerable to static liquefaction. A sample with a smaller initial void ratio or larger rolling friction coefficient will reinforce the stability of the structure and reduce the likelihood of failure. For microscopic properties, the evolution of coordination numbers, contact forces, force chains and the anisotropies of the assemblies are explored and discussed. Rotational resistance helps increase the shear stress of the granular material, and the microscopic parameters indicate that the assembly has a strong anisotropy and a stable structure to resist the increasing loading.

Natural soils are one of the most inherently variables in the ground. Although the significance of inherent soil variability in relation to reliable predictions of consolidation rates of soil deposits has long been realized, there have been few studies that addressed the issue of soil variability for the problem of ground improvement by prefabricated vertical drains. Despite showing valuable insights into the impact of soil spatial variability on soil consolidation by prefabricated vertical drains, available stochastic works on this subject are based on a single-drain (or unit cell) analyses. However, how the idealized unit cell solution can be a supplement to the complex multi-drain systems for spatially variable soils has never been addressed in the literature. In this study, a rigorous stochastic finite elements modeling approach that allows the true nature of soil spatial variability to be considered in a reliable and quantifiable manner, both for the single-drain and multi-drain systems, is presented. The feasibility of performing an analysis based on the unit cell concept as compared with the multi-drain analysis is assessed in a probabilistic context. It is shown that with proper input statistics representative of a particular domain of interest, both the single-drain and multi-drain analyses yield almost identical results.

A novel, simplified approach is presented in order to compute variations of grading in granular assemblies during confined comminution under quasi-static compression. The method is based on a population balance equation and requires a breakage probability, considered here as a probabilistic phenomenon that takes into account the particle strength and the loading condition of individual grains. Under basic assumptions, a simple breakage probability can be defined in order to get a valuable result for engineering applications and powder technology. The size effect in the strength of individual particles is introduced according to Weibull's theory. The particle loading and the cushioning effect in the granular packing are accounted for by considering the orientations of the contact forces obtained from 3D discrete element method simulations of highly polydisperse materials. The method proposed could have a value for engineering purposes in powder technology and geomechanics and gives a general framework for further research developments based on population balance.

Quantitative assessment of the risk of submarine landslides is an essential part of the design process for offshore oil and gas developments in deep water, beyond the continental shelf. Landslides may be triggered by a reduction in shear strength of subsea sediments over a given zone, caused for example by seismic activity. Simple criteria are then needed to identify critical conditions whereby the zone of weakness could grow catastrophically to cause a landslide. A number of such criteria have been developed over the last decade, based either on ideas drawn from fracture mechanics, or considering the equilibrium of the initial weakened zone and adjacent process zones of gradually softening material. Accounting for the history of the weak zone initiation is critical for derivation of reliable propagation criteria, in particular considering dynamic effects arising from accumulating kinetic energy of the failing material, which will allow the failure to propagate from a smaller initial zone of weakened sediments. Criteria are developed here for planar conditions, taking full account of such dynamic effects, which are shown to be capable of reducing the critical length of the softened zone by 20% or more compared with criteria based on static conditions. A numerical approach is used to solve the governing dynamic equations for the sliding material, the results from which justify assumptions that allow analytical criteria to be developed for the case where the initial softening occurs instantaneously. The effect of more gradual softening is also explored.

In the absence of initial cracks, the material behavior is limited by its strength, usually defined in homogeneous conditions (of stress and strain). Beyond this limit, in quasi-brittle case, cracks may propagate and the material behavior tends to be well described by fracture mechanics. Discrete element approaches show consistent results dealing with this transition during rupture. However, the calibration of the parameters of the numerical models (i.e., stiffness, strength, and toughness) may be quite complex and sometimes only approximative. Based on a brittle rupture criterion, we analyze the biaxial response of uncracked samples. Thus, tensile and compressive strengths are analytically identified and become direct parameters of our discrete model. Furthermore, a physically reliable crack initiation (and subsequent propagation) is shown to be induced during rupture and verified by the simulation of three-point bending and diametral compression tests. Copyright © 2016 John Wiley & Sons, Ltd.

In the present study, we have developed a numerical method which can simulate the dynamic behaviour of a seabed ground during gas production from methane hydrate-bearing sediments. The proposed method can describe the chemo-thermo-mechanical-seismic coupled behaviours, such as phase changes from hydrates to water and gas, temperature changes and ground deformation related to the flow of pore fluids during earthquakes. In the first part of the present study, the governing equations for the proposed method and its discretization are presented. Then, numerical analyses are performed for hydrate-bearing sediments in order to investigate the dynamic behaviour during gas production. The geological conditions and the material parameters are determined using the data of the seabed ground at Daini-Atsumi knoll, Eastern Nankai Trough, Japan, where the first offshore production test of methane hydrates was conducted. A predicted earthquake at the site is used in the analyses.

Regarding the seismic response to the earthquake which occur during gas production process, the wave profiles of horizontal acceleration and horizontal velocity were not extensively affected by the gas production. Hydrate dissociation behaviour is sensitive to changes in the pore pressure during earthquakes. Methane hydrate dissociation temporarily became active in some areas because of the main motion of the earthquake, then methane hydrate dissociation brought about an increase in the average pressure of the fluids during the earthquake. And, it was this increase in average pore pressure that finally caused the methane hydrate dissociation to cease during the earthquake.

Compaction and associated fluid flow are fundamental processes in sedimentary basin deformation. Purely mechanical compaction originates mainly from pore fluid expulsion and rearrangement of solid particles during burial, while chemo-mechanical compaction results from Intergranular Pressure-Solution (IPS) and represents a major mechanism of deformation in sedimentary basins during diagenesis. The aim of the present contribution is to provide a comprehensive 3D framework for constitutive and numerical modeling of purely mechanical and chemo-mechanical compaction in sedimentary basins. Extending the concepts that have been previously proposed for the modeling of purely mechanical compaction in finite poroplasticity, deformation by IPS is addressed herein by means of additional viscoplastic terms in the state equations of the porous material. The finite element model integrates the poroplastic and poroviscoplastic components of deformation at large strains. The corresponding implementation allows for numerical simulation of sediments accretion/erosion periods by progressive activation/deactivation of the gravity forces within a fictitious closed material system. Validation of the numerical approach is assessed by means of comparison with closed-form solutions derived in the context of a simplified compaction model. The last part of the paper presents the results of numerical basin simulation performed in one dimensional setting, demonstrating the ability of the modeling to capture the main features in elastoplastic and viscoplastic compaction.

This paper presents a generalized, rigorous and simple large strain solution for the undrained expansion of a vertical cylindrical cavity in critical state soils using a rate-based plasticity formulation: the initial stress field is taken as anisotropic, that is with horizontal stresses that differ from the vertical stress, and the soil is assumed to satisfy any two-invariant constitutive model from the critical state (Cam-clay) family; no simplifying assumption is made during the mathematical derivation; calculating the effective stresses around the cavity requires the solution of a nonlinear equation by means of the Newton–Raphson method in combination with quadrature. Cavity expansion curves and stress distributions in the soil are then presented for different critical state models (including the modified Cam-clay model). The solution derived can be useful for estimating the instantaneous response of saturated low-permeability soils around piles and self-boring pressuremeters and can serve as trustworthy benchmark for numerical analysis codes. Copyright © 2016 John Wiley & Sons, Ltd.

The kinematic approach in combination with numerical simulation is used to examine the effect of pore water pressure on tunnel face stability. Pore water pressure distribution obtained by numerical calculations using *FLAC*^{3D} is used to interpolate the pore water pressure on a 3D rotational collapse mechanism. Comparisons are made to check the present approach against other solutions, showing that the present approach improves the existing upper bound solutions. Results obtained indicate that critical effective face pressure increases with water table elevation. Several normalized charts are also presented for quick evaluation of tunnel face stability. At the end of the paper, the influence of anisotropic permeability on tunnel face stability is also discussed, showing that the isotropic model leads to an overestimation of the necessary tunnel face pressure for anisotropic soils. Copyright © 2016 John Wiley & Sons, Ltd.

Crack growth in hot brittle rocks, driven by thermal cooling, was simulated using a coupled two-dimensional discrete element and heat transport model that explicitly includes the random initiation and subsequent propagation of interacting cracks. The model clearly predicts that a quasi-hierarchical array of subparallel cracks, oriented along the direction of the temperature gradient, is formed under small to moderately large thermally generated strain load conditions. The simulation results also demonstrate that, after an initial transient, thermal cracks propagate in a stable fashion with a velocity that scales with ~ *t*^{− 1/2}. However, under large thermal strain loads, a more complicated geometry composed of cracks that curve and coalesce develops during the later stages of crack growth. Copyright © 2016 John Wiley & Sons, Ltd.

In this paper, a fully coupled model is developed for numerical modeling of hydraulic fracturing in partially saturated weak porous formations using the extended finite element method, which provides an effective means to simulate the coupled hydro-mechanical processes occurring during hydraulic fracturing. The developed model is for short fractures where plane strain assumptions are valid. The propagation of the hydraulic fracture is governed by the cohesive crack model, which accounts for crack closure and reopening. The developed model allows for fluid flow within the open part of the crack and crack face contact resulting from fracture closure. To prevent the unphysical crack face interpenetration during the closing mode, the crack face contact or self-contact condition is enforced using the penalty method. Along the open part of the crack, the leakage flux through the crack faces is obtained directly as a part of the solution without introducing any simplifying assumption. If the crack undergoes the closing mode, zero leakage flux condition is imposed along the contact zone. An application of the developed model is shown in numerical modeling of pump-in/shut-in test. It is illustrated that the developed model is able to capture the salient features bottomhole pressure/time records exhibit and can extract the confining stress perpendicular to the direction of the hydraulic fracture propagation from the fracture closure pressure. Copyright © 2016 John Wiley & Sons, Ltd.

This paper develops a three-layer model and elastic solutions to capture nonlinear response of rigid, passive piles in sliding soil. Elastic solutions are obtained for an equivalent force per unit length *p _{s}* of the soil movement. They are repeated for a series of linearly increasing

- On-pile pressure in rotationally restrained, sliding layer reduces by a factor α, which resembles the
*p*-multiplier for a laterally loaded, capped pile, but for its increase with vertical loading (embankment surcharge), and stiffness of underlying stiff layer: α = 0.25 and 0.6 for a shallow, translating and rotating piles, respectively; α = 0.33–0.5 and 0.8–1.3 for a slide overlying a stiff layer concerning a uniform and a linearly increasing pressure, respectively; and α = 0.5–0.72 for moving clay under embankment loading. - Ultimate state is well defined using the ratio of passive earth pressure coefficient over that of active earth pressure. The subgrade modulus for a large soil movement may be scaled from model tests.
- The normalised rotational stiffness is equal to 0.1–0.15 for the capped piles, which increases the pile displacement with depth.

The three-layer model solutions well predict nonlinear response of capped piles subjected to passive loading, which may be used for pertinent design. Copyright © 2016 John Wiley & Sons, Ltd.

Hydraulic fracturing (HF) of underground formations has widely been used in different fields of engineering. Despite the technological advances in techniques of *in situ* HF, the industry uses semi-analytical tools to design HF treatment. This is due to the complex interaction among various mechanisms involved in this process, so that for thorough simulations of HF operations a fully coupled numerical model is required.

In this study, using element-free Galerkin (EFG) mesh-less method, a new formulation for numerical modeling of hydraulic fracture propagation in porous media is developed. This numerical approach, which is based on the simultaneous solution of equilibrium and continuity equations, considers the hydro-mechanical coupling between the crack and its surrounding porous medium. Therefore, the developed EFG model is capable of simulating fluid leak-off and fluid lag phenomena.

To create the discrete equation system, the Galerkin technique is applied, and the essential boundary conditions are imposed via penalty method. Then, the resultant constrained integral equations are discretized in space using EFG shape functions. For temporal discretization, a fully implicit scheme is employed. The final set of algebraic equations that forms a non-linear equation system is solved using the direct iterative procedure.

Modeling of cracks is performed on the basis of linear elastic fracture mechanics, and for this purpose, the so-called diffraction method is employed. For verification of the model, a number of problems are solved. According to the obtained results, the developed EFG computer program can successfully be applied for simulating the complex process of hydraulic fracture propagation in porous media. Copyright © 2016 John Wiley & Sons, Ltd.

This paper presents a hypoplastic constitutive model for the viscous behavior of frozen soil. The model is composed of a ‘solid’ part and a ‘fluid’ part. The solid part is based on the extended hypoplastic model, and the fluid part is dependent on the second time derivative of strain. The performance of the model is demonstrated by simulating some uniaxial compression tests at different strain rates. Moreover, the model can describe in a unified way the three stages of typical creep tests, that is, primary, secondary, and tertiary stage. Copyright © 2016 John Wiley & Sons, Ltd.

Concrete-faced rockfill dam (CFRD) is a popular alternative to traditional dam types in the last two decades. The modelling of CFRD involves complex multi-body contact and strong geometry and material nonlinearities. We present a numerical approach for the modelling of CFRDs in this paper. Based on the dual-mortar finite element method, the presented approach considers different parts of rockfill and all concrete slabs as independent deformable continuum. The multi-body contacts are modelled using Lagrange multipliers with a weak form segment-to-segment contact strategy. To alleviate instability induced by strong geometry nonlinearity in the slab–slab contact, we propose a mixed type of constraints for the tangential contact. A general transformation scheme is introduced to simplify the implementation of contact constraints. Three-dimensional analysis of Tianshengqiao-1 CFRD is performed. The nonlinear and time-dependent deformation of the rockfill is considered. We study the influence of the rockfill deformation on the reliability of the concrete face. Three major concerns of the face, that is, the axial compression, the slab–slab separation and the face-rockfill separation, are discussed in detail. The numerical results are compared with data from *in-situ* observation. Copyright © 2016 John Wiley & Sons, Ltd.

Modeling of progressive development of zones of large inelastic shear deformation (shear band) that results from strain-softening behavior of sensitive clays could explain the failure mechanisms of large landslides. Because of toe erosion, a shear band can be initiated and propagated upward (inward) from the river bank. On the other hand, upslope surcharge loading could generate shear bands that might propagate down towards the river bank. In the present study, upward and downward propagation of shear bands and failure of sensitive clay slopes are modeled using the Coupled Eulerian Lagrangian approach in Abaqus finite element (FE) software. It is shown that the formation and propagation of shear bands are significantly influenced by kinematic constraints that change with displacements of the soil masses, and therefore the propagation of an existing shear band might be stopped and new shear bands could be formed. The main advantages of the present FE modeling are: (i) extremely large strains in the shear bands can be successfully simulated without numerical issues, (ii) a priori definition of shearing zones is not required to tackle severe strains; instead, the FE program automatically identifies the critical locations for shear band formation and propagation. Toe erosion could significantly increase the slope failure potential because of upslope surcharge loading. FE analyses with a thick and thin sensitive clay layers show that the global failure could occur at lower surcharge loads in the former as compared to the latter cases. Copyright © 2016 John Wiley & Sons, Ltd.

This paper uses Biot's poroelasticity approach to examine the consolidation behaviour of a rigid foundation with a frictionless base in contact with a poroelastic halfspace. The mathematical development of the mixed boundary value problem involves a set of dual integral equations in the Laplace transform domain which cannot be conveniently solved by employing conventional procedures. In this paper, a numerical solution is developed using a scheme where the contact normal stress is approximated by a discretized equivalent. The influence of limiting drainage boundary conditions at the entire surface of the halfspace on the degree of consolidation of the rigid circular foundation is investigated. The results obtained in this study are compared with the corresponding results given in the literature. Copyright © 2016 John Wiley & Sons, Ltd.

In the practice of geotechnical engineering, the case of a ring footing carrying a set of concentrated point loads is a common problem. At times, the induced vertical and angular displacements for the ring footing need to be evaluated at a relatively precise level. By making use of the governing set of equations derived for the case of a general curved beam, expressions that can be easily implemented in modern computing software are derived for the vertical and angular displacements of a ring footing of rectangular cross section as functions of the radial position. The loading case considered is a vertical point load, and the soil is modelled as elastic. Estimates of the displacements have been shown for a common range of practical applications. The behaviour for a set of concentrated loads may be evaluated using the derived equations through direct superposition. Nonlinear finite element analysis is used to evaluate the vertical deflection and angular twist of the ring foundation. Numerical analysis performed for three ring foundations with different radii and cross sections is reported to validate the accuracy of the derived analytical solution. Copyright © 2016 John Wiley & Sons, Ltd.

A new analytical proof is presented for steady-state seepage in recharged heterogeneous unconfined aquifers. The paper also presents a detailed procedure and important rules for performing correctly numerical studies of unsaturated seepage. Once a numerical solution is calibrated with field data, using a set of spatially distributed values for hydraulic conductivity *K* and effective infiltration *EI*, any new numerical analysis with a set of *αK* and *αEI* values, where *α* is a constant, yields an equally good calibration. However, if the effective porosities of each layer are unchanged, the groundwater velocities are multiplied by *α*, whereas the travel times are divided by *α*, which may help to select *α* in order to match known travel time data. This is a clear example of multiple solutions to an inverse problem. The paper underlines the role and the need to finely mesh unsaturated zones and also contacts between layers to reach the asymptotic convergence range, as it was carried out to verify the proof and as it should be completed to study any seepage problem. A few consequences of the new analytical proof and the rigorous procedure are shown with examples. Copyright © 2016 John Wiley & Sons, Ltd.

By virtue of a pair of scalar potentials for the displacement of the solid skeleton and the pore fluid pressure field of a saturated poroelastic medium, an alternative solution method to the Helmholtz decomposition is developed for the wave propagation problems in the framework of Biot's theory. As an application, a comprehensive solution for three-dimensional response of an isotropic poroelastic half-space with a partially permeable hydraulic free surface under an arbitrarily distributed time-harmonic internal force field and fluid sources is developed. The Green's functions for the poroelastic fields, corresponding to point, ring, and disk loads, are reduced to semi-infinite complex-valued integrals that can be evaluated numerically by an appropriate quadrature scheme. Analytical and numerical comparisons are made with existing elastic and poroelastic solutions to illustrate the quality and features of the solution. Copyright © 2016 John Wiley & Sons, Ltd.

The main purpose of the paper is to present a relatively simple, yet realistic, constitutive model for simulations of structured sensitive clays. The proposed constitutive model can simulate 1-D and isotropic consolidation, and drained and undrained shear response of sensitive structured clay.

The proposed sensitive bounding surface model is based on concepts from the modified Cam clay model (Roscoe and Burland, 1968) and bounding surface plasticity (Dafalias and Herrmann, 1982), with the addition of a simple degradation law. The key material parameters are *M*, *λ*, *κ*, and *ν* from the modified Cam clay framework, *h* from the bounding surface framework to model a smoothed elasto-plastic transition, and *ω _{v}*,

The model has separate parameters to model destructuration caused by volumetric strain and deviatoric strain. The model is capable of modeling unusual behavior of strain softening during 1-D compression (i.e., a reduction of effective stress as void ratio decreases). A good match between test results and the model simulation is demonstrated. Copyright © 2016 John Wiley & Sons, Ltd.

No abstract is available for this article.

]]>The mechanical behavior of granular materials is characterized by strong nonlinearity and irreversibility. These properties have been differently described by a variety of constitutive models. To test any constitutive model, experimental data relative to the nature of the incremental stress–strain response of the material is desirable. However, this type of laboratory data is scarce because of being expensive and difficult to obtain. The discrete element method has been used several times as an alternative to obtain incremental responses of granular materials. Crushable grains add one extra source of irreversibility to granular materials. Crushability has been variously incorporated into different constitutive models. Again, it will be helpful to obtain incremental responses of crushable granular materials to test these models, but the experimental difficulties are increased. Making use of a recently introduced crushing model for discrete element simulation, this paper presents a new procedure to obtain incremental responses in discrete analogs of granular crushable materials. The parallel probe approach, previously used for uncrushable discrete analogs, is here extended to account for the presence of crushable grains. The contribution of grain crushing to the incremental irreversible strain is identified and separately measured. Robustness of the proposed method is examined in detail, paying particular attention to aspects such as dynamic instability or crushing localization. The proposed procedure is later applied to map incremental responses of a discrete analog of Fontainebleau sand on the triaxial plane. The effect of stress ratio and granular state on plastic flow characteristics is highlighted. Copyright © 2016 John Wiley & Sons, Ltd.

In this article, we evaluate geomechanics of fluid injection from a fully penetrating vertical well into an unconsolidated formation confined with stiff seal rocks. The coupled behavior of an isotropic, homogeneous sand layer is studied under injection pressures that are high enough to induce plasticity yet not fracturing. Propagation of the significant influence zone surrounding the injection borehole, quantified by the extent of the plastic domain in the elasto-plastic model, is examined for the first time. First, a new fully coupled axisymmetric numerical model is developed. A comprehensive assessment is performed on pore pressures, stresses/strains, and failure planes during the entire transient period of an injection cycle. Results anticipate existence of five distinctive zones in terms of plasticity state: liquefaction at wellbore; two inner plastic domains surrounding the wellbore, where failure occurs along two planes and major principal stress is in vertical direction; remaining of the plastic domain, where formation fails along one plane and major principal stress is in radial direction; and a non-plastic region. Failure mechanism at the wellbore is found to be shear followed by liquefaction. Next, a novel methodology is proposed based on which new weakly coupled poro-elasto-plastic analytical solutions are derived for all three stress/strain components. Unlike previous studies, extension of the plastic zone is obtained as a function of injection pressure, incorporating plasticity effects on the subsequent elastic domain. Solutions, proven to be a good approximation of numerical simulations, offer a huge advantage as the run time of coupled numerical simulations is considerably long. Copyright © 2016 John Wiley & Sons, Ltd.

One-dimensional mathematical models for vapor-phase volatile organic compound (VOC) diffusion through composite cover barriers are presented. An analytical solution to the model was obtained by the method of separation of variables. The results obtained by the proposed solution agree well with those obtained by a numerical analysis. Based on the proposed analytical model, the VOC breakthrough curves of five different composite covers are compared. The effects of degree of saturation of geosynthetic clay liner (GCL) or compacted clay liner (CCL) on VOC migration in the composite covers are then presented. Results show that the composite cover barriers provide much better diffusion barriers for VOC than the single CCL. The top surface steady-state flux for a composite barrier, consisting of a 1.5 mm geomembrane (GM) and a 20 cm CCL, can be 8.3 times lower than that for a 30 cm CCL. The surface steady-state flux for the case with (1.5 mm GM + 6 mm GCL) was found to be 2.3 times lower than that for the case with (1.5 mm GM + 20 cm CCL). The degree of saturation *S _{r}* of the CCL has a great influence on VOC migration in composite covers when

This paper presents an analytical-numerical approach to obtain the distribution of stresses and deformations around a reinforced tunnel. The increase in the radial stress of the reinforced tunnel, based on the performance of a bolt, is modeled by a function, which its maximum value is in the vicinity of the bolt periphery and it exponentially decreases in the far distance from the bolt. On the basis of this approach, the shear stiffness between the bolt and the rock mass and the shear stress distribution around the bolt within the rock mass are also analytically obtained. The results are compared with those obtained by the assumption of ‘uniform increase of radial stress’ method, which is made by the previous studies. The analyses show when the bolts' spacing is large, the safety factor must be increased if the ‘uniform increase of radial stress’ method is used for the design.

Finally, a procedure is introduced to calculate the non-equal deformation of the rock mass between the bolts at any radius that can be useful to compute the bending moment in shotcrete layer in New Austrian Tunnelling Method (NATM) approach. Copyright © 2016 John Wiley & Sons, Ltd.

An investigation is made to present analytical solutions provided by a Winkler model approach for analysis of piled rafts with nodular pile subjected to vertical loads in nonhomogeneous soils. The vertical stiffness coefficient along a piled raft with the nodular pile in nonhomogeneous soils is derived from the displacement given by the Mindlin solution for elastic continuum analysis. The vertical stiffness coefficients for the bases of the raft and the nodular part in the nodular pile in a soil are expressed by the Muki solution for the 3-D elastic analysis. The relationship between settlement and vertical load on the pile base is presented considering the Mindlin solution and the equivalent thickness in the equivalent elastic method. The interaction factor between the shaft of the nodular pile and the soil is expressed taking into account the Mindlin solution and the equivalent elastic modulus. The relationship between settlement and vertical load for a piled raft with the nodular pile in nonhomogeneous soils is obtained by using the recurrence equation of influence factors of the pile for each layer. The percentage of each load carried by both nodular pile and raft subjected to vertical load is represented through the vertical influence factors proposed here. Comparison of the results calculated by the present method for piled rafts with nodular piles in nonhomogeneous soils has shown good agreement with those obtained from the finite element method and a field test. Copyright © 2016 John Wiley & Sons, Ltd.

This note presents a new method to derive closed-form expressions describing the horizontal response of an end-bearing pile in viscoelastic soil subjected to harmonic loads at its head. The soil surrounding the pile is assumed as a linearly viscoelastic layer. The propagation of waves in the soil and pile is treated mathematically by three-dimensional and one-dimensional theories, respectively. Unlike previous studies of the problem, the formulation presented allows the governing equations of the soil to be solved directly, eliminating the need to introduce potential functions. Accordingly, the dynamic response of the pile is obtained by means of the initial parameter method, invoking the requirement for continuity at the pile–soil interface. It is demonstrated that the derived compact solution matches exactly an existing solution that utilises potential functions to formulate the problem. Copyright © 2016 John Wiley & Sons, Ltd.