Hans Muhlhaus

A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials

We review the methods available for large deformation simulations of geomaterials before presenting a Lagrangian integration point finite element method designed specifically to tackle this problem. In our Ellipsis code, the problem domain is represented by an Eulerian mesh and an embedded set of Lagrangian integration points or particles. Unknown variables are computed at the mesh nodes and the Lagrangian particles carry history variables during the deformation process. This method is ideally suited to model fluid-like behavior of continuum solids which are frequently encountered in geological contexts. We present benchmark examples taken from the geomechanics area.

Longevity and stability of cratonic lithosphere: Insights from numerical simulations of coupled mantle convection and continental tectonics

[1] The physical conditions required to provide for the tectonic stability of cratonic crust and for the relative longevity of deep cratonic lithosphere within a dynamic, convecting mantle are explored through a suite of numerical simulations. The simulations allow chemically distinct continents to reside within the upper thermal boundary layer of a thermally convecting mantle layer. A rheologic formulation, which models both brittle and ductile behavior, is incorporated to allow for plate-like behavior and the associated subduction of oceanic lithosphere. Several mechanisms that may stabilize cratons are considered. The two most often invoked mechanisms, chemical buoyancy and/or high viscosity of cratonic root material, are found to be relatively ineffective if cratons come into contact with subduction zones. High root viscosity can provide for stability and longevity but only within a thick root limit in which the thickness of chemically distinct, high-viscosity cratonic lithosphere exceeds the thickness of old oceanic lithosphere by at least a factor of 2. This end-member implies a very thick mechanical lithosphere for cratons. A high brittle yield stress for cratonic lithosphere as a whole, relative to oceanic lithosphere, is found to be an effective and robust means for providing stability and lithospheric longevity. This mode does not require exceedingly deep strength within cratons. A high yield stress for only the crustal or mantle component of the cratonic lithosphere is found to be less effective as detachment zones can then form at the crust-mantle interface which decreases the longevity potential of cratonic roots. The degree of yield stress variations between cratonic and oceanic lithosphere required for stability and longevity can be decreased if cratons are bordered by continental lithosphere that has a relatively low yield stress, i.e., mobile belts. Simulations that combine all the mechanisms can lead to crustal stability and deep root longevity for model cratons over several mantle overturn times, but the dominant stabilizing factor remains a relatively high brittle yield stress for cratonic lithosphere.

Instability, softening and localization of deformation

Abstract Material strain softening is commonly taken as a necessary and sufficient condition for localization in deforming rocks. However, there is a wide range of experimental and theoretical information which shows that localization can occur in sands, brittle rocks and ductile metals under strain-hardening conditions. This paper aims to bring these two contrasted views together. Three separate criteria are necessary in order to understand localization behaviour. The first involves the stability of the deforming system. The second determines whether a deforming system will undergo bifurcation so as to cease deforming in a homogeneous mode and instead deform in an inhomogeneous mode such as barrelling or localization. The stability and bifurcation criteria are independent of each other since barrelling is a stable mode whereas localization is unstable. The third criterion establishes if the unstable bifurcation mode is one of localization or of some other kind. Localization may arise from the presence of vertices on the yield surface (as in the case of pressure insensitive, rate dependent metals and in brittle rocks due to the development of preferred microfractures for slip) or from the constitutive relation being such that the plastic strain-rate vector is not normal to the yield surface (as in the cases of pressure sensitive, dilatant rocks, of materials deforming by crystal-plastic processes involving dislocation cross-slip and/or climb, and of visco-plastic materials in which voids are forming due to diffusive processes). It is important to distinguish between material and system softening (or hardening) behaviour. The theory for a kinematically unconstrained shortening experiment (that is, rigid, frictionless platens) indicates that localization can occur in strain-hardening materials but the system must strain-soften from then on; that is, localization occurs at peak stress for the system even though the material may continue to harden (or soften). However, the addition of kinematic constraints (such as friction at elastic platens, a constraint to deform in plane strain or at constant volume) means that localization may occur in a system that is monotonically strain hardening. Shear zones in naturally deformed rocks show ample evidence of dilatant behaviour in that evidence for the passage of large volumes of fluid during localization is common as is the development of dilatant vein systems. As such, since shear zones are strongly constrained by the elastic and (limited) plastic response of the relatively undeformed rocks surrounding the shear zones, strain-hardening behaviour of the system is to be expected as the norm, even if the rocks within the shear zones are undergoing material strain-softening.

Finite element modelling of temperature gradient driven rock alteration and mineralization in porous rock masses

We present finite element simulations of temperature gradient driven rock alteration and mineralization in fluid saturated porous rock masses. In particular, we explore the significance of production/annihilation terms in the mass balance equations and the dependence of the spatial patterns of rock alteration upon the ratio of the roll over time of large scale convection cells to the relaxation time of the chemical reactions. Special concepts such as the gradient reaction criterion or rock alteration index (RAI) are discussed in light of the present, more general theory. In order to validate the finite element simulation, we derive an analytical solution for the rock alteration index of a benchmark problem on a two-dimensional rectangular domain. Since the geometry and boundary conditions of the benchmark problem can be easily and exactly modelled, the analytical solution is also useful for validating other numerical methods, such as the finite difference method and the boundary element method, when they are used to deal with this kind of problem. Finally, the potential of the theory is illustrated by means of finite element studies related to coupled flow problems in materially homogeneous and inhomogeneous porous rock masses.

8 – Continuum Models for Layered and Blocky Rock

Publisher Summary This chapter presents the framework of Cosserat continuum theory (CT) to model the behavior of granular, layered, and blocky rock. The materials of geomechanics are distinguished from other engineering materials such as steel and concrete, primarily by their visible inhomogeneity. The mechanical modeling of such materials is complex but is needed for the numerical and analytical prediction or back-analysis of forces and displacements within rock bodies around engineering structures. Nonhomogeneous materials can be modeled by conventional continuum theories only if the characteristic fabric length of the material is vanishingly small, as compared to some characteristic structural length. The chapter discusses the most important kinematic and static relationships for plane, infinitesimal deformations of Cosserats continuum. The absence of any characteristic length for the localized strain in conventional continuum models is the reason for the critical dependence of finite element calculations in related boundary value problems on the employed finite element grid. Analytical solutions of boundary value problems of Cosserat theory require approximately the same mathematical effort as corresponding boundary value problems of conventional continuum mechanics.

Generalised homogenisation procedures for granular materials

Engineering materials are generally non-homogeneous, yet standard continuum descriptions of such materials are admissible, provided that the size of the non-homogeneities is much smaller than the characteristic length of the deformation pattern. If this is not the case, either the individual non-homogeneities have to be described explicitly or the range of applicability of the continuum concept is extended by including additional variables or degrees of freedom. In the paper the discrete nature of granular materials is modelled in the simplest possible way by means of finite-difference equations. The difference equations may be homogenised in two ways: the simplest approach is to replace the finite differences by the corresponding Taylor expansions. This leads to a Cosserat continuum theory. A more sophisticated strategy is to homogenise the equations by means of a discrete Fourier transformation. The result is a Kunin-type non-local theory. In the following these theories are analysed by considering a model consisting of independent periodic 1D chains of solid spheres connected by shear translational and rotational springs. It is found that the Cosserat theory offers a healthy balance between accuracy and simplicity. Kunin’s integral homogenisation theory leads to a non-local Cosserat continuum description that yields an exact solution, but does not offer any real simplification in the solution of the model equations as compared to the original discrete system. The rotational degree of freedom affects the phenomenology of wave propagation considerably. When the rotation is suppressed, only one type of wave, viz. a shear wave, exists. When the restriction on particle rotation is relaxed, the velocity of this wave decreases and another, high velocity wave arises.

Finite element analysis of steady-state natural convection problems in fluid-saturated porous media heated from below

In this paper, a progressive asymptotic approach procedure is presented for solving the steady-state Horton-Rogers-Lapwood problem in a fluid-saturated porous medium. The Horton-Rogers-Lapwood problem possesses a bifurcation and, therefore, makes the direct use of conventional finite element methods difficult. Even if the Rayleigh number is high enough to drive the occurrence of natural convection in a fluid-saturated porous medium, the conventional methods will often produce a trivial non-convective solution. This difficulty can be overcome using the progressive asymptotic approach procedure associated with the finite element method. The method considers a series of modified Horton-Rogers-Lapwood problems in which gravity is assumed to tilt a small angle away from vertical. The main idea behind the progressive asymptotic approach procedure is that through solving a sequence of such modified problems with decreasing tilt, an accurate non-zero velocity solution to the Horton-Rogers-Lapwood problem can be obtained. This solution provides a very good initial prediction for the solution to the original Horton-Rogers-Lapwood problem so that the non-zero velocity solution can be successfully obtained when the tilted angle is set to zero. Comparison of numerical solutions with analytical ones to a benchmark problem of any rectangular geometry has demonstrated the usefulness of the present progressive asymptotic approach procedure. Finally, the procedure has been used to investigate the effect of basin shapes on natural convection of pore-fluid in a porous medium

The role of mobile belts for the longevity of deep cratonic lithosphere: The Crumple Zone Model

Many Archean cratons are surrounded by Proterozoic mobile belts that have experienced episodes of tectonic re-activation over their lifetimes. This suggests that mobile belt lithosphere may be associated with long lived, inherited weakness. It is proposed that the proximity of this weakness can increase the longevity of deep Archean lithosphere by buffering Archean cratons from mantle derived stresses. The physical plausibility of this idea is explored through numerical simulations of mantle convection that include continents and allow for material rheologies that model the combined brittle and ductile behavior of the lithosphere. Within the simulations, the longevity of deep cratonic lithosphere does increase if it is buffered by mobile belts that can fail at relatively low stress levels.

Computer simulation of single-layer buckling☆

Abstract This paper aims to clarify the influence that initial perturbations have in controlling the buckling process. The questions are: how is the Biot—Ramberg dominant wavelength modified by the presence of finite initial perturbations? How is the shape of the resultant folds influenced by the initial geometry? In answering these questions we also revisit many of the results already embedded in the literature for viscous materials but contrast the behaviour of these materials with those of strongly pressure-dependent elastic—plastic materials; this paper represents the first time that the buckling behaviour of such materials has been reported. For layers with a series of initial small perturbations, the current results confirm that fold wavelength and growth rate are controlled by competence contrast ( R ). Wavelength selection also occurs in a layer involving perfectly-sinusoidal small perturbations, resulting in a dominant wavelength different from the input one. For layers with an isolated initial perturbation, both R and initial perturbation geometries influence buckling. If the width of the initial perturbation is smaller than a critical width related to R , significant growth of the initial perturbation is possible. When the width of the initial perturbation is larger than the critical width, the simple growth of the perturbation is possible only for early stages. It then splits into two or more secondary perturbations according to R . These perturbations can all grow into finite folds in elastic—viscous models, but only some of them can do so in elastic—plastic models.

Theoretical and numerical analyses of convective instability in porous media with upward throughflow

Exact analytical solutions have been obtained for a hydrothermal system consisting of a horizontal porous layer with upward throughflow. The boundary conditions considered are constant temperature, constant pressure at the top, and constant vertical temperature gradient, constant Darcy velocity at the bottom of the layer. After deriving the exact analytical solutions, we examine the stability of the solutions using linear stability theory and the Galerkin method. It has been found that the exact solutions for such a hydrothermal system become unstable when the Rayleigh number of the system is equal to or greater than the corresponding critical Rayleigh number. For small and moderate Peclet numbers (Pe less than or equal to 6), an increase in upward throughflow destabilizes the convective flow in the horizontal layer. To confirm these findings, the finite element method with the progressive asymptotic approach procedure is used to compute the convective cells in such a hydrothermal system. Copyright (C) 1999 John Wiley & Sons, Ltd.