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.


Table 1: System of SI units for thermal problems
Length m m m cm
Density kg/m3 103kg/m3 106kg/m3 106g/cm3
Temperature K K K K
Time s s s s
Specific Heat J/(kg k) 10-3J/(kg K) 10-6J/(kg K) 10-6cal/(g K)
Thermal Conductivity W/(mk) W/(mk) W/(mk) (cal/s)/cm2K4
Convective Heat-Trans. Coefficient W/(m2K) W/(m2K) W/(m2K) (cal/s)/cm2K
Radiative Heat-Trans. Coefficient W/(m2K4) W/(m2K4) W/(m2K4) (cal/s)/cm2K4
Flux Strength W/m2 W/m2 W/m2 (cal/s)/cm2
Source Strength W/m3 W/m3 W/m3 (cal/s)/cm3
Decay Constant s-1 s-1 s-1 s-1

Table 2: System of Imperial units for thermal problems
Length ft in
Density slugs/ft3 snails/in3
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)/(ft2 R) (Btu/hr)/(in2 R)
Radiative Heat-Trans. Coefficient (Btu/hr)/(ft2 R4) (Btu/hr)/(in2 R4)
Flux Strength (Btu/hr)/ft2 (Btu/hr)/in2
Source Strength (Btu/hr)/ft3 (Btu/hr)/in3
Decay Constant hr-1 hr-1

\[\begin{split}\mbox{~where~~} \begin{array}{lcl} 1\mbox{K} & = & 1.8 \mbox{R}; \\ 1\mbox{J} & = & 0.239 \mbox{ cal} = 9.48 \times 10^{-4} \mbox{ Btu}; \\ 1\mbox{J/kg K} & = & 2.39 \times 10^{-4} \mbox{ btu/1b R}; \\ 1\mbox{W} & = & 1 \mbox{ J/s} = 0.239 \mbox{ cal/s} = 3.412 \mbox{ Btu/hr}; \\ 1\mbox{W/m K} & = & 0.578 \mbox{ Btu/(ft/hr R); and} \\ \mbox{W/m}^2\mbox{K} & = & 0.176 \mbox{ Btu/ft}^2\mbox{hr R.} \end{array}\end{split}\]

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:

\[\begin{split}\begin{array}{lcl} \mbox{Temp(}^\circ \mbox{C)} & = & \frac{5}{9} \ast \mbox{[Temp(}^\circ \mbox{F) - 32]}; \\ \mbox{Temp(}^\circ \mbox{F)} & = & \mbox{[1.8} \ast \mbox{Temp(}^\circ \mbox{C)] + 32}; \\ \mbox{Temp(}^\circ \mbox{C)} & = & \mbox{Temp(K) - 273.15; and} \\ \mbox{Temp(}^\circ \mbox{F)} & = & \mbox{Temp(R) - 459.67.} \end{array}\end{split}\]

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

(1)\[- \frac{\partial q_i}{\partial x_i} + q_v = \rho C_v \frac{\partial T}{\partial t}\]

where \(q_i\) is the heat-flux vector [W/m2], \(q_v\) is the volumetric heat-source intensity or power density [W/m3], \(\rho\) is the mass density [kg/m3], \(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

(2)\[q_i = -k_{ij} \frac{\partial T}{\partial x_j}\]

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

(3)\[\Delta R = \alpha R \Delta T\]

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

(4)\[\Delta \overline{F}_{ n} = - \bar{k}_n A \Delta U_n = - \bar{k}_n A (\bar \alpha \bar{L} \Delta T)\]

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\):

(5)\[\alpha = \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/m3], 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

(7)\[\left< {\partial q_i \over \partial x_i} \right> = {1 \over V} \int_V{\partial q_i \over \partial x_i}\,dV\]

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

(8)\[\int_V{\partial q_i \over \partial x_i}\,dV = \int_S{q_i n_i}\,dS\]

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

(9)\[\int_S{q_i n_i}\,dS = \sum_{p=1}^N{q_i^{(p)} n_i^{(p)} \Delta S^{(p)} } = \sum_{p=1}^N{Q^{(p)} }\]

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

(10)\[\left< {\partial q_i \over \partial x_i} \right> = {1 \over V} \sum_{p=1}^N{Q^{(p)}}\]

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

(11)\[-\sum_{p=1}^N{Q^{(p)}} + Q_v = m C_v {\partial T \over \partial t}\]

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

(12)\[Q = -{\Delta T \over \eta L}\]

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

(13)\[{\partial T \over \partial t} = {1 \over m C_v}\left[ -\sum_{p=1}^N{Q^{(p)}} + Q_v \right] = {1 \over m C_v} \tilde Q\]

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

(14)\[\begin{split}\begin{array}{lcl} T_{< t + \Delta t >} &=& T_{< t >} + \Delta T \\ &= & T_{< t >} + \Delta t \left[\frac{1}{m C_v} \tilde Q \right]_{< t >} \end{array}\end{split}\]

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

(15)\[T_{< t + \Delta t >} = T_0 \left[ 1 - {\Delta t \over m C_v} \sum_{p=1}^N{1 \over \eta^{(p)} L^{(p)}} \right]\]

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

(16)\[\Delta t < {m C_v \over \sum{1 \over \eta L}}\]

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

(17)\[\left< q_i \right> = {1 \over V} \int_V{q_i\,dV}\]

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\):

(18)\[\left< q_i \right> = {1 \over V} \sum_{p=1}^M{q_i^{(p)} V^{(p)}} = {1 \over V} \sum_{p=1}^M{q_i^{(p)} A^{(p)} l^{(p)}}\]

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

(19)\[q_i = {Q n_i \over A}\]

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

(20)\[q_i = -{\Delta T n_i \over \eta L A}\]

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

(21)\[\Delta T = L n_j {\partial T \over \partial x_j}\]

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

(22)\[\left< q_i \right> = - \left[ {1 \over V} \sum_{p=1}^M{l^{(p)} n_i^{(p)} n_j^{(p)} \over \eta^{(p)}} \right] {\partial T \over \partial x_j}\]

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

(23)\[k_{ij} = {1 \over V} \sum_{p=1}^M{l^{(p)} n_i^{(p)} n_j^{(p)} \over \eta^{(p)}}\]

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

(24)\[k_{ij} = \left( { 1-n \over \sum\limits_{N_b} V^{(b)} } \right) { \sum\limits_{N_p}{ l^{(p)} n_i^{(p)} n_j^{(p)} \over \eta^{(p)}} }\]

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:

(25)\[\begin{split}k = \begin{cases} {1\over 3}\left( k_{11} + k_{22} + k_{33} \right) = {1 \over 3 V \eta} \sum_{p=1}^M{l^{(p)} \left( n_1^2 + n_2^2 + n_3^2 \right)} = {1 \over 3 V \eta} \sum_{p=1}^M{l^{(p)}} \qquad {\rm (3D)} \\ {1\over 2}\left( k_{11} + k_{22} \right) = {1 \over 2 V \eta} \sum_{p=1}^M{l^{(p)} \left( n_1^2 + n_2^2\right)} = {1 \over 2 V \eta} \sum_{p=1}^M{l^{(p)}} \qquad {\rm (2D)} \end{cases}\end{split}\]

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

(26)\[\begin{split} \eta &= {1 \over 3 V k} \sum_{p=1}^M{l^{(p)}} {\rm (3D)} \\ \eta &= {1 \over 2 V k} \sum_{p=1}^M{l^{(p)}} {\rm (2D)}\end{split}\]

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

(27)\[\begin{split}{|k_{12}| \over k} &= {2 \sum{l |n_1 n_2|} \over \sum{l}} \\ {\left| k_{11} - k_{22} \right| \over k} &= {2 \sum{l|n_1^2-n_2^2|} \over \sum{l}} \\ {|k_{13}| \over k} &= {2 \sum{l |n_1 n_3|} \over \sum{l}} \\ {|k_{23}| \over k} &= {2 \sum{l |n_2 n_3|} \over \sum{l}} \\ {\left| k_{11} - k_{33} \right| \over k} &= {2 \sum{l|n_1^2-n_3^2|} \over \sum{l}} \\ {\left| k_{22} - k_{33} \right| \over k} &= {2 \sum{l|n_2^2-n_3^2|} \over \sum{l}}\end{split}\]

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

(28)\[\begin{split}\begin{cases} \eta = { 1 \over 3k } \left( { 1-n \over \sum\limits_{N_b} V^{(b)} } \right) { \sum\limits_{N_p} l^{(p)} } {\rm (3D)} \\ \eta = { 1 \over 2k } \left( { 1-n \over \sum\limits_{N_b} V^{(b)} } \right) { \sum\limits_{N_p} l^{(p)} } {\rm (2D)} \end{cases}\end{split}\]

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

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

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])

(30)\[\Delta t_{ther} = { L_c^2 / \kappa }\]

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

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

For rock and soil, \(\kappa\) is on the order of 10-6 m2/s, \(\rho\) is on the order of 103 kg/m3, and \(K + (4/3) G\) is on the order of 1010 N/m2. Using these orders of magnitude in Equation (31), we see that the ratio of thermal-to-mechanical timescales is on the order of 109 \(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.