# PFC Thermal Formulation

## Mathematical Model Description

The variables involved in heat conduction in a continuum are the temperature and the heat-flux vector. We wish to compute the time evolution of these two fields throughout the material. These variables are related by the continuity equation and transport laws derived from Fourier’s law of heat conduction. Substitution of Fourier’s law into the continuity equation yields the differential equation of heat conduction, which may be solved for particular geometries and properties, given specific boundary and initial conditions.

Conventions and Systems of Units

Vector and tensor quantities are expressed using indicial notation with respect to a fixed right-handed rectangular Cartesian coordinate system. The indices range over the set {1, 2, 3} for the three-dimensional case, and the set {1, 2} for the two-dimensional case. The Einstein summation convention is employed; thus, the repetition of an index in a term denotes a summation with respect to that index over its range. Vertical braces denote the magnitude of a vector or the absolute value of a scalar. A dot over a variable indicates a derivative with respect to time (e.g., \(\dot x_i = \partial x_i/\partial t\)). SI units are used to illustrate parameters and dimensions of variables. All thermal quantities must be given in a consistent set of units. No conversions are performed by the program. The following two tables present examples of consistent sets of units for thermal quantities.

Length | m | m | m | cm |

Density | kg/m^{3} |
10^{3}kg/m^{3} |
10^{6}kg/m^{3} |
10^{6}g/cm^{3} |

Temperature | K | K | K | K |

Time | s | s | s | s |

Specific Heat | J/(kg k) | 10^{-3}J/(kg K) |
10^{-6}J/(kg K) |
10^{-6}cal/(g K) |

Thermal Conductivity | W/(mk) | W/(mk) | W/(mk) | (cal/s)/cm^{2}K^{4} |

Convective Heat-Trans. Coefficient | W/(m^{2}K) |
W/(m^{2}K) |
W/(m^{2}K) |
(cal/s)/cm^{2}K |

Radiative Heat-Trans. Coefficient | W/(m^{2}K^{4}) |
W/(m^{2}K^{4}) |
W/(m^{2}K^{4}) |
(cal/s)/cm^{2}K^{4} |

Flux Strength | W/m^{2} |
W/m^{2} |
W/m^{2} |
(cal/s)/cm^{2} |

Source Strength | W/m^{3} |
W/m^{3} |
W/m^{3} |
(cal/s)/cm^{3} |

Decay Constant | s^{-1} |
s^{-1} |
s^{-1} |
s^{-1} |

Length | ft | in |

Density | slugs/ft^{3} |
snails/in^{3} |

Temperature | R | R |

Time | hr | hr |

Specific Heat | (32.17)^{-1} Btu/(1b R) |
(32.17)^{-1} Btu/(1b R) |

Thermal Conductivity | (Btu/hr)/(ft R) | (Btu/hr)/(in R) |

Convective Heat-Trans. Coefficient | (Btu/hr)/(ft^{2} R) |
(Btu/hr)/(in^{2} R) |

Radiative Heat-Trans. Coefficient | (Btu/hr)/(ft^{2} R^{4}) |
(Btu/hr)/(in^{2} R^{4}) |

Flux Strength | (Btu/hr)/ft^{2} |
(Btu/hr)/in^{2} |

Source Strength | (Btu/hr)/ft^{3} |
(Btu/hr)/in^{3} |

Decay Constant | hr^{-1} |
hr^{-1} |

Temperatures may be expressed in the Kelvin, Rankine, Celsius, and Fahrenheit scales. The Kelvin temperature scale uses the unit “Kelvin” (symbol K), the Rankine temperature scale uses the unit “Rankine” (symbol R), the Celsius temperature scale uses the unit “degree Celsius” (symbol °C), and the Fahrenheit temperature scale uses the unit “degree Fahrenheit” (symbol °F). The temperature scales are related by:

Governing Equations

By assuming that strain changes play a negligible role in influencing the temperature—a valid assumption for quasi-static mechanical problems involving solids and liquids—the heat-conduction equation for a continuum is given by

where \(q_i\) is the heat-flux vector [W/m^{2}], \(q_v\) is the volumetric heat-source intensity or power
density [W/m^{3}], \(\rho\) is the mass density [kg/m^{3}], \(C_v\) is the specific heat at constant
volume [J/kg°C], and \(T\) is the temperature [°C]. Note that for nearly all solids and liquids, the specific heats
at constant pressure and at constant volume are essentially equal. Consequently, \(C_p\) and \(C_v\) can be used
interchangeably.

Fourier’s law for a continuum defines the relation between the heat-flux vector and the temperature gradient as

where \(k_{ij}\) is the thermal-conductivity tensor [W/m°C].

Mechanical Coupling: Thermal Strains

Thermal strains are produced in the PFC material by accounting for the thermal expansion of the particles. The linear parallel bond contact model also accounts for the thermal expansion of the bonding material that joins them (but note that among the built-in contact models, only the linear parallel bond model includes thermal expansion). The thermal expansion is applied as follows.

Given a temperature change of \(\Delta T\), we change each particle radius, \(R\), such that

where \(\alpha\) is the coefficient of linear thermal expansion associated with the particle.

If a bonded linear parallel bond is present at the mechanical contact associated with a thermal contact, then we account for the expansion of the bond material by assuming that only the normal component of the force vector carried by the bond, \(\Delta \overline{F}_{ n}\), will be affected by the temperature change. We envision an isotropic expansion of the bond material, which effectively changes the bond length, \(\bar L\). This is modeled by changing the normal component of the bond force vector as

where \(\bar k_n\) is the bond normal stiffness, \(A\) is the area of the bond cross-section, \(\bar \alpha\) is the expansion coefficient of the bond material (taken equal to the average value of the expansion coefficients of the particles at the two ends of the pipe associated with the bond), \(\bar L\) is the bond length (taken equal to the distance between the centroids of the two particles at the ends of the pipe associated with the bond), and \(\Delta T\) is the temperature increment (taken equal to the average temperature change of the two particles at the ends of the pipe associated with the bond).

The value of \(\alpha\) in Equation (3) is a micro-property (see section below) and can be set equal to the coefficient of linear thermal expansion of a three-dimensional continuum, denoted by \(\alpha_t\):

Thermal Micro-properties

The following four micro-properties are used by the PFC thermal logic. The values assigned will differ depending on whether the PFC model is intended to represent a granular material (as an unbonded assembly of particles) or a solid continuum (as a bonded assembly of particles).

Density, \(\rho\) [kg/m^{3}], of each ball: When representing a granular material, the value of \(\rho\) should be set equal to the density of each grain. When representing a solid continuum, the value of \(\rho\) should be chosen such that the PFC material has the same total thermal mass in a fixed volume of material as does the physical material. This condition is satisfied by setting(6)\[\rho = \rho_t/(1-n)\]where \(\rho_t\) is the density of the solid continuum, and \(n\) is the average porosity of the PFC material. For densely packed and bonded PFC assemblies, the value of \(n\) is fixed by the initial microstructure and is approximately 0.16/0.40. The value of \(\rho\) is set with the

`ball attribute density`

command.

Specific heat at constant volume, \(C_v\) [J/kg°C], of each ball: This micro-property can be set equal to the value for a solid continuum. The value of \(C_v\) is set with the`ball thermal attribute sheat`

command.

Coefficient of linear thermal expansion, \(\alpha\) [1/°C], of each ball: When representing a granular material, the value of \(\alpha\) should correspond with that of each grain. This micro-property can be set equal to the value for a solid continuum, \(\alpha_t\), in the case of PFC. The value of \(\alpha\) is set with the`ball thermal attribute thexp`

command.

Thermal resistance per unit length, \(\eta\) [°C/W m], of each thermal contact: By scanning the pipe network of a given PFC thermal material, one can compute the value of \(\eta\) that, if assigned to all pipes, will produce a thermally isotropic material with a mean macroscopic conductivity, k [W/m°C] (see below). Alternatively, the value of \(\eta\) can be set directly with the contact.thermal.property.thres command.

Boundary and Initial Conditions

Initial conditions correspond to a given temperature field. Temperatures are assigned to balls, clumps, and walls with the
`ball thermal attribute temperature`

,
`clump thermal attribute temperature`

, or `wall thermal attribute temperature`

commands, respectively.
A ball or a clump may act as a heat-generating source if its applied thermal power is specified with the
`ball thermal attribute appliedpower`

or `clump thermal attribute appliedpower`

command, respectively.
Note that walls only act as a fixed temperature boundary condition.

In general, boundary conditions require a description of the local boundary in terms of an outwardly directed unit normal vector. However, in the PFC material, there is no explicit representation of the boundary (unlike in a continuum code, where zone surfaces make up the boundary). This problem is overcome by specifying the power being input to a boundary surface instead of the flux. (To convert the input power to a flux (power per cross-sectional area), divide the total power by the area normal to the direction of heat flow.)

Note that actual physical boundaries of the particle assembly are adiabatic by default (i.e., no heat will flow out of the
boundary of the
pipe network). Convective boundary conditions are not presently supported, but they could be added in the same way as in
*FLAC3D* ([Itasca2013]), by specifying the temperature of the surrounding fluid and the convective heat-transfer coefficient.

## Numerical Formulation

The PFC thermal material is represented as a network of heat reservoirs (associated with each particle) and thermal contacts
(associated with the mechanical contacts). Heat flow occurs via conduction in the active thermal contacts that connect the
reservoirs. Each ball
and clump represents a heat reservoir. Associated with each reservoir is a temperature, a mass, a volume, the specific heat,
\(C_v\), and coefficient of linear thermal expansion, \(\alpha\). The temperature is computed during the thermal
simulation, and the mass and volume are determined from the ball or clump density, size, and assumed shape. A thermal contact is
associated with each contact, but only some of these thermal contacts are active. By default, a thermal contact is active if the
host mechanical contact
is active (the criterion for a thermal contact to be active can be altered with the contact.thermal.activate and
contact.thermal.inhibit commands).
A thermal contact joins two reservoirs, and heat flow only occurs through thermal contacts. Associated with each thermal
contact is a power, \(Q\),
and a thermal contact model, which is responsible for updating the power \(Q\) based on the temperature increments in
the contacting pieces. There are, for now,
only two built-in thermal contact models: the *thermal null* contact model, which does not update the thermal contact power,
and the *thermal pipe* contact model.
The thermal pipe contact model defines three modifiable properties: a thermal resistance, \(\eta\), a coefficient of
linear thermal expansion, \(\bar \alpha\), and a
temperature, \(T\). Each thermal contact has a finite length that is equal to the distance between the centroids of the
two connected pieces.
Thermally induced strains are produced in the PFC material by modifying particle size and the force carried in each bond of
linear parallel bond models,
to account for heating of both particles and bonding material.

Numerical Discretization

The thermal material is discretized into a network of heat reservoirs and thermal pipes. Heat flow occurs via conduction in the active thermal contacts that connect the reservoirs. Radiative and convective heat transfer are not included in the present formulation.

We begin by establishing the heat-conduction equation for a single reservoir. Associate a control volume, \(V\), with the reservoir. (Typically, the sum of all reservoir control volumes equals the total material volume.) The heat outflow per unit volume is given by the divergence of \(q_i\). The average value of the divergence of \(q_i\) in the reservoir is defined by

The volume integral can be replaced with a surface integral by application of the Gauss divergence theorem to the reservoir:

where \(n_i\) is the outward unit normal vector of the surface, \(S\), of the control volume. If we restrict heat to flow only in the \(N\) thermal pipes associated with the reservoir, then the surface integral can be replaced by the summation

where the superscript (\(p\)) designates the value of the associated variable for pipe \(p\), and \(Q^{(p)}\) is the power in pipe \(p\) that is flowing out of the reservoir. (Note that \(q_i^{(p)} \Delta S^{(p)} = Q^{(p)} n_i^{(p)}\).)

Substitution of Equations (8) and (9) into Equation (7) yields

The heat-conduction equation for a single reservoir is found by substituting Equation (10) into Equation (1) to yield

where \(Q_v\) is the heat-source intensity, \(m\) is the thermal mass, and \(C_v\) is the specific heat at constant volume.

If we regard each pipe as a one-dimensional object with a thermal resistance per unit length of \(\eta\), then the power in a pipe is given by

where \(\Delta T\) is the temperature difference between the two reservoirs on each end of the pipe, and \(L\) is the pipe length.

We complete the discretization of Equation (11) by expressing the time derivative using forward finite differences. First, we rewrite the equation as

where \(\tilde Q\) is referred to as the out-of-balance power. Starting with an initial temperature field, the power in each pipe is updated via Equation (1). Then, reservoir temperatures are updated (provided the temperature is not fixed) using the forward finite-difference expression

where \(\Delta t\) is the thermal timestep.

Numerical stability of this explicit scheme can only be achieved if the timestep remains below a limiting value. To derive the stability criterion, we consider the situation in which a single reservoir is given a temperature perturbation, \(T_0\), from a zero initial state. After one thermal timestep, the new reservoir temperature is

To prevent alternating signs of temperature as the computation is repeated for successive \(\Delta t\) (which means that the temperature is oscillating about the initial state), the coefficient of \(T_0\) in the above relation must be positive. Such a requirement is satisfied if

In practice, it is found that the thermal calculation is stable if the thermal timestep is set equal to the minimum value (over all reservoirs) of Equation (16).

Relation between Thermal Resistance and Thermal Conductivity

We describe a procedure to estimate the conductivity tensor, \(k_{ij}\), of the PFC thermal material. Then, using this procedure, we compute the value of \(\eta\) that, if assigned to all pipes, will produce a thermally isotropic material with a mean macroscopic conductivity \(k\).

The average heat flux in a volume, \(V\), of material is defined by

where \(q_i\) is the heat-flux vector acting throughout the volume. By restricting heat to flow only in the thermal pipes, the integral can be replaced by a sum over all \(M\) pipes contained within \(V\):

where \(A^{(p)}\) is the cross-sectional area of pipe \(p\), and \(l^{(p)}\) is the length of pipe \(p\) contained within \(V\). The heat flux in a pipe is given by

where \(Q\) is the power in the pipe, and \(n_i\) is the unit vector directed along the pipe. After substitution of Equation (1), the heat flux becomes

where \(\Delta T\) is the temperature difference across the pipe. If the mean temperature gradient in the material, \({\partial T \over \partial x_j}\), is assumed to apply at the micro-level within each pipe, then

Substituting Equations (20) and (21) into Equation (18) yields

By comparing this expression with Fourier’s law in Equation (2), we find that the thermal conductivity tensor of the PFC thermal material is given by

Given a PFC thermal material, the values of the thermal conductivity tensor can be computed using Equation (23). However, definition of the volume in that expression is problematic because of particles that intersect the measurement region. Therefore, the following procedure is employed to evaluate \(k_{ij}\) over a region defined by a measurement sphere. We associate each reservoir with a PFC particle, and note that in a statistically uniform assembly, the volume associated with each particle is \({V^{(b)}/(1-n)}\), where \(n\) is the porosity and \(V^{(b)}\) is the particle volume. Substituting into Equation (23) yields

where the summations are taken over the \(N_b\) balls with centroids contained within the measurement sphere, and the \(N_p\) active thermal pipes of these balls; and:

- \(n\)
- is the porosity within the measurement sphere;
- \(V^{(b)}\)
- is the volume of ball \((b)\);
- \(l^{(p)}\)
- is the length of thermal pipe \((p)\) associated with balls in the measurement sphere. (A ball is considered to be in the sphere if its centroid lies in the sphere. If both balls used by the pipe are in the sphere, then \(l = L\) (the pipe length). If only one ball is in the sphere, then \(l = LV_1/(V_1+V_2)\), where \(V_1\) and \(V_2\) are the volumes of the balls inside and outside of the sphere, respectively.);
- \(n_i^{(p)}\)
- is the unit normal vector directed along pipe \((p)\); and
- \(\eta^{(p)}\)
- is the thermal resistance of pipe \((p)\).

In a thermally isotropic material, the thermal conductivity tensor can be expressed in terms of a single parameter, \(k\), such that \(k_{ij} = k\delta_{ij}\). For such a material, assuming uniform \(\eta\) for all pipes:

which can be rewritten to express \(\eta\) in terms of \(k\) as

Equation (26) gives the value of \(\eta\) that, if assigned to all pipes, will produce a thermally isotropic material with a mean macroscopic conductivity, \(k\). The summation in Equation (26) can be evaluated with a single scan through the particle assembly. A check that the particle assembly is suitable for representing a thermally isotropic material is to confirm that the ratios

are small (The four last ratio above only apply to the 3D case). When applying Equation (26) to a region defined by a measurement sphere, we use the expression

where the terms are defined in Equation (24).

## Solving Thermal Problems

PFC supports both thermal-only and coupled thermal-mechanical analysis. The thermal calculation mode is enabled by
issuing the `model configure thermal`

command. This allocates memory required to implement the thermal logic. After issuing
this command, the thermal calculation mode is on, and the mechanical calculation mode is on. The calculation modes can be
modified at any time during a PFC analysis by issuing the `model thermal`

and `model mechanical`

commands with the
on/off keywords. When the thermal calculation mode is on, the thermal timestep will be calculated automatically,
but a user-defined timestep can be selected using the `model mechanical timestep`

thermal.fix, or a maximum timestep limit can be
set with
the `model mechanical timestep`

thermal.max command. The thermal micro-properties must be specified for all particles,
and initial and boundary conditions must be assigned.

Thermal-Only Analysis

The command set thermal on mechanical off is used to request a thermal-only calculation. The `model cycle`

command may
then be given to execute a given number of thermal steps. Alternatively, the `model solve`

command may be given to continue
cycling until either a particular thermal time has been reached (with the solve.thermal.time keyword) or a steady-state
thermal solution has
been obtained (the solve.thermal.ratio-average keyword). The thermal timestep and total thermal time are printed to the screen
when `model cycle`

or `model solve`

is given.

Thermal-Mechanical Analysis

When performing a thermal-mechanical analysis, one must be aware of the different timescales associated with the thermal and mechanical processes. Because the critical timestep for each process corresponds to the time needed for information to propagate from one particle to the next, these values can be used to compare timescales.

In a dense system, the timestep associated with the mechanical process, which corresponds with the time for a \(p\)-wave to propagate a distance \(L_c\) through the medium, can be estimated as

where \(\rho\) is the mass density, \(K\) is the bulk modulus, \(G\) is the shear modulus, and \(L_c\) is a characteristic length. The timestep associated with the thermal process, which corresponds with the time for the diffusion front to propagate a distance \(L_c\) through the medium, can be estimated as ([Itasca2013])

where \(\kappa\) is the thermal diffusivity (\(k / \rho C_v\)) and \(L_c\) is the characteristic length (volume of solid divided by surface area exchanging heat). If we set both characteristic lengths equal to one another, then the ratio of thermal-to-mechanical timescales may be expressed as

For rock and soil, \(\kappa\) is on the order of 10^{-6} m^{2}/s, \(\rho\) is on the order of
10^{3} kg/m^{3}, and \(K + (4/3) G\) is on the order of 10^{10} N/m^{2}. Using these orders of
magnitude in Equation (31), we see that the ratio of thermal-to-mechanical timescales is on the order of 10^{9}
\(L_c\). This ratio is very large, even if \(L_c\) is 1 mm. In practice, mechanical effects can be assumed to occur
instantaneously when compared to diffusion effects. This assumption is adopted in PFC, where one or more mechanical
sub-steps may be taken after each thermal step, and the thermal time is not incremented during the mechanical sub-steps.

In most modeling situations, the initial mechanical conditions 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 is to be investigated, for instance, at a certain thermal time, a thermal-only calculation (set thermal on mechanical off) may be performed until the desired time is reached. The thermal calculation may then be switched off (set thermal off) and the mechanical calculation switched on (set mechanical on), and the system may be cycled to mechanical equilibrium before the response is analyzed.

For nonlinear mechanical models (i.e., when bond breakages occur or when the contact fabric changes under load), 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 adjust according to the different timescales involved. 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.
Stopping conditions for which mechanical and thermal equilibrium will be assumed can be specified. The command
`model mechanical substep`

specifies the maximum number of mechanical sub-steps to be performed for each thermal
step. If the `model solve`

command is given, then, after each thermal step, sufficient mechanical sub-steps (with a maximum
of \(m\)) will be taken until mechanical equilibrium has been achieved. (If a `model cycle`

\(n\) command is given,
then \(n\) thermal steps will be taken, with \(m\) mechanical sub-steps being taken after each thermal step.)

References

[Itasca2013] | (1, 2) Itasca Consulting Group Inc. FLAC3D (Fast Lagrangian Analysis of Continua in 3 Dimensions), Version 5.01. Minneapolis, Minnesota: ICG, 2013. |

Was this helpful? ... | FLAC3D © 2019, Itasca | Updated: Feb 25, 2024 |