Continuity statement for the wetting liquid phase in a porous medium

This section defines the continuity statement for the liquid flowing through the mesh in an Abaqus porous model. The medium may be wholly or partially saturated with this liquid. When the medium is only partially saturated, the remainder of the voids is filled with another fluid. An example is a geotechnical problem, with soil containing water and air: continuity is written for the water phase

This page discusses:

See Also
In Other Guides
Coupled Pore Fluid Diffusion and Stress Analysis
About Pore Fluid Flow Properties

ProductsAbaqus/Standard

The wetting liquid can attach to and, thus, be trapped by certain solid particles in the medium: this volume of trapped liquid attached to solid particles forms a “gel.”

A porous medium is modeled approximately in Abaqus by attaching the finite element mesh to the solid phase. Liquid can flow through this mesh. A continuity equation is, therefore, required for the liquid, equating the rate of increase in liquid mass stored at a point to the rate of mass of liquid flowing into the point within the time increment. This continuity statement is written in a variational form as a basis for finite element approximation. The liquid flow is described by introducing Darcy's law or, alternatively, Forchheimer's law. The continuity equation is satisfied approximately in the finite element model by using excess wetting liquid pressure as the nodal variable (degree of freedom 8), interpolated over the elements. The equation is integrated in time by using the backward Euler approximation. The total derivative of this integrated variational statement of continuity with respect to the nodal variables is required for the Newton iterations used to solve the nonlinear, coupled, equilibrium and continuity equations. This expression is also derived in this section.

Consider a volume containing a fixed amount of solid matter. In the current configuration this volume occupies space V with surface S. In the reference configuration it occupied space V0. Wetting liquid can flow through this volume: at any time the volume of such “free” liquid (liquid that can flow if driven by pressure) is written Vw. Wetting liquid can also become trapped in the volume, by absorption into the gel. The volume of such trapped liquid is written Vt.

The total mass of wetting liquid in the control volume is

Vρw[dVw+dVt]=Vρw(nw+nt)dV,

where ρw is the mass density of the liquid.

The time rate of change of this mass of wetting liquid is

ddt(Vρw(nw+nt)dV)=V1Jddt(Jρw(nw+nt))dV.

The mass of wetting liquid crossing the surface and entering the volume per unit time is

-SρwnwnvwdS,

where vw is the average velocity of the wetting liquid relative to the solid phase (the seepage velocity) and n is the outward normal to S.

Equating the addition of liquid mass across the surface S to the rate of change of liquid mass within the volume V gives the wetting liquid mass continuity equation

V1Jddt(Jρw(nw+nt))dV=-SρwnwnvwdS.

Using the divergence theorem and because the volume is arbitrary, this provides the pointwise equation

1Jddt(Jρw(nw+nt))+x(ρwnwvw)=0.

The equivalent weak form is

Vδuw1Jddt(Jρw(nw+nt))dV+Vδuwx(ρwnwvw)dV=0,

where δuw is an arbitrary, continuous, variational field. This statement can also be written on the reference volume:

V0δuwddt(Jρw(nw+nt))dV0+V0δuwJx(ρwnwvw)dV0=0.

In Abaqus/Standard this continuity statement is integrated approximately in time by the backward Euler formula, giving

V0δuw[(Jρw(nw+nt))t+Δt-(Jρw(nw+nt))t]dV0+ΔtV0δuw[Jx(ρwnwvw)]t+ΔtdV0=0,

which, over the current volume, is

Vδuw[(ρw(nw+nt))t+Δt-1Jt+Δt(Jρw(nw+nt))t]dV+ΔtVδuw[x(ρwnwvw)]t+ΔtdV=0.

We now drop the subscript t+Δt by adopting the convention that any quantity not explicitly associated with a point in time is taken at t+Δt.

The divergence theorem allows the equation to be rewritten as

V[δuw(ρwρw0(nw+nt)-1J(ρwρw0(J(nw+nt))t)-Δtρwρw0nwδuwxvw]dV+ΔtSδuwρwρw0nwnvwdS=0,

where—for convenience—we have normalized the equation by the density of the liquid in the reference configuration, ρw0.

Since nw=sn, this is the same as

(1)V[δuw(ρwρw0(sn+nt)-1J(ρwρw0J(sn+nt))t)-Δtρwρw0snδuwxvw]dV+ΔtSδuwρwρw0snnvwdS=0.