Mathematical Model Description

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\). (\(f\) can be 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 this section for conversions to other systems of units).

The following dimensionless numbers are useful in the characterization of transient heat conduction:

Characteristic length

(1)\[L_c = {{volume \ of \ solid}\over{surface \ area \ exchanging \ heat}}\]

Thermal diffusivity

(2)\[\kappa = {{k}\over{\rho\ C_v}}\]

where:

\(k\)

is the thermal conductivity;

\(\rho\)

is the density; and

\(C_v\)

is the specific heat at constant volume.

Characteristic time

(3)\[t_c = {{L_c^2}\over{\kappa}}\]

Conduction

The variables involved in heat conduction in FLAC3D are the temperature and the three components of the heat flux. These variables are related through the energy-balance equation and transport laws derived from Fourier’s law of heat conduction. Substitution of Fourier’s law into the energy-balance equation yields the differential equation of heat conduction, which may be solved for particular geometries and properties, given specific boundary and initial conditions. Thermal volumetric-strains are introduced into the incremental mechanical and fluid constitutive laws to account for thermal-stress and thermal-pore pressure couplings.

Energy-Balance Equation

The differential expression of the energy balance has the form

(4)\[-q_{i,i} + q_v = {{\partial\zeta}\over{\partial t}}\]

where \(q_i\) is the heat-flux vector in \([W/m^2]\), \(q_v\) is the volumetric heat-source intensity in \([W/m^3]\), and \(\zeta\) is the heat stored per unit volume in \([J/m^3]\). In general, temperature changes may be caused by changes in both energy storage and volumetric strain, \(\epsilon\), and the thermal constitutive law relating those parameters can be expressed as

(5)\[{{\partial T}\over{\partial t}} = M_{th} \left( {{\partial \zeta}\over{\partial t}} - \beta_{th} {{\partial \epsilon}\over{\partial t}} \right)\]

where \(M_{th}\) and \(\beta_{th}\) are material constants, and \(T\) is temperature.

In FLAC3D, we consider a particular case of this law for which \(\beta_{th}\) = 0 and \(M_{th} = 1/(\rho C_v)\). \(\rho\) is the mass density of the medium in \([kg/m^3]\), and \(C_v\) is the specific heat at constant volume in \([J/kg\ ^{\circ}C]\). The hypothesis here is that strain changes play a negligible role in influencing the temperature, a valid assumption for quasi-static mechanical problems involving solids and liquids. Accordingly,

(6)\[{{\partial \zeta}\over{\partial t}} = \rho\ C_v {{\partial T}\over{\partial t}}\]

Substitution of Equation (6) in Equation (4) yields the energy-balance equation

(7)\[-q_{i,i} +q_v = \rho C_v {{\partial T}\over{\partial t}}\]

Note that, for nearly all solids and liquids, the specific heats at constant pressure and at constant volume are essentially equal. Consequently, \(C_v\) and \(C_p\) can be used interchangeably.

Transport Law

The basic law that defines the relation between the heat-flux vector and the temperature gradient is Fourier’s law. For a stationary, homogeneous, isotropic solid, this constitutive law is given in the form

(8)\[q_i = - k T_{,i}\]

where \(T\) is the temperature [\(^{\circ}\)C], and \(k\) is the thermal conductivity in \([W/m ^{\circ}C]\).

Isotropic Heat Conduct

For isotropic conduct, the thermal conductivity matrix has the form:

(9)\[\begin{split}\boldsymbol{k} = k \cdot \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\end{split}\]

For this case, one conductivity scalar parameter is sufficient.

An alphabetical list of properties of the isotropic heat conduction thermal model is given here.

Anisotropic Heat Conduct

For anisotropic conduct, the thermal conductivity matrix has the form:

(10)\[\begin{split}\boldsymbol{k} = \begin{bmatrix} k_{xx} & 0 & 0 \\ 0 & k_{yy} & 0 \\ 0 & 0 & k_{zz} \end{bmatrix}\end{split}\]

For this case, three conductivity parameters corresponding to different directions are required.

In general,

(11)\[\begin{split}\boldsymbol{k} = \begin{bmatrix} k_{xx} & k_{xy} & k_{xz} \\ k_{xy} & k_{yy} & k_{yz} \\ k_{xz} & k_{yz} & k_{zz} \end{bmatrix}\end{split}\]

For this case, six conductivity parameters are required.

An alphabetical list of properties of the anisotropic heat conduction thermal model is given here.

Boundary and Initial Conditions

Substitution of Equation (8) for \(q_i\) in Equation (7) yields the differential equation for heat conduction. Initial conditions correspond to a given temperature field. Boundary conditions are generally expressed in terms of temperature or the component of the heat-flux vector normal to the boundary. In this version of FLAC3D, four types of conditions are considered, corresponding to: (1) given temperature; (2) given component of the flux normal to the boundary; (3) convective boundaries; and (4) insulated (adiabatic) boundaries.

In FLAC3D, a convective boundary condition has the form

(12)\[q_n = h (T - T_e)\]

where \(q_n\) is the component of the flux normal to the boundary in the direction of the exterior normal, \(h\) is the convective heat-transfer coefficient \([W/m^2$ $^{\circ}C]\), \(T\) is the temperature of the boundary surface, and \(T_e\) is the temperature of the surrounding fluid [\(^{\circ}\)C]. Note that in the numerical formulation used in FLAC3D, boundaries are adiabatic by default.

Mechanical Coupling — Thermal Strains

Solution of thermal-stress problems requires reformulation of the incremental stress-strain relations, which is accomplished by subtracting the portion due to temperature change from the total strain increment. Because free thermal expansion results in no angular distortion in an isotropic material, the shearing-strain increments are unaffected. The thermal-strain increments associated with the free expansion corresponding to temperature increment \(\Delta T\) have the form

(13)\[\Delta \epsilon_{ij} = \alpha_t \Delta T \delta_{ij}\]

where \(\alpha_t\ [1/^{\circ}C]\) is the coefficient of linear thermal expansion, and \(\delta_{ij}\) is the Kronecker delta.

The differential equation of motion and the rate of strain-velocity relations are based upon mechanical considerations, and are unchanged for thermomechanics.

Fluid Coupling — Thermally Induced Pore Pressures

The heat transfer may be coupled to the groundwater calculation by taking into account pore-pressure change caused by thermal expansion of the saturated matrix. The constitutive law has the form (see this this equation in FLAC3D Fluid-Thermal-Mechanical-Formulation — Mathematical Description)

(14)\[{{\partial p}\over{\partial t}} = M \left( {{\partial \zeta}\over{\partial t}} - \alpha {{\partial \epsilon}\over{\partial t}} + \beta {{\partial T}\over{\partial t}} \right)\]

where \(\zeta\) is the variation of fluid content, \(M\) is the Biot modulus, \(\alpha\) is the Biot coefficient, \(\epsilon\) is the volumetric strain, and \(\beta\) is the undrained volumetric thermal expansion. In FLAC3D, the pore-pressure increment \(\Delta p_{ther}\), caused by a temperature change \(\Delta T\), has the form

(15)\[\Delta p_{ther} =\ M \beta \Delta T\]

\(\beta\ [1/^{\circ}C]\) is defined as the pore-pressure variation divided by (\(3 \alpha_t\ M\)) per unit temperature change in an undrained constrained test (no deformation). For an ideal porous material, this coefficient is related to the volumetric thermal-expansion coefficients for the grains, \(3 \alpha_t\), and the fluid, \(3 \alpha_f\), through the formula

(16)\[\beta = 3[\alpha_t(\alpha - n) + \alpha_f n]\]

where \(n\) = porosity.

For an incompressible solid constituent, (\(\alpha\) = 1 and \(M = {K_f} / {n}\)), the pore-pressure increment becomes

(17)\[\Delta p_{ther} =\ {K_f \over n} \beta \Delta T\]

where \(K_f\) = fluid bulk modulus.

The corresponding change in total stress is taken into consideration in the constitutive laws through the relation

(18)\[\sigma_{ij} =\ - \Delta p_{ther} \delta_{ij}\]

Advection

The mechanisms of convective heat transfer in porous media considered in the FLAC3D implementation are forced convection, in which the heat is carried by the fluid motion, and free convection, which accounts for fluid motion caused by density differences due to temperature variations.

The working assumptions are those of local thermal equilibrium (flow at low Reynolds number is considered) and small density variations (Boussinesq approximation applies). The basic hydrothermal equations for saturated fluid flow processes are provided below.

Energy Balance for Convective-Diffusive Heat Transport

(19)\[c^T {{\partial T} \over {\partial t}} + \nabla \cdot {\bf q}^T + \rho_0 c_w {\bf q}_w \cdot \nabla T - q^T_v = 0\]

where \(T\) is temperature, \({\bf q}^T\) is thermal flux, \({\bf q}_w\) is fluid specific discharge, \(q^T_v\) is volumetric heat source intensity, \(\rho_0\) and \(c_w\) are reference density and specific heat of the fluid, respectively, and \(c^T\) is the effective specific heat. The effective specific heat is defined as

(20)\[c^T = \rho_d C_v + nS \rho_0 c_w\]

where \(\rho_d\) and \(C_v\) are solid matrix bulk density and bulk specific heat, respectively, \(n\) is porosity, and \(S\) is saturation.

Fluid Mass Balance (Slightly Compressible Fluid)

(21)\[{{\partial P} \over {\partial t}} = M \biggl(- \nabla \cdot {\bf q}_w - \alpha {{\partial \epsilon} \over {\partial t}} + \beta {{\partial T} \over {\partial t}} \biggr)\]

where \(M\) is the Biot modulus, \(\epsilon\) is the volumetric strain, and \(\beta\) is the volumetric thermal expansion of the porous matrix. Note that the term \(n \beta_f\) (\(\beta_f\) is volumetric thermal expansion of the fluid), which enters in the expression of the coefficient \(\beta\), is usually neglected in the Boussinesq approximation, and that the last two terms in Equation (21) are only present if mechanical coupling is considered.

Transport Laws

Heat Transport (Fourier Law)

(22)\[{\bf q}^T = - k^T \nabla T\]

where \(k^T\) is the effective thermal conductivity, which is isotropic in the advection formulation. The effective thermal conductivity is defined in terms of the solid and fluid thermal conductivities, \(k^T_s\) and \(k^T_w\), by the equation

(23)\[k^T = k^T_s + n S k^T_w*\]

Fluid Transport (Darcy Law)

(24)\[{\bf q}_w = -{\bf k} \nabla (P - \rho_w {\bf g} \cdot {\bf x})\]

where \({\bf k}\) is the fluid mobility coefficient (intrinsic permeability, \(\boldsymbol{\kappa}\), over dynamic viscosity, \(\mu_w\)), and \(\rho_w\) is the fluid density. In this equation, fluid density is related to temperature changes by the linear equation

(25)\[\rho_w = \rho_0 [1 - \beta_f (T - T_0)]\]

where \(\beta_f\) is the volumetric thermal expansion of the fluid, and \(T_0\) is the reference temperature.

In these equations, the fluid density variations due to temperature changes are significant only in their generation of buoyancy forces (i.e., the Boussinesq approximation).

Thermal-Mechanical-Pore Pressure Coupling

Coupling to the mechanical constitutive equations is done via “effective” normal strain rates (i.e., total-minus-thermal normal strain-rate components), and pore pressure changes, which affect effective stresses. The advective logic may be used in combination with any of the mechanical constitutive models available for FLAC3D. In the example of an elastic material in small strain, the mechanical constitutive relations are

(26)\[{{\partial \sigma_{ij}} \over {\partial t}} + \alpha {{\partial P} \over {\partial t}} \delta_{ij} = 2G \biggl({{\partial \epsilon_{ij}} \over {\partial t}} - \alpha_t {{\partial T} \over {\partial t}} \delta_{ij} \biggr) + \biggl(K - {2 \over 3} G \biggr) \cdot \biggl({{\partial \epsilon_{kk}} \over {\partial t}} - 3 \alpha_t {{\partial T} \over {\partial t}} \biggr) \delta_{ij}\]

where \(\sigma_{ij}\) and \(\epsilon_{ij}\) are total stresses and strains, \(\alpha\) is the Biot coefficient, \(K\) and \(G\) are bulk and shear moduli, \(\alpha_t\) is linear thermal-expansion coefficient of the solid matrix, and \(\delta_{ij}\) is the Kronecker delta.

Initial and Boundary Conditions

In a coupled simulation, initial and boundary conditions must be provided for each of the processes involved. In principle, any combination of the fluid and thermal boundary conditions documented in the FLAC3D manual can be used in an advection simulation: Dirichlet and Neuman boundary conditions can be assigned, and point and volume sources can be specified for both fluid and thermal calculations.

An alphabetical list of properties of the advection-conduction thermal model is given here.