Mathematical Model Description

Conventions and Definitions

As a notation convention, the symbol ai denotes component i of the vector [a] in a Cartesian system of reference axes; Aij is component (i,j) of tensor [A]. Also, f,i is used to represent the partial derivative of f with respect to xi. (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

(1)Lc=volume of solidsurface area exchanging heat

Thermal diffusivity

(2)κ=kρ Cv
where: k is the thermal conductivity;
  ρ is the density; and
  Cv is the specific heat at constant volume.

Characteristic time

(3)tc=L2cκ

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 energy-balance equation and transport laws derived from Fourier’s law of heat conduction. Substitution of Fourier’s law into the energy-balance equation yields the differential equation of heat conduction, which may be solved for particular geometries and properties, given specific boundary and initial conditions. Thermal volumetric-strains are introduced into the incremental mechanical and fluid constitutive laws to account for thermal-stress and thermal-pore pressure couplings.

Energy-Balance Equation

The differential expression of the energy balance has the form

(4)qi,i+qv=ζt

where qi is the heat-flux vector in [W/m2], qv is the volumetric heat-source intensity in [W/m3], and ζ is the heat stored per unit volume in [J/m3]. In general, temperature changes may be caused by changes in both energy storage and volumetric strain, ϵ, and the thermal constitutive law relating those parameters can be expressed as

(5)Tt=Mth(ζtβthϵt)

where Mth and βth are material constants, and T is temperature.

In FLAC3D, we consider a particular case of this law for which βth = 0 and Mth=1/(ρCv). ρ is the mass density of the medium in [kg/m3], and Cv is the specific heat at constant volume in [J/kg C]. The hypothesis here is that strain changes play a negligible role in influencing the temperature, a valid assumption for quasi-static mechanical problems involving solids and liquids. Accordingly,

(6)ζt=ρ CvTt

Substitution of Equation (6) in Equation (4) yields the energy-balance equation

(7)qi,i+qv=ρCvTt

Note that, for nearly all solids and liquids, the specific heats at constant pressure and at constant volume are essentially equal. Consequently, Cv and Cp can be used interchangeably.

Transport Law

The basic law that defines the relation between the heat-flux vector and the temperature gradient is Fourier’s law. For a stationary, homogeneous, isotropic solid, this constitutive law is given in the form

(8)qi=kT,i

where T is the temperature [C], and k is the thermal conductivity in [W/mC].

Isotropic Heat Conduct

For isotropic conduct, the thermal conductivity matrix has the form:

(9)k=k[100010001]

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:

(10)k=[kxx000kyy000kzz]

For this case, three conductivity parameters corresponding to different directions are required.

In general,

(11)k=[kxxkxykxzkxykyykyzkxzkyzkzz]

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 qi 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 heat-flux 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

(12)qn=h(TTe)

where qn is the component of the flux normal to the boundary in the direction of the exterior normal, h is the convective heat-transfer coefficient [W/m2$$C], T is the temperature of the boundary surface, and Te is the temperature of the surrounding fluid [C]. Note that in the numerical formulation used in FLAC3D, boundaries are adiabatic by default.

Mechanical Coupling — Thermal Strains

Solution of thermal-stress problems requires reformulation of the incremental stress-strain 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 shearing-strain increments are unaffected. The thermal-strain increments associated with the free expansion corresponding to temperature increment ΔT have the form

(13)Δϵij=αtΔTδij

where αt [1/C] is the coefficient of linear thermal expansion, and δij is the Kronecker delta.

The differential equation of motion and the rate of strain-velocity 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 pore-pressure change caused by thermal expansion of the saturated matrix. The constitutive law has the form (see this this equation in FLAC3D Fluid-Thermal-Mechanical-Formulation — Mathematical Description)

(14)pt=M(ζtαϵt+βTt)

where ζ is the variation of fluid content, M is the Biot modulus, α is the Biot coefficient, ϵ is the volumetric strain, and β is the undrained volumetric thermal expansion. In FLAC3D, the pore-pressure increment Δpther, caused by a temperature change ΔT, has the form

(15)Δpther= MβΔT

β [1/C] is defined as the pore-pressure variation divided by (3α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 thermal-expansion coefficients for the grains, 3αt, and the fluid, 3αf, through the formula

(16)β=3[αt(αn)+αfn]

where n = porosity.

For an incompressible solid constituent, (α = 1 and M=Kf/n), the pore-pressure increment becomes

(17)Δpther= KfnβΔT

where Kf = fluid bulk modulus.

The corresponding change in total stress is taken into consideration in the constitutive laws through the relation

(18)σij= Δptherδij

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 Convective-Diffusive Heat Transport

(19)cTTt+qT+ρ0cwqwTqTv=0

where T is temperature, qT is thermal flux, qw is fluid specific discharge, qTv is volumetric heat source intensity, ρ0 and cw are reference density and specific heat of the fluid, respectively, and cT is the effective specific heat. The effective specific heat is defined as

(20)cT=ρdCv+nSρ0cw

where ρd and Cv are solid matrix bulk density and bulk specific heat, respectively, n is porosity, and S is saturation.

Fluid Mass Balance (Slightly Compressible Fluid)

(21)Pt=M(qwαϵt+βTt)

where M is the Biot modulus, ϵ is the volumetric strain, and β is the volumetric thermal expansion of the porous matrix. Note that the term nβf (βf is volumetric thermal expansion of the fluid), which enters in the expression of the coefficient β, 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)

(22)qT=kTT

where kT 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, kTs and kTw, by the equation

(23)kT=kTs+nSkTw

Fluid Transport (Darcy Law)

(24)qw=k(Pρwgx)

where k is the fluid mobility coefficient (intrinsic permeability, κ, over dynamic viscosity, μw), and ρw is the fluid density. In this equation, fluid density is related to temperature changes by the linear equation

(25)ρw=ρ0[1βf(TT0)]

where βf is the volumetric thermal expansion of the fluid, and T0 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).

Thermal-Mechanical-Pore Pressure Coupling

Coupling to the mechanical constitutive equations is done via “effective” normal strain rates (i.e., total-minus-thermal normal strain-rate 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

(26)σijt+αPtδij=2G(ϵijtαtTtδij)+(K23G)(ϵkkt3αtTt)δij

where σij and ϵij are total stresses and strains, α is the Biot coefficient, K and G are bulk and shear moduli, αt is linear thermal expansion coefficient of the solid matrix, and δ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 advection-conduction thermal model is given here.