Fracture flow using a MultiApp approach: Equations

Background

PorousFlow can be used to simulate flow through fractured porous media in the case when the fracture network is so complicated it cannot be incorporated into the porous-media mesh. The fundamental premise is the fractures can be considered as lower dimensional entities within the higher-dimensional porous media. This page is part of a set that describes a MOOSE MultiApp approach to simulating such models:

The heat equation in fracture-matrix systems

Uncoupled equations

The 3D heat equation is

(1) Here:

  • is the temperature, with SI units K;

  • is the volumetric heat capacity of the medium (which is the product of specific heat capacity and density), with SI units J.m.K;

  • is time and is the vector of spatial derivatives;

  • is the thermal conductivity, with SI units J.s.m.K.

This equation (multiplied by the test functions) gets integrated over the elemental volume in the finite-element scheme. The fracture is actually a 3D object, but it is so thin that temperature (and and ) do not vary appreciably over its thickness (aperture). Therefore, the integral over the thickness may be done explicitly, and the fracture modelled as a 2D object, with heat equation reading

(2)

where is the fracture thickness (aperture), and are derivatives transverse to the fracture. This equation is integrated (by MOOSE) over the fracture-element 2D area in the finite-element method.

commentnote

The preceding equation implies all the Kernels for heat and mass flow in the fracture system must be multiplied by the fracture aperture . That is, in the MOOSE input file governing the 2D physics, the coefficient of the CoefTimeDerivative Kernel is , and the coefficient of the AnisotropicDiffusion is . If tensorial quantities such as or the permeability tensor are anisotropic, they need to be expressed in 3D space.

Coupling

In order to specify the coupling between the fracture and matrix system precisely, let us introduce some notation. Suppose the position of the fracture is given by a level-set function, (where , and are the Cartesian coordinates): that is, the fracture is at points where , but not at points where . For example, if the fracture occupies the plane, then would be a suitable function. The one-dimensional Dirac delta function, , has SI units m, and is zero everywhere except on the fracture. When is integrated over a direction normal to the fracture, the result is . When is integrated over a volume containing the fracture, the result is , where is the area of the fracture in the volume.

The approach explored in this page uses two temperature variables, and , which are the temperature in the matrix and fracture, respectively. is defined throughout the entire matrix (including the embedded fracture), while is defined on the fracture only. These are assumed to obey:

(3)

In these equations, is the heat-transfer coefficient between the matrix and the fracture, with SI units J.m.s.K. If then the heat-transfer term takes heat energy from the system and applies it to the system. The heat-transfer coefficient plays a central role in this page and is discussed in detail below.

When the second (fracture) equation is integrated over a (2D) portion of the fracture of area , the heat-energy transferred to the matrix is , where the indicates the average over . This is equal to the heat-energy transferred according to the first equation, when it is integrated over a volume containing the same portion of fracture (remember ). Therefore, heat energy is conserved in this system.

The heat-transfer coefficient

Heat-transfer coefficients have been used by engineers to accurately account for the complicated processes that occur at the interface between two materials. If the two materials have temperature and , respectively, and are connected by an area then the rate of heat transfer between them is , with SI units J.s. Heat-transfer coefficients have been estimated for many different situations (see, for instance Wikipedia).

To motivate a numerical value for , assume that the fracture width is tiny compared with any relevant length scales in the matrix (such as the finite-element sizes). While the temperature distribution in the immediate vicinity of the fracture will be governed by complicated processes, the point of is to hide that complexity. Assume there is a "skin" around the fracture in which the complicated processes occur. Denote the temperature just outside this skin by , so the heat transfer through the skin is

(4) where is the fracture temperature, is the heat-transfer coefficient of the skin, with SI units J.m.s.K, and is the heat flux, with SI units J.s.m. A discussion on the numerical value of is found below.

Just as the fracture may be considered two-dimensional (implemented by including its aperture in the fracture Kernels and using 2D finite elements), the size of the skin is also ignored. Thereby, the matrix "sees" the fracture as an object of (not ) sitting within it. Consider Figure 1, where an object of temperature sits between two finite element nodes, of temperature and .

Figure 1: An object of temperature sits a distance from a finite-element node of temperature , and a distance from a node of temperature .

The flow from the fictitious "skin" object to each matrix finite-element nodes is governed by the heat equation. In the linear limit, Eq. (1) implies

(5)

Assuming the sum of these equals , that is : Using this equality allows to be written in terms of the other temperature values, and substituting the expression into Eq. (4) yields the heat-loss from the fracture: Notice that the prefactors are exactly what a linear-lagrange element in MOOSE would prescribe to a DiracKernel sitting at the fracture position. Hence, the effective heat-transfer coefficient that would be used in a MOOSE model of this situation is

commentnote

Notice a key assumptions made in the presentation above — the linear approximation of Eq. (5), where it is assumed that the matrix element sizes are small enough to resolve the physics of interest. This means that phenomena associated with short-time, small-scale physics won't be resolvable in models with large elements, when using the heat-transfer coefficients recommended in this page. This is because that in various real-life scenarios there may be no heat flow between the "0" position and the skin even if (contradicting Eq. (5) ) because, for instance, changes of take some time to be felt by the "0" position. The solution is to use smaller elements to resolve this short-time, small-scale phenomena. The small fracture network model provides an example.

The same analysis may be performed in the 2D-3D situation. While Figure 1 is only a 1D picture, flow from the fracture to the matrix only occurs in the normal direction, so it also represents the 2D-3D situation. With an arbitrary-shaped element containing a portion of a fracture, the heat flows to each node, and hence , depend on the shape of the element. However, following the approach of the Peaceman borehole Peaceman borehole, the following rule-of-thumb is suggested for MOOSE simulations:

(6)

where depends on the matrix element surrounding the fracture:

  • is the component of the matrix thermal conductivity in the normal direction to the fracture;

  • is the average of the distances between the fracture and the nodes that lie on its right side;

  • is the average of the distances between the fracture and the nodes that lie on its left side (opposite the right side).

The result of Eq. (6) depends on and . In most settings, it is appropriate to assume that since this corresponds to making a shift of the fracture position by an amount less than the finite-element size. Since the accuracy of the finite-element scheme is governed by the element size, such small shifts introduce errors that are smaller than the finite-element error. This means Eq. (6) reads

(7)

The heat-transfer through the skin is likely to exceed heat conduction through the skin, since the heat transfer involves faster processes such as convection. Therefore, , where is the skin size and the numerical quantity . This yields

(8)

where the final limit is for , which is likely to be reasonably correct in most simulations.

Finally, consider the case which has very large fracture elements compared with matrix elements. Then it could happen that , where is now the fracture area modelled by one finite-element fracture node, and is the volume modelled by one finite-element matrix node. Then the single fracture node can apply a lot of heat to the "small" matrix node, which will likely cause oscillations or even numerical instability if is too large, in the MultiApp approach.

Mass flow

A "mass-transfer coefficient" may be used to model fluid mass transfer between the fracture and matrix. The equations for all fluid phases are structurally identical, but just have different numerical values for viscosity, relative permeability, etc, so in this section consider just one phase. Define the potential (with SI units Pa) by

(9)

Here

  • is the porepressure (SI units Pa);

  • is the fluid density (SI units kg.m);

  • is the gravitational acceleration (SI units m.s or Pa.m.kg);

  • is the coordinate direction pointing "upwards" (i.e. in the opposite direction to gravity) with SI units m. If gravity acted in a different direction then would be replaced by another coordinate direction.

In the Darcy setting modelled by PorousFlow, the equivalent of the heat transfer found in Eq. (3) is

and the equivalent of Eq. (4) is

Here

  • is the mass flux through the skin, with SI units kg.s.m.

  • is the mass-transfer coefficient of the skin, with SI units kg.s.m.Pa.

The heat-flow arguments presented above may now be followed using instead of , and in place of , to yield a suggestion for the mass-transfer coefficient in the 2D-3D situation:

(10)

where depends on the matrix element surrounding the fracture in the following way

  • is the density of the fluid phase (SI units kg.m).

  • is the component of the matrix permeability tensor in the normal direction to the fracture (SI units m).

  • is the relative permeability of the fluid phase (dimensionless).

  • is the dynamic viscosity of the fluid phase (SI units Pa.s).

  • is the average of the distances between the fracture and the nodes that lie on its right side.

  • is the average of the distances between the fracture and the nodes that lie on its left side (opposite the right side).

This appears in the mass-transfer: , where is given in Eq. (9). The arguments that led to Eq. (7) and Eq. (8) lead in this case to

(11)

There are numerical subtleties in these expressions that may not impact many simulations, but modellers should be aware of the following.

  • seen by the fracture is the matrix nodal values interpolated to the fracture node. There are clearly many ways of doing this interpolation, eg, could be interpolated from the fracture nodes, or could be evaluated at the fracture node given its position, porepressure and temperature.

  • Similar remarks hold for .

  • Stability of the numerics will be improved by evaluating at the upwind position. That is, if fluid is flowing from the fracture to the matrix, this term would best be evaluated using of the skin, while for the reverse flow the evaluation should be done using the nodal matrix values. However, it is usually more convenient to fix the evaluation to use the elemental-averaged values of the matrix (using a constant, monomial AuxVariable for the matrix).

Heat advection

Heat advection follows the same principle as mass transfer. The heat transfer by advection through the skin:

where is the fluid enthalpy, and the other terms are given above. In the finite-element setting , and

(12)

The arguments that led to Eq. (7) and Eq. (8) lead in this case to

(13)

As in the mass-transfer case, the numerics will be impacted by where these terms are evaluated.

Combinations of transfer

In many simulations, heat transfer by both conduction and convection will be active, in addition to mass transfer. The transfers are simply added together. For instance