In this paper, we used a theoretical model for the variation of Eulerian porosity, which takes into account the adsorption process known to be the main mechanism of production or sequestration of gas in many reservoir of coal. This process is classically modeled using Langmuir's isotherm. After implementation in Code_Aster, a fully coupled thermo-hydro-mechanical analysis code for structures calculations, we used numerical simulations to investigate the influence of coal's hydro-mechanical properties (Biot's coefficient, bulk modulus), Langmuir's adsorption parameters, and the initial liquid pressure in rock mass during CO_{2} injection in coal. These simulations showed that the increase in the values of Langmuir's parameters and Biot's coefficient promotes a reduction in porosity because of the adsorption process when the gas pressure increases. Low values of bulk modulus increase the positive effect (i.e., increase) of hydro-mechanical coupling on the porosity evolution. The presence of high initial liquid pressure in the rock mass prevents the progression of injected gas pressure when CO_{2} dissolution in water is taken into account. Copyright © 2014 John Wiley & Sons, Ltd.

Numerical modeling has now become an indispensable tool for investigating the fundamental mechanisms of toxic nonaqueous phase liquid (NAPL) removal from contaminated groundwater systems. Because the domain of a contaminated groundwater system may involve irregular shapes in geometry, it is necessary to use general quadrilateral elements, in which two neighbor sides are no longer perpendicular to each other. This can cause numerical errors on the computational simulation results due to mesh discretization effect. After the dimensionless governing equations of NAPL dissolution problems are briefly described, the propagation theory of the mesh discretization error associated with a NAPL dissolution system is first presented for a rectangular domain and then extended to a trapezoidal domain. This leads to the establishment of the finger-amplitude growing theory that is associated with both the corner effect that takes place just at the entrance of the flow in a trapezoidal domain and the mesh discretization effect that occurs in the whole NAPL dissolution system of the trapezoidal domain. This theory can be used to make the approximate error estimation of the corresponding computational simulation results. The related theoretical analysis and numerical results have demonstrated the following: (1) both the corner effect and the mesh discretization effect can be quantitatively viewed as a kind of small perturbation, which can grow in unstable NAPL dissolution systems, so that they can have some considerable effects on the computational results of such systems; (2) the proposed finger-amplitude growing theory associated with the corner effect at the entrance of a trapezoidal domain is useful for correctly explaining why the finger at either the top or bottom boundary grows much faster than that within the interior of the trapezoidal domain; (3) the proposed finger-amplitude growing theory associated with the mesh discretization error in the NAPL dissolution system of a trapezoidal domain can be used for quantitatively assessing the correctness of computational simulations of NAPL dissolution front instability problems in trapezoidal domains, so that we can ensure that the computational simulation results are controlled by the physics of the NAPL dissolution system, rather than by the numerical artifacts. Copyright © 2014 John Wiley & Sons, Ltd.

Water retention in compacted clays is dominated by multi-modal pore size distribution which evolves during hydro-mechanical paths depending on water content and stress history. A description of the evolutionary fabric has been recently introduced in models for water retention, but mostly on a heuristic base. Here, a possible systematic approach to account for evolving pore size distribution is presented, and its implications in models for water retention are discussed. The approach relies on quantitative information derived from mercury intrusion porosimetry data. The information is exploited to introduce physically based evolution laws for the parameters of water retention models. These laws allow tracking simultaneously the evolution of the aggregated fabric and the consequent hydraulic state of compacted clays. The influence of clay microstructure, mechanical constraints and water content changes on the water retention properties is highlighted and quantified from experimental data on different compacted soils with different activity of the clayey fraction. The framework is discussed with reference to a widespread water retention model and validated against experimental data on a Sicilian scaly clay compacted to different dry densities and subjected to a number of hydro-mechanical paths. Copyright © 2014 John Wiley & Sons, Ltd.

In this work, we investigate the main pumping parameters that influence a fluid-driven fracture in cohesive poroelastic and poroelastoplastic weak formations. These parameters include the fluid viscosity and the injection rate. The first parameter dominates in the mapping of the propagation regimes from toughness to viscosity, whereas the second parameter controls the storage to leak-off dominated regime through diffusion.

The fracture is driven in weak permeable porous formation by injecting an incompressible viscous fluid at the fracture inlet assuming that the fracture propagates under plane strain conditions. Fluid flow in the fracture is modeled by lubrication theory. Pore fluid movement in the porous formation is based on the Darcy law. The coupling follows the Biot theory, whereas the irreversible rock deformation is modeled with the Mohr–Coulomb yield criterion with associative flow rule. Fracture propagation criterion is based on the cohesive zone approach. Leak-off is also considered. The investigation is performed numerically with the FEM to obtain the fracture opening, length, and propagation pressure versus time.

We demonstrate that pumping parameters influence the fracture geometry and fluid pressures in weak formations through the viscous fluid flow and the diffusion process that create back stresses and large plastic zones as the fracture propagates. It is also shown that the product of the propagation velocity and fluid viscosity, *µv* that appears in the scaling controls the magnitude of the plastic zones and influences the net pressure and fracture geometry. These findings may explain partially the discrepancies in net pressures between field measurements and conventional model predictions for the case of weak porous formation. Copyright © 2014 John Wiley & Sons, Ltd.

The failure process of brittle rock submitted to a compression state of stress with different confining pressures is investigated in this paper based on discrete element method (DEM) simulations. In the DEM model, the rock sample is represented by bonding rigid particles at their contact points. The numerical model is first calibrated by comparing the macroscopic response with the macroscopic response of Beishan granite obtained from laboratory tests. After the validation of numerical model in terms of macroscopic responses, the failure process of the DEM model under unconfined and confined compression is studied in micro-scale in detail. The contact force network and its relation to the development of micro-cracks and evolution of major fractures are studied. Confining pressure will prohibit the development of tensile cracks and hence alter the failure patterns. An in-depth analysis of micro-scale response is carried out, including the orientation distribution and probability density of stress acting on parallel bonds, the effect of particle size heterogeneity on bond breakage and the evolution of fabric tensor and coordination number of parallel bond. The proposed micromechanical analysis will allow us to extract innovative features emerged from the stresses and crack evolution in brittle rock failure process. Copyright © 2014 John Wiley & Sons, Ltd.

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

The aim of this study is to investigate the effect of pre-existing, or structural, cracks on dynamic fragmentation of granite. Because of the complex behavior of rock materials, a continuum approach is employed relying upon a plasticity model with yield surface locus as a quadratic function of the mean pressure in the principal stress space coupled with an anisotropic damage model. In particular, Bohus granite rock is investigated, and the material parameters are chosen based on previous experiments. The equation of motion is discretized using a finite element approach, and the explicit time integration method is employed. The pre-existing cracks are introduced in the model by considering sets of elements with negligible tensile strength that leads to their immediate failure when loaded in tension even though they still carry compressive loads as crack closure occurs because of compressive stresses. Previously performed edge-on impact tests are reconsidered here to validate the numerical model. Percussive drilling is simulated, and the influence of the presence of pre-existing cracks is studied. The results from the analysis with different crack lengths and orientations are compared in terms of penetration stiffness and fracture pattern. It is shown that pre-existing cracks in all investigated cases facilitate the drilling process. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents a new method to derive the analytical solution for the vertical impedance of an end-bearing pile in viscoelastic soil. The soil is assumed as a homogeneous and isotropic layer, and the pile is considered as a one-dimensional Euler rod. Considering both the vertical and radial displacements of soil and soil–pile coupled vibration, the governing equations of the soil and pile are established. The volumetric strain of soil is obtained by transformation on the equations of soil and variable separation method. Then the vertical and radial displacements of soil are obtained accordingly. The displacement response and impedance function of pile are derived based on the continuity assumption of the displacement and stress between the pile and soil. The solution is verified by being compared with an existing solution obtained by introducing potential functions. Furthermore, a comparison with two other simplified solutions is conducted. Numerical examples are presented to analyze the vibration characteristics of the pile. Copyright © 2014 John Wiley & Sons, Ltd.

A new data-mining approach is presented for modelling of the stress–strain and volume change behaviour of unsaturated soils considering temperature effects. The proposed approach is based on the evolutionary polynomial regression (EPR), which unlike some other data-mining techniques, generates a transparent and structured representation of the behaviour of systems directly from raw experimental (or field) data. The proposed methodology can operate on large quantities of data in order to capture nonlinear and complex relationships between contributing variables. The developed models allow the user to gain a clear insight into the behaviour of the system. Unsaturated triaxial test data from the literature were used for development and verification of EPR models. The developed models were also used (in a coupled manner) to produce the entire stress path of triaxial tests. Comparison of the EPR model predictions with the experimental data revealed the robustness and capability of the proposed methodology in capturing and reproducing the constitutive thermomechanical behaviour of unsaturated soils. More importantly, the capability of the developed models in accurately generalizing the predictions to unseen data cases was illustrated. The results of a sensitivity analysis showed that the models developed from data are able to capture and represent the physical aspects of the unsaturated soil behaviour accurately. The merits and advantages of the proposed methodology are also discussed. Copyright © 2014 John Wiley & Sons, Ltd.

The hydraulic fracturing propagation regimes in the plane strain model are uniformly investigated using a numerical method based on the finite element method. The regimes range from toughness-dominated cases to viscosity-dominated cases, covering zero leak-off situations and small leak-off situations. Unlike the asymptotic solutions, the numerical method is independent of the energy dissipation regimes and fluid storage regimes. The numerical method pays no special attention to the fracture tip, and it simulates fracture tip behaviors by increasing the number of functions in a natural and uniform manner. The numerical method is verified by comparing its results with the asymptotic solutions. The effect of the model sizes on the numerical method is discussed along with the robustness of the numerical method. Copyright © 2014 John Wiley & Sons, Ltd.

The interaction between twin-parallel tunnels affects the tunnelling-induced ground deformation, which may endanger the nearby structures. In this paper, an analytical solution is presented for problems in determining displacements and stresses around deforming twin-parallel tunnels in an elastic half plane, on the basis of complex variable theory. As an example, a uniform radial displacement was assumed as the boundary condition for each of the two tunnels. Special attention was paid to the effects of tunnel depth and spacing between the two tunnels on the surface movement to gain deep insight into the effect of the interaction between twin-parallel tunnels using the proposed analytical approach. It is revealed that the influence of twin tunnel interaction on surface movements diminishes with both the increase of the tunnel depth and the spacing between the two tunnels. The presented analytical solution manifests that, similar to most of the existing numerical results, the principle of superposition can be applied to determine ground deformation of twin-parallel tunnels with a certain large depth and spacing; otherwise, the interaction effect between the two tunnels should be taken into account for predicting reliable ground movement. Copyright © 2014 John Wiley & Sons, Ltd.

The paper is motivated by the long-term safety analysis of the CO_{2} geological storage. We present a methodology for the assessment of the geomechanical impact of progressive rock dissolution. The method is based on the use of X-ray tomography and the numerical dissolution technique. The influence of evolution of the microstructure on the macroscopic properties of the rock is analysed by using periodic homogenization method. The numerical computations show progressive degradation of all components of the stiffness (orthotropic) tensor. Moreover, the evolution of associated mass transfer properties (as tortuosity and conductivity tensors), by using the periodic homogenization method, is also calculated. The correlation between the mechanical parameters and the transfer properties during the dissolution process is presented. The results show that the highest increase of the hydraulic conductivity (in direction *Y*) is not associated with the highest decrease of Young modulus in this direction. Moreover, the highest decrease of Young modulus (in the direction *X*) is not associated with percolation in this direction. Finally, an incremental law to calculate settlement, in case of a rock with evolving microstructure, is proposed. The solution of the macroscopic settlement problem under constant stress and drained conditions showed that the geomechanical effects of the rock dissolution are rather limited. Copyright © 2014 John Wiley & Sons, Ltd.

In this work, the macroscopic Richards equation for moisture transport is established in unsaturated porous media using periodic homogenization. By performing dimensional analysis on microscopic equations of moisture transfer, dimensional numbers characterizing moisture transport appear. The application of the asymptotic homogenization leads to the classical Richards equation, which is justified rigorously this way. Moreover, we obtain an accurate definition of the homogenized diffusion tensor of moisture involving the geometric properties of the microstructure and known transport properties of the material. A different behavior for the transport of water vapor between hygroscopic and super-hygroscopic region is revealed. Finally, a simple 2D example where an analytical solution exists is addressed. Copyright © 2014 John Wiley & Sons, Ltd.

Within the framework of our discontinuous deformation analysis for rock failure algorithm, this paper presents a two-dimensional coupled hydromechanical discontinuum model for simulating the rock hydraulic fracturing process. In the proposed approach, based on the generated joint network, the calculation of fluid mechanics is performed first to obtain the seepage pressure near the tips of existing cracks, and then the fluid pressure is treated as linearly distributed loads on corresponding block boundaries. The contribution of the hydraulic pressure to the initiation/propagation of the cracks is considered by adding the components of these blocks into the force matrix of the global equilibrium equation. Finally, failure criteria are applied at the crack tips to determine the occurrence of cracking events. Several verification examples are simulated, and the results show that this newly proposed numerical model can simulate the hydraulic fracturing process correctly and effectively. Although the numerical and experimental verifications focus on one unique preexisting crack, because of the capability of discontinuous deformation analysis in simulating block-like structures, the proposed approach is capable of modeling rock hydraulic fracturing processes. Copyright © 2014 John Wiley & Sons, Ltd.

The present study pertains to the development of a foundation model for predicting the behavior of geosynthetic reinforcement railway track system rested on soft clay subgrade. The ballast and sub-ballast layers have been idealized by Pasternak shear layer. The geosynthetic layer is represented by a stretched rough elastic membrane. Burger model has been used to characterize the soft clay subgrade. Numerical solutions have been obtained by adopting the finite difference scheme combined with non-dimensioning the governing equations of the proposed model. The results confirm that the present model is quite capable of predicting the time-dependent settlement response of geosynthetic reinforcement railway track system placed on soft clay subgrade. The surface settlement profile and mobilized tensile load of geosynthetics has been evaluated by considering variation in the wheel load, sleeper width, thickness of ballast and sub-ballast layers and shear modulus of ballast and sub-ballast layers. It has been observed that an increase in the sleeper width by 24% results in the reduction in central settlement and mobilized tensile load by 6.5% and 20.1%, respectively. It was found that with a 50% increase in the thickness of the ballast layer, the central settlement has decreased by 7.3% and the mobilized tension at the zone of maximum curvature has increased by 24.6%. However, with an increase in the thickness of the sub-ballast layer, a considerable reduction in both central settlement and the mobilization of tension on geosynthetic has been noticed. The pattern of variation of settlement and mobilized tension for an increase in the shear modulus of ballast and sub-ballast material was found to be almost similar. Copyright © 2014 John Wiley & Sons, Ltd.

In this work, we describe a meshless numerical method based on local collocation with RBFs for the solution of the poroelasticity equation. The RBF finite collocation approach forms a series of overlapping nodal stencils, over which an RBF collocation is performed. These local collocation systems enforce the governing PDE operator throughout their interior, with the intersystem communication occurring via the collocation of field variables at the stencil periphery. The method does not rely on a generalised finite differencing approach, whereby the governing partial differential operator is reconstructed at the global level to drive the solution of the PDE. Instead, the PDE governing and boundary operators are enforced directly within the local RBF collocation systems, and the sparse global assembly is formed by reconstructing the value of the field variables at the centrepoint of the local stencils. In this way, the solution of the PDE is driven entirely by the local RBF collocation, and the method more closely resembles the approach of the full-domain RBF collocation method. By formulating the problem in this fashion, high rates of convergence may be attained without the computational cost and numerical ill-conditioning issues that are associated with the full-domain RBF collocation approach.

An analytical solution is formulated for a 2D poroelastic fluid injection scenario and is used to verify the proposed implementation of the method. Highly accurate solutions are produced, and convergence rates in excess of sixth order are observed for each field variable (i.e. pressure and displacement) and field-variable derivative (i.e. pressure gradients and stresses). The stress and displacement fields resulting from the solution of the poroelasticity equation are then used to describe the formation and propagation of microfractures and microfissures, which may form in the presence of large shear strain, in terms of a continuous damage variable which modifies the mechanical and hydraulic properties of the porous medium. The formation of such hydromechanical damage, and the resulting increase in hydraulic conductivity, is investigated for a pressurised injection into sandstone. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents a numerical scheme for fluid-particle coupled discrete element method (DEM), which is based on poro-elasticity. The motion of the particles is resolved by means of DEM. While within the proposition of Darcian regime, the fluid is assumed as a continuum phase on a Eulerian mesh, and the continuity equation on the fluid mesh for a compressible fluid is solved using the FEM. Analytical solutions of traditional soil mechanics examples, such as the isotropic compression and one-dimensional upward seepage flow, were used to validate the proposed algorithm quantitatively. The numerical results showed very good agreement with the analytical solutions, which show the correctness of this algorithm. Sensitivity studies on the effect of some influential factors of the coupling scheme such as pore fluid bulk modulus, volumetric strain calculation, and fluid mesh size were performed to display the accuracy, efficiency, and robustness of the numerical algorithm. It is revealed that the pore fluid bulk modulus is a critical parameter that can affect the accuracy of the results. Because of the iterative coupling scheme of these algorithms, high value of fluid bulk modulus can result in instability and consequently reduction in the maximum possible time-step. Furthermore, the increase of the fluid mesh size reduces the accuracy of the calculated pore pressure. This study enhances our current understanding of the capacity of fluid-particle coupled DEM to simulate the mechanical behavior of saturated granular materials. Copyright © 2014 John Wiley & Sons, Ltd.

Irregularly shaped (IRS) particles widely exist in many engineering and industrial fields. The macro physical and mechanical properties of the particle system are governed by the interaction between the particles in the system. The interaction between IRS particles is more complicated because of their complex geometric shape with extremely irregular and co-existed concave and convex surfaces. These particles may interlock each other, making the sliding and friction of IRS particles more complex than that of particles with regular shape. In order to study the interaction of IRS particles more efficiently, a refined method of constructing discrete element model based on computed tomography scanning of IRS particles is proposed. Three parameters were introduced to control the accuracy and the number of packing spheres. Subsequently, the inertia tensor of the IRS particle model was optimized. Finally, laboratory and numerical open bottom cylinder tests were carried out to verify the refined modeling method. The influence of particle shape, particle position, and mesoscopic friction coefficient on the interaction of particles was also simulated. It is noteworthy that with the increase of mesoscopic friction coefficient, the fluidity of IRS particle assembly decreases, and intermittent limit equilibrium state may appear. Copyright © 2014 John Wiley & Sons, Ltd.

Geomechanical models are often used to predict the impact on land surface of fluid withdrawal from deep reservoirs, as well as investigating measures for mitigation. The ability to accurately simulate surface displacements, however, is often impaired by limited information on the geomechanical parameters characterizing the geological formations of interest. In this study, we employ an ensemble smoother, a data assimilation algorithm, to provide improved estimates of reservoir parameters through assimilation of measurements of both horizontal and vertical surface displacement into geomechanical model results. The method leverages the demonstrated potential of remote sensing techniques developed in the last decade to provide accurate displacement data for large areas of the land surface.

For evaluation purposes, the methodology is applied to the case of a disk-shaped reservoir embedded in a homogeneous, isotropic, and linearly elastic half space, subject to a uniform change in fluid pressure. Multiple sources of uncertainty are investigated, including the radius, *R*, the thickness, *h*, and the depth, *c*, of the reservoir; the pore pressure change, *Δp*; porous medium's vertical uniaxial compressibility, *c _{M}*, and Poisson's ratio,

A methodology is developed in SPH framework to analyze the behavior of preexisting multiple intersecting discontinuities or joints in rock material. The procedure does not require any additional unknowns to represent discontinuities and to capture velocity jump across them. Instead, a discontinuity is represented by a set of joint particles placed along the discontinuity plane, in which relative velocity and traction vector is evaluated, obeying the Mohr–Coulomb friction law with zero tension constrain. For failure of continuous rock material, the Drucker–Prager yield criterion with tensile cracking is employed in the elastic-plastic constitutive model. Free-sip, no-slip, and symmetric boundary conditions are also implemented in SPH framework for proper representation of physical system.

The paper analyzes behavior of a rock sample having a discontinuity plane under uniaxial loading and compares velocity and stress with a theoretical solution derived considering effective vertical stiffness of the joint planes. The efficacy of the proposed method is successfully demonstrated by solving another two problems of jointed rock mass under uniaxial and gravitational loading conditions.Copyright © 2014 John Wiley & Sons, Ltd.

An approximate solution for the effects of high strain rates, and gradual strength degradation, on the penetration resistance of penetrometers can be obtained by combining the strain-path method with the classical upper bound theorem. The stream path calculations require the integration of the material constitutive equation along the streamlines. Unless the geometry is simple so that the integration can be evaluated analytically, numerical procedures are required to backtrack streamlines. The strain at any location is calculated by finding the streamline that passes through the given point and integrating the strain rate along that streamline from its inlet boundary. Thus, the calculations can be complicated, and errors can be accumulated during the calculation procedure.

This paper presents an efficient approach for evaluating cumulative strains around penetrometers without the need to backtrack individual streamlines. In this approach, the strain components are treated as field variables. The global solution is obtained using the streamline upwind Petrov–Galerkin method. The new method together with an Eulerian-based finite element formulation was used to study the cone penetration test and evaluate the effect of strain softening on the cone resistance. Copyright © 2014 John Wiley & Sons, Ltd.

The strength parameter *m _{i}* in the Hoek–Brown strength criterion is empirical and was developed by trial and error. To better understand the fundamental relationship between

It is well known that for a sufficiently high seepage velocity, the governing flow law of porous media is nonlinear (*J*. *Computers* & *Fluids* 2010; **39**: 2069–2077). However, this fact has not been considered in the studies of soil-pore fluid interaction and in conventional soil mechanics. In the present paper, a fully explicit dynamic finite element method is developed for nonlinear Darcy law. The governing equations are expressed for saturated porous media based on the extension of the Biot (*J*. *Appl*. *Phys*. 1941; **12**: 155–164) formulation. The elastoplastic behavior of soil under earthquake loading is simulated using a generalized plasticity theory that is composed of a yield surface along with non-associated flow rule. Numerical simulations of porous media subjected to horizontal and vertical components of ground motion excitations with different permeability coefficients are carried out; while computed maximum pore water pressure is specially taken into consideration to make the difference between Darcy and non-Darcy flow regimes tangible. Finally, the effect of non-Darcy flow on the evaluated liquefaction potential of sand in comparison to conventional Darcy law is examined. Copyright © 2014 John Wiley & Sons, Ltd.

Slope stability analysis of soil with a weak layer sandwiched between two strong layers is considered as a complex geotechnical problem. In this problem, the objective function is non-convex and discontinuous with the presence of multiple strong local minima. Classical optimization techniques fail to converge to a valid solution unless a proper initial trial is adopted. Even though many new optimization algorithms have emerged, they have not been applied to geotechnical problems yet. In the present study, some recent swarm intelligence algorithms are adopted for some complicated example of slope stability problems and benchmarked with the traditional particle swarm optimization algorithm. From the results, it seems the levy flight krill herd algorithm is the most efficient method over proposed algorithms for this kind of problem. Copyright © 2014 John Wiley & Sons, Ltd.

This paper addresses the geotechnical engineering problem of evaluating the ultimate bearing capacity of a strip foundation resting upon a reinforced soil, by means of the yield design homogenization approach. The analysis is notably focused on the determination of the macroscopic strength criterion of such reinforced soils, where both constituents are purely cohesive, which can be conveniently expressed through the notion of anisotropic cohesion. A comprehensive comparison is made between the classical configuration of reinforcing columns and the more original one of orthogonal reinforcing trenches. Among the most outstanding results of the analysis is the conclusion that the cross trench configuration is notably more efficient in terms of load bearing capacity than the reinforcement by columns, notably when significantly inclined loading is concerned. Copyright © 2014 John Wiley & Sons, Ltd.

A four-node plane parametric element AQGβ6-I is constructed on the basis of the quadrilateral area coordinate, the generalized conforming principle and the projection technique with a penalty factor *β* within an interval of 0–1. When *β* = 0, the element has excellent bending performance. When *β* = 1, the element can pass patch test strictly; its performance is as good as many famous elements. When *β* value is between 0 and 1, such as *β* = 0.5, the element can arrive at a compromise between (relatively) low sensitivity to mesh distortion and perfect convergence. The work provides an illuminating method to alleviate a difficult problem in finite element modelling using the four-node quadrilateral element, which can pass the strict patch test, but has poor performance in bending dominated problem; on the contrary, it has excellent performance in bending dominated problem but cannot pass the strong patch test.

The AQGβ6-I with the convergence formulation (*β* = 1) is then applied to coupled solid-deformation/fluid-flow simulation for porous geomaterials. The computational examples are carried out to demonstrate that the AQGβ6-I (*β* = 1) element is not only stable, reliable and efficient but also of high accuracy. The present study provides a good applicable element for finite element simulations of solid-deformation/fluid-flow for porous geomaterials. Copyright © 2014 John Wiley & Sons, Ltd.

Large underground caverns are increasingly being considered for the construction of industrial facilities and transportation infrastructure in order to optimize the use of surface land in large urban cities. Due to the geological constraints underground, it is sometimes necessary to construct a cavern close to an existing cavern. Pillars serve as an underground support element in twin caverns, without which it is difficult to sustain the weight of the overburden materials. If the strength of a pillar is exceeded, it will fail, and the load that it carried will be transferred and thereby contribute to the collapse of the twin caverns. The lack of confinement in slender pillars also contributes to the complete collapse of pillars at relatively low stress magnitudes. From a design point of view, understanding the pillar failure mechanism and the interaction effect between twin caverns is essential. This paper presents a numerical investigation on the influence of various design parameters on twin cavern interaction. Pillar performances with respect to the changes to the maximum principal stress in the pillar, the peak vertical stress, and the peak principal stress difference in the pillar core are studied in order to examine the failure mechanisms and to identify situations in which there is significant cavern interaction and overlap of the plastic zones. Copyright © 2014 John Wiley & Sons, Ltd.

Conceived as a potential alternative to the classical design methods employed for analyzing the stability of underground structures driven in jointed rocks, the homogenization approach stems from the heuristic idea that, from a macroscopic point of view, a rock mass cut by a network of joints may be perceived as a homogenized continuum. The strength properties of the latter can be theoretically obtained from the failure conditions of its individual constituents: rock matrix and joint interfaces.

At the material level, the limit analysis reasoning is used in the context of homogenization to formulate the homogenized strength criterion of a jointed rock mass in the particular situation of a single set of parallel joints. As it could be expected, the obtained closed-form expressions show the strength anisotropy induced by joint preferential orientation. The support functions (*π* functions) associated with the homogenized strength criterion are also determined in both plane strain and three-dimensional cases. This criterion is then applied to the investigation of stability analysis of a tunnel excavated in a jointed rock mass. Upper bounds estimated of the stability factor are derived from the implementation of the kinematic approach directly on the homogenized underground structure. Finally, the approach is applied to analyze and discuss the collapse of the Pinheiros subway station (São Paulo, Brazil). Copyright © 2014 John Wiley & Sons, Ltd.

A transversely isotropic multi-layered half-space, with axis of material symmetry perpendicular to the free surface, supports a flexible either annular or solid circle foundation. The contact area of the foundation and the half-space is considered to be both frictionless and tensionless. The foundation is assumed to be affected by a vertical static axisymmetric load. Detailed analysis of the interaction of these two systems with different thickness of layers is the target of this paper. With the use of ring load Green's functions for both the foundation and the continuum half-space, an integral equation accompanied with some inequalities is introduced to model the complex BVP. With the incorporation of ring-shape FEM, we are capable of capturing both regular and singular solution smoothly. The validity of the combination of the analytical and numerical method is proved with comparing the results of this paper with a number of benchmark cases of both linear and nonlinear interaction of circular and annular foundation with half-space. Some new illustrations are presented to portray the aspect of the anisotropy and layering of the half-space. Copyright © 2014 John Wiley & Sons, Ltd.

This paper develops a multiple response surfaces approach to approximate the limit state function for slope failure by second-order polynomial functions, to incorporate the variation of the most probable slip surfaces, and to evaluate the slope failure probability *p _{f}*. The proposed methodology was illustrated through a cohesive soil slope example. It is shown that the

The variation of number of the most probable slip surfaces is studied at different scale of fluctuation (*λ*) values. It is found that when full correlation assumed for each of random fields (i.e., spatial variability is ignored), the number of the most probable slip surfaces is equal to the number of random fields (in this study, it is 3). When the spatial variability grows significantly, the number of the most probable slip surfaces or number of multiple response surfaces firstly increases evidently to a higher value and then varies slightly. In addition, the contribution of a specific most probable slip surface varies dramatically at different spatial variability level, and therefore, the variation of the most probable slip surfaces should be accounted for in the reliability analysis. The multiple response surfaces approach developed in this paper provides a limit equilibrium method and MCS-based means to incorporate such a variation of the most probable slip surfaces in slope reliability analysis. Copyright © 2014 John Wiley & Sons, Ltd.

Failure in geotechnical engineering is often related to tension-induced cracking in geomaterials. In this paper, a coupled meshless method and FEM is developed to analyze the problem of three-dimensional cracking. The radial point interpolation method (RPIM) is used to model cracks in the smeared crack framework with an isotropic damage model. The identification of the meshless region is based on the stress state computed by FEM, and the adaptive coupling of RPIM and FEM is achieved by a direct algorithm. Mesh-bias dependency, which poses difficulties in FEM-based cracking simulations, is circumvented by a crack tracking algorithm. The performance of our scheme is demonstrated by two numerical examples, that is, the four-point bending test on concrete beam and the surface cracks caused by tunnel excavation. Copyright © 2014 John Wiley & Sons, Ltd.

Damage in the form of cracks is predicted to assess the susceptibility of a tunnel to failure due to a blast. The material-point method is used in conjunction with a decohesive failure model as the basis for the numerical simulations. The assumption of a cylindrical charge as the source for the blast allows the restriction of plane strain and two-dimensional analyses. In the simulation, a further restriction of a single pressure pulse is used as the source of stress waves that are reflected and refracted after reaching the free surface of the tunnel wall. Three critical zones of significant cracking in the vicinity of a tunnel are identified as potential contributors to tunnel failure. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents the applications of the differential evolution (DE) algorithm in back analysis of soil parameters for deep excavation problems. A computer code, named Python-based DE, is developed and incorporated into the commercial finite element software ABAQUS, with a parallel computing technique to run an FE analysis for all trail vectors of one generation in DE in multiple cores of a cluster, which dramatically reduces the computational time. A synthetic case and a well-instrumented real case, that is, the Taipei National Enterprise Center (TNEC) project, are used to demonstrate the capability of the proposed back-analysis procedure. Results show that multiple soil parameters are well identified by back analysis using a DE optimization algorithm for highly nonlinear problems. For the synthetic excavation case, the back-analyzed parameters are basically identical to the input parameters that are used to generate synthetic response of wall deflection. For the TNEC case with a total of nine parameters to be back analyzed, the relative errors of wall deflection for the last three stages are 2.2, 1.1, and 1.0%, respectively. Robustness of the back-estimated parameters is further illustrated by a forward prediction. The wall deflection in the subsequent stages can be satisfactorily predicted using the back-analyzed soil parameters at early stages. Copyright © 2014 John Wiley & Sons, Ltd.

This paper presents a general coupling extended multiscale FEM (GCEMs) for solving the coupling problem of elasto-plastic consolidation of heterogeneous saturated porous media. In the GCEMs, the numerical multiscale base functions for the solid skeleton and fluid phase of the coupling system are all constructed on the basis of the equivalent stiffness matrix of the unit cell, which not only contain the interaction between the solid and fluid phases but also consider the time effect. Furthermore, in order to improve the computational accuracy for two-dimensional problems, a multi-node coarse element strategy for the GCEMs is proposed, and a two-scale iteration algorithm for the elasto-plastic consolidation analysis is developed. Some one-dimensional and two-dimensional homogeneous and heterogeneous numerical examples are carried out to validate the proposed method through the comparison with the coupling multiscale FEM and standard FEM. Numerical results show that the newly developed GCEMs can almost preserve the same convergent property as the standard FEM and also possesses the advantages of high computational efficiency. In addition, the GCEMs can be easily applied to other coupling multifield and multiphase transient problems. Copyright © 2014 John Wiley & Sons, Ltd.

The scaled boundary finite-element method, a semi-analytical computational scheme primarily developed for dynamic stiffness of unbounded domains, is applied to the analysis of unsteady seepage flow problems. This method is based on the finite-element technology and gains the advantages of the boundary element method as well. Only boundary of the domain is discretized, no fundamental solution is required and singularity problems can be modeled rigorously. Anisotropic and non-homogeneous materials satisfying similarity are modeled with no additional efforts. In this study, firstly, formulation of the method for the transient seepage flow problems is derived followed by its solution procedures. The accuracy, simplicity and applicability of the method are demonstrated via four numerical examples of transient seepage flow – three of them are available in the literature. Homogenous, non-homogenous, isotropic and anisotropic material properties are considered to show the versatility of the technique. Excellent agreement with the finite-element method is observed. The method out-performs the finite-element method in modeling singularity points. Copyright © 2014 John Wiley & Sons, Ltd.

While methane hydrates (MHs) can be present in various forms in deep seabeds or permafrost regions, this paper deals with MH-bearing sediments (MHBS) where the MH has formed bonds between sand grains. A bond model based on experimentally validated contact laws for cemented granules is introduced to describe the mechanical behavior of the MH bonds. The model parameters were derived from measured values of temperature, water pressure and MH density. Bond width and thickness adopted for each bond of the MHBS were selected based on the degree of MH saturation. The model was implemented into a 2D distinct element method code. A series of numerical biaxial standard compression tests were carried out for various degrees of MH saturation. A comparison with available experimental data shows that the model can effectively capture the essential features of the mechanical behavior of MHBS for a wide range of levels of hydrate saturation under drained and undrained conditions. In addition, the analyses presented here shed light on the following: (1) the relationship between level of cementation and debonding mechanisms taking place at the microscopic level and the observed macro-mechanical behavior of MHBS and (2) the relationship between spatial distribution of bond breakages and contact force chains with the observed strength, dilatancy and deformability of the samples. Copyright © 2014 John Wiley & Sons, Ltd.

By incorporating the nonlinear variation of a soil's compressibility and permeability during the process of consolidation, an analytical solution for the radial consolidation of vertical drains has been developed for a general time-variable loading. The general solution was verified for the cases of instantaneous loading and ramp loading. Detailed solutions were further derived for two special loading schemes: multistage loading and preloading–unloading–reloading. The nonlinear consolidation behavior of a vertical drain subjected to these two types of loading schemes was then investigated by a parametric study. The results show that the loading rate, the ratio of the compressibility index to the permeability index (*C*_{c}/*C*_{k}), and the initial stress state have a significant influence on the consolidation rate. A smaller value of *C*_{c}/*C*_{k}, a larger initial stress, or a fast loading rate always leads to a rapid consolidation rate. During the unloading period, a negative excess pore water pressure may occur, and a slower unloading rate may reduce this negative value. Copyright © 2014 John Wiley & Sons, Ltd.

The present study extends the numerical manifold method to include the hydro-mechanical model to investigate the effect of water flow in fractures on the stability of rock structures, particularly slopes. The proposed flow model is verified by a simple 2-D flow problem in a homogeneous aquifer. Combining the water flow model with the earlier developed fracture evolution technique, the entire failure process of the rock slope due to a heavy rain is simulated. The results illustrate that the developed numerical manifold method can not only determine the trigger factor of the crack initiation but also model the failure processes related to crack initiation, propagation, block formation, detachment and sliding due to the water effect successfully. Copyright © 2014 John Wiley & Sons, Ltd.

A series of micromechanical tests were conducted to investigate the bond failure criterion of bonded granules considering the effect of bond thickness, with the aim of enhancing the bond contact model used in the distinct element simulations of cemented geomaterials. The granules were idealized in a two-dimensional context as one pair of aluminum rods bonded by resin epoxy or cement. The mechanical responses of nearly 500 rod pairs were tested under different loading paths to attain the yield loads of bonded granules at variable bond thickness. This study leads to a generic bond failure criterion incorporating the effect of the bond thickness. The results show that the bond compressive resistance largely decreases with increasing bond thickness owing to the presence of the confinement at the bond-particle interface. The strength envelopes obtained from the combined shear compression tests and combined torsion compression tests have identical functional form, and they decrease in size with increasing bond thickness but remain unchanged in shape. Given the same cementation material, the generic bond strength envelope in a three-dimensional contact force space under different loading paths remains the same in shape but shrinks with the increase of bond thickness. Copyright © 2014 John Wiley & Sons, Ltd.

The paper presents closed-form solutions for stress and displacement influence functions for stress discontinuity (SD) and displacement discontinuity (DD) elements, for a two-dimensional plane-strain elastic, transversely anisotropic medium. The solutions for SD elements are based on Kelvin's problem and for DD elements on the concept of dipoles. Stress and displacement influence functions are derived for the following elements: constant SD, linear SD, constant DD, linear DD, square root DD, parabolic DD, constant DD surface, and linear DD surface elements. The formulations are incorporated into FROCK, a hybridized boundary element method code, and are validated by providing comparisons between the results from FROCK and the finite element code ABAQUS. A limited parametric analysis shows the effects of slight anisotropy on the stress field around the tip of a crack and of the orientation of the crack with respect to the axes of elastic symmetry. Copyright © 2014 John Wiley & Sons, Ltd.

A quasi-static homogeneous drained triaxial compression test on cohesionless sand under constant lateral pressure was simulated using a three-dimensional discrete element method. Grains were modelled by means of particle clusters composed of rigid spheres or spheres with contact moments imitating irregular particle shapes. Attention was paid to the effect of initial void ratio and grain shape mixture on the shear strength, volume changes, force chains, kinetic, elastic and dissipated energies. In addition, the effect of the mean grain size, grain size distribution, grain size range, specimen size and roughness and stiffness of boundaries was numerically analysed in initially dense sand. Some numerical results were compared with available experimental results. Copyright © 2014 John Wiley & Sons, Ltd.

Within the framework of soil–pile interaction, a novel displacement scheme for the transverse kinematic response of single piles to vertically propagating S waves is proposed on the basis of the modified Vlasov foundation model. The displacement model contains a displacement function along the pile axis and an attenuation function along the radial direction. The governing equations and boundary conditions of the two undetermined functions are obtained in a coupled form by using Hamilton's principle. An iterative algorithm is adopted to decouple and solve the two unknown functions. In light of the governing equation of the pile kinematics, a mechanical model is proposed to evaluate the present method on a physical basis considering material damping. The coefficient of the equivalent Winkler spring is derived explicitly as function of the displacement decay parameter *γ* and soil Poisson's ratio. A parametric study is performed to investigate the effects of the soil–pile system properties on the kinematic response of single piles. The results show that the dimensionless pile length controls the transverse kinematics of piles. In terms of the theory of beams on elastic foundation, the classification limits of the dimensionless pile length may be *π* ∕ 4 and *π*, respectively. Copyright © 2014 John Wiley & Sons, Ltd.

Recent study indicates that the response of rigid passive piles is dominated by elastic pile–soil interaction and may be estimated using theory for lateral piles. The difference lies in that passive piles normally are associated with a large scatter of the ratio of maximum bending moment over maximum shear force and induce a limiting pressure that is ~1/3 that on laterally loaded piles. This disparity prompts this study.

This paper proposes pressure-based pile–soil models and develops their associated solutions to capture response of rigid piles subjected to soil movement. The impact of soil movement was encapsulated into a power-law distributed loading over a sliding depth, and load transfer model was adopted to mimic the pile–soil interaction. The solutions are presented in explicit expressions and can be readily obtained. They are capable of capturing responses of model piles in a sliding soil owing to the impact of sliding depth and relative strength between sliding and stable layer on limiting force prior to ultimate state. In comparison with available solutions for ultimate state, this study reveals the 1/3 limiting pressure (of the active piles) on passive piles was induced by elastic interaction. The current models employing distributed pressure for moving soil are more pertinent to passive piles (rather than plastic soil flow). An example calculation against instrumented model piles is provided, which demonstrates the accuracy of the current solutions for design slope stabilising piles. Copyright © 2014 John Wiley & Sons, Ltd.

An analysis method for transient groundwater flow during slug tests performed in vertical cutoff walls is presented. The analytical solution for evaluating hydraulic conductivity of vertical cutoff walls is derived by applying the method of images to the previously developed analytical solution that is exclusively applicable to an infinite aquifer. Two distinct boundary conditions are considered to account for the configuration of the vertical cutoff wall: the wall-soil formation interfaces with or without the existence of filter cakes, that is, constant-head boundary and no-flux boundary conditions. A series of type curves is constructed from the analytical solution and compared with those of a partially penetrated well within an aquifer. The constant-head boundary condition provides faster hydraulic head recovery than the aquifer case. On the other hand, the no-flux boundary condition leads to a delayed hydraulic head recovery. The greater the shape factor and well offset from the center of the cutoff wall, and the smaller the width of the cutoff wall, the greater the effect of the boundary condition observed in the type curves. This result shows the significance of considering proper boundary conditions at the vertical cutoff wall in analyzing slug tests. Copyright © 2014 John Wiley & Sons, Ltd.

Mechanical trapping (or straining) of fine particles is a key mechanism in many filtration systems. For example, the performance of rapid sand filters depends in part on mechanical trapping of larger fine particles, while relying on adsorptive processes to trap very small fine particles and microbes. The ability to trap these particles is directly related to the construction of the packed bed used for filtration in this system. Thus, the ability to model the effect of the inner structure of the packed bed can lead to more efficient design for improved filtration. Because of its significant efficiency, gravitational sphere packing is employed in this work to simulate a bed of mono-sized randomly packed spheres. The simulated bed provides a way to visualize the pore network and estimate the pore size distribution associated with the void space between particles. Furthermore, by subsequently introducing fine particles into the bed, we evaluate the mass-rate of fine particles passing through and possibly saturating the packed bed. Results show that fine particles between 15% and 25% of the coarse particle size can be physically strained within the randomly packed bed. These results differ significantly from the results obtained assuming a periodically spaced bed. The technique therefore provides an efficient yet accurate alternative for understanding how fine particles pass through a coarse particulate medium. Copyright © 2014 John Wiley & Sons, Ltd.

In the present work, it is attempted to derive a macroscopic constitutive law for the elastic deformation of granular materials, based on microscopic considerations. For the sake of simplicity, the solution is restricted to two dimensions, that is, a random assembly of infinitely extended cylinders. After examining pairwise contact interactions, the elastic energy rate of the assembly is retrieved in a discrete form. Introducing the probability density function of the contact orientations, the continuum form of the elastic energy density rate is evaluated as a function of generalized strains and curvatures, and their rates. The stresses and couple stresses result as dual variables to the generalized strain and curvature rates. Some properties of the resulting model are discussed, examples are presented and conclusions are drawn.Copyright © 2014 John Wiley & Sons, Ltd.

In this paper, a series of multimaterial benchmark problems in saturated and partially saturated two-phase and three-phase deforming porous media are addressed. To solve the process of fluid flow in partially saturated porous media, a fully coupled three-phase formulation is developed on the basis of available experimental relations for updating saturation and permeabilities during the analysis. The well-known element free Galerkin mesh-free method is adopted. The partition of unity property of MLS shape functions allows for the field variables to be extrinsically enriched by appropriate functions that introduce existing discontinuities in the solution field. Enrichment of the main unknowns including solid displacement, water phase pressure, and gas phase pressure are accounted for, and a suitable enrichment strategy for different discontinuity types are discussed. In the case of weak discontinuity, the enrichment technique previously used by Krongauz and Belytschko [*Int. J. Numer. Meth. Engng.*, 1998; 41:1215–1233] is selected. As these functions possess discontinuity in their first derivatives, they can be used for modeling material interfaces, generating only minor oscillations in derivative fields (strain and pressure gradients for multiphase porous media), as opposed to unenriched and constrained mesh-free methods. Different problems of multimaterial poro-elasticity including fully saturated, partially saturated one, and two-phase flows under the assumption of fully coupled extended formulation of Biot are examined. As a further development, problems involved with both material interface and impermeable discontinuities, where no fluid exchange is permitted across the discontinuity, are considered and numerically discussed. Copyright © 2014 John Wiley & Sons, Ltd.

A series of laboratory experiments and numerical simulations are conducted to explore the characteristics of mixtures composed of sand and rubber particles of the same median diameter. The mixtures are prepared with different volumetric sand fractions (sf = V_{sand}/V_{total}). The experiment focuses on assessing the strain level on the characteristics of the mixture with the volume fraction of each component. Numerical simulations using the discrete element method are performed to obtain insight into the microscale behavior and internal mechanism of the mixtures. The experimental results show that the behavior of the mixtures is dependent on the relative sand and rubber particles composition with variation in the strain levels. The numerical simulation reveals the effect of the soft rubber particle inclusion in the mixture on the micromechanical parameters. In low sand fraction mixtures, a high shear stress along the contact is mobilized, and the stress state is driven to a more anisotropic condition because of the relatively high particle friction angle of the rubber. The rubber particles play different roles with the strain level in the mixture, including increasing the coordination number and controlling plasticity of the mixture in a small strain, preventing buckling of the force chain in an intermediate strain, and leading to contractive behavior in a large strain. Copyright © 2014 John Wiley & Sons, Ltd.

This paper deals with analytical and numerical modelling of the internal stress generated in argillaceous rocks during humidification/desiccation processes, which is an essential issue for damage study. This local stress field arises from two mechanisms: (i) complex interactions between free swelling/shrinking clay matrix and non-strained inclusions of carbonate and quartz and (ii) a self-restraint effect induced by the moisture gradient during the transient moisture exchange process. The inclusion–matrix interaction is investigated in different cases. Firstly, the analytical solution of the stress around a cylindrical inclusion embedded in an infinite swelling matrix is derived: The inclusion would suffer tension (compression) under humidification (desiccation), and the resulting cracking patterns are discussed. Then, the problem of two inclusions with different distances in an infinite swelling matrix is considered, and it is shown that the local stress around an inclusion will be perturbed and amplified by neighbouring inclusions. Finally, an inclusion outcropping at the free surface of a swelling matrix is modelled as to investigate the effect of free surface: The inclusion–matrix interface undergoes shear stresses of which the maximum is found at the free surface. In addition to the inclusion–matrix interaction, the self-restraint effect is investigated: The induced stress is maximal at the beginning of humidification/desiccation processes and vanishes gradually with time. The quantity of the self-restraint stress is strongly controlled by the hydric loading rate. Copyright © 2014 John Wiley & Sons, Ltd.

The paper examines ion (chloride) transport equations in porous media (concrete) integrated over a representative elementary volume, that is to say, averaging over the macroscopic level the phenomena that occur really at the pore scale. There are three basic variables to be used: chloride concentration, moisture and temperature. The diffusion process is examined, in addition to other phenomena such as convection (the motion of dissolved substances caused by flow of water in a pore solution of partially saturated media) or chloride binding (the capacity of free chloride of being chemically bound, particularly with C_{3}A to form Friedel salts). Contrary to other approaches, such effects are not considered by means of apparent diffusion coefficients but by developing the complete set of time-dependent equations for both the chloride concentration within the pore solution and the moisture content within the pore space.

Once the general model is described, the system of equations can be solved numerically by means of a two-dimensional finite element formulation. The main objective is to reproduce results of experimental tests by means of a priori parameter estimation, according to the characteristics of materials and external environment conditions, thereby superseding the well-known best fit a posteriori through Fick's second equation.

While the introduction of hygrometric conditions and convection phenomena appears to be of high significance, other factors like temperature, surface concentration, chloride binding or equivalent hydration time are analysed too. The proposed model can reproduce bidimensional complex geometries, for example, cracked concrete cover, as well as variable surface condition. An application case is developed through a realistic model of the geometry of a crack. Copyright © 2014 John Wiley & Sons, Ltd.

A numerical scheme for the computation of the permeability of complex microstructures is presented. As a darcean counterpart of the FFT-based scheme in elasticity, the method is designed to be directly coupled with 3D imaging techniques of porous samples, without meshing or definition of an equivalent pore network. The method relies on the variational principle of Hashin and Shtrikman, which ensures a rigorous upper bound status to the estimated permeabilities and provides an energetically consistent rule for heterogeneous voxels comprising both the solid and fluid phases. Copyright © 2014 John Wiley & Sons, Ltd.

An analytical solution is developed in this paper to investigate the dynamic response of a large-diameter end-bearing pipe pile subjected to torsional loading in viscoelastic saturated soil. The wave propagation in saturated soil and pile are simulated by Biot's two-phased linear theory and one-dimensional elastic theory, respectively. The dynamic equilibrium equations of the outer soil, inner soil, and pile are established. The solutions for the outer and inner soils in frequency domain are obtained by Laplace transform technique and the separation of variables method. Then, the dynamic response of the pile is obtained on the basis of the perfect contacts between the pile and the outer soil as well as the inner soil. The results in this paper are compared with that of a solid pile in elastic saturated soil to verify the validity of the solution. Furthermore, the solution in this paper is compared with the classic plane strain solution to verify the solution further and check the accuracy of the plane strain solution. Numerical results are presented to analyze the vibration characteristics and illustrate the effect of the soil parameters and the geometry size of the pile on the complex impedance and velocity admittance of the pile head. Finally, the displacement of the soil at different depth and frequency is analyzed. Copyright © 2014 John Wiley & Sons, Ltd.

In porous media, chemical species that dissolve in pore water can be transported via diffusion mechanisms or advective fluxes, close to or far away from where precipitation occurs. In the case of a high-level radioactive waste disposal system, compacted bentonite is used in a buffer material in an engineering barrier system to minimize the amount of specific nuclides that breach into the surrounding host rock. To minimize breaching, it is very important to understand the transport mechanism of multiple chemical species in porous media. In the following research, we introduced FEM analysis methods using the results of the molecular dynamics simulation and homogenization analysis (MD/HA) method. First, the diffusion coefficients of ions (Cl^{−}, I^{−}, and Na^{+}) in different water layers of Na-beidellite were calculated using the MD/HA procedure under various dry density (1.2, 1.6, and 2.0 Mg/m^{3}) and temperature (293, 323, and 363 K) conditions. Next, using FEM analysis that used the MD/HA results as input parameters, the diffusion behaviors of ions in porous media were calculated. The results indicate that the diffusion coefficients of the interlayer water in Na-beidellite are different from the diffusion coefficients under dry density conditions. Further, the concentration profiles (*C _{t}/C*