Discretized equilibrium statement for a porous medium

Abaqus/Standard provides an equation to determine the discretized equilibrium in a porous medium.

Related Topics
In Other Guides
Coupled pore fluid diffusion and stress analysis
Geostatic stress state

ProductsAbaqus/Standard

Equilibrium is expressed by writing the principle of virtual work for the volume under consideration in its current configuration at time t:

Vσ:δεdV=StδvdS+Vf^δvdV,

where δv is a virtual velocity field, δε=defsym(δv/x) is the virtual rate of deformation, σ is the true (Cauchy) stress, t are surface tractions per unit area, and f^ are body forces per unit volume.

For our system f^ will often include the weight of the wetting liquid,

fw=(sn+nt)ρwg,

where ρw is the density of the wetting liquid and g is the gravitational acceleration, which we assume to be constant and in a constant direction (so that, for example, the formulation cannot be applied directly to a centrifuge experiment unless the model in the machine is small enough that g can be treated as constant). For simplicity we consider this loading explicitly so that any other gravitational term in f^ is associated only with the weight of the dry porous medium. Thus, we write the virtual work equation as

(1)Vσ:δεdV=StδvdS+VfδvdV+V(sn+nt)ρwgδvdV,

where f are all body forces except the weight of the wetting liquid.

In a finite element model equilibrium is approximated as a finite set of equations by introducing interpolation functions. The notation used to indicate such discretization are those quantities with uppercase superscripts (for example, vN), which represent nodal variables, with the summation convention adopted for the superscripts. The interpolation is assumed to be based on material coordinates in the material skeleton (a “Lagrangian” formulation).

For simplicity, in this section we consider only the case where the problem has no internal constraints—such as incompressibility—and the discretization is made entirely by approximating equilibrium: this results in the displacement (or stiffness) method. Mixed formulation (“hybrid”) elements are available for porous medium analysis with Abaqus/Standard, but consideration of such formulations does not require any important extension of the development at this stage.

The virtual velocity field is interpolated by

δv=NNδvN,

where NN(Si) are interpolation functions defined with respect to material coordinates, Si.

The virtual rate of deformation is interpolated as

δε=βNδvN,

where, in the simplest case,

βN=sym(δNNx),

although more general forms are used in some of the elements in Abaqus.

The virtual work equation is thus discretized as

δvNVβN:σdV=δvN[SNNtdS+VNNfdV+V(sn+nt)ρwNNgdV],

where the δvN are assumed to be independent.

The term conjugate to δvN on the left-hand side of this equation is referred to subsequently as the internal force array, IN:

IN=defVβN:σdV.

Likewise, the external force array, PN, is taken from the right-hand side:

PN=defSNNtdS+VNNfdV+V(sn+nt)ρwNNgdV

(PN includes any d'Alembert forces).

Choosing each δvN to be nonzero in turn expresses equilibrium as a balance of internal and external forces:

(2)IN-PN=0.

These discretized equilibrium equations, together with the continuity equation discussed in Continuity statement for the wetting liquid phase in a porous medium, define the state of the porous medium. The equilibrium equations are written at the end of a time increment when implicit integration is used and, for all but the simplest cases, they are nonlinear. Newton's method is often used for their solution. Also, small, linear perturbations of the system are sometimes of interest (an example is the small vibration problem). These considerations imply a need for the Jacobian matrix of the system, which defines the variation of each term in the equations with respect to the basic variables of the discretized problem, which—for this case—are the nodal positions, xN (or, equivalently, the displacements xN-XN), and the nodal wetting liquid pressure values, uwN. Symbolically we write such a variation of a term, f say, as df, meaning

df=deffxNdxN+fuwPduP.

From the variation of discretized equilibrium, Equation 2, the term dPN gives rise to the mass matrix (for the d'Alembert forces) and the “load stiffness matrix” in the Jacobian. The load stiffness matrix is discussed in Elements and Loading and Constraints for particular load types. The load stiffness term associated with the weight of the wetting liquid is

-V1Jd[J(sn+nt)ρw]NNgdV,

where

J=def|dVdV0|

is the ratio of volume in the current configuration to volume in the reference configuration.

The term dIN is

dIN=dVβN:σdV=V[1JβN:d(Jσ)+σ:dβN]dV.

The first term includes d(Jσ), which is the variation of stress caused by variations in nodal positions and pore liquid pressure values. In a continuum sense (that is, before the spatial discretization of the solution variables) this term is defined by the effective stress principle and by the constitutive assumptions used for the material and is discussed in more detail below. Introducing the spatial discretization into the second term provides a contribution to the initial stress matrix.

Since the effective stress, σ¯, is generally stored as components associated with spatial directions, the rotation of the material during an increment must be included in the formulation. This issue is discussed in detail in Rate of deformation and strain increment, Stress rates, State storage, and Solid element formulation. For the purpose of the present development we assume that the variation of stress is

d(Jσ)=d(J(1-nt)σ¯)-d(Jntp¯tI)+J(dΩσ+σdΩT)-d(Jχuw)I,

where d(Jσ¯) is the variation in effective stress associated with constitutive response in the material (that is, caused by variations in the strain or other state variables) and dΩ=defasym(dx/x) is the spin of the material. Using this assumption, the Jacobian contribution from stress in the porous medium is

(3)dIN=V[1JβN:{d(J(1-nt)σ¯)-d(Jntp¯tI)}+σ:(dβN+2βNdΩ)-βN:I(χ+uwdχduw)duw-βN:IχuwI:dε]dV,

where dε=defsym(dx/x) is the strain rate (the “rate of deformation”) so that

dJJ=trace(dε)=I:dε.