Computational Geomechanics Research Articles

Constitutive modelling in computational geomechanics – 61st Rankine lecture, British Geotechnical Association, 2023

Constitutive models are an essential part of computational modelling in geotechnics; they are at the heart of almost all theoretical predictions of geotechnical structures. How the stress–strain (and perhaps time) response of soil (and rock) is represented in these mathematical models is usually the key to successful prediction of the behaviour of geotechnical structures. However, the important details of these models, particularly the idealisations that are made, may be poorly or incompletely understood, or ignored, sometimes at significant cost to the unwary analyst. Indeed, the capabilities and the shortcomings of these models, especially the more advanced models, are not always easy to ascertain. In some cases, determination of the values of the input parameters is not straightforward. Consequently, it may be difficult to determine which model to select for a particular task. This paper explores some of the more important developments in the constitutive modelling of soils and addresses some of these issues of potential concern. The need for such models and the various attributes and capabilities that the commonly used models possess are reviewed. Also discussed is the issue of matching a particular model to the geotechnical problem at hand, which model attributes are required and why. The intention is to place emphasis on the physical basis of these models, rather than explore their mathematical complexity in detail. Some of the constitutive models encoded in the software packages used routinely in geotechnical practice are reviewed, and discussion is also provided on their specific limitations. Examples of practical applications, involving the solution of boundary and initial value problems, are described to illustrate both the advantages and some of the limitations of both commonly used and highly advanced constitutive models.

Strength and stiffness parameters for hardening soil model of rockfill materials

In recent years, an increasing number of rockfill dams have been constructed in many developing countries around the world. One difficulty in the engineering design and construction of rockfill dams is the characterisation of the rockfill materials. Thanks to advancements in computational geomechanics, a finite element analysis is usually employed in geotechnical engineering design. To perform such analyses, a proper constitutive model with relevant parameters that reflect the stress–strain-strength characteristics of rockfill materials are required. This study presents a procedure to determine the strength and stiffness parameters of rockfill materials that can be used with the Hardening Soil (HS) model to characterise the behaviour of rockfill materials. Based on triaxial and oedometer test results, a set of HS model parameters can be determined. To implement and validate the procedure, six rockfill materials from the Nam Ngum 2 Dam construction project in Laos PDR are employed in this study. The influences of stress level and particle breakage on the stress–strain/volumetric-strain behaviours of rockfill materials are discussed. Finally, the procedure to determine the material parameters for rockfill materials presented in this study can also be applied to other rockfill materials.

Physics-based reduced order model for computational geomechanics

Numerical models are an essential tool in stability analysis, design, and construction for geotechnical engineering. Yet, numerical modeling is too time-consuming in practical engineering. In this study, a physics-based reduced order model (ROM) was developed to approximate the displacement and stress field of geotechnical structure by combining Latin hypercube sampling (LHS), a numerical method and proper orthogonal decomposition (POD). The set of design variables were constructed using LHS. A numerical model was used to generate the snapshots based on the design variables. POD was used to compute the eigenvalues and eigenvectors of a spatial Gram matrix, which was constructed based on snapshots. The first K eigenfunction vectors were determined based on eigenvalues and eigenvectors of the spatial Gram matrix. The interpolation matrix of elements was computed using a radial basin function (RBF), and then the vector of an element was determined by solving the penalized linear systems. To determine the new design variables, the coefficient of ROM was determined based on the RBF and the vector of elements, and the unknown field variables were predicted based on the ROM. The ROM was illustrated and verified for a circular tunnel. The results showed that the predicted displacement and stress field were in excellent agreement with both the analytical and numerical solutions. The physics-based ROM characterized well the deformation and failure mechanism of the surrounding rock mass and can be used to replace a numerical model for back analysis, optimal design, and uncertainty analysis of geotechnical engineering, thereby eliminating costly repetitive computations.

Mitigation of ancient coal mining hazards to overhead line equipment structures

Subsidence and collapse of structures supporting overhead line equipment (OLE) due to legacy coal mining can be catastrophic during both the construction and operational phases, causing damage to construction plant, injuries and fatalities to the construction and maintenance workforce, failure of the track and disruption to the operational railway. When applied to linear railway schemes spanning considerable distances, industry guidelines for appraising and mitigating these mining hazards are overly conservative with a disproportionate cost. On phase 4 of the North West Electrification Project between Manchester Victoria and Euxton Junction, an innovative risk-based method was devised and implemented to identify and mitigate mining hazards to OLE structures. Formulation of these value-engineering solutions prevented programme drift and led to significant cost savings. Details of this methodology, including a detailed desk study, computational geomechanics, a well-targeted ground investigation, void remediation and real-time foundation monitoring, are discussed in this paper.

Application of the smoothed point interpolation methods in computational geomechanics: A comparative study

In this study, computational efficiency and performance of the smoothed point interpolation methods (SPIMs) in coupled problems of geomechanics are examined and compared with those of the standard linear finite element method (FEM). In the SPIMs, the problem domain is first discretised using a triangular background mesh, and the smoothing domains are then constructed using the cells of the background mesh. Different approaches for construction of the smoothing domains result in different SPIMs including Edge-based SPIM (ESPIM), cell-based SPIM (CSPIM), and node-based SPIM (NSPIM), each featuring different computational characteristics. For each type of the SPIM, only the simplest and most widely used supporting node selection technique is adopted in this study. The field function is approximated using shape functions constructed adopting either the polynomial point interpolation method (PIM), or radial point interpolation method (RPIM). The governing equations are derived from the equilibrium equation of the solid phase, and momentum and mass balance equations for the flow phase. The equations are then discretised in space using the generalized smoothed Galerkin (GS-Galerkin) method, and in time using the three-point time discretisation technique. The performance of the SPIMs in coupled geomechanical problems are then thoroughly assessed and compared to those of the linear FEM through extensive numerical analyses of four benchmark problems in saturated porous media covering both linear and non-linear material behaviours, and the conclusions drawn by analysing the results.

Code verification in computational geomechanics: Method of manufactured solutions (MMS)

Code verification is the process of ensuring that there are no coding mistakes or algorithm deficiencies in a scientific computing code. In this paper, code verification in computational geomechanics based on the method of manufactured solutions (MMS) is discussed. The application of this technique, which is the most rigorous method among all code verification techniques, is presented through verification of an in-house geomechanics code based on the Edge-based Smoothed Point Interpolation Method (ESPIM) for spatial, and three-point scheme for temporal discretisation of the governing equations. The MMS, as a generic approach to construct all-inclusive and complex test solutions for code verification, is explained. Two manufactured solutions along with an order of accuracy study are employed in several different setups to obtain the numerical and formal order of accuracies of the code. Results from the mesh and time step refinement studies confirm the verification of our code, and demonstrate the application of the MMS and order of convergence study in coupled problems in geomechanics.

Finite Element Modelling of Coupled Fluid-Flow and Geomechanical Aspects for the Sustainable Exploitation of Reservoirs: The Case Study of the Cavone Reservoir

The study of the seismogenic mechanical effects induced by oil & gas activities is a socially impacting issue of environmental engineering as well as a challenging task in computational geomechanics. It requires the solution of a coupled problem governed by poroelastic and fluid flow equations in a faulted domain in the presence of in situ stress fields. As a viable alternative to state-of-the-art academical computational models, the present study contributes a simplified methodology based on a commercial Finite Element multiphysics software. The focus is on the evaluation of the link between the oil & gas activities of the Cavone oilfield reservoir, located in North Italy and adjacent to the Mirandola fault, and the recent seismic sequence that struck Emilia in May 2012. An operational coupled fluid-geomechanical procedure is developed where the Cavone reservoir is subjected to the typical in situ stresses, and the nearby Mirandola fault is modelled as an impervious thin layer.

Hybrid Subsea CCS System-CO2 Storage in Offshore Ultra Deep-Water Salt Caverns

Computational geomechanics is the field of computational mechanics where geo-engineering systems are analyzed using the prinCurrently there is a demand for CCS of large quantities of CO₂ associated with CH₄ in the pre-salt offshore oil fields in Brazil. The pre-salt reservoirs have as caprock 2000 meters of continuous rock salt. Rock salt is a special geomaterial.

Current Trends in Computational Geomechanics and Geotechnical Engineering

Current Trends in Computational Geomechanics and Geotechnical Engineering

Application of mesh-free smoothed particle hydrodynamics (SPH) for study of soil behavior

The finite element method (FEM), discrete element method (DEM), and Discontinuous deformation analysis (DDA) are among the standard numerical techniques applied in computational geo-mechanics. However, in some cases there no possibility for modelling by traditional finite analytical techniques or other mesh-based techniques. The solution presented in the current study as a completely Lagrangian and mesh-free technique is smoothed particle hydrodynamics (SPH). This method was basically applied for simulation of fluid flow by dividing the fluid into several particles. However, several researchers attempted to simulate soil-water interaction, landslides, and failure of soil by SPH method. In fact, this method is able to deal with behavior and interaction of different states of materials (liquid and solid) and multiphase soil models and their large deformations. Soil indicates different behaviors when interacting with water, structure, instrumentations, or different layers. Thus, study into these interactions using the mesh based grids has been facilitated by mesh-less SPH technique in this work. It has been revealed that the fast development, computational sophistication, and emerge of mesh-less particle modeling techniques offer solutions for problems which are not modeled by the traditional mesh-based techniques. Also it has been found that the smoothed particle hydrodynamic provides advanced techniques for simulation of soil materials as compared to the current traditional numerical methods. Besides, findings indicate that the advantages of applying this method are its high power, simplicity of concept, relative simplicity in combination of modern physics, and particularly its potential in study of large deformations and failures.

Special Issue on Experimental and Computational Geomechanics for Unsaturated Soils

Special Issue on Experimental and Computational Geomechanics for Unsaturated Soils

On Application of the Maximum Entropy Meshless Method for Large Deformation Analysis of Geotechnical Problems

Traditional grid-based numerical techniques such as the Finite Element Method (FEM) are known to suffer when large deformations of the continuum are encountered. As such, there has been limited success using this class of methods to solve many of the complex problems encountered in computational geomechanics. The potential of Meshfree techniques for addressing this perceived deficiency has been recognised. This study presents a robust Maximum Entropy Meshless (MEM) method for the analysis of problems involving geometrical nonlinearity in computational geomechanics. The method is validated via simulation of an undrained layer of soil under a rigid and rough strip footing undergoing large deformations and its merit is demonstrated through a comparison of the results with those obtained via the FEM.

Pulse Fracturing in Shale Reservoirs: Geomechanical Aspects, Ductile/Brittle Transition, and Field Implications

Summary This paper presents a brief review of misconceptions on industry-standard brittleness/ductility definitions; geomechanical aspects and numerical evaluation of pulse fracturing by use of an advanced constitutive model implemented within the ANSYS® Autodyn® (ANSYS 2014); and fracture-network patterns because of pulse loading in shale that show ductile/brittle transition. In shale gas, one of the primary goals is to create extensive fracture networks that can remain open during production. Field experience has shown that not all shale formations respond to hydraulic fracturing effectively. It is important to identify and accurately design alternative fracturing techniques that would overcome some of the limitations. Pulse-fracturing rates and peak loads can be customized to lie between hydraulic and explosive fracturing. This technique has the potential to shatter shale, in particular by triggering a ductile/brittle transition at an optimized pulse rate. To date, operational considerations of pulse fracturing for success remain qualitative. Recent advances in computational geomechanics help us quantify the effect of key operational parameters for field applications. Further to this, advanced constitutive models implemented for these analyses have the benefit of simulating ductile/brittle transition, if the stress state and loading conditions dictate that the material should. This study on pulse fracturing shows that for a certain combination of reservoir, geomechanical, and pulse-loading parameters, induced fractures can propagate in multiple directions. This phenomenon might promote a self-propping mechanism for a network of fractures. At the end of this paper, favorable conditions when the pulse-fracturing technique would work and key parameters that trigger ductile/brittle transition are summarized and presented.

Application of a GPU-accelerated hybrid preconditioned conjugate gradient approach for large 3D problems in computational geomechanics

This paper presents a hybrid preconditioning technique for Conjugate Gradient method and discusses its parallel implementation on Graphic Processing Unit (GPU) for solving large sparse linear systems arising from application of interior point methods to conic optimization problems in the context of nonlinear Finite Element Limit Analysis (FELA) for computational Geomechanics. For large 3D problems, the use of direct solvers in general becomes prohibitively expensive due to exponentially growing memory requirements and computational time. Besides, the so-called saddle-point systems resulting from use of optimization framework is not an exemption. On the other hand, although preconditioned iterative methods have moderate storage requirements and therefore can be applied to much larger problems than direct methods, they usually exhibit high number of iterations to reach convergence. In present paper, we show that this problem can be effectively tackled using the proposed hybrid preconditioner along with an elaborate implementation on GPU. Furthermore, numerical results verify the robustness and efficiency of the proposed technique.

An ALE method for penetration into sand utilizing optimization-based mesh motion

The numerical simulation of penetration into sand is one of the most challenging problems in computational geomechanics. The paper presents an arbitrary Lagrangian–Eulerian (ALE) finite element method for plane and axisymmetric quasi-static penetration into sand which overcomes the problems associated with the classical approaches. An operator-split is applied which breaks up solution of the governing equations over a time step into a Lagrangian step, a mesh motion step, and a transport step. A unique feature of the ALE method is an advanced hypoplastic rate constitutive equation to realistically predict stress and density changes within the material even at large deformations. In addition, an efficient optimization-based algorithm has been implemented to smooth out the non-convexly distorted mesh regions that occur below a penetrator. Applications to shallow penetration and pile penetration are given which make use of the developments.

An adaptive preconditioned conjugate gradient approach for extremely large optimisation problems arising in computational geomechanics

The feasibility of preconditioned conjugate gradient (PCG) for solving extremely large sparse linear systems arising in three-dimensional finite element limit analysis of geotechnical stability problems is studied. With problem size increasing beyond few millions of variables, the use of direct solvers becomes prohibitive (both in storage and computational times) even on modern computers. This motivates the study on possible employment of iterative methods to solve systems of equations resulting from interior point method-based limit analysis formulations. In this paper, an adaptive preconditioning technique based on incomplete Cholesky (IC) preconditioner, as well as a novel implementation scheme for preconditioned conjugate gradient method, is proposed to effectively deal with highly ill-conditioned systems, such as those which occur at last iterations of IPM methods. Numerical tests on some ultra-large scale sample geotechnical problems, along with comparison with other potential solution methods, verify the robustness and efficiency of the proposed approach.

Parallel preconditioned conjugate gradient method for large sparse and highly ill-conditioned systems arising in computational geomechanics

The efficiency of parallel preconditioned conjugate gradient PCG algorithm for solving large sparse linear systems arising from application of interior point methods to conic optimisation problems in the context of nonlinear finite element limit analysis FELA for computational geomechanics is studied. For large 3D problems, the use of direct solvers in general becomes prohibitively expensive owing to exponentially growing memory requirements and computational time. And the so-called saddle-point systems resulting from use of optimisation framework is not an exemption. On the other hand, although preconditioned iterative methods have moderate storage requirements and therefore can be applied to much larger problems than direct methods, they usually exhibit high number of iterations to reach convergence. In the present paper, we show that this problem can be effectively tackled using efficient variants of sparse approximate inverse preconditioners along with an elaborate parallel implementation on multicore CPUs and significant improvements can be achieved by parallel implementation on graphic processing unit GPU. Furthermore, the efficiency of our proposed implementation is verified by the presented numerical results.

Wave-induced dynamic response of saturated multi-layer porous media: Analytical solutions and validity regions of various formulations in non-dimensional parametric space

In this paper, dynamic response of saturated-layered porous media under harmonic waves is evaluated through a semi-analytical solution. The coupled differential equations governing the dynamics of saturated or nearly saturated porous media such as soils containing all the inertial terms of solid and fluid phases are presented for a multi-layer system. Possible simplifications of the equations which are called formulations are introduced based upon the presence of inertial terms associated with the phases. The semi-analytical solutions to the response of multiple layers for all the formulations are presented in terms of pore water pressure and stress variations considering a set of non-dimensional parameters and their respective ratios. Validity of the formulations is presented in a non-dimensional parametric space. The maximum discrepancies in the pore pressure response of the formulations leading to validity regions are illustrated for typical dynamic problems. Subsequently, the effects of layering and drainage conditions on these regions are also presented. The proposed semi-analytical solution may be served as a benchmark one for validating the coupled numerical solutions, which can be used to deal with real scientific and geo-engineering problems in the emerging field of computational geomechanics.

Large deformation analysis of geomechanics problems by a combined rh-adaptive finite element method

Finite element analysis of large deformation problems is a major challenge in computational geomechanics, due largely to the severe mesh distortion that may occur after updating the spatial configuration of the nodal points using a conventional Updated-Lagrangian approach. There are two alternative and reasonably well-known strategies to tackle this issue of mesh distortion, viz., the r-adaptive and h-adaptive methods. The r-adaptive finite element method has been designed to eliminate possible mesh distortion by changing and optimising the locations of the nodal points without modifying the overall topology of the mesh adopted to solve a given problem. In order to obtain an accurate solution by this method a relatively fine mesh is required right from the beginning of the analysis, which potentially increases the overall analysis time. On the other hand, the h-adaptive finite element method improves the accuracy of the solution by gradually decreasing the size of the elements based on an error assessment method. However, this approach may leave the very small elements in the mesh vulnerable to distortion. To eliminate the individual drawbacks of these adaptive methods, while preserving the accuracy of the solution, a combined rh-adaptive finite element method has been developed and is presented in this paper for the analysis of sophisticated problems of geomechanics that involve large deformations and changing boundary conditions. The proposed method is designed to improve the accuracy of the solution obtained using the h-refinement strategy while successfully avoiding the mesh distortion by the use of r-adaptive refinement. It is shown that this new combination can significantly increase the efficiency of adaptive finite element methods.

A comparative study of error assessment techniques for dynamic contact problems of geomechanics

Numerical analysis of dynamic large deformation problems is one of the most challenging and sophisticated tasks in computational geomechanics. The complexity is mainly due to nonlinear soil behaviour, large deformations accompanied by severe mesh distortion, changing boundary conditions due to contact, and time-dependent behaviour. A typical example of such problems is the dynamic penetration of an object into a layer of soil due to its initial kinetic energy. This paper introduces an h-adaptive finite element method to tackle penetration as well as indentation problems of geomechanics in the presence of inertia forces. The h-adaptive finite element procedure automatically changes and optimises the density of the finite element mesh in a region to obtain a more accurate solution as well as to eliminate or reduce possible mesh distortion. A key component of an h-adaptive strategy is the error assessment technique. Although several methods exist for estimating the error in a finite element domain, the advantages as well as the capability of such methods in many geotechnical problems, particularly those involved with dynamic forces and changing boundary conditions, are not clearly understood. This paper describes a comparative study between three alternative error estimation techniques, including those based on the energy norm, the Green–Lagrange strain, and the plastic dissipation. The specific problems investigated involve the static and dynamic analysis of a strip footing on an undrained layer of soil, as well as the penetration of an object into a layer of sand. For these dynamic contact problems of geomechanics, the numerical results clearly show that the error estimator based on the Green–Lagrange strain outperforms the other two methods.