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

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.

Input File Usage

Use the following option to define the initial gap opening directly:

SECTION CONTROLS, INITIAL GAP OPENING

Abaqus/CAE Usage

Initial gap opening is not supported in Abaqus/CAE.

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.

Input File Usage

Use the following option to define fully saturated isotropic permeability:

PERMEABILITY, TYPE=ISOTROPIC, SPECIFIC= $γ_{w}$

Abaqus/CAE Usage

Property module: material editor: OtherPore FluidPermeability: Type: Isotropic, Specific weight of wetting liquid:  $γ_{w}$

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}$ .

Input File Usage

Use the following option to define the tangential flow in a Newtonian fluid:

GAP FLOW, TYPE=NEWTONIAN, KMAX

Abaqus/CAE Usage

Property module: material editor: OtherPore FluidGap Flow: Type: Newtonian: Toggle on Maximum Permeability and enter the value of  $k_{\max}$

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.

Input File Usage

GAP FLOW, TYPE=POWER LAW

Abaqus/CAE Usage

Property module: material editor: OtherPore FluidGap Flow: Type: Power law

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⁷.

Input File Usage

GAP FLOW, TYPE=BINGHAM PLASTIC, PTAU= $P_{τ}$

Abaqus/CAE Usage

The Bingham plastic fluid model is not supported in Abaqus/CAE.

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⁷.

Input File Usage

GAP FLOW, TYPE=HERSCHEL-BULKLEY, PTAU= $P_{τ}$

Abaqus/CAE Usage

The Herschel-Bulkley fluid model is not supported in Abaqus/CAE.

Multiple Fluid Flow

Abaqus allows the use of multiple fluids in an analysis, up to a maximum of four fluid types. This capability is useful in situations where the spatial distribution of the different fluid types is known (precomputed) at all times during the analysis. In other words, the type of fluid at each node in the domain is known at all times during the analysis.

The definition of multiple fluids in an analysis involves assigning a unique integer value that acts as an identifier for each fluid type and associating each fluid identifier with a specific predefined field variable used to specify the spatial distribution of the fluid as a function of time. You can utilize the fluid identifier to specify the viscosity of each fluid. Only the fluid viscosity can be different for different fluids. All other properties (for example, density, thermal, etc.) are assumed to be same for all the fluids. The field variable associated with a fluid type allows you to predefine the fluid type at a node as a function of time using amplitude definitions or as functions of both time and space utilizing user subroutine UFIELD or USDFLD.

You specify a value for each field variable, associated with a fluid identifier, at each node. It is recommended that you specify a value for the field variable ( $F V$ ) such that $0 \leq F V \leq 1.0$ . Abaqus assigns the fluid type at an integration point to be the one for which the field variable has the maximum value at that point. If all field variables have the same numerical value at an integration point, Abaqus assumes the first fluid (fluid identifier equal to 1) to be active.

Input File Usage

*GAP FLOW, FLUID ID=Fluid Id number, DEFINITION=NEWTONIAN or POWER LAW or HERSCHEL-BULKLEY or BINGHAM PLASTIC
FIELD, VARIABLE=Field variable number
FLUID ID FIELD
Fluid Id number, Field variable number

Abaqus/CAE Usage

Specifying multiple fluid flow is not supported in Abaqus/CAE.

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.

Input File Usage

FLUID LEAKOFF

Abaqus/CAE Usage

Property module: material editor: OtherPore FluidFluid Leakoff: Type: Coefficients

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.

Input File Usage

FLUID LEAKOFF, DEPENDENCIES

Abaqus/CAE Usage

Property module: material editor: OtherPore FluidFluid Leakoff: Type: Coefficients: Toggle on Use temperature-dependent data and select the number of 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.

Input File Usage

FLUID LEAKOFF, USER

Abaqus/CAE Usage

Property module: material editor: OtherPore FluidFluid Leakoff: Type: User

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 - α}{α}} .

Input File Usage

Use the following option to specify the density of the pore fluid:

DENSITY, PORE FLUID

Abaqus/CAE Usage

Property module: material editor: GeneralDensity

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.

Input File Usage

Use the following option to define the initial gap directly:

INITIAL CONDITIONS, TYPE=INITIAL GAP
element number or element set, omit values for D

Abaqus/CAE Usage

Initial gap definition is not supported in Abaqus/CAE.

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.

Input File Usage

Use the following option to assign initial damage values:

INITIAL CONDITIONS, TYPE=INITIAL GAP
element number or element set, D at each integration point

Abaqus/CAE Usage

Initial gap definition is not supported in Abaqus/CAE.

Additional Considerations

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

Cohesive Offset

A cohesive offset (see Cohesive Offset)refers to a value of the separation in the normal direction, $δ_{n} = δ_{n}^{o f f}$ , below which the resistance to additional closure $(δ_{n} > 0, {\dot{δ}}_{n} < 0)$ of the cohesive surfaces increases significantly. This capability is useful in preventing the closure of a fracture when changes in external loading conditions might otherwise result in closure. A practical application of this capability occurs during hydraulic fracturing with slurry, where solid proppant particles (the solid part of the slurry) help ‘prop open’ the fracture when external hydraulic pressure loads are removed.

Unsymmetric Matrix Storage and Solution

The pore pressure cohesive element matrices are unsymmetric; therefore, unsymmetric matrix storage and solution might 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.
FLVF: When multiple fluids are defined, the components of FLVF represent the volume fraction of each fluid type. In particular, FLVF_i refers to the volume fraction of fluid type i.