Defining the Constitutive Response of Fluid Transitioning from Darcy Flow to Poiseuille Flow

The cohesive element fluid flow model:

is typically used in geotechnical applications, where fluid flow continuity within the cohesive element and through the interface must be maintained;
supports the transition from Darcy flow to Poiseuille flow (gap flow) as damage in the element initiates and evolves;
enables modeling of an additional resistance layer on the surface of the cohesive element to model fluid leakoff into the formation;
enables fluid pressure on the cohesive element surface to contribute to its mechanical behavior, which enables the modeling of hydraulically driven fracture;
can be used only in conjunction with traction-separation behavior;
supports fluid flow continuity between intersecting layers of cohesive pore pressure elements; and
enables gravity-induced fluid flux modeling.

This page discusses:

Defining Pore Fluid Flow Properties
Tangential Flow
Normal Flow across Gap Surfaces
Flow Flux Induced by Gravity
Transition from Darcy Flow to Poiseuille Flow
Additional Considerations
Output

Defining Pore Fluid Flow Properties

The fluid constitutive response consists of the following:

Tangential flow along the cohesive element midplane, which can be modeled as either Darcy or Poiseuille flow; and
Normal flow (also referred to as leakoff) across the cohesive element, which can reflect resistance due to caking or fouling effects.

You assign the tangential and normal flow properties separately.

The flow patterns of the pore fluid in the element are shown in Figure 1. The fluid is assumed to be incompressible, and the formulation is based on a statement of flow continuity that considers tangential and normal flow and the rate of opening of the cohesive element.

Tangential Flow

Tangential flow in a cohesive element will transition from Darcy flow to Poiseuille flow as damage in the element initiates and evolves. The transition is designed to approximate the changing nature of fluid flow through an initially undamaged porous material (Darcy flow) to flow in a crack (Poiseuille flow) as the material is damaged. You must specify the fluid constitutive response for both types of flow.

Gap Opening

The tangential flow equations for the cohesive element are solved within a gap along the length of the element. The gap opening, $d$ , is defined as

d = t_{c u r r} - t_{o r i g} + g_{i n i t} = \hat{d} + g_{i n i t},

where $t_{c u r r}$ and $t_{o r i g}$ are the current and original cohesive element geometrical thicknesses, respectively; $g_{i n i t}$ is the initial gap opening, which has a default value of 0.002; and $\hat{d}$ represents the physical crack opening that is used for Poiseuille flow once the element is damaged. A schematic illustration is shown in Figure 2. $g_{i n i t}$ is not a physical quantity. It is used by Abaqus/Standard to ensure that the flow equations can be solved robustly when the physical gap is closed (i.e., $\hat{d} = 0$ ). As $\hat{d}$ increases, the effect of $g_{i n i t}$ on the flow equations is diminished, as described in Transition from Darcy Flow to Poiseuille Flow.

Darcy Flow

Darcy flow defines a simple relationship between the volumetric flow rate of a fluid and the fluid pressure gradient in a porous material. The relationship is defined by the expression

q d = - \frac{k d}{γ_{w}} \nabla p,

where $k$ is the permeability, $\nabla p$ is the pressure gradient along the cohesive element, $d$ is the gap opening, and $γ_{w}$ is the fluid specific weight.

Poiseuille Flow

In Abaqus/Standard Poiseuille flow within cohesive elements refers to the steady viscous flow between two parallel plates. For this flow, you can specify either a Newtonian fluid or a power law fluid.

Newtonian Fluid

In the case of a Newtonian fluid the volume flow rate density vector is given by the expression

q d = - k_{t} \nabla p,

where $k_{t}$ is the tangential permeability (the resistance to the fluid flow), $\nabla p$ is the pressure gradient along the cohesive element, and $d$ is the gap opening.

Abaqus defines the tangential permeability, or the resistance to flow, according to Reynold's equation:

k_{t} = \frac{d^{3}}{12 μ},

where $μ$ is the fluid viscosity and $d$ is the gap opening. You can also specify an upper limit on the value of $k_{t}$ .

Power Law Fluid

In the case of a power law fluid the constitutive relation is defined as

τ = K {\dot{γ}}^{α},

where $τ$ is the shear stress, $\dot{γ}$ is the shear strain rate, $K$ is the fluid consistency, and $α$ is the power law coefficient. Abaqus defines the tangential volume flow rate density as

q d = - (\frac{2 α}{1 + 2 α}) {(\frac{1}{K})}^{\frac{1}{α}} {(\frac{d}{2})}^{\frac{1 + 2 α}{α}} {∥ \nabla p ∥}^{\frac{1 - α}{α}} \nabla p,

where $d$ is the gap opening.

Bingham Plastic Fluid

In the case of a Bingham plastic fluid the volume flow rate density vector is given by the expression

q d = - (\frac{1}{3}) (\frac{1}{K}) {(\frac{d}{2})}^{3} \nabla p {(1 - \frac{2 τ_{0}}{2 d ‖ \nabla p ‖})}^{2} (2 + \frac{2 τ_{0}}{d ‖ \nabla p ‖}),

where $K$ is the fluid consistency, $τ_{0}$ is the yield stress, and $d$ is the gap opening. The unyielded fluid is modeled as a Newtonian fluid with viscosity equal to $P_{τ} \times K$ , where $P_{τ}$ has a default value of 10⁷.

Herschel-Bulkley Fluid

In the case of a Herschel-Bulkley fluid the volume flow rate density vector is given by the expression

q d = - 2 [(\frac{α}{1 + α}) (\frac{α}{1 + 2 α})] {(\frac{1}{K})}^{\frac{1}{α}} {(\frac{d}{2})}^{\frac{1 + 2 α}{α}} {‖ \nabla p ‖}^{\frac{1 - α}{α}} {(1 - \frac{2 τ_{0}}{d ‖ \nabla p ‖})}^{\frac{1 + α}{α}} (\frac{1}{α} + 1 + \frac{2 τ_{0}}{d ‖ \nabla p ‖}) \nabla p,

where $K$ is the fluid consistency, $α$ is the power law coefficient, $τ_{0}$ is the yield stress, and $d$ is the gap opening. The unyielded fluid is modeled as a Newtonian fluid with viscosity equal to $P_{τ} \times K$ , where $P_{τ}$ has a default value of 10⁷.

Normal Flow across Gap Surfaces

You can permit normal flow by defining fluid leak-off coefficients for the pore fluid material. These coefficients define a pressure-flow relationship between the cohesive element's middle nodes and its adjacent surface nodes. The fluid leak-off coefficients can be interpreted as the permeability of a finite layer of material on the cohesive element surfaces, as shown in Figure 3.

The normal flow is defined as

q_{t} = c_{t} (p_{i} - p_{t})

and

q_{b} = c_{b} (p_{i} - p_{b}),

where $q_{t}$ and $q_{b}$ are the flow rates into the top and bottom surfaces, respectively; $c_{t}$ and $c_{b}$ are the fluid leak-off coefficients at the top and bottom element surfaces, respectively; $p_{i}$ is the midface pressure; and $p_{t}$ and $p_{b}$ are the pore pressures on the top and bottom surfaces, respectively.

Defining Leak-off Coefficients as a Function of Temperature and Field Variables

You can optionally define leak-off coefficients as functions of temperature and field variables.

Defining Leak-off Coefficients in a User Subroutine

User subroutine UFLUIDLEAKOFF can also be used to define more complex leak-off behavior, including the ability to define a time accumulated resistance, or fouling, through the use of solution-dependent state variables.

Flow Flux Induced by Gravity

In the presence of a distributed gravity load the tangential flow rate density vector is given by the expression

q d = - k_{t} (- ρ g_{t}),

where $k_{t}$ is the tangential permeability as defined above, $g_{t}$ is the projection of the gravity vector onto the midsurface of the cohesive element, and $ρ$ is the pore fluid density. For Darcy flow,

k_{t} = \frac{k d}{γ_{w}} .

For Poiseuille flow, in the case of a Newtonian fluid,

k_{t} = \frac{d^{3}}{12 μ} .

In the case of a power law fluid,

k_{t} = (\frac{2 α}{1 + 2 α}) {(\frac{1}{K})}^{\frac{1}{α}} {(\frac{d}{2})}^{\frac{1 + 2 α}{α}} {∥ \nabla p ∥}^{\frac{1 - α}{α}} .

Transition from Darcy Flow to Poiseuille Flow

For a Newtonian fluid the transition from Darcy flow to Poiseuille flow as a function of the damage variable, $D$ , is described by the expression

q d = - ((1 - D \hat{F} (\hat{d})) \frac{k}{γ_{w}} g_{i n i t} + D \hat{F} (\hat{d}) (\frac{{(\hat{d})}^{3}}{12 μ})) (\nabla p - ρ g_{t}),

\hat{F} (\hat{d}) = {\begin{matrix} 0 & \hat{d} < 0 \\ \frac{\hat{d}}{g_{i n i t}} & 0 \leq \hat{d} \leq g_{i n i t} \\ 1 & g_{i n i t} < \hat{d} \end{matrix} .

The above relationship also supports a transition from Poiseuille flow back to Darcy flow as the physical gap, $\hat{d}$ , in a damaged element closes. The flow transition equation for a power law fluid is obtained similarly.

Initially Open Elements

You can define an initial gap to identify elements that are fully damaged; that is, $D = 1.0$ at the integration points of the elements.

Assigning Initial Damage Values

You can define an initial gap to identify elements and assign $D$ directly to the integration points. If you assign an initial damage variable to any of the integration points but not all of them, a value of $D = 0.0$ is assigned to the integration points to which you did not assign a value.

If an element set is used, you must ensure that all elements within the set have the proper uniform order of integration points.

Additional Considerations

Your use of cohesive element fluid properties and your property values can impact your solution in some cases.

Unsymmetric Matrix Storage and Solution

The pore pressure cohesive element matrices are unsymmetric; therefore, unsymmetric matrix storage and solution may be needed to improve convergence (see Matrix Storage and Solution Scheme in Abaqus/Standard).

Large Coefficient Values

You must make sure that the tangential permeability or fluid leak-off coefficients are not excessively large. If either coefficient is many orders of magnitude higher than the permeability in the adjacent continuum elements, matrix conditioning problems may occur, leading to solver singularities and unreliable results.

Meshing Requirement at Intersections of Cohesive Elements

When different layers of cohesive pore pressure elements intersect, a common midsurface node must be shared by all elements to support fluid flow continuity. Figure 4 shows a two-dimensional mesh example of intersecting elements. Elements 10, 20, 30, and 40 share the same middle node, 100, at the intersecting point.

Pore Fluid Contact Property Consideration

To maintain pore pressure continuity between contact surfaces, you can define a pore fluid contact property (see Pore Fluid Contact Properties). If pore pressure degrees of freedom exist on both sides of a contact interface, the pore fluid contact property will enforce some pore pressure continuity. The default behavior is infinite permeability, which is equivalent to no resistance to the fluid flux between the two contact surfaces. If the top and bottom surfaces of the pore pressure cohesive elements are used to define the contact surface, the fluid flux should be controlled by the cohesive element definition. Therefore, the pore fluid contact property definition should be deactivated. You can remove its effects by setting the contact permeability to be zero and setting the cutoff gap fill to be zero.

Output

The following output variables are available when flow is enabled in pore pressure cohesive elements:

GFVR: Gap fluid volume rate.
PFOPEN: Fracture opening.
LEAKVRT: Leak-off flow rate at element top.
ALEAKVRT: Accumulated leak-off flow volume at element top.
LEAKVRB: Leak-off flow rate at element bottom.
ALEAKVRB: Accumulated leak-off flow volume at element bottom.
FLDVEL: Material point fluid velocity.