Solving Thermal-Only and Coupled-Thermal Problems

FLAC3D has the ability to perform thermal-only and coupled thermal fluid-flow analyses. In addition, both types of simulation can be coupled to a mechanical analysis. In all cases, the first step in the command procedure is to issue a model configure thermal command so that extra memory can be assigned for the thermal calculation. Then, a choice must be made between explicit and implicit thermal algorithms. By default, the explicit algorithm will be selected, but the implicit mode of calculation may be activated or deactivated at any stage of the thermal-only or coupled calculation using the zone thermal implicit on or zone thermal implicit off command. In the explicit mode, the thermal timestep will be calculated automatically, but a smaller timestep can be selected using the model thermal timestep command. The magnitude of the timestep must be specified by the user in the implicit mode. This is done by issuing a model thermal timestep command. The thermal model and properties must be specified for any zones that conduct heat. There are four thermal models available in this version of the code (see zone thermal cmodel): null, isotropic, anisotropic and advection-conduction, as well as a hydration thermal model called hydration. The null model is used to exclude part of the grid from the thermal calculation. (The thermal domain is defined by the assembly of zones with a non-null thermal model.) Note that flux boundary conditions, for instance, are assigned through a range to the boundaries of the thermal domain. Also, a zone to which a null mechanical model has been assigned will only be ignored in a thermal calculation, provided a null model is also assigned to the zone. The isotropic and anisotropic models are used for heat conduction analyses in which thermal conductivity is isotropic and anisotropic, respectively. Note that anisotropic can be used for isotropic conductivity analyses; however, it is computationally less efficient than the isotropic model for this case. The advection-conduction model is used for conduction-advection calculations; thermal conductivity is assumed to be isotropic for this model. Note that any combination of thermal constitutive models is allowed in a conduction-advection analysis; however, the transport of heat by advection will only be considered in zones to which the advection-conduction model has been assigned. Material density must also be specified. Initial and boundary conditions are assigned to complete the thermal problem setup.

Thermal Conduction-Only Analysis

The command zone mechanical active off (and zone fluid active off, if model configure fluid is specified) is used to request a thermal-only calculation. The model step (or model cycle) command may then be given to execute a given number of thermal steps. To stop the calculation when a particular thermal time is reached, a model solve thermal time-total command may be issued.

The thermal properties, isotropic thermal conductivity and specific heat, are assigned with the zone thermal property command followed by the keyword conductivity and specific-heat, respectively.

A steady-state thermal solution can be obtained by simply issuing the model solve thermal command. By default, the thermal calculation will be performed until the ratio of the unbalanced heat flux at a gridpoint divided by the applied heat flux magnitude falls below the value of 10^-5 for all gridpoints in the model. This limiting ratio produces a steady-state heat flux condition for most thermal analyses. The ratio can be adjusted with the model solve thermal ratio command.

The thermal timestep and total thermal time are printed to the screen when model solve thermal cycles or model solve thermal time is given. The timestep and other thermal calculation-mode information can be obtained with the zone thermal list command.

Thermal Advection-Only Analysis

The logic for advection-conduction relies on sequential stepping of fluid-only and thermal-only calculations; it is activated in zones to which the thermal constitutive model advection-conduction has been assigned. The constitutive model advection-conduction (see zone thermal property) is responsible for calculation of thermal flux and temperature gradient from gridpoint temperatures for the zones. The model also stores the thermal properties of conductivity, specific heat, thermal expansion and reference temperature for the zone, as well as fluid density (for buoyancy calculation) and fluid-specific discharge. Conductivity is assumed to be isotropic for this model. The spatial distribution can be arbitrary for all properties. Finally, both explicit and implicit thermal and fluid calculations can be performed with the advection-conduction logic.

In most instances, the fluid configuration will be selected when using the advection logic. Free advection can only be activated (in zones containing the advection-conduction model), provided this configuration is adopted. Forced convection can be used with or without the fluid configuration. In the latter case, the (known) fluid specific discharge is assigned as a model property.

In the fluid configuration (model configure fluid), the user can assign several properties (see zone thermal property):

conductivity	matrix thermal conductivity
temperature-reference	reference temperature
expansion-fluid	coefficient of volumetric thermal expansion for the fluid
conductivity-fluid	fluid thermal conductivity
specific-heat-fluid	fluid specific heat
specific-heat	matrix specific heat
expansion	matrix linear thermal expansion coefficient

In this case, the properties conductivity-effective, specific-heat-effective, specific-discharge-x, specific-discharge-y and specific-discharge-z are calculated automatically by the code and should only be accessed for reading. (If values are assigned by the user, they will be ignored.)

If the fluid configuration is not selected, then the properties conductivity-effective, specific-heat-effective, specific-discharge-x, specific-discharge-y and specific-discharge-z must be assigned values; the properties specific-heat, conductivity and conductivity-fluid are not used. In this case, the properties to be assigned by the user are (see zone thermal property):

conductivity-effective	effective thermal conductivity
specific-heat-effective	effective specific heat over matrix density
specific-discharge-x	\(x\)-component of specific discharge
specific-discharge-y	\(y\)-component of specific discharge
specific-discharge-z	\(z\)-component of specific discharge
temperature-reference	reference temperature
specific-heat-fluid	fluid specific heat
expansion	matrix linear thermal expansion coefficient

In both configuration cases, the advective term in the heat equation will be activated only if both fluid density (zone water density) and fluid specific heat (zone thermal property specific-heat-fluid) are nonzero.

Steady-State Conduction Solution

The initial conditions for an advection simulation often correspond to a steady-state conduction solution. These conditions can be established by initializing pore pressures and temperatures directly, if they are known, or by letting the code do the calculation. In the latter case, the procedure is to perform the thermal calculation first, using a zero value of reference fluid density to prevent addition of the advection term in the energy balance equation (using zone water density = 0, with model fluid active off and model thermal active on), and then perform the fluid calculation (with model fluid active on and model thermal active off) after resetting the fluid reference density to a realistic value.

Forced Advection and Free Advection Simulations

The property expansion-fluid (see zone thermal property) can be set to zero to run simulations in which free advection is neglected. To model free advection, and neglect forced advection, the fluid reference density can be set to zero during the thermal calculation, and reset to its actual value for the fluid flow calculation.

Synchronization of Fluid and Thermal Times

Note that the modules must be run independently in this implementation of the advection logic: fluid flow must be turned off during thermal stepping, and thermal should be off during the flow calculation.

It is also the user’s responsibility to synchronize fluid and thermal times as the simulation proceeds. To do this for explicit calculations, the following procedure can be followed. First, the default timesteps for fluid and thermal calculation are determined by performing one “preview” calculation step of each individual module, and recording the timestep value displayed on the screen. The default fluid (thermal) timestep is a function of grid size, fluid (thermal) properties of mobility, fluid bulk modulus and porosity (effective conduction, specific heat). Second, fluid and thermal timesteps are selected such that: (a) one is an integer multiple of the other; and (b) each timestep is smaller than or equal to its default value. (The model fluid timestep or model thermal timestep command is used for this.) Finally, the simulation is run in an alternating sequence of thermal and fluid-flow steps: one step of the module with the largest timestep; and enough steps of the other module to close the time gap. (An example of a FISH function devised to do this can be found in the Horton-Rogers-Lapwood Convection Problem.) In addition, if mechanical coupling is considered, enough mechanical steps should be included in this sequence to keep the system in quasi-static equilibrium.

Thermal-Mechanical Analysis

The thermal option can be combined with the mechanical calculation to perform a thermal-mechanical analysis with FLAC3D. All of the features of the thermal calculation, including transient and steady-state heat transfer, and thermal solution by either explicit or implicit algorithm, are available in a thermal-mechanical calculation.

The thermal-mechanical coupling is provided by the influence of temperature change on the volumetric change of a zone (see this equation). The linear thermal expansion coefficient is assigned via the zone thermal property expansion command.

Thermal-mechanical analysis may be performed with any of the built-in mechanical models in FLAC3D. Note, however, that the thermal model is not made null automatically for zones that are made null mechanically. The zone thermal cmodel assign null command must be given for thermal null zones.

Consider the transport of heat by conduction: there are two issues that must be taken into consideration when performing a thermal-mechanical analysis. The first one is concerned with the different time scales associated with the thermal and mechanical processes. Because timesteps correspond to the time needed for the information to propagate from one node to the next, their formulation can be used to compare time scales.

Typically, the timesteps associated with the mechanical process are of the form

(1)\[\Delta t_{mech} = \sqrt{{\rho}\over{K + (4/3) G}} \ L_c\]

where:	\(K\)	=	bulk modulus;
	\(G\)	=	shear modulus;
	\(\rho\)	=	density; and
	\(L_c\)	=	characteristic length.

Hence, the ratio of thermal-conduction to mechanical timestep may be expressed as (see this equation, \(h\) = 0 and \(m\) = 1)

(2)\[{{\Delta t_{ther}}\over{\Delta t_{mech}}} = \sqrt{{K + 4/3 G}\over{\rho}} \ {{L_c}\over {\kappa}}\]

where \(\kappa\) is the thermal diffusivity. For non-metals, this property is of the order of 10^-6 m²/s, at most. For rock and soil, \(\rho\) is of the order of 10³ kg/m³, while \(K + (4/3) G\) is approximately 10¹⁰ N/m². Using those orders of magnitude in Equation (2), we see that the ratio of thermal-to-mechanical time scales is of the order 10⁵ \(L_c\). This ratio remains very large, even if \(L_c\) is of the order of 1 mm. In practice, mechanical effects can be assumed to occur instantaneously when compared to diffusion effects. This is also the approach adopted in FLAC3D, where no time is associated with any of the mechanical steps taken in association with the thermal steps.

The second issue concerns the thermal-mechanical coupling in FLAC3D. This coupling occurs only in one direction: temperature changes cause thermal strains that influence the stresses to occur; while the thermal calculation is unaffected by the mechanical changes taking place.

Typically, the initial mechanical conditions in a simulation correspond to a state of equilibrium that must first be achieved before the coupled analysis is started. If the medium is elastic, and the thermal-mechanical response must be investigated at a certain thermal time without analyzing the transients, the thermal and mechanical simulations can be carried out in succession. First, a thermal-only calculation is performed until the desired time (model thermal active on and model mechanical active off). Following this, the system may be cycled to mechanical equilibrium, using the model thermal active off and model mechanical active on commands, before the response is analyzed.

For nonlinear mechanical models (i.e., when plasticity is involved), the thermal changes must be communicated to the mechanical module at closer time intervals to respect the path-dependency of the system. Typically, at small dimensionless thermal time, a certain number of mechanical steps must be taken for each thermal step to allow the system to reach quasi-static equilibrium. At large dimensionless thermal time, if the system approaches thermal equilibrium, several thermal timesteps may be taken without significantly disturbing the mechanical state of the medium. A corresponding numerical simulation may be controlled manually by alternating between thermal-only and mechanical-only modes.

Such a tedious task may be avoided by using the model solve command in combination with appropriate settings: model thermal active on and model mechanical active on make both thermal and mechanical processes active. model solve mechanical ratio will set a limit to the out-of-balance force ratio under which mechanical “equilibrium” will be assumed. model mechanical substep n and model mechanical slave on declares the mechanical module as the slave which must perform n sub-cycles (or fewer, if equilibrium is detected — see previous setting) for each step taken by the master. model thermal substep n declares the thermal module as the master. By default, the thermal module will be the master with m = 100 if the mechanical module is declared as slave.

If, for each thermal timestep, the mechanical module needs only one sub-step to reach equilibrium, then the number of thermal sub-steps (which are now slave sub-steps to the mechanical module) is doubled, but it never exceeds the number m. The system reverts to the original setting of sub-steps whenever the succession of one mechanical timestep for each thermal group of sub-steps is broken.

The following information is printed to the screen while cycling with the model solve command in this mode:

• total step number (including thermal and mechanical sub-steps)
• total sub-step number for master module at previous step
• total sub-step number for slave module at previous step
• current module being calculated
• current sub-step number
• current maximum mechanical unbalanced force ratio
• current maximum thermal unbalanced force ratio
• current timestep for master module

This information can be monitored to ensure that the calculation is proceeding as intended.

In a third approach, the model step command may be used while both mechanical and thermal modules are on (model thermal active on and model mechanical active on). In this option, one mechanical step will be taken for each thermal step. Here, thermal steps are assumed to be so small that one mechanical step is enough to re-equilibrate the system mechanically after each thermal step is taken. Input Instructions for Thermal Analysis may be consulted for a complete list of available command options.

Thermal Coupling to Pore Pressure

FLAC3D can also be used to solve coupled thermal-fluid flow problems in which temperature change affects fluid pore pressures. This coupled response can also be coupled to the mechanical stress calculation.

The coupling is induced by the influence of temperature change on volumetric change of the saturated matrix volume. The undrained volumetric thermal expansion coefficient, \(\beta\), is prescribed as an input property assigned by the zone fluid property undrained-thermal-coefficient.

In order to apply thermal-fluid flow coupling, the FLAC3D grid must be configured for both the thermal and fluid flow options (model configure fluid and model configure thermal). If advection of heat is not considered, the thermal and fluid flow calculation modes can be run independently by suppressing one mode or the other (i.e., zone thermal active on and zone fluid active off or zone thermal active off and zone fluid active on), or the modes can be run concurrently (zone thermal active on and zone fluid active on). In the latter case, the thermal and fluid flow timesteps are synchronized with each other, based upon the smaller of the two. If advection is considered, the fluid and thermal calculations must be performed is sequence. In addition, fluid and thermal time scales must be synchronized by the user, following the procedure described in Synchronization of Fluid and Thermal Times.