Soil properties are indispensable input parameters in geotechnical design and analysis. In engineering practice, particularly for projects with relatively small or medium sizes, soil properties are often not measured directly, but estimated from geotechnical design charts using results of some commonly used laboratory or *in situ* tests. For example, effective friction angle ϕ′ of soil is frequently estimated using standard penetration test (SPT) *N* values and design charts relating SPT *N* values to ϕ′. Note that directly measured ϕ′ data are generally not available when (and probably why) the use of design charts is needed. Because design charts are usually developed from past observation data, on either empirical or semi-theoretical basis, uncertainty is unavoidably involved in the design charts. This situation leads to two important questions in engineering practice: how good or reliable are the soil properties estimated in a specific site when using the design charts? (or how to measure the performance of the design charts in a specific site?); and how to incorporate rationally the model uncertainty when estimating soil properties using the design charts? This paper aims to address these two questions by developing a Bayesian statistical approach. In this paper, the second question is firstly addressed (i.e., soil properties are probabilistically characterized by rationally incorporating the model uncertainty in the design chart). Then, based on the characterization results obtained, an index is proposed to evaluate the site-specific performance of design charts (i.e., to address the first question). Equations are derived for the proposed approach, and the proposed approach is illustrated using both real and simulated SPT data. Copyright © 2016 John Wiley & Sons, Ltd.

An efficient and robust algorithm to numerically detect material instability or bifurcation is of great importance to the understanding and simulation of material failure in computational and applied mechanics. In this work, an intelligence optimizer, termed the particle swarm optimization, is introduced to the numerical solution of material bifurcation problem consisting of finding the bifurcation time as well as the corresponding bifurcation directions. The detection of material bifurcation is approached as a constrained minimization problem where the determinant of the acoustic tensor is minimized. The performance of the particle swarm optimization method is tested through numerical bifurcation analysis on both small and finite deformation material models in computational inelasticity with increasing complexity. Compared with conventional numerical approaches to detect material bifurcation, the proposed method demonstrates superior performance in terms of both computational efficiency and robustness. Copyright © 2016 John Wiley & Sons, Ltd.

Since the early work of Athey (1930), there have been many attempts to describe the various nonlinear behaviors of rocks and soils in terms of functionals having only a few parameters, while nevertheless being able to fit the complicated available data with satisfactory accuracy. Such approaches have not been universally applied however, and the present analyses are intended to draw attention to the possibility of using such nonlinear fitting methods on old as well as new data sets. In particular, some special emphasis is placed here on re-examining the well-known laboratory data of Coyner (1984) on rocks in light of such modeling tools, and we find that the nonlinear approach again has several clear advantages – especially in terms of reducing the number of variables needed to describe the observed behavior of both bulk modulus and porosity of rocks undergoing large changes in pressure. Copyright © 2016 John Wiley & Sons, Ltd.

In order to evaluate analytically the ITZ volume fraction (*f _{ITZ}*) in concrete, a three phase model is proposed for the random concrete microstructure using the Voronoï tessellation. Within this model, the ITZ local thickness is a statistical variable depending on the local paste thickness available between each couple of neighbouring aggregates. The

Evaluating the induced subsidence is a critical step in multi-seam longwall mining. Numerical modelling can be a cost-effective approach to this problem. Numerical evaluation of longwall mining-induced subsidence is much more complicated when more than one seam is to be extracted. Only a few research works have dealt with this problem. This paper discusses the essential requirements of a robust numerical modelling approach to simulation of multi-seam longwall mining-induced subsidence. In light of these requirements, the previous works on this topic are critically reviewed. A simple yet robust FEM-based modelling approach is also proposed that is capable of simulating caving process, rock mass deterioration and subsidence around multi-seam excavations. The effectiveness of this approach in comparison with two other conventional FEM approaches is demonstrated through numerical examples of two different multi-seam mining configurations. Results show that the proposed numerical modelling approach is the only robust method, which is capable of simulating multi-seam subsidence in both demonstrated cases. Copyright © 2016 John Wiley & Sons, Ltd.

This paper presents a new simplified method, based on Hypothesis B, for calculating the consolidation settlements of double soil layers exhibiting creep. In the new simplified Hypothesis B method, different stress–strain states including over-consolidation and normal consolidation states can be considered with the help of the ‘equivalent time’ concept. Zhu and Yin method and US Navy method are adopted to calculate the average degree of consolidation for a double soil layer profile. This new simplified Hypothesis B method is then used to calculate the consolidation settlements of double soil layers, which have two different total thicknesses of soil layer (4 m and 8 m) and three different OCR values (Over-Consolidation Ratio, *OCR* = 1, 1.5 and 2). The accuracy and verification of this new simplified method are examined by comparing the calculated results with simulation results from a fully coupled finite element (FE) program using a soft soil creep model. Four cases of double layer soil profiles are analyzed. Hypothesis A method with US Navy method for the average degree of consolidation has also been used to for calculating consolidation settlements of the same cases. For *Case I*(4*m*) and *Case III*(8*m*), it is found that curves of the new simplified Hypothesis B method using both Zhu and Yin method and US Navy method are very close to the results from FE simulations with the *relative errors* within 8.5%. For *Case II*(4*m*) and *Case IV*(8*m*), it is found that curves of the new simplified Hypothesis B method using Zhu and Yin method agree better with results from FE simulations with the *relative errors* within 11.7% than curves of the new simplified Hypothesis B method adopting US Navy method with the *relative error* up to 36.1%. Curves of Hypothesis A method adopting US Navy method have the *relative error* up to 55.0% among all four cases. In overall, the new simplified Hypothesis B method is suitable for calculation of consolidation settlements of double soil layers exhibiting creep, in which, Zhu and Yin method is recommended to obtain the average degree of consolidation. Copyright © 2016 John Wiley & Sons, Ltd.

For civil engineering structures with a tightness role, structural permeability is a key issue. In this context, this paper presents a new proposition of a numerical modelling of leakage rate through a cracked concrete structure undergoing mode I cracking. The mechanical state of the material, considered in the framework of continuum mechanics based on finite element modelling, is described by means of the stress-based nonlocal damage model which takes into account the stress state and provides realistic local mechanical fields. A semi-discrete method based on the strong discontinuity approach to estimate crack opening is then considered in the post-treatment phase. Using a Poiseuille's like relation, the coupling between the mechanical state of the material and its dry gas conductivity is performed.

For validation purposes, an original experimental campaign is conducted on a dry concrete disc loaded in a splitting setup. During the loading, gas conductivity and digital image correlation analysis are performed. The comparison with the 3D experimental mechanical global response highlights the performance of the mechanical model. The comparison between crack openings measured by digital image correlation and estimated by the strong discontinuity method shows a good agreement.

Finally, the results of the semi-discrete approach coupled with the gas conductivity compared with experimental data show a good estimation of the structural conductivity. Consequently, if the mechanical problem is well modelled at the global scale, then the proposed approach provides good estimation of gas conductivity. Copyright © 2016 John Wiley & Sons, Ltd.

The paper proposes a stress-driven integration strategy for Perzyna-type viscoplastic constitutive models, which leads also to a convenient algorithm for viscoplastic relaxation schemes. A generalized trapezoidal rule for the strain increment, combined with a linearized form of the yield function and flow rules, leads to a plasticity-like compliance operator that can be explicitly inverted to give an algorithmic tangent stiffness tensor also denoted as the m-AGC tangent operator. This operator is combined with the stress-prescribed integration scheme, to obtain a natural error indicator that can be used as a convergence criterion of the intra-step iterations (in physical viscoplasticity), or to a variable time-step size in viscoplastic relaxation schemes based on a single linear calculation per time step. The proposed schemes have been implemented for an existing zero-thickness interface constitutive model. Some numerical application examples are presented to illustrate the advantages of the new schemes proposed. Copyright © 2016 John Wiley & Sons, Ltd.

The microstructure of rock was numerically reproduced by a polygonal grain-based model, and its mechanical behavior was examined by performing the uniaxial compression test and Brazilian tests via the Universal Distinct Element Code. The numerical results of the model demonstrated good agreement with the experimental results obtained with rock specimens in terms of the stress–strain behavior, strength characteristics, and brittle fracture phenomenon. An encouraging result is that the grain-based model-Universal Distinct Element Code model can reproduce a low ratio of tensile to compressive strength of 1/20 to 1/10 without the need for an additional process. This finding is ascribed to the fact that the geometrical features of polygons can effectively capture the effects of angularity, finite rotation, and interlocking of grains that exist in reality. A numerical methodology to monitor the evolution of micro-cracks was developed, which enabled us to examine the progressive process of the failure and distinguish the contribution of tensile cracking to the process from that of shear cracking. From the observations of the micro-cracking process in reference to the stress–strain relation, crack initiation stress, and crack damage stress, it can be concluded that the failure process of the model closely resembles the microscopic observations of rock. We also carried out a parametric study to examine the relationships between the microscopic properties and the macroscopic behavior of the model. Depending on the micro-properties, the model exhibited a variety of responses to the external load in terms of the strength and deformation characteristics, the evolution of micro-cracks, and the post-peak behavior. Copyright © 2016 John Wiley & Sons, Ltd.

The present study investigates propagation of a cohesive crack in non-isothermal unsaturated porous medium under mode I conditions. Basic points of skeleton deformation, moisture, and heat transfer for unsaturated porous medium are presented. Boundary conditions on the crack surface that consist of mechanical interaction of the crack and the porous medium, water, and heat flows through the crack are taken into consideration. For spatial discretization, the extended finite element method is used. This method uses enriched shape functions in addition to ordinary shape functions for approximation of displacement, pressure, and temperature fields. The Heaviside step function and the distance function are exploited as enrichment functions for representing the crack surfaces displacement and the discontinuous vertical gradients of the pressure and temperature fields along the crack, respectively. For temporal discretization, backward finite difference scheme is applied. Problems solved from the literature show the validity of the model as well as the dependency of structural response on the material properties and loading. Copyright © 2016 John Wiley & Sons, Ltd.

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

The ultimate bearing capacity problem of column-reinforced foundations under inclined loading is investigated within the framework of static and kinematic approaches of yield design theory. The configuration of a native soft clayey soil reinforced by either a group of purely cohesive columns (lime-column technique) or a group of purely frictional columns (stone-column technique) is analyzed under plane strain conditions. First, lower bound estimates are derived for the ultimate bearing capacity by considering statically admissible piecewise linear stress distributions that comply with the local strength conditions of the constitutive materials. The problem is then handled by means of the yield design kinematic approach of limit analysis through the implementation of several failure mechanisms, allowing the formulation of upper bound estimates for the ultimate bearing capacity. A series of finite element limit load solutions obtained from numerical elastoplastic simulations suggests that the predictions derived from the kinematic approach appear to be more accurate than the estimates obtained from the static approach. Comparison with available results obtained in the context of yield design homogenization demonstrates the accuracy of the proposed direct analysis, which may therefore be viewed as complementary approach to homogenization-based approaches when a small number of columns is involved. Copyright © 2016 John Wiley & Sons, Ltd.

The theory of fractional calculus has been successfully applied to model the triaxial behaviour of soils under static loading conditions. However, limited work has been carried out in using the fractional calculus to describe the cyclic behaviour of granular soils. In this paper, a fractional order constitutive model for granular soils under drained cyclic loading is proposed by incorporating the concept of fractional rate for strain accumulation. The fractional rate for strain accumulation is obtained from the analysis of the experimental data by utilizing the fractional calculus. Comparison between the test results and model predictions is presented. The key feature of the proposed model is that it can reasonably characterise the cyclic deformation of granular soils under both low and high loading cycles. Copyright © 2016 John Wiley & Sons, Ltd.

Because of a multitude of steep slopes being constructed adjacent to roadways, there is greater concern of landslide occurrence, particularly in instances where poor geomaterials are present. Installation of piles along the slope is one commonly adopted method. This paper presents the assessment of the stability of a rock slope with stabilizing piles based on kinematic analysis. The pile effect is introduced with a resultant lateral force and a moment. Upper bound solutions of the pile's lateral force are derived with a log-spiral rotational failure mechanism. The slope performance based on the bearing capacity of surcharge loading is also discussed with consideration of pore water pressure. In order to substantiate the derived theoretical solutions, numerical analysis with optimization technique is carried out. Results demonstrate that rock materials with high quality are conducive to ensure slope stability. Reduced lateral force on the pile is produced with lower rock weight, slope height, and surcharge loading. Finally, the safety factor and stability coefficient are discussed to complete the evaluation of the slope stability. Copyright © 2016 John Wiley & Sons, Ltd.

The storage of petroleum products above ground surface has many constraints and limitations. A viable alternative is to excavate large underground spaces in rock to provide a safer way for oil storage. Soft rock formations such as salt domes provide suitable conditions from environmental and operational aspects. The potential for high volume storage and low permeability are among advantages of oil storage in caverns excavated in salt rocks. The complicated shape of oil storage caverns, complex behavior of salt rock, and boundary conditions associated with large underground openings are major challenges in the design of salt caverns excavated for oil storage purposes. In this study, the deformation mechanism and stability of salt caverns were investigated. A comprehensive 3D numerical study was carried out to investigate the effects of cavern size and depth, salt rock deformation modulus, and ground in-situ stress regime on the behavior of large salt caverns. The stress field and deformation mechanisms were studied numerically aiming at shedding lights into the design aspects of salt caverns for oil storage. The analysis results show that the cavern safety factor is reduced as a function of cavern depth and storage volume. Also, with decrease in *k* (ratio of horizontal to vertical in-situ stress), the stability of salt caverns will increase; however, with increase in salt rock young modulus, the sensitivity of cavern stability to *k* ratio is reduced. Copyright © 2016 John Wiley & Sons, Ltd.

An elasto-viscoplastic constitutive model for asphaltic materials is presented within the context of bounding surface plasticity theory, taking into account the effects of the stress state, void binder degree of saturation, temperature and strain rate on the material behaviour. A stress state dependent non-linear elasticity model is introduced to represent time-independent recoverable portion of the deformation. The consistent visco-plasticity framework is utilised to capture the rate-dependent, non-recoverable strain components. The material parameters introduced in the model are identified, and their determination from conventional laboratory tests is discussed. The capability of the model to reproduce experimentally observed response of asphaltic materials is demonstrated through numerical simulations of several laboratory test data from the literature. Copyright © 2016 John Wiley & Sons, Ltd.

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

This paper describes a fully coupled finite element/finite volume approach for simulating field-scale hydraulically driven fractures in three dimensions, using massively parallel computing platforms. The proposed method is capable of capturing realistic representations of local heterogeneities, layering and natural fracture networks in a reservoir. A detailed description of the numerical implementation is provided, along with numerical studies comparing the model with both analytical solutions and experimental results. The results demonstrate the effectiveness of the proposed method for modeling large-scale problems involving hydraulically driven fractures in three dimensions. © 2016 The Authors. International Journal for Numerical and Analytical Methods in Geomechanics published by John Wiley & Sons Ltd.

This paper integrates random field simulation of soil spatial variability with numerical modeling of coupled flow and deformation to investigate consolidation in spatially random unsaturated soil. The spatial variability of soil properties is simulated using the covariance matrix decomposition method. The random soil properties are imported into an interactive multiphysics software COMSOL to solve the governing partial differential equations. The effects of the spatial variability of Young's modulus and saturated permeability together with unsaturated hydraulic parameters on the dissipation of excess pore water pressure and settlement are investigated using an example of consolidation in a saturated-unsaturated soil column because of loading. It is found that the surface settlement and the pore water pressure profile during the process of consolidation are significantly affected by the spatially varying Young's modulus. The mean value of the settlement of the spatially random soil is more than 100% greater than that of the deterministic case, and the surface settlement is subject to large uncertainty, which implies that consolidation settlement is difficult to predict accurately based on the conventional deterministic approach. The uncertainty of the settlement increases with the scale of fluctuation because of the averaging effect of spatial variability. The effects of spatial variability of saturated permeability *k _{sat}* and air entry parameters are much less significant than that of elastic modulus. The spatial variability of air entry value parameters affects the uncertainties of settlement and excess pore pressure mostly in the unsaturated zone. Copyright © 2016 John Wiley & Sons, Ltd.

This paper focuses on the sensitivity analysis for coupled thermo–hydro–mechanical problems employing both local and global sensitivity methods. A derivative-based method is used in the local sensitivity approach, whereas the random balance designs method is used for the global sensitivity analysis. The main goal is to investigate the effect of uncertainties in the constitutive parameters on the results from nonlinear coupled thermo–hydro–mechanical analyses of unsaturated soil behavior whose modeling generally involves large sets of constitutive relations. Knowing the parameter sensitivity allows to qualitatively assess the validity of the results obtained by computational simulations of high-risk situations, for example, emerging nuclear waste repositories. Copyright © 2016 John Wiley & Sons, Ltd.

Elastic lateral dynamic impedance functions are defined as the ratio of the lateral dynamic force/moment to the corresponding lateral displacement/rotation at the top ending of a foundation at very small strains. Elastic lateral dynamic impedance functions have a defining influence on the natural frequencies of offshore wind turbines supported on cylindrical shell type foundations, such as suction caissons, bucket foundations, and monopiles. This paper considers the coupled horizontal and rocking vibration of a cylindrical shell type foundation embedded in a fully saturated poroelastic seabed in contact with a seawater half-space. The formulation of the coupled seawater–shell–seabed vibration problem is simplified by treating the shell as a rigid one. The rigid shell vibration problem is approached by the integral equation method using ring-load Green's functions for a layered seawater-seabed half-space. By considering the boundary conditions at the shell–soil interface, the shell vibration problem is reduced to Fredholm integral equations. Through an analysis of the corresponding Cauchy singular equations, the intrinsic singular characteristics of the problem are rendered explicit. With the singularities incorporated into the solution representation, an effective numerical method involving Gauss–Chebyshev method is developed for the governing Fredholm equations. Selected numerical results for the dynamic contact load distributions, displacements of the shell, and lateral dynamic impedance functions are examined for different shell length–radius ratio, poroelastic materials, and frequencies of excitation. Copyright © 2016 John Wiley & Sons, Ltd.

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. © 2016 The Authors. International *Journal for Numerical and Analytical Methods in Geomechanics* published by 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.

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.

No abstract is available for this article.

]]>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.

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.

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.

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.

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.

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.

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.