Numerical Formulation

Conduction

The energy-balance equation, together with Fourier’s law, is solved in FLAC3D using a finite-difference approach based on a medium discretization into zones composed of two overlays of constant heat-flux tetrahedra. The numerical scheme rests on a nodal formulation of the energy-balance equation. This formulation is equivalent to the mechanical formulation (presented in Theoretical Background that leads to the nodal form of Newton’s law. It is obtained by substituting the temperature, heat-flux vector and temperature gradient for velocity vector, stress tensor and strain-rate tensors, respectively. The resulting system of ordinary differential equations is solved using two modes of discretization in time, corresponding to explicit and implicit formulations. The principal results are summarized below.

Finite-Difference Approximation to Space Derivatives

By convention, the nodes of a tetrahedron are referred to locally by a number from 1 to 4, face \(n\) is opposite node ;math:\(n\), and the superscript \((f)\) relates to the value of the associated variable on face \(f\).

A linear temperature variation is assumed within a tetrahedron; the temperature gradient, expressed in terms of nodal values of the temperature by application of the Gauss divergence theorem, can be written as

(1)\[T_{,j} = - {{1}\over{3V}} \sum_{l=1}^4 T^l\ n_j^{(l)}\ S^{(l)}\]

where \([n]^{(l)}\) is the exterior unit vector normal to face \(l\), \(S\) is the face surface area, and \(V\) is the tetrahedron volume.

Nodal Formulation of the Energy-Balance Equation

The energy-balance equation can be expressed as

(2)\[q_{i,i} + b^* = 0\]

where

(3)\[b^* = \rho C_v {{\partial T}\over{\partial t}} - q_v\]

is the equivalent of the instantaneous “body forces” used in the mechanical node formulation. First consider a single tetrahedron. Using this analogy, the nodal heat \(Q_e^n [W]\), \(n\) = 1,4, in equilibrium with the tetrahedron heat flux and body forces, may be expressed as

(4)\[Q_e^n = Q_t^n - {{q_v V}\over{4}} + m^n\ C_v^n {{dT}\over{dt}}^n\]

where

(5)\[Q_t^n = {{q_i n_i^{(n)} S^{(n)}}\over 3}\]

and

(6)\[m^n = {{\rho V}\over{4}}\]

In principle, the nodal form of the energy-balance equation is established by requiring that, at each global node, the sum of equivalent nodal heat (\(-Q_e^n\)) from all tetrahedra meeting at the node (divided by two for two overlays) and nodal contribution (\(Q_w^n\)) of applied boundary fluxes and sources be zero.

The components of the tetrahedron heat-flux vector in Equation (5) are related to the temperature gradient by means of the transport law (see Fourier’s law). In turn, the components of the temperature gradient can be expressed in terms of the tetrahedron’s nodal temperatures using Equation (1) .

In order to save computer time, a zone formulation is adopted in FLAC3D. At each node, \(n\), of a particular zone, the sum, \(Q_z^n\), of contributions (Equation (5)) from all tetrahedra belonging to the zone and meeting at the node is formed and divided by two for two overlays. Local zone matrices \([M]\) that relate nodal zone values, \(Q_z^n\), to nodal zone temperatures, \(T^n\) (\(n\) = 1,8) are assembled. Because the matrices are symmetrical, 36 components are computed. These components are saved at the beginning of the computation and updated every 10 steps in large-strain mode. By definition of zone matrices, we have

(7)\[Q_z^n = M_{nj} T^j\]

where \([T]\) is the local vector of nodal zone temperatures.

In turn, global nodal values \(Q_T^n\) are obtained by superposition of zone contributions. Taking some liberty with the notation, we write for each global node \(n\),

(8)\[Q_T^n = C_{nj} T^j\]

where \([C]\) is the global matrix, and \([T]\) is, in this context, the global vector of nodal temperatures.

Again considering a single tetrahedron, we write

(9)\[- \sum Q_e^n + \sum Q_w^n = 0\]

where, for simplicity of notation, the \(\sum\) sign is used to represent summation of the contributions at global node \(n\) of all zones meeting at that node. (A zone contribution consists of contributions of all tetrahedra involved in the zone, divided by two for two overlays.) Using Equation (4) in Equation (9), we obtain, after some manipulations,

(10)\[{{d T^n}\over{d t}} = - {{1} \over {\sum[m C_p]^n}} \left[Q_T^n + \sum Q_{app}^n \right]\]

where \(Q_T^n\) is a function of the nodal temperatures defined in Equation (8), and \(\sum Q_{app}^n\) is the known contribution of applied volume sources, boundary fluxes and point sources, defined as

(11)\[\sum Q_{app}^n = - \sum {\left[{{q_v V}\over{4}} + Q_w\right]}^n\]

Equation (10) is the nodal form of the energy-balance equation at node \(n\); the right-hand side term, \(Q_T^n + \sum Q_{app}^n\), is referred to as the “out-of-balance heat.” One such expression holds at each global node involved in the discretization. Together they form a system of ordinary differential equations that is solved in FLAC3D using both explicit and implicit finite-difference schemes. The domain of application of each scheme is discussed below.

Explicit Finite-Difference Formulation

In the explicit formulation, the temperature at a node is assumed to vary linearly over a time interval, \(\Delta t\): the derivative in the left-hand side of Equation (10) is expressed using forward finite differences; and the out-of-balance heat is evaluated at time \(t\). Starting with an initial temperature field, nodal temperatures at incremental time values are updated, provided the temperature is not fixed, using the expression

(12)\[T_{<t+\Delta t>}^n = T_{<t>}^n + \Delta T_{<t>}^n\]

where

(13)\[\Delta T_{<t>}^n = {\chi}^n \left[ Q_{T <t>}^n + \sum Q_{app<t>}^n \right]\]

and

(14)\[{\chi}^n = - {{\Delta t}\over{\sum [m C_p]^n}}\]

Numerical stability of the explicit scheme can only be achieved if the timestep remains below a limiting value.

Stability Criterion

To derive the stability criterion, we consider the situation in which a node, \(n\), in an assembly of zones is given a temperature perturbation, \(T_0\), from a zero initial state. Using Equation (8), we obtain

(15)\[Q_T^n = C_{nn}\ T_0\]

If node \(n\) belongs to a convective boundary, we have

(16)\[\sum Q^n_{app} = D_{nn}\ T_o\]

where \(D_{nn}\) is used to represent the temperature coefficient in the global convective term at node \(n\).

After one thermal timestep, the new temperature at node \(n\) is (see Equations (12) through (14))

(17)\[T_{<\Delta t>}^n = T_0 \left[ 1 - {{\Delta t}\over{\sum [m C_p]^n}} (C_{nn} + D_{nn}) \right]\]

To prevent alternating signs of temperatures as the computation is repeated for successive \(\Delta t\), the coefficient of \(T_0\) in the preceding relation must be positive. Such a requirement implies that

(18)\[\Delta t < {{\sum [m C_p]^n}\over{C_{nn} + D_{nn}}}\]

The value of the timestep used in FLAC3D is the minimum nodal value of the right-hand side in Equation (18), multiplied by the safety factor of 0.8.

To assess the influence of the parameters involved, it is useful to keep in mind the following representation of the critical timestep. If \(L_c\) is the smallest tetrahedron characteristic length, we may write an expression of the form

(19)\[\Delta t_{cr} = {{1}\over{m}} {\left [ {{\kappa}\over{L_c^2}} + {{h}\over{\rho\ C_v L_c}} \right ]}^{-1}\]

where \(\Delta t_{cr}\) is the critical timestep, \(\kappa\), \(h\), \(\rho\) and \(C_v\) are defined previously in Mathematical Model Description, and \(m\) is a constant, larger than unity, which depends on the medium geometrical discretization (e.g., \(m\) = 6 for a regular discretization in cubes — see Karlekar and Desmond 1982).

The critical timestep in Equation (19) corresponds to a measure of the characteristic time needed for the diffusion “front” to propagate throughout the tetrahedron. To estimate the time needed for the front to propagate throughout a particular domain, a similar expression can be used, provided \(L_c\) is interpreted as the characteristic length of the domain under consideration. Consider the case where no convection occurs. By taking the ratio of characteristic times for the domain and the tetrahedron, we see that the number of steps needed to model the propagation of the diffusion process throughout that domain is proportional to the ratio of square characteristic lengths for the domain and the tetrahedron. That number may be so large that the use of the explicit method alone could become prohibitive.

On the other hand, the advantage of this first-order method is that the calculated timestep is usually small enough to follow nodal temperature fluctuations accurately.

Implicit Finite-Difference Formulation

The requirement that \(\Delta t\) should be restricted in size to ensure stability sometimes results in an extremely small timestep (of the order of a fraction of a second), especially when transient conduction in multiple layers is involved. The implicit formulation eliminates this restriction, but it involves solving simultaneous equations at each timestep.

The implicit formulation in FLAC3D uses the Crank-Nicolson method, in which the temperature at a node is assumed to vary quadratically over the time interval \(\Delta t\). In this method, the derivative \(dT/dt\) in Equation (10) is expressed using a central-difference formulation corresponding to a half-timestep, while the out-of-balance heat is evaluated by taking average values at \(t\) and \(t + \Delta t\). In this formulation, we have

(20)\[T^n_{<t + \Delta t>} = T_{<t>}^n + \Delta T_{<t>}^n\]

where

(21)\[\Delta T_{<t>}^n = {\chi}^n \left[{{1}\over{2}} \left(Q_{T <t + \Delta t>}^n + Q_{T <t>}^n \right) + \sum \overline{Q^n_{app}}\right]\]

and

(22)\[\sum \overline{Q^n_{app}} = {1 \over 2} \left( \sum Q^n_{app<t + \Delta t>} + \sum Q^n_{app<t>} \right)\]

From the definition in Equation (8), we can write

(23)\[Q_{T<t>}^n = C_{nj} T_{<t>}^j\]

and

(24)\[Q_{T<t + \Delta t>}^n = C_{nj} T_{<t + \Delta t>}^j\]

Using (12), this last expression becomes

(25)\[Q_{T<t + \Delta t>}^n = Q_{T<t>}^n + C_{nj} \Delta T_{<t>}^j\]

After substitution of Equation (25) into Equation (21), we obtain (using Equation (23))

(26)\[\Delta T_{<t>}^n = {\chi}^n \left[C_{nj} T_{<t>}^j + {{1}\over{2}} C_{nj} \Delta T_{<t>}^j + \sum \overline{Q_{app}^n}\right]\]

Finally, regrouping terms,

(27)\[\left[\delta_{nj} - {{\chi}^n\over{2}} C_{nj}\right] \Delta T_{<t>}^j = {\chi}^n \left[C_{nj} T_{<t>}^j + \sum \overline{Q_{app}^n}\right]\]

For simplicity of notation, we define the known matrix \([A]\) and vector \([b_{<t>}]\) as

(28)\[A_{nj} = \delta_{nj} - {{\chi}^n\over{2}} C_{nj}\]

and

(29)\[b_{n<t>} = {\chi}^n \left[C_{nj} T_{<t>}^j + \sum \overline{Q_{app}^n}\right]\]

With this convention, we can write Equation (27) in the form

(30)\[A_{nj} \Delta T_{<t>}^j = b_{n<t>}\]

Such an equation can be written for each of the nodes involved in the grid. The resulting system of equations is solved in FLAC3D using Jacobi’s iterative method. In this approach, temperature increments at iteration \(r + 1\) and node \(n\) are calculated using the recurrence relation

(31)\[\Delta T_{<t>}^{n <r+1>} = \Delta T_{<t>}^{n <r>} + {{1}\over{A_{nn}}} \left[-A_{nj} \Delta T_{<t>}^{j <r>} + b_{n<t>}\right]\]

where Einstein’s notation convention applies to index \(j\) only. sing the definition by Equation (28) of \([A]\), Equation (31) takes the form

(32)\[\Delta T_{<t>}^{n <r+1>} = \Delta T_{<t>}^{n <r>} + {{1}\over{1 - {{\chi}^n\over{2}} C_{nn}}} \left[{{\chi}^n\over{2}} C_{nj} \Delta T_{<t>}^{j <r>} - \Delta T_{<t>}^{n <r>} + b_{n<t>}\right]\]

The initial approximation is chosen such that

(33)\[\Delta T_{<t>}^{n <0>} = 0\]

Hence, the first approximation \(\Delta T_{<t>}^{n <1>}\) is calculated using the same formula as that used in the explicit scheme divided by \(\left( 1 - {{\chi^n}\over 2} C_{nn}\right)\). (Note that in FLAC3D, temperature-dependent boundary conditions (contained in \(b_n\)) are updated in the implicit iterative procedure.)

In FLAC3D, a minimum of 3 and a maximum of 500 iterations are considered, and the criterion for detection of convergence has the form

(34)\[\max\limits_n \ \ \big\vert\ \ \Delta T^{n<r+1>}_{<t>} - \Delta T^{n<r>}_{<t>} \ \ \big\vert\ \ < 10^{-2}\ \left( \max\limits_n \ \ \big\vert\ \ \Delta T^{n<r>}_{<t>} \ \ \big\vert \right)\]

The magnitude of the timestep must be selected in relation to both convergence and accuracy of the implicit scheme. Although the Crank-Nicolson method is stable for all positive values of \(\Delta t\) (for no convection), the convergence of Jacobi’s method is not unconditionally guaranteed unless the matrix \([A]\) is strictly diagonally dominant:

(35)\[\sum\limits_{{j=1}\atop{j\ne k}}^n \ \ \bigg\vert\ \ A_{kj} \ \ \bigg\vert\ \ < \ \ \bigg\vert\ \ A_{kk} \ \ \bigg\vert\]

for \(1 \le i \le n\) (see Dahlquist and Bjorck 1974). According to the definition by Equation (28) of \(A_{nj}\) and Equation (14) of \(\chi^n\), this sufficient condition is always fulfilled for sufficiently small values of \(\Delta t\). If convergence of Jacobi’s method is not achieved, an error message is issued. It is then necessary to either reduce the magnitude of the implicit timestep, or use the explicit method. (The explicit timestep may be used as a lower bound.)

Also, although the implicit method is second-order accurate, some insight may be needed in order to select the appropriate timestep. Indeed, its value must remain small compared to the wavelength of any nodal temperature fluctuation. Typically, the explicit method is used earlier in the run or in its perturbed stages, while the implicit method is preferred for the remainder of the simulation. (Alternatively, the implicit method could be used with the explicit timestep value for extra accuracy.)

Computation time and computer memory are two factors that must be taken into consideration when selecting the implicit approach in FLAC3D. In the implicit method, a set of equations must be solved at each timestep, requiring a minimum of three iterations. The amount of calculation required for one iterative step is approximately equal to that needed for one timestep in the explicit scheme. Also, intermediate values must be stored in the iterative procedure, requiring extra memory to be allocated in comparison with the explicit scheme. Those inconveniences, however, can be more than offset by the much larger timestep generally permitted by the implicit method, or by the gain in accuracy allowed.

Thermal-Stress Coupling

The heat transfer may be coupled to thermal-stress calculations at any time during a transient simulation. The coupling occurs in one direction only (i.e., the temperature may result in stress changes, but mechanical changes in the body resulting from force application do not result in temperature change). This restriction is not believed to be of great significance here, since the energy changes for quasi-static mechanical problems are usually negligible. The stress change in a triangular zone is given by (from this equation):

(36)\[\Delta \sigma _{ij}\ =\ -3K\ \alpha_t\ \Delta T \delta _{ij}\]

The above assumes a constant temperature in each triangular zone that is interpolated from the surrounding gridpoints. This stress is added to the zone stress state prior to application of the constitutive law.

Mechanical properties can also be made a function of temperature change by accessing temperature and property values via FISH.

Thermal Pore-Pressure Coupling

Pore-pressure change resulting from temperature change is found by adding a contribution to the thermal-mechanical term in the calculation for nodal pore pressure (this equation in Fluid-Mechanical Interaction).

The corresponding change in total stress is calculated using the expression (from this equation, this equation, and this equation in Fluid-Mechanical Interaction):

(37)\[\Delta \sigma_{ij} =\ - \alpha \ M \beta \ \Delta T\ \delta _{ij}\]

which becomes, for incompressible solid grains (\(\alpha\) = 1 and \(M = {{K_f} \over {n}}\)),

(38)\[\Delta \sigma_{ij} =\ - {K_f \over n} \beta \ \Delta T\ \delta _{ij}\]

Groundwater properties can also be made a function of temperature change by accessing temperature and property values via FISH.

Advection

The development for convection-conduction is an extension to the explicit, finite-difference fluid and thermal-conduction implementations described earlier. It consists of adding an advection term in the thermal energy balance, as formulated in this equation, to account for forced convection by the fluid. In addition, the logic for natural convection is considered by means of the dependency (this equation) of fluid density on temperature in the transport equation (this equation).

The numerical implementation of the thermal-conduction part of the advection logic is done using the finite-difference formulation described in Conduction. The algorithm is modified to account for convection by adding a heat contribution to the out-of-balance flux at the nodes, before new temperatures are evaluated for the thermal step. The convection contribution is calculated on a zone basis. First, the heat source term \(\rho_0 c_w q_w \cdot \nabla T V\) is computed for each zone to which the convection model is assigned, at each thermal timestep. In this formula, \(q_w\) is fluid-specific discharge for the zone, \(\nabla T\) is temperature gradient, and \(V\) is zone volume. The temperature gradient is evaluated from nodal temperatures in each tetrahedra, using the Gauss divergence theorem, and an average is taken for the zone. The heat source term is then distributed to the nodes and added to their out-of-balance flux. Finally, the local fluid density used in the transport equation is updated in terms of zone temperature, according to the law in this equation, to account for natural convection.

The advection logic couples fluid and thermal calculations, and relies on sequential stepping of fluid and thermal modules. It is activated by means of the thermal constitutive model, zone thermal cmodel assign advection-conduction. The advection term is added during the thermal calculation, and requires the input of fluid-specific discharge. The discharge can be provided manually, or by the fluid flow calculation. The temperature field is used to evaluate the buoyancy term during the fluid flow calculation.

Stability and Accuracy

The stability of the numerical scheme for solving the energy balance equation (this equation) is governed by two numbers: the Peclet number for the grid, \(P_e\); and the grid Courant number, \(C_r\). For uniform specific heat, \(c^T\), the energy balance equation (this equation) can be written in the form

(39)\[{{\partial T} \over {\partial t}} = \nabla \cdot (D \nabla T) - {\bf v} \cdot \nabla T\]

where, by definition,

(40)\[D = {{k^T} \over {c^T}}\]

and

(41)\[{\bf v} = {{\rho_0 c_w} \over c^T} {\bf q}_w\]

The Peclet number gives a measure of the relative effects of advection and diffusion, and is defined as

(42)\[P_e = {{\|{\bf v}\|} \over {D / \Delta L}}\]

where \(\Delta L\) is the grid size. The Courant number is a measure of the advective distance covered by a particle during a timestep \(\Delta t^T\) over \(\Delta L\):

(43)\[C_r = {{\|{\bf v}\| \Delta t^T} \over {\Delta L}}\]

In the one-dimensional case, for instance, the classical constraints are \(P_e \leq 2\) and \(C_r \leq 1\) (see, for example, Perrochet and Bérod 1993). In the current implementation, the explicit fluid and thermal timesteps are calculated based on fluid and thermal diffusivities, and are identical to those adopted in the fluid and conduction logic. The code thus sets no intrinsic limit on the grid Peclet number. However, the accuracy is likely to be different for different number values, and influenced by grid size and timestep.

Warning: Note that the explicit timestep adopted by FLAC3D is not unconditionally stable for the advection logic.