Mathematical Model Description
Conventions and Definitions
As a notation convention, the symbol \(a_i\) denotes component \(i\) of the vector \([a]\) in a Cartesian system of reference axes; \(A_{ij}\) is component \((i,j)\) of tensor \([A]\). Also, \(f_{,i}\) is used to represent the partial derivative of \(f\) with respect to \(x_i\). (\(f\) can be a scalar variable or a vector component.)
The Einstein summation convention applies only on indices \(i\), \(j\) and \(k\), which take the values 1, 2, 3 for components that involve spatial dimensions. The indices may take any values when used in matrix equations.
SI units are used to illustrate parameters and dimensions of variables (see this section for conversions to other systems of units).
The following dimensionless numbers are useful in the characterization of transient heat conduction:
Characteristic length
Thermal diffusivity
where: 
\(k\) 
is the thermal conductivity; 
\(\rho\) 
is the density; and 

\(C_v\) 
is the specific heat at constant volume. 
Characteristic time
Conduction
The variables involved in heat conduction in FLAC3D are the temperature and the three components of the heat flux. These variables are related through the energybalance equation and transport laws derived from Fourier’s law of heat conduction. Substitution of Fourier’s law into the energybalance equation yields the differential equation of heat conduction, which may be solved for particular geometries and properties, given specific boundary and initial conditions. Thermal volumetricstrains are introduced into the incremental mechanical and fluid constitutive laws to account for thermalstress and thermalpore pressure couplings.
EnergyBalance Equation
The differential expression of the energy balance has the form
where \(q_i\) is the heatflux vector in \([W/m^2]\), \(q_v\) is the volumetric heatsource intensity in \([W/m^3]\), and \(\zeta\) is the heat stored per unit volume in \([J/m^3]\). In general, temperature changes may be caused by changes in both energy storage and volumetric strain, \(\epsilon\), and the thermal constitutive law relating those parameters can be expressed as
where \(M_{th}\) and \(\beta_{th}\) are material constants, and \(T\) is temperature.
In FLAC3D, we consider a particular case of this law for which \(\beta_{th}\) = 0 and \(M_{th} = 1/(\rho C_v)\). \(\rho\) is the mass density of the medium in \([kg/m^3]\), and \(C_v\) is the specific heat at constant volume in \([J/kg\ ^{\circ}C]\). The hypothesis here is that strain changes play a negligible role in influencing the temperature, a valid assumption for quasistatic mechanical problems involving solids and liquids. Accordingly,
Substitution of Equation (6) in Equation (4) yields the energybalance equation
Note that, for nearly all solids and liquids, the specific heats at constant pressure and at constant volume are essentially equal. Consequently, \(C_v\) and \(C_p\) can be used interchangeably.
Transport Law
The basic law that defines the relation between the heatflux vector and the temperature gradient is Fourier’s law. For a stationary, homogeneous, isotropic solid, this constitutive law is given in the form
where \(T\) is the temperature [\(^{\circ}\)C], and \(k\) is the thermal conductivity in \([W/m ^{\circ}C]\).
Isotropic Heat Conduct
For isotropic conduct, the thermal conductivity matrix has the form:
For this case, one conductivity scalar parameter is sufficient.
An alphabetical list of properties of the isotropic heat conduction thermal model is given here.
Anisotropic Heat Conduct
For anisotropic conduct, the thermal conductivity matrix has the form:
For this case, three conductivity parameters corresponding to different directions are required.
In general,
For this case, six conductivity parameters are required.
An alphabetical list of properties of the anisotropic heat conduction thermal model is given here.
Boundary and Initial Conditions
Substitution of Equation (8) for \(q_i\) in Equation (7) yields the differential equation for heat conduction. Initial conditions correspond to a given temperature field. Boundary conditions are generally expressed in terms of temperature or the component of the heatflux vector normal to the boundary. In this version of FLAC3D, four types of conditions are considered, corresponding to: (1) given temperature; (2) given component of the flux normal to the boundary; (3) convective boundaries; and (4) insulated (adiabatic) boundaries.
In FLAC3D, a convective boundary condition has the form
where \(q_n\) is the component of the flux normal to the boundary in the direction of the exterior normal, \(h\) is the convective heattransfer coefficient \([W/m^2$ $^{\circ}C]\), \(T\) is the temperature of the boundary surface, and \(T_e\) is the temperature of the surrounding fluid [\(^{\circ}\)C]. Note that in the numerical formulation used in FLAC3D, boundaries are adiabatic by default.
Mechanical Coupling — Thermal Strains
Solution of thermalstress problems requires reformulation of the incremental stressstrain relations, which is accomplished by subtracting the portion due to temperature change from the total strain increment. Because free thermal expansion results in no angular distortion in an isotropic material, the shearingstrain increments are unaffected. The thermalstrain increments associated with the free expansion corresponding to temperature increment \(\Delta T\) have the form
where \(\alpha_t\ [1/^{\circ}C]\) is the coefficient of linear thermal expansion, and \(\delta_{ij}\) is the Kronecker delta.
The differential equation of motion and the rate of strainvelocity relations are based upon mechanical considerations, and are unchanged for thermomechanics.
Fluid Coupling — Thermally Induced Pore Pressures
The heat transfer may be coupled to the groundwater calculation by taking into account porepressure change caused by thermal expansion of the saturated matrix. The constitutive law has the form (see this this equation in FLAC3D FluidThermalMechanicalFormulation — Mathematical Description)
where \(\zeta\) is the variation of fluid content, \(M\) is the Biot modulus, \(\alpha\) is the Biot coefficient, \(\epsilon\) is the volumetric strain, and \(\beta\) is the undrained volumetric thermal expansion. In FLAC3D, the porepressure increment \(\Delta p_{ther}\), caused by a temperature change \(\Delta T\), has the form
\(\beta\ [1/^{\circ}C]\) is defined as the porepressure variation divided by (\(3 \alpha_t\ M\)) per unit temperature change in an undrained constrained test (no deformation). For an ideal porous material, this coefficient is related to the volumetric thermalexpansion coefficients for the grains, \(3 \alpha_t\), and the fluid, \(3 \alpha_f\), through the formula
where \(n\) = porosity.
For an incompressible solid constituent, (\(\alpha\) = 1 and \(M = {K_f} / {n}\)), the porepressure increment becomes
where \(K_f\) = fluid bulk modulus.
The corresponding change in total stress is taken into consideration in the constitutive laws through the relation
Advection
The mechanisms of convective heat transfer in porous media considered in the FLAC3D implementation are forced convection, in which the heat is carried by the fluid motion, and free convection, which accounts for fluid motion caused by density differences due to temperature variations.
The working assumptions are those of local thermal equilibrium (flow at low Reynolds number is considered) and small density variations (Boussinesq approximation applies). The basic hydrothermal equations for saturated fluid flow processes are provided below.
Energy Balance for ConvectiveDiffusive Heat Transport
where \(T\) is temperature, \({\bf q}^T\) is thermal flux, \({\bf q}_w\) is fluid specific discharge, \(q^T_v\) is volumetric heat source intensity, \(\rho_0\) and \(c_w\) are reference density and specific heat of the fluid, respectively, and \(c^T\) is the effective specific heat. The effective specific heat is defined as
where \(\rho_d\) and \(C_v\) are solid matrix bulk density and bulk specific heat, respectively, \(n\) is porosity, and \(S\) is saturation.
Fluid Mass Balance (Slightly Compressible Fluid)
where \(M\) is the Biot modulus, \(\epsilon\) is the volumetric strain, and \(\beta\) is the volumetric thermal expansion of the porous matrix. Note that the term \(n \beta_f\) (\(\beta_f\) is volumetric thermal expansion of the fluid), which enters in the expression of the coefficient \(\beta\), is usually neglected in the Boussinesq approximation, and that the last two terms in Equation (21) are only present if mechanical coupling is considered.
Transport Laws
Heat Transport (Fourier Law)
where \(k^T\) is the effective thermal conductivity, which is isotropic in the advection formulation. The effective thermal conductivity is defined in terms of the solid and fluid thermal conductivities, \(k^T_s\) and \(k^T_w\), by the equation
Fluid Transport (Darcy Law)
where \({\bf k}\) is the fluid mobility coefficient (intrinsic permeability, \(\boldsymbol{\kappa}\), over dynamic viscosity, \(\mu_w\)), and \(\rho_w\) is the fluid density. In this equation, fluid density is related to temperature changes by the linear equation
where \(\beta_f\) is the volumetric thermal expansion of the fluid, and \(T_0\) is the reference temperature.
In these equations, the fluid density variations due to temperature changes are significant only in their generation of buoyancy forces (i.e., the Boussinesq approximation).
ThermalMechanicalPore Pressure Coupling
Coupling to the mechanical constitutive equations is done via “effective” normal strain rates (i.e., totalminusthermal normal strainrate components), and pore pressure changes, which affect effective stresses. The advective logic may be used in combination with any of the mechanical constitutive models available for FLAC3D. In the example of an elastic material in small strain, the mechanical constitutive relations are
where \(\sigma_{ij}\) and \(\epsilon_{ij}\) are total stresses and strains, \(\alpha\) is the Biot coefficient, \(K\) and \(G\) are bulk and shear moduli, \(\alpha_t\) is linear thermalexpansion coefficient of the solid matrix, and \(\delta_{ij}\) is the Kronecker delta.
Initial and Boundary Conditions
In a coupled simulation, initial and boundary conditions must be provided for each of the processes involved. In principle, any combination of the fluid and thermal boundary conditions documented in the FLAC3D manual can be used in an advection simulation: Dirichlet and Neuman boundary conditions can be assigned, and point and volume sources can be specified for both fluid and thermal calculations.
An alphabetical list of properties of the advectionconduction thermal model is given here.
Was this helpful? ...  Itasca Software © 2022, Itasca  Updated: Mar 09, 2023 