# 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

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)

where \(\kappa\) is the thermal diffusivity.
For non-metals, this property is of the order of 10^{-6} m^{2}/s, at most.
For rock and soil, \(\rho\) is of the order of 10^{3} kg/m^{3},
while \(K + (4/3) G\) is approximately 10^{10} N/m^{2}.
Using those orders of magnitude in Equation (2),
we see that the ratio of thermal-to-mechanical time scales is of the order 10^{5} \(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.

Was this helpful? ... | FLAC3D © 2019, Itasca | Updated: Oct 08, 2020 |