Modelling large-deforming fluid-saturated porous media using an Eulerian incremental formulation

Modelling large-deforming fluid-saturated porous media using an Eulerian incremental formulation

Eduard Rohan Vladimír Lukeš European Centre of Excellence, NTIS – New Technologies for Information Society Faculty of Applied Sciences, University of West Bohemia, Univerzitnií 22, 30614 Pilsen, Czech Republic

The paper deals with modelling fluid saturated porous media subject to large deformation. An Eulerian incremental formulation is derived using the problem imposed in the spatial configuration in terms of the equilibrium equation and the mass conservation. Perturbation of the hyperelastic porous medium is described by the Biot model which involves poroelastic coefficients and the permeability governing the Darcy flow. Using the material derivative with respect to a convection velocity field we obtain the rate formulation which allows for linearization of the residuum function. For a given time discretization with backward finite difference approximation of the time derivatives, two incremental problems are obtained which constitute the predictor and corrector steps of the implicit time-integration scheme. Conforming mixed finite element approximation in space is used. Validation of the numerical model implemented in the SfePy code is reported for an isotropic medium with a hyperelastic solid phase. The proposed linearization scheme is motivated by the two-scale homogenization which will provide the local material poroelastic coefficients involved in the incremental formulation.

large deformation, fluid saturated hyperelastic porous media, updated Lagrangian formulation, Biot model, finite element method

1 Introduction

Within the small strain theory, behaviour of the fluid-saturated porous materials (FSPM) was described by M.A. Biot in Biot (1941) and further developed in Biot (1955). It has been proved (see e.g. Auriault et al. (1990); Auriault and Boutin (1992) cf. Rohan et al. (2015)) that the poroelastic coefficients obtained by M.A. Biot and D.G. Willis in Biot and Willis (1957) hold also for quasistatic situations, cf. Brown et al. (2014). In the dynamic case the inertia forces cannot be neglected, however, for small fluid pressure gradients leading to slow interstitial flows in pores under slow deformation and within the limits of the small strain theory, the effective poroelasticity coefficients of the material describing the macroscopic behaviour can be obtained by homogenization of the static fluid-structure interaction at the level of micropores, see also Rohan et al. (2013).

Behaviour of the FSPM at finite strains has been studied within the mixture theories Bowen (1976, 1982) and also using the macroscopic phenomenological approach Biot (1977); Coussy et al. (1998). The former approach has been pursued in works which led to the theory of porous media (TPM) De Boer (1998); Diebels and Ehlers (1996); Ehlers and Bluhm (2002), including elastoplasticity Borja and Alarcón (1995), or extended for multi-physics problems, see e.g. Markert (2008). Micromechanical approaches based on a volume averaging were considered in de Buhan et al. (1998); Dormieux et al. (2002). The issue of compressibility for the hyperelastic solid phase was considered in Gajo (2010), cf. De Boer and Bluhm (1999). A very thorough revision of the viscous flows in hyperelastic media with inertia effects was given in Chapelle and Moireau (2014). There is a vast body of the literature devoted to the numerical modelling of porous media under large deformation, namely in the context of geomechanical applications, thus, taking into account plasticity of the solid phase Carter et al. (1979); Prevost (1982); Nazem et al. (2008). Although the fluid retention in unsaturated media is so far not well understood in the context of deterministic models, a number of publications involve this phenomenon in computational models Meroi et al. (1995); Song and Borja (2013). Concerning the problem formulation, three classical approaches can be used. Apart of the total Lagrangian (TL) formulations employed in the works related to the TPM De Boer (1998); Diebels and Ehlers (1996); Ehlers and Bluhm (2002), the updated Lagrangian (UL) formulation has been considered in Meroi et al. (1995); Li et al. (2004); Shahbodagh-Khan et al. (2015). As an alternative, the ALE formulation provides some advantages, like the prevention of the FE mesh distorsions, cf. Nazem et al. (2008), cf. Brown et al. (2014). Transient dynamic responses of the hyperelastic porous medium has been studied in Li et al. (2004), within the framework of the UL formulation with the linearization developed in Borja and Alarcón (1995). In Shahbodagh-Khan et al. (2015) the model has been enriched by the hydraulic hysteresis in unsaturated media.

It is worth to recall the generally accepted open problem of “missing one more equation” in the phenomenological theories of the FSPM. This deficiency can be circumvented by upscaling. For this, in particular, the homogenization of periodic media can be used, although several other treatments have been proposed, see e.g. Gray and Schrefler (2007) where the Biot’s theory was recovered within the linearization concept. The conception of homogenization applied to upscale the Biot continuum in the framework of the updated Lagrangian formulation has been proposed Rohan et al. (2006), cf. Rohan (2006); recently in Brown et al. (2014), where the ALE formulation was employed and some simplification were suggested to arrive at a computationally tractable two-scale problem. In Rohan and Lukeš (2015) we described a nonlinear model of the Biot-type poroelasticity to treat situations when the deformation has a significant influence on the permeability tensor controlling the seepage flow and on the other poroelastic coefficients. Under the small deformation assumption and using the first order gradient theory of the continuum, the homogenized constitutive laws were modified to account for their dependence on the macroscopic variables involved in the upscaled problem. For this, the sensitivity analysis known from the shape optimization was adopted.

In this paper we propose a consistent Eulerian incremental formulation for the FSPM which is intended as the basis for multi-scale computational analysis of FSPM undergoing large deformation. In the spatial configuration, we formulate the dynamic equilibrium equation using the concept of the effective stress involved in the Biot model. The fluid redistribution is governed by the Darcy flow model with approximate inertia effects, and by the mass conservation where the pore inflation is controlled through the Biot stress coupling coefficients. The compressibility of the solid and fluid phases are respected by the Biot modulus. The solid skeleton is represented by the neo-Hookean hyperelastic model employed in Li et al. (2004), although any finite strain energy function could be considered. First, the rate formulation, i.e. the Eulerian formulation, is derived by the differentiation of the residual function with respect to time using the conception of the material derivative consistent with the objective rate principle. The time discretization leads to the updated Lagrangian formulation whereby the out-of-balance terms are related to the residual function associated with the actual reference configuration. Although, in this paper, we do not describe behaviour at the pore level, the proposed formulation is coherent with the homogenization procedure which can provide the desired effective material parameters at the deformed configuration. In relationships with existing works Rohan et al. (2006); Rohan (2006) and Brown et al. (2014), this is an important aspect which justifies the proposed formulation as an improved modelling framework for a more accurate multi-scale computational analysis of the porous media.

The paper is organized, as follows. In Section 2 we present a model of the hyperelastic FSPM described in the spatial configuration and formulate the nonlinear problem. The inertia effects related to both phases are accounted for in Definition 2.2, although in Definition 2.2 we consider the simplified dynamics. Section 3 constitutes the main part of the paper. There we introduce the rate Eulerian formulation using the above mentioned differentiation of the residual function. Then we explain the time discretization which leads to the finite time step incremental formulation suitable for the numerical computation using the finite element (FE) method. In this procedure, we neglect the inertia terms associated with the fluid seepage in the skeleton. Definitions 3.3.1 and 3.3.2 provide two approximate incremental formulations related the predictor and corrector steps, respectively. In Section 4, we provide a numerical illustration and validation of the incremental formulation using 2D examples considered in the work Li et al. (2004), namely the 1D compression test and the consolidation of a partially loaded layer. For the validation we use the comparison between responses obtained by the proposed model with those obtained by linear models. Besides inspection of the displacement fields, we show the influence of the time step length on the fluid content, strain energy and the dissipation.

Basic notations

Through the paper we shall adhere to the following notation.The spatial position in the medium is specified through the coordinates with respect to a Cartesian reference frame. The Einstein summation convention is used which stipulates implicitly that repeated indices are summed over. For any two vectors , the inner product is , for any two 2nd order tensors the trace of is . The gradient with respect to coordinate is denoted by . We shall use gradients of the vector fields; the symmetric gradient of a vector function u, i.e. the small strain tensor is where the transpose operator is denoted by the superscript . The components of a vector u will be denoted by for , thus . To understand the tensorial notation, in this context, , so that , where . By we denote the closure of a bounded domain . n is a unit normal vector defined on a boundary , oriented outwards of . The real number set is denoted by . Other notations are introduced through the text.

2 Model of large deforming porous medium

Behaviour of the FSPM is governed by the equilibrium equations and the mass conservation governing the Darcy flow in the large deforming porous material. Initial configuration associated with material coordinates , and with domain is mapped to the spatial (current) configuration at time associated with spatial coordinates and with domain . Thus, where is continuously differentiable. By we denote the deformation gradient given by .

2.1 Governing equations – Biot model

We shall first introduce the governing equations for the problem of fluid diffusion through hyperelastic porous skeleton. These involve the Cauchy stress tensor , the displacement field , the bulk pressure111We adhere to the standard sign convention: the positive pressure means the negative volumetric stress induced by compression. and the perfusion velocity which describes motion of the fluid relative to the solid phase. The constitutive laws for reads as


where is the pore fluid pressure, expresses the relative volume change of the skeleton, is the left deformation tensor, is the Biot coupling coefficient. Further, by “eff” we refer to the effective (hyperelastic) stress, which is related to a strain energy function. Here we consider a neo-Hookean-type material model employed in Li et al. (2004), which is parameterized by and the Lamé constant ; for a small strain approximation, expresses the shear modulus.

The inertia of the solid and fluid phases is expressed by the skeleton acceleration and by the relative acceleration of the fluid


where are the solid and fluid phase velocities, respectively222This expression for is merely approximate since the dot means the material derivative with respect to the skeleton, for details see e.g. . Chapelle and Moireau (2014).. It is now possible to establish the Darcy law extended for dynamic flows in the moving porous skeleton. According to Coussy Coussy (2004)


where is the hydraulic permeability associated with the well known Darcy law, is the volume force acting on the fluid and a is an additional acceleration accounting for the tortuosity effect. The fluid acceleration can be approximated using the one of the skeleton and using the relative acceleration introduced in (3) multiplied by a correction factor , see also Coussy (2004).

Thus, the balance of forces (equilibrium equation) and the volume conservation read as


where is the Biot compressibility, g is the volume force acting on the mixture, is the mean density.

We remark, that the system (1)-(6), first introduced by M. Biot, see e.g. Coussy (2004); de Boer (2000), is based on the mixture theory. In such approach the microstructure is represented by volume fractions, whereas any other features of the microstructural geometry and topology can only be regarded by values of the constitutive coefficients. Although, in many realistic situations the discussed model cannot be recovered by a more rigorous homogenization treatment (which reflects genuinely the geometry of the microstructure and takes in to account internal friction induced by the fluid viscosity) — this phenomenological model can still be used to treat those media for which a detailed description of the microstructure would be too complex and practically impossible.

2.2 Nonlinear problem formulation

The differential equations (5)-(6) which hold in , are supplemented by the boundary conditions on which is decomposed into boundary segments according to the following splits:


In general, we may consider


Since the first and the second time derivatives are involved in (4)-(6), the initial conditions should involve , and also . It is worth noting that these imposed in must be consistent with the boundary conditions on at time .. On , the impermeable boundary will be considered for the sake of simplicity, i.e. .

In order to establish the weak formulation of the problem defined in Section 2.2, we need the admissibility sets for displacements and for pressures admissibility sets (we do not specify functional spaces, assuming enough of regularity for all involved unknown functions)


and the associated spaces for the test functions


We consider the force equilibrium equation and the fluid content conservation equation satisfied in the weak sense. Further we assume that the material properties ensure the following properties:


whereby K and B are positive definite.

Remark 1. Incompressibility issue. When an incompressible solid is considered, the stress decomposition in the solid phase reads , thus, in (1). This is the Terzaghi effective stress principle confirmed by both the micromechanical and homogenization approaches, see e.g. de Buhan et al. (1998) and Rohan et al. (2015). It is worth noting, that the fluid compressibility has no influence on the stress decomposition. In general, in (6) reflects the bulk compressibility of both the phases. For an incompressible solid, . Assuming the solid incompressibility, expresses the porosity change, thus, the fluid content variation due to the skeleton deformation. Therefore, the strain energy of the hyperelastic solid should be independent of . However, in this paper we employ the same neo-Hookean material used in Li et al. (2004) where the strain energy depends on , although in the mass conservation the solid phase is considered as incompressible. We admit this inconsistency to be able to validate the formulation proposed in the present paper using the results obtained in Li et al. (2004). It is worth to remark, that the issue of compressible components has been discussed in a detail in Gajo (2010).

We can now establish the weak formulation of the dynamic flow-deformation problem which respects the inertia terms induced by the mixture and also by the relative fluid flow in the skeleton, cf. Chapelle and Moireau (2014).

Definition 1. (Weak solution for the three field formulation) The weak solution of the problem introduced in Section 2.2 is the triplet which for all satisfies




The time discretization and the numerical solutions will be considered only for situations when reduced dynamics provides a sufficient physical approximation. By neglecting effects of the accelerations , the following two field formulation can be obtained.

Definition 2. (Weak solution of the two filed formulation — reduced dynamics) The weak solution of the problem introduced in Section 2.2 is the couple which for all satisfies




In both the Definitions 2.2 and 2.2, the domain , such that . It reveals the implicit definition of the solution being defined in the domain which itself depends on the solution. To avoid such a complication, the three classical approaches, i.e. the TL, UL, or ALE formulations can be used. We use the incremental UL (updated Lagrangian) formulation where the recent domain serves as the reference configuration, see e.g. Crisfield (1997); Simo and Hughes (1998) for details.

Besides the kinematics associated with the strain filed, also the solution-dependent volume fraction of the fluid and the fluid density contribute to the nonlinearity of the model equations,


where is the reference pressure and the last expression of can be approximated by .Nevertheless, below we confine our treatment to situations, where can be neglected.

The material coefficients , B and K are defined using the following formulas:


where is the bulk modulus of the fluid, is a constant permeability tensor and .

3 Consistent incremental formulation

In this section we derive the consistent incremental formulation for the model introduced above. It will be done in two steps. By time differentiation with a given perturbation velocity field introduced in (18), associated with the configuration transformation, first we derive the rate of the residual function which allows to define an approximated problem for computing increments of the displacement and pressure field associated with a time increment . As a next step, the time discretization is considered. A suitable approximation of the time derivatives occurring in the rate formulation leads to a linear subproblem which can be solved for the increments and, consequently, to updated state variables of the porous medium.

We shall need the perturbation velocity field defined in which satisfies . This allows us to introduce the perturbed position and the perturbed domain


where is the perturbation time. Using the convection velocity , we can express the associated material derivative which will be used to derive the rate form of the residuum in (15).

Remark 2. Material derivative and convection field . The material derivative is introduced by virtue of the perturbation (18). It is worth noting that this differentiation is used in the sensitivity analysis of functionals associated with optimal shape problem where the partial differential constraints are considered. In this context, the convection velocity field called “the design velocity” is established to link perturbations of the designed boundary with perturbations of the material point in the domain. The following formulae can be derived easily, see e.g. Haslinger and Neittaanmaki (1988); Haug et al. (1986).


The incremental formulation which is defined below is based on (14) rewritten for time , whereby the residual function can be approximated by , as follows. Let and assume the associated configuration (the primary state fields and all dependent quantities at ) can be estimated, while we assume that is known333In fact, an approximation is known rather than an exact residual which should be zero. In general, denoting by the test functions


where is the increment is the differential of the residual functional evaluated at time with respect to the state variation and configuration perturbation associated with the time increment . For , (20) provides a first order approximation of using the tangent defined at . In general, a secant-based approximation can be established when an estimate of is available at . This approach is followed in this work; in Section 3.3 we introduce a predictor-corrector algorithm within the time discretization scheme applied to (20).

Definition 3. (Incremental formulation) Given fields and defined in and , respectively, the new state at time satisfies




pointwise at any . The updated domain is obtained using updated coordinates, .

This definition can be applied to introduce a time integration numerical scheme based on the predictor-corrector steps: after a time discretization, by putting , (21) yields the predictor problem for computing an approximate increment which can be used to introduce an intermediate state at . Then (21) provides a “corrected approximate increment”.

3.1 The Lie derivative of the Cauchy stress

To derive the rate form of the force-equilibrium equation we need the Lie derivative of the Cauchy stress which is related to the 2nd Piola-Kirchhoff stress S by


recalling the notation . Since the perturbation of deformation gradient F is , the Lie derivative of the Cauchy stress and of the Kirchhoff stress are related to the material derivatives and by the following expressions




Application of the material derivative with respect to the convection field , see Remark 3 to differentiate the first integral in (15) yields (we abbreviate )


where is defined using the constitutive relationship. In the poroelastic medium, the stress is decomposed into its effective and volumetric parts: . The Lie derivative of the effective part of the Cauchy stress can be expressed in terms of the tangential stiffness tensor and linear velocity strain (with respect to the spatial coordinates, i.e. where )


Above is the tangential elasticity tensor which can be identified for a given strain energy function. The Lie derivative of the volumetric part can be expressed,


and substituted into (26) to express the terms involving , thus


3.2 Rate formulation in the Euler configuration

The rate form is introduced to derive a consistent linearization of the nonlinear problem. We shall need the material derivative of the configuration-dependent variables. The perturbation defined by (18) induces perturbation of , and ,


The following differentials of K and hold, see (17):


Obviously, , since is constant.

In general we consider non-conservative loading by boundary tractions acting on . The classical result of the sensitivity analysis in shape optimization yields the material derivative of the virtual power induced by the traction forces h,


where n is the unit normal vector to and is the 2-multiple of the mean curvature of the deformed surface.

3.2.1 Incremental operators

The following expressions are obtained from (26)-(30). We consider a reference state and its perturbation associated with the convection . Differentiation of yields


Differentiation of yields


In both (33) and (34) we need to express w and its perturbation which can be obtained pointwise in ,


Using (33)-(35) substituted in equation (21), we obtain the problem for computing increments which are related to a finite time increment . As the convection field is associated with and since also the rates , and must be expressed in terms of , (21) is not a linear problem which could be solved directly. To be able to solve this evolutionary problem in time, a suitable time discretization must be introduced, as will be explained in the next section.

3.3 Time discretization of the incremental subproblem

In the present section we establish a semidiscretized formulation which leads to linear subproblems defined at each time level within the framework of the updated reference configurations and time integration scheme based on the rate formulation treated in the preceding section. Obviously, time discretization of problem (21) is not unique. Implicit Runge-Kutta methods are introduced by means of intermediate time stages, see e.g. Park and Lee (2009). In computational mechanics, the Newmark- integration schemes which involve only two consecutive time levels became very popular Li et al. (2004); Shahbodagh-Khan et al. (2015). When solving a nonlinear problem with dissipation, the naturally arising but seldom considered question concerns with the influence of the time step length and the number of iterations employed to solve the nonlinear problem by the Newton-Raphson method on the dissipated energy.

Here we consider a multistep approach based on combination of forward and backward finite differences to discretize the time derivatives and the convection ; optionally a two-level computational scheme of the “leap frog” type, Park and Lee (2009) can be defined which uses intermediate time levels. For this we introduce a uniform time discretization of the time interval upon introducing the time levels , , such that . By we refer to an approximated value of . Further, we define and consider the following approximation of time derivatives at ,


Below we propose a predictor-corrector time-integration scheme, whereby the corrector uses the rate of the residual at an interpolated time level. The need for the 3rd time derivative in the rate form equations (in fact is involved) apparently induces the 4 step numerical scheme, see (36), however, as will be shown later, the predictor step circumvents this complication so that 3 consecutive steps are involved only. In order to get a linear subproblem arising from (33)-(35), the convection velocity field involved in these expressions can be approximated by the backward difference, i.e. , while the time derivatives are approximated using (36) with . For the corrector step, will be approximated using an estimate of computed by the predictor.

In the rest of this paper, we shall confine to the simplified dynamics of the porous medium which disregards the inertia induced by the relative fluid-solid motion. The time discretization of (15) yields the following approximation


for all .

The discretization of the incremental operators (33)-(35) evaluated at time with yields expressions where we use the following notation and ; all other coefficients are evaluated according to the configuration . From (33) we get


Eq. (35) is multiplied by , so that in (15) is replaced by . The discretization results in


where the differential of the seepage velocity is given by the following approximation:


Above can be substituted by .

In order to simplify notation, we shall introduce the following tensors and material coefficients depending on the solution at time and on other variables associated with time :


Further, we shall employ an abbreviated notation:


Two incremental linearized subproblems will be formulated which provide the solutions of the problem introduced in the Definition 3. In this context, we emphasize the different meaning of u and which now denote the increments and . By and we denote sets of admissible increments and respectively. For a given time discretization with a time step , the sets and are defined according to the Dirichlet boundary conditions considered according to (36). Also the tensors and are defined in , and the volume load is given at .

3.3.1 The predictor step

The predictor step is based on the approximation (20) with evaluated at


where is approximated by , thus, recalling (42). For all quantities associated with the configuration at we use the following simplified notation: ; all other coefficients are evaluated according to the configuration .

By virtue of Definition 3, the following is to be solved: Find , such that