FLAC3D Theory and Background • Fluid-Mechanical Interaction

# FLAC3D Fluid-Thermal-Mechanical-Formulation — Mathematical Description

In this section, we consider fluid-mechanical coupling to the heat-conduction logic (see fluid-thermal coupling via the convection logic). Most engineering analyses involving simultaneous deformation and fluid/thermal diffusion mechanisms are carried out using uncoupling techniques (e.g., in this section). In addition to those calculation modes, FLAC3D provides for the option of coupled fluid-thermal-mechanical analysis (i.e., in which the mechanical response of a porous material can be studied under transient fluid flow and/or thermal conditions). Although this section is mainly concerned with the modeling of deformation-fluid diffusion problems, the general equations governing the fluid-thermal-mechanical response in FLAC3D are presented here for completeness.

The formulation of coupled deformation-diffusion processes in FLAC3D is done within the framework of the quasi-static Biot theory, and can be applied to problems involving single-phase Darcy flow in a porous medium. Various types of fluids, including water and oil, can be represented with this model.

The variables involved in the description of fluid flow through porous media are pore pressure, saturation and the three components of the specific discharge vector. These variables are related through the fluid mass-balance equation, Darcy’s law for fluid transport, a constitutive equation specifying the fluid response to changes in pore pressure, saturation, volumetric strains and temperature, and an equation of state relating pore pressure to saturation in the unsaturated range. Pore pressure and temperature influences are involved in the mechanical constitutive laws to complete the thermal fluid-flow mechanical coupling.

Assuming the volumetric strain and temperature rates are known, substitution of the mass balance equation into the fluid constitutive relation, using Darcy’s law, yields a differential equation in terms of pore pressure and saturation that may be solved for particular geometries, properties, boundary and initial conditions.

## Conventions and Definitions

As a notation convention, the symbol $$a_i$$ denotes component $$i$$ of the vector $$\{ a \}$$ in a Cartesian system of reference axes; $$A_{ij}$$ is component $$(i,j)$$ of tensor $$[A]$$. Also, $$f_{,i}$$ is used to represent the partial derivative of $$f$$ with respect to $$x_i$$. (The symbol $$f$$ denotes either a scalar variable or a vector component.)

The Einstein summation convention applies only on indices $$i$$, $$j$$ and $$k$$, which take the values 1, 2, 3 for components that involve spatial dimensions. The indices may take any values when used in matrix equations.

SI units are used to illustrate parameters and dimensions of variables. See System of Units for conversions to other systems of units.

The following reference quantities are useful in the characterization of transient fluid flow.

Characteristic length:

(1)$L_c = {{volume \ of \ flow \ domain}\over{surface \ area \ of \ flow \ domain}}$

Fluid diffusivity:

(2)$c = {{k} \over {1 / M + {\alpha}^2 / {\alpha}_1}}$

where $$k$$ is mobility coefficient (FLAC3D permeability), $$M$$ is the Biot modulus, $$\alpha$$ is the Biot coefficient, $${\alpha}_1 = K + 4/3 G$$, $$K$$ is the drained bulk modulus, and $$G$$ is the shear modulus of the porous material.

The Biot coefficient takes into account the grain compressibility for the porous material. If $$\alpha$$ is equal to unity, the grains are considered to be incompressible and the Biot modulus $$M$$ is equal to $$K_f / n$$, where $$K_f$$ is fluid bulk modulus and $$n$$ is porosity. The fluid diffusivity becomes

(3)$c = {{k} \over {n / K_f + 1 / {\alpha}_1}}$

Note that for flow-only calculations (rigid material), the fluid diffusivity is

(4)$c = {k M}$

See Properties and Units for Fluid-Flow Analysis for additional discussion on the relations between these properties.

Characteristic time:

(5)$t_c = {{L_c^2}\over{c}}$

## Governing Differential Equations

The differential equations describing the fluid-thermal-mechanical response of a porous material, for heat transported by conduction in the material domain, are as follows.

Transport Law

The fluid transport is described by Darcy’s law:

(6)$q_i = - k_{il} \ \hat{k}(s) \left[ p - {\rho}_f x_j g_j \right] _{,l}$

where $$q_i$$ is the specific discharge vector, $$p$$ is fluid pore pressure, $$k$$ is the tensor of absolute mobility coefficients (FLAC3D permeability tensor) of the medium, $$\hat{k}(s)$$ is the relative mobility coefficient, which is a function of fluid saturation, $$s$$, $${\rho}_f$$ is the fluid density, and $$g_i$$, $$i$$ = 1,3 are the three components of the gravity vector.

For future reference, the quantity $$\phi = (p - {\rho}_f x_j g_j) / ({\rho}_f g)$$ (where $$g$$ is the modulus of the gravity vector) is defined as the head, and $${\rho}_f g \phi$$ as the pressure head.

Heat flow by conduction is described by Fourier’s law of heat transport:

(7)$q_i^T = - k^T_{ij} T_{,j}$

where $$q_i^T$$ is the heat flux vector, $$T$$ is temperature and $$k^T$$ is the thermal conductivity tensor. Note that temperature, flux, convective and adiabatic boundary conditions are considered in the FLAC3D formulation.

Balance Laws

For small deformations, the fluid mass balance may be expressed as

(8)$-q_{i,i} + q_v = {{\partial\zeta}\over{\partial t}}$

where $$q_v$$ is the volumetric fluid source intensity in [1/sec], and $$\zeta$$ is the variation of fluid content or variation of fluid volume per unit volume of porous material due to diffusive fluid mass transport, as introduced by Biot (1956).

The thermal energy balance is expressed as

(9)$-q_{i,i}^T + q_v^T = {{\partial \zeta^T} \over {\partial t}}$

where $$\zeta^T$$ is the heat stored per unit volume of porous material, and $$q_v^T$$ is the volumetric heat source intensity.

The balance of momentum has the form

(10)$\sigma_{ij,j} + \rho g_i = \rho {{d v_i} \over{d t}}$

where $$\rho = \rho_d + n s \rho_w$$ is the bulk density, and $$\rho_d$$ and $$\rho_w$$ are the densities of the dry matrix and the fluid, respectively, $$n$$ is porosity, and $$s$$ is saturation.

Constitutive Laws

Changes in the variation of fluid content are related to changes in pore pressure, $$p$$, saturation, $$s$$, mechanical volumetric strains, $${\epsilon}$$, and temperature, $$T$$. The response equation for the pore fluid is formulated as

(11)${1 \over M} {{\partial p} \over {\partial t}} + {n \over s} {{\partial s} \over {\partial t}}= {1 \over s} {{\partial \zeta} \over {\partial t}} - \alpha {{\partial \epsilon} \over {\partial t}} + \beta {{\partial T} \over {\partial t}}$

where $$M$$ is Biot modulus [N/m2], $$n$$ is the porosity, $$\alpha$$ is Biot coefficient and $$\beta$$ is the undrained thermal coefficient [1/°C], which takes into account the fluid and grain thermal expansions.

In the FLAC3D formulation, the influence of capillary pressure is neglected and the fluid pressure is taken as zero when saturation is less than one. The relative fluid mobility is related to saturation by a cubic law, which is zero for zero saturation and one for full saturation:

(12)$\hat{k}(s) = s^2(3 - 2s)$

Fluid flow in the unsaturated zone is thus solely driven by gravity. While the influence of gravity is not required to saturate an initially dry medium, gravity drives the process of desaturation. In this case, some level of residual saturation is present because the apparent permeability, $$k \ \hat{k}(s)$$, goes to zero as the saturation approaches zero.

The thermal constitutive law is expressed as

(13)${1 \over M^T} {{\partial T} \over {\partial t}} = {{\partial \zeta^T} \over {\partial t}}$

where $$M^T = 1/(\rho C_v)$$, $$\rho$$ is the mass density of the medium, and $$C_v$$ is the specific heat at constant volume.

Note that for nearly all solids and liquids, the specific heats at constant pressure and constant volume are essentially equal. Consequently, $$C_v$$ and $$C_{\rho}$$ can be used interchangeably.

The constitutive response for the porous solid has the form (see notation conventions of Theoretical Background)

(14)$\check \sigma_{ij} + \alpha {{\partial p} \over {\partial t}} \delta_{ij} = H(\sigma_{ij}, \xi_{ij} - \xi_{ij}^T, \kappa)$

where the quantity on the left-hand side of the equation is Biot effective stress, $$\check\sigma_{ij}$$ is the co-rotational stress rate, $$H$$ is the functional form of the constitutive law, $$\kappa$$ is a history parameter, $$\delta_{ij}$$ is the Kronecker delta, and $$\xi_{ij}$$ is the strain rate.

The thermal strain rates are expressed as (see Thermal Strains)

(15)$\xi_{ij}^T = {\alpha}_t {{\partial T} \over {\partial t}} \delta_{ij}$

where $$\alpha_t$$ is the coefficient of linear thermal expansion.

In particular, the elastic relations that relate effective stresses to strains are (small strain)

(16)$\sigma_{ij} - \sigma_{ij}^o + \alpha (p - p^o) \delta_{ij} = 2G (\epsilon_{ij} - \epsilon_{ij}^T) + (K - {2 \over 3} G) (\epsilon_{kk} - \epsilon_{kk}^T)$

where the superscript $$^o$$ refers to the initial state, $$\epsilon_{ij}$$ is the strain, and $$K$$ and $$G$$ are the bulk and shear moduli of the drained elastic solid.

Biot and Terzaghi Effective Stress

The expression for Biot effective stress, $$\sigma_{ij} + \alpha p \delta_{ij}$$, may be derived from compliance principles (see Effect of the Biot Coefficient); as such, it characterizes the deformability of the solid. On the other hand, Terzaghi effective stress, $$\sigma_{ij} + p \delta_{ij}$$, measures the stress level sustained by the solid matrix; it is still being used to detect failure in a plastic material (or potential failure in an elastic solid).

Compatibility Equations

The relation between strain rate and velocity gradient is

(17)$\xi_{ij} = {1 \over 2} \left[v_{i,j} + v_{j,i} \right]$

## Fluid Flow Boundary and Initial Conditions in FLAC3D

Initial conditions correspond to known pressure and saturation. Those two values must be consistent with the FLAC3D formulation: pore pressure must be zero if saturation is less than one, and vice-versa. (By default, in FLAC3D, initial pore pressure is zero and saturation is one.)

Four types of boundary conditions are considered, corresponding to (1) given pore pressure, (2) given component of the specific discharge normal to the boundary, (3) leaky boundaries, and (4) impermeable boundaries.

In FLAC3D, boundaries are impermeable by default.

A leaky boundary condition has the form

(18)$q_n = h (p - p_e)$

where $$q_n$$ is the component of the specific discharge normal to the boundary in the direction of the exterior normal, $$h$$ is the leakage coefficient in [m3/N-sec], $$p$$ is the pore pressure at the boundary surface, and $$p_e$$ is the pore pressure in the region to or from which leakage is assumed to occur.

Additional combined boundary conditions, such as uniform pressure for given total influx over a boundary segment, can be imposed using FISH.

Note that saturation cannot be imposed as a boundary condition (this is different from FLAC). Also, to keep the model fully saturated, the fluid tension may be set to a large negative number.