# Analytical Thermal Formulation

## Introduction

The main purpose of this option is to model heat sources embedded in the material. Point heat sources are placed either individually, or in lines and grids, to represent point, line or plane sources. The stress changes caused by temperature changes (due to thermal expansion) are equilibrated by using the standard mechanical FLAC3D code to take a “snapshot” of the transient stresses at different times.

Note that the analytical thermal option in FLAC3D does not provide a general-purpose heat transfer program. Its purpose is specifically oriented to solving design problems associated with nuclear waste disposal.

There are several advantages to this method of calculating the temperatures:

1. The infinite thermal boundary is automatically incorporated.
2. The calculation procedure is extremely quick, compared with finite element or other numerical methods.
3. The calculation of the temperature at any time is independent of the temperature at previous times, so that it takes the same amount of time to calculate the solution regardless of the time at which the solution is required.
4. The mechanical boundary conditions can be applied correctly, because the standard mechanical logic in FLAC3D is used for the mechanical part of the calculation.
5. Inhomogeneous and anisotropic mechanical properties can be used.

There are two limitations:

1. The material is assumed to be thermally homogeneous and isotropic with constant properties.
2. Only a few restricted thermal boundary conditions can be applied (i.e., adiabatic and isothermal planes).

In spite of the limitations of this method of solving thermal-mechanical problems, the thermal coupling in FLAC3D is still a powerful aid when solving problems with heat sources, such as nuclear waste isolation problems.

## Theory

The basic theoretical components of thermomechanical computations in FLAC3D are the determination of the temperature distribution due to the heat sources, and the calculation of the stresses induced by these temperature changes. Calculation of stress changes due to thermal effects is described in Thermal-Stress Coupling. The section below describes the analytical calculation if temperatures.

### The Temperature Due to a Distribution of Heat Sources

The simplest possible heat source to consider is the instantaneous point source in an infinite medium. The temperature distribution resulting from such a source is given by the simple expression (Carslaw and Jaeger 1959)

$T = {Q \thinspace \exp \bigl({-r^2 \over 4Kt} \bigr) \over {{8(\pi K t)}^{3/2} \rho \thinspace C_p}}$

where:

$$Q$$ = quantity of heat in source;

$$\rho$$ = density;

$$C_p$$ = specific heat;

$$K$$ = thermal diffusivity = $$k/ \rho C_p$$;

$$r$$ = distance between source and observation point;

$$t$$ = time after source released; and

$$k$$ = thermal conductivity.

The temperature at a point due to any spatial and temporal distribution of heat sources can be found by considering the distribution as a sum of point sources. In the limit, as the number of sources becomes large, and the strength of each is reduced so that the total heat output is constant, the sum becomes an integral over time, space (a line, surface or volume) or both.

For example, the temperature due to a time-dependent heat source of strength $$Q(t)$$ can be written as

(1)$T = {1 \over {8(\pi K)^{3/2}} \thinspace \rho \thinspace C_p} \ \int_o^t \ {{Q(t') \thinspace e^{-r^2} / 4K (t - t')} \over {(t - t')^{3/2}}} \ dt'$

For the case that the heat source is constant (i.e., $$Q(t)$$ = $$H(t)$$, where $$H(t)$$ is the Heaviside step function, Equation (1) becomes (Carslaw and Jaeger 1959))

(2)$T = {1 \over 4 \pi Kr} \ \hbox{erfc}\ {r \over {(4Kt)^{1/2}}}$

where ferfc ($$x$$) = complementary error function.

The temperature due to an exponentially decaying source of the form $$Q(t) = Q_e \exp(-At)$$ is given by (Crouch and McClain 1978)

(3)$T = {{Q_e \over {4\pi Kr}}\ \exp \bigl({{-r^2} \over {4Kt}} \bigr) \cdot Re \biggl[ w (At)^{1/2} + {{ir} \over {(4Kt)^{1/2}}} \biggr]}$

where:

$$w(z) = {{2iz} \over {\pi}}\ \int_0^{\inf}\ {e^{-t^2} \over {z^2 - t^2}}\ dt,\ Im\{z\} > 0$$

$$Re$${$$z$$} = real part of $$z$$

$$Re$${$$z$$} = imaginary part of $$z$$; and

$$i$$ = (-1)1/2.

Line, plane and volume sources can be written as integrals of a point source over a distance, an area or a volume, respectively. In this way, completely arbitrary distributions of heat sources can be represented by analytical expressions. In FLAC3D, the integration is carried out numerically, so that line sources are represented as a series of point sources along a line segment, and plane sources are represented by a grid of point sources. Thus, to determine the temperature at any point at any instant in time, a sum of terms must be computed, where each term is of the form of the right-hand side of Equation (2) or Equation (3).

### Boundary Conditions

Because the point source solution used in FLAC3D applies to an infinite domain, only a limited number of boundary conditions can be specified. Orthogonal isothermal planes and planes of symmetry can be specified. Symmetric planes are effected by applying an image source of equal magnitude on the opposite side of the plane of symmetry, while isothermal planes require an image source of opposite sign (a heat sink). The effect of having two equal sources is that each produces a flux across the plane, but the fluxes are in opposite directions so that the net flux is zero. Thus, symmetry planes are equivalent to adiabatic boundaries. If the sources are of opposite sign, on the other hand, the net flux is from the heat source to the heat sink. In this case, the two sources each contribute to the temperature at any point on the plane, but because the two contributions are equal in magnitude but of opposite sign, the effect is a net temperature change of zero (i.e., an isothermal condition). Figure 1 illustrates the application of adiabatic and isothermal conditions on different faces, and shows the heat source combinations required to achieve these conditions. Each boundary plane specified doubles the amount of calculation required, so that if three different planes are defined, eight times as much calculation is needed. Figure 1: Adiabatic and isothermal boundaries and their FLAC3D representations.

## Application in FLAC3D

Heat source characteristics such as number of components, decay constants and proportion of each component are input by the user. In addition, the location and initial strengths of arrays of point sources are input. Symmetric and isothermal planes may be specified, as may the initial temperature and the time at which the solution is desired. Once all these quantities have been input, the temperatures are calculated at all gridpoints. Stresses are not calculated at this stage, but when mechanical cycling is performed, the average temperature change for each zone is used to calculate a stress change in each zone. The stress state thus produced will not be in equilibrium, but is treated as an initial condition for the mechanical calculation. The standard FLAC3D mechanical calculation is used to reach equilibrium.

## Input Commands

The commands required for solving thermal-mechanical problems with FLAC3D are given below. The commands are presented in groups corresponding to their function. It is assumed that the reader is familiar with |flac3d| and the mechanical calculations it performs, as well as the input commands.

## FISH Functions

FISH functions that may be used to help with analytical thermal analyses are shown below

## Verification

First, the temperature distribution due to a single point source is compared with an analytical solution. Then, to test the thermal logic thoroughly, different combinations of point, line and grid sources are examined. These results should be expected to agree extremely well with the analytical solutions, since FLAC3D is essentially just evaluating the analytical solutions. These first examples merely confirm that the coding is correct, and illustrate the use of the program.

Finally, the correctness of the mechanical solution is demonstrated by comparing the stress state due to a point source with that of a point source in an infinite medium. This is a true test of the thermomechanical coupling, since the FLAC3D mechanical logic is used to equilibrate the stresses. This test is particularly rigorous because the analytical solution contains a singularity at the source itself, producing steep gradients around the source.

### A Single Non-decaying Point Heat Source

Consider a 30 kW heat source, located at the origin, in an infinite medium with the following properties:

 density 1500 kg/m3 specific heat 500 J/kg/C0 thermal conductivity 2.4 W/m/C0

From these properties, a thermal diffusivity of 3.2 × 10−6 m2/s is calculated.

The following data file is used to calculate the temperature distribution due to this source. The accuracy of thermal calculations in FLAC3D is independent of the mechanical grid, but the temperatures are calculated at gridpoints, so a grid from which the calculated temperature distribution can be extracted is specified.

PointSource.dat

model new
model config thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,2,0 point 3 0,0,10 size 1 1 5
zone thermal analytical prop diffus=3.2e-6 cond=2.4
zone thermal analytical source point 0,0,0 strength 30000 time-start 0

fish def get_temp(stab)
;
; Put temperatures into input table
;
loop foreach local gp gp.list
end_if
end_loop
end

program call "solutions.fis"
zone thermal analytical solve time=2e6
[get_temp('2e6 s - FLAC3D')]
[point_source_line(2e6,'2e6 s - Analytical',30000)]

zone thermal analytical solve time=2e7
[get_temp('2e7 s - FLAC3D')]
[point_source_line(2e7,'2e7 s - Analytical',30000)]

zone thermal analytical solve time=2e11
[get_temp('2e11 s - FLAC3D')]
[point_source_line(2e11,'2e11 s - Analytical',30000)]



Figure 2 shows the results of the FLAC3D calculations compared with the analytical solution, which is obtained using Equation (2). FLAC3D uses an analytical expression for the temperature calculations, so it would be expected that the temperatures would agree exactly. The slight differences are due to the methods of calculating the complementary error function, ferfc($$x$$). This function is not readily available in computer languages, but must be approximated by polynomials. The polynomial used to evaluate the analytical solution differs from that used in FLAC3D. The difference is less than one percent everywhere, except at the last few points for a time of 2e6 seconds, where the temperatures are so small that the difference is not significant. Note that temperatures are only shown at a distance of 2 m from the source because, if placed closer than this, the singularity at the source location would dominate the figure. Figure 2: Temperature distribution around a constant heat source.

### Superposition of Several Non-decaying Sources

This data file illustrates the superposition of four different sources, with different strengths and starting times. Once again, the solution is compared with the analytic solution (below), and, again, the FLAC3D results are in excellent agreement.

SeveralSources.dat

model new
model config thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,2,0 point 3 0,0,20 size 1 1 10

fish def sources
global pos1 = vector(2,2,10)
global pos2 = vector(1,2,10)
global pos3 = vector(2,1,10)
global pos4 = vector(2,2,5)
global str1 = 7000
global str2 = 10000
global str3 = 4000
global str4 = 9000
global ttime1 = 1e7
global ttime2 = 1e7
global ttime3 = 2e7
global ttime4 = 3e7
end
[sources]

zone thermal analytical source point [pos1] strength [str1] time-start [ttime1] group 's1'
zone thermal analytical source point [pos2] strength [str2] time-start [ttime2] group 's2'
zone thermal analytical source point [pos3] strength [str3] time-start [ttime3] group 's1'
zone thermal analytical source point [pos4] strength [str4] time-start [ttime4] group 's2'
;
; use two different source specifications to test summing
zone thermal analytical source component c-num 1 frac 0.2 c-num 2 frac 0.8 range group 's1'
zone thermal analytical source component c-num 1 frac 0.4 c-num 2 frac 0.6 range group 's2'

zone thermal analytical prop diffus=3.2e-6 cond=2.4

; t less than start for all sources
zone thermal analytical solve time=8e6

; t = 2e6 for 2 sources
zone thermal analytical solve time=1.2e7
program call "solutions.fis"
program call "four_sources.fis"
[four_sources(1.2e7)]

; t=2e7 for two, 1e7 for 1, 0 for the other
zone thermal analytical solve time=3e7
[four_sources(3e7)]



Solution for several sources at 1.2e7s

     Position       FLAC3D   Analytical          %error
(   0,   0,  0)     0.217415     0.217415    -4.76841e-09
(   2,   0,  0)     0.234689     0.234689    -4.75302e-09
(   0,   2,  0)     0.263318     0.263319    -4.72827e-09
(   0,   0,  2)      1.29244      1.29244    -4.35278e-09
(   2,   2,  0)     0.284381     0.284381    -4.71239e-09
(   0,   2,  2)      1.59242      1.59242    -4.29714e-09
(   2,   0,  2)      1.40406      1.40406    -4.33167e-09
(   2,   2,  2)      1.73186      1.73186    -4.27511e-09
(   0,   0,  4)      6.03594      6.03594    -3.89927e-09
(   0,   2,  4)       7.6861       7.6861     -3.8184e-09
(   2,   0,  4)      6.63536      6.63536    -3.86921e-09
(   2,   2,  4)      8.47482      8.47482    -3.78632e-09
(   0,   0,  6)      22.5854      22.5854    -3.42289e-09
(   0,   2,  6)      31.0016      31.0016    -3.29656e-09
(   2,   0,  6)      25.4547      25.4547    -3.37809e-09
(   2,   2,  6)       35.328       35.328    -3.24632e-09
(   0,   0,  8)      66.7079      66.7079    -2.97905e-09
(   0,   2,  8)      114.102      114.102    -2.75705e-09
(   2,   0,  8)      79.6213      79.6213    -2.91026e-09
(   2,   2,  8)      145.745      145.745    -2.66811e-09
(   0,   0, 10)      114.102      114.102    -2.75705e-09
(   0,   2, 10)      325.441      325.441    -2.38054e-09
(   2,   0, 10)      145.745      145.745    -2.66811e-09
(   0,   0, 12)      66.7079      66.7079    -2.97905e-09
(   0,   2, 12)      114.102      114.102    -2.75705e-09
(   2,   0, 12)      79.6213      79.6213    -2.91026e-09
(   2,   2, 12)      145.745      145.745    -2.66811e-09
(   0,   0, 14)      22.5854      22.5854    -3.42289e-09
(   0,   2, 14)      31.0016      31.0016    -3.29656e-09
(   2,   0, 14)      25.4547      25.4547    -3.37809e-09
(   2,   2, 14)       35.328       35.328    -3.24632e-09
(   0,   0, 16)      6.03594      6.03594    -3.89927e-09
(   0,   2, 16)       7.6861       7.6861     -3.8184e-09
(   2,   0, 16)      6.63536      6.63536    -3.86921e-09
(   2,   2, 16)      8.47482      8.47482    -3.78632e-09
(   0,   0, 18)      1.29244      1.29244    -4.35278e-09
(   0,   2, 18)      1.59242      1.59242    -4.29714e-09
(   2,   0, 18)      1.40406      1.40406    -4.33167e-09
(   2,   2, 18)      1.73186      1.73186    -4.27511e-09
(   0,   0, 20)     0.217415     0.217415    -4.76841e-09
(   0,   2, 20)     0.263318     0.263319    -4.72827e-09
(   2,   0, 20)     0.234689     0.234689    -4.75302e-09
(   2,   2, 20)     0.284381     0.284381    -4.71239e-09


Solution for several sources at 3e7s

     Position       FLAC3D   Analytical          %error
(   0,   0,  0)      22.4082      22.4082    -3.05422e-09
(   2,   0,  0)      22.9336      22.9336     -3.0465e-09
(   0,   2,  0)      23.3051      23.3051    -3.03619e-09
(   0,   0,  2)      35.6482      35.6482     -2.8699e-09
(   2,   2,  0)      23.8552      23.8552    -3.02842e-09
(   0,   2,  2)      37.4742      37.4742    -2.84756e-09
(   2,   0,  2)       36.751       36.751    -2.85948e-09
(   2,   2,  2)      38.6496      38.6496    -2.83706e-09
(   0,   0,  4)      57.8748      57.8748    -2.68276e-09
(   0,   2,  4)      62.1213      62.1213     -2.6537e-09
(   2,   0,  4)      60.4811      60.4811    -2.66801e-09
(   2,   2,  4)      65.0042      65.0042    -2.63881e-09
(   0,   0,  6)      97.5049      97.5049    -2.49894e-09
(   0,   2,  6)      109.892      109.892    -2.45788e-09
(   2,   0,  6)      105.045      105.045    -2.47655e-09
(   2,   2,  6)      119.036      119.036    -2.43522e-09
(   0,   0,  8)      169.583      169.583    -2.33686e-09
(   0,   2,  8)      222.197      222.197    -2.26957e-09
(   2,   0,  8)      199.528      199.528    -2.29828e-09
(   2,   2,  8)      271.074      271.074    -2.23097e-09
(   0,   0, 10)      237.168      237.168    -2.25835e-09
(   0,   2, 10)      454.242      454.242    -2.14416e-09
(   2,   0, 10)      344.254      344.254     -2.1885e-09
(   0,   0, 12)      169.583      169.583    -2.33686e-09
(   0,   2, 12)      222.197      222.197    -2.26957e-09
(   2,   0, 12)      199.528      199.528    -2.29828e-09
(   2,   2, 12)      271.074      271.074    -2.23097e-09
(   0,   0, 14)      97.5049      97.5049    -2.49894e-09
(   0,   2, 14)      109.892      109.892    -2.45788e-09
(   2,   0, 14)      105.045      105.045    -2.47655e-09
(   2,   2, 14)      119.036      119.036    -2.43522e-09
(   0,   0, 16)      57.8748      57.8748    -2.68276e-09
(   0,   2, 16)      62.1213      62.1213     -2.6537e-09
(   2,   0, 16)      60.4811      60.4811    -2.66801e-09
(   2,   2, 16)      65.0042      65.0042    -2.63881e-09
(   0,   0, 18)      35.6482      35.6482     -2.8699e-09
(   0,   2, 18)      37.4742      37.4742    -2.84756e-09
(   2,   0, 18)       36.751       36.751    -2.85948e-09
(   2,   2, 18)      38.6496      38.6496    -2.83706e-09
(   0,   0, 20)      22.4082      22.4082    -3.05422e-09
(   0,   2, 20)      23.3051      23.3051    -3.03619e-09
(   2,   0, 20)      22.9336      22.9336     -3.0465e-09
(   2,   2, 20)      23.8552      23.8552    -3.02842e-09


### Boundary Conditions

A simple model with a single source at coordinates (1,1,1) is run with a symmetry boundary at the z=0 plane. The result is compared to the same model with no symmetry plane and a “mirror image” source at (1,1,-1). The temperature results are identical.

The same model is run with an isothermal boundary on the x=0 plane and this is compared to a model with a negative source at (-1,1,1). Again, the results are the same. The data files for these two tests are shown below.

symmetry.dat

model new
model config thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,2,0 point 3 0,0,10 size 1 1 5
zone thermal analytical prop diffus=3.2e-6 cond=2.4
zone thermal analytical source point 1,1,1 strength 30000

zone thermal analytical sym-planes z

zone thermal analytical solve time=2e6

[file.all('sym.txt','text') = gp.temp(::gp.list)]

====================================================

; this should give the same answer

model new
model config thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,2,0 point 3 0,0,10 size 1 1 5
zone thermal analytical prop diffus=3.2e-6 cond=2.4

zone thermal analytical source point 1,1,1 strength 30000
zone thermal analytical source point 1,1,-1 strength 30000

zone thermal analytical solve time=2e6

[file.all('sym_two_sources.txt','text') = gp.temp(::gp.list)]


isothermal.dat

model new
model config thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,2,0 point 3 0,0,10 size 1 1 5
zone thermal analytical prop diffus=3.2e-6 cond=2.4
zone thermal analytical source point 1,1,1 strength 30000

zone thermal analytical iso-planes x

zone thermal analytical solve time=2e6

[file.all('iso.txt','text') = gp.temp(::gp.list)]

====================================================

; this should give the same answer

model new
model config thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,2,0 point 3 0,0,10 size 1 1 5
zone thermal analytical prop diffus=3.2e-6 cond=2.4

zone thermal analytical source point 1,1,1 strength 30000
zone thermal analytical source point -1,1,1 strength -30000

zone thermal analytical solve time=2e6

[file.all('iso_two_sources.txt','text') = gp.temp(::gp.list)]


### Initial Conditions

The thermal sources created with the zone thermal analytical source command add temperature contributions on top of any existing temperatures in the model. To test this, the simple one-source model was rerun, but with an initial temperature of 100 degrees applied throughout a model. Alos, the source strength was decreased by an order of magnitude to make the effect of the initial temoerature more obvious.

The result is compared to the analytical soluion plus 100 degrees. The model is behaving correctly as shown below.

initial.dat

model new
model config thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,2,0 point 3 0,0,10 size 1 1 5

zone thermal cmodel assign iso
zone gridpoint ini temperature 100
zone thermal analytical prop diffus=3.2e-6 cond=2.4
zone thermal analytical source point 0,0,0 strength 3000 time-start 0

fish def get_temp(itab)
;
; Put temperatures into input table
;
loop foreach local gp gp.list
end_if
end_loop
end

program call "solutions.fis"

; solve to 2e6 s
zone thermal analytical solve time=2e6
[get_temp('2e6 s - FLAC3D')]

[point_source_line(2e6,'2e6 s - Analytical',3000)]

;add 100 degreees to analytical solution
local tab_ = table.find(stab)
loop local i (1, table.size(tab_))
table.y(tab_,i) += 100.0
end_loop
end

; solve to 2e7 s
zone thermal analytical solve time=2e7
[get_temp('2e7 s - FLAC3D')]
[point_source_line(2e7,'2e7 s - Analytical',3000)] Figure 3: Temperature distribution around a constant heat source with a 100 degree initial temperature applied.

### Lines and Grids of Sources

To verify the operation of the commands to generate lines and grids of sources, one problem is solved with three different sets of data. The first file uses points only, the second uses lines, and the third uses grids. Note how the use of the grid and line commands drastically reduces the amount of input required. It can easily be verified that these three data sets generate the same temperature distribution.

Data Files

points.dat

;
; collection of thermal point sources
;
model new
model config thermal
zone create brick p 0 0,0,0 p 1 2,0,0 p 2 0,10,0 p 3 0,0,2 size 1 5 1

zone thermal analytical prop diffus=3.2e-6 cond=2.4

; make points
fish def make_points
loop local i (1,3)
local x_ = i*0.5
loop local j (1,3)
local y_ = j*0.5
loop local k (1,3)
local z_ = k*0.5
command
zone thermal analytical source point [x_] [y_] [z_] strength 3000 time-start 1e7
zone thermal analytical source point [x_] [y_] [z_] strength 4000 time-start 2e7
zone thermal analytical source point [x_] [y_] [z_] strength 5000 time-start 0
end_command
end_loop
end_loop
end_loop
end
[make_points]

; solve
zone thermal analytical solve time=2.2e7

; dump temperatures to file
[file.all('points.txt','text') = gp.temp(::gp.list)]


lines.dat

;
; collection of thermal line sources
;
model new
model config thermal
zone create brick p 0 0,0,0 p 1 2,0,0 p 2 0,10,0 p 3 0,0,2 size 1 5 1

zone thermal analytical prop diffus=3.2e-6 cond=2.4

zone thermal analytical source line point-1 0.5 0.5 0.5 point-2 1.5 0.5 0.5 ...
num=3 stre 3000   t-start 1e7
zone thermal analytical source line point-1 0.5 0.5 0.5 point-2 1.5 0.5 0.5 ...
num=3 stre 4000   t-start 2e7
zone thermal analytical source line point-1 0.5 0.5 0.5 point-2 1.5 0.5 0.5 ...
num=3 stre 5000   t-start 0
;
zone thermal analytical source line point-1 0.5 0.5 1.0 point-2 1.5 0.5 1.0 ...
num=3 stre 3000   t-start 1e7
zone thermal analytical source line point-1 0.5 0.5 1.0 point-2 1.5 0.5 1.0 ...
num=3 stre 4000   t-start 2e7
zone thermal analytical source line point-1 0.5 0.5 1.0 point-2 1.5 0.5 1.0 ...
num=3 stre 5000   t-start 0
;
zone thermal analytical source line point-1 0.5 0.5 1.5 point-2 1.5 0.5 1.5 ...
num=3 stre 3000   t-start 1e7
zone thermal analytical source line point-1 0.5 0.5 1.5 point-2 1.5 0.5 1.5 ...
num=3 stre 4000   t-start 2e7
zone thermal analytical source line point-1 0.5 0.5 1.5 point-2 1.5 0.5 1.5 ...
num=3 stre 5000   t-start 0
;
;
zone thermal analytical source line point-1 0.5 1.0 0.5 point-2 0.5 1.0 1.5 ...
num=3 stre 3000   t-start 1e7
zone thermal analytical source line point-1 0.5 1.0 0.5 point-2 0.5 1.0 1.5 ...
num=3 stre 4000   t-start 2e7
zone thermal analytical source line point-1 0.5 1.0 0.5 point-2 0.5 1.0 1.5 ...
num=3 stre 5000   t-start 0
;
zone thermal analytical source line point-1 1.0 1.0 0.5 point-2 1.0 1.0 1.5 ...
num=3 stre 3000   t-start 1e7
zone thermal analytical source line point-1 1.0 1.0 0.5 point-2 1.0 1.0 1.5 ...
num=3 stre 4000   t-start 2e7
zone thermal analytical source line point-1 1.0 1.0 0.5 point-2 1.0 1.0 1.5 ...
num=3 stre 5000   t-start 0
;
zone thermal analytical source line point-1 1.5 1.0 0.5 point-2 1.5 1.0 1.5 ...
num=3 stre 3000   t-start 1e7
zone thermal analytical source line point-1 1.5 1.0 0.5 point-2 1.5 1.0 1.5 ...
num=3 stre 4000   t-start 2e7
zone thermal analytical source line point-1 1.5 1.0 0.5 point-2 1.5 1.0 1.5 ...
num=3 stre 5000   t-start 0
;
zone thermal analytical source line point-1 0.5 1.5 0.5 point-2 0.5 1.5 1.5 ...
num=3 stre 3000   t-start 1e7
zone thermal analytical source line point-1 0.5 1.5 0.5 point-2 0.5 1.5 1.5 ...
num=3 stre 4000   t-start 2e7
zone thermal analytical source line point-1 0.5 1.5 0.5 point-2 0.5 1.5 1.5 ...
num=3 stre 5000   t-start 0
;
zone thermal analytical source line point-1 1.0 1.5 0.5 point-2 1.0 1.5 1.5 ...
num=3 stre 3000   t-start 1e7
zone thermal analytical source line point-1 1.0 1.5 0.5 point-2 1.0 1.5 1.5 ...
num=3 stre 4000   t-start 2e7
zone thermal analytical source line point-1 1.0 1.5 0.5 point-2 1.0 1.5 1.5 ...
num=3 stre 5000   t-start 0
;
zone thermal analytical source line point-1 1.5 1.5 0.5 point-2 1.5 1.5 1.5 ...
num=3 stre 3000 t-start 1e7
zone thermal analytical source line point-1 1.5 1.5 0.5 point-2 1.5 1.5 1.5 ...
num=3 stre 4000 t-start 2e7
zone thermal analytical source line point-1 1.5 1.5 0.5 point-2 1.5 1.5 1.5 ...
num=3 stre 5000 t-start 0

; solve
zone thermal analytical solve time=2.2e7

; dump temperatures to file
[file.all('lines.txt','text') = gp.temp(::gp.list)]


grids.dat

;
; collection of thermal grid sources
;
model new
model config thermal
zone create brick p 0 0,0,0 p 1 2,0,0 p 2 0,10,0 p 3 0,0,2 size 1 5 1

zone thermal analytical prop diffus=3.2e-6 cond=2.4

zone thermal analytical source grid point-1 (0.5 0.5 0.5) point-2 (0.5 0.5 1.5) ...
point-3 (1.5 0.5 1.5) n-12=3 n-23=3 ...
stre 3000 t-start 1e7
zone thermal analytical source grid point-1 (0.5 0.5 0.5) point-2 (0.5 0.5 1.5) ...
point-3 (1.5 0.5 1.5) n-12=3 n-23=3 ...
stre 4000 t-start 2e7
zone thermal analytical source grid point-1 (0.5 0.5 0.5) point-2 (0.5 0.5 1.5) ...
point-3 (1.5 0.5 1.5) n-12=3 n-23=3 ...
stre 5000 t-start 0

zone thermal analytical source grid point-1 (0.5 1.0 0.5) point-2 (0.5 1.0 1.5) ...
point-3 (0.5 1.5 1.5) n-12=3 n-23=2 ...
stre 3000 t-start 1e7
zone thermal analytical source grid point-1 (0.5 1.0 0.5) point-2 (0.5 1.0 1.5) ...
point-3 (0.5 1.5 1.5) n-12=3 n-23=2 ...
stre 4000 t-start 2e7
zone thermal analytical source grid point-1 (0.5 1.0 0.5) point-2 (0.5 1.0 1.5) ...
point-3 (0.5 1.5 1.5) n-12=3 n-23=2 ...
stre 5000 t-start 0
;
zone thermal analytical source grid point-1 (1.0 1.0 0.5) point-2 (1.0 1.0 1.5) ...
point-3 (1.0 1.5 1.5) n-12=3 n-23=2 ...
stre 3000 t-start 1e7
zone thermal analytical source grid point-1 (1.0 1.0 0.5) point-2 (1.0 1.0 1.5) ...
point-3 (1.0 1.5 1.5) n-12=3 n-23=2 ...
stre 4000 t-start 2e7
zone thermal analytical source grid point-1 (1.0 1.0 0.5 point-2 )(1.0 1.0 1.5) ...
point-3 (1.0 1.5 1.5) n-12=3 n-23=2 ...
stre 5000 t-start 0
;
zone thermal analytical source grid point-1 (1.5 1.0 0.5) point-2 (1.5 1.0 1.5) ...
point-3 (1.5 1.5 1.5) n-12=3 n-23=2 ...
stre 3000 t-start 1e7
zone thermal analytical source grid point-1 (1.5 1.0 0.5) point-2 (1.5 1.0 1.5) ...
point-3 (1.5 1.5 1.5) n-12=3 n-23=2 ...
stre 4000 t-start 2e7
zone thermal analytical source grid point-1 (1.5 1.0 0.5) point-2 (1.5 1.0 1.5) ...
point-3 (1.5 1.5 1.5) n-12=3 n-23=2 ...
stre 5000 t-start 0

; solve
zone thermal analytical solve time=2.2e7

; dump temperatures to file
[file.all('grids.txt','text') = gp.temp(::gp.list)]


### Decaying Sources

The following data file demonstrates the solution for a decaying heat source.

decaying-source.dat

;
; thermal simulation of a decaying point sources
;
model new
model config thermal

zone create brick p 0 0,0,0 p 1 2,0,0 p 2 0,10,0 p 3 0,0,2 size 1 5 1

zone thermal analytical source point 0 0 0 stre 62.8318e3 t-start 0

zone thermal analytical source components com-num 1 frac 0.2 decay 1.5 ...
comp-num 2 frac 0.8 decay 1.5

zone thermal analytical prop diffus=.166667 cond 5

zone thermal analytical solve time=6

fish def get_temp(stab)
;
; Put temperatures into input table
;
loop foreach local gp gp.list
end_if
end_loop
end

[get_temp('FLAC3D')] ;  FLAC3D results at 6 s

table '3DEC' import '3dec.txt'


The resulting temperature distribution is shown in Figure 4. The solution was checked against 3DEC which in turn was checked againat Equation (3), using tabulated values of $$w(z)$$ found in Abramowitz and Stegun (1972). FLAC3D agrees extremely well with the analytical solution. Figure 4: Temperature distribution around a decaying heat source.

### Thermomechanical Effects

The primary purpose of the thermal logic in FLAC3D is to examine thermomechanical effects, not just thermal effects. To test this capability, the response of an infinite medium to an instantaneous point source is determined.

Nowacki (1962) gives the solution for an instantaneous point source in an infinite medium as

$u_i = \phi,_i$
$\sigma_{ij} = 2G \bigl[\phi,_{ij} - \thinspace\delta_{ij} \phi,_{kk}\bigr]$

where:

$$\phi(R,t) = {mQ \over {4 \pi R}} \ \hbox{erf} \ \bigl[{K \over {(\beta^{1/2}}} \bigr]$$;

$$R$$ = distance from source;

$$t$$ = time;

$$m$$ = $${\alpha (1 + \nu)} \over {(1 - \nu)}$$ ;

$$\nu$$ = Poisson’s ratio;

$$G$$ = shear modulus;

$$\beta$$ = 4 $$Kt$$;

$$K$$ = thermal diffusivity;

$$Q$$ = “source strength”; = $$Q_o / \rho\ C_p$$;

$$Qo$$ = heat released ($$J$$);

$$\rho$$ = density;

$$Cp$$ = specific heat; and

$$\alpha$$ = thermal-expansion coefficient.

From these equations, we can obtain the stresses and displacements along the $$z$$-axis ($$x = y$$ = 0):

$\sigma_{xx} = \sigma_{yy} = {mG \over {2\pi z^3}} \biggl[-V_1 + {2z \over {(\pi)^{1/2}}} \ V_2 + {{4\ z^3} \over {(\pi)^{1/2}}} \ V_3 \biggr]$
$\sigma_{zz} = {mG \over {2\pi z^3}} \biggl[2V_1 - {{4z} \over {(\pi)^{1/2}}} \ V_2 \biggr]$
$\sigma_{xy} = \sigma_{xz} = \sigma_{yz} = 0$
$u_x = u_y = 0$
$u_z = {m \over {2\pi z}} \biggl[{V_2 \over {(\pi)^{1/2}}} - {V_1 \over {2z}} \biggr]$

where:

$$V_1(z,t) = Q\ \hbox{erf} \ \bigl[{z \over {(\beta)^{1/2}}} \bigr]$$;

$$V_2(z,t) = \bigl[{Q \over {(\beta)^{1/2}}} \ \hbox{exp}(-z^2/\beta) \bigr]$$; and

$$V_3(z,t) = V_2(z,t)/\beta$$

The following material properties and source parameters are used to solve this problem with FLAC3D:

 thermal conductivity 1 W/m/°C density 2000 kg/m3 specific heat 500 J/kg/°C linear thermal-expansion coefficient 5 × 10−3 /°C Young’s modulus 1.01 GPa Poisson’s ratio 0.2

The total heat input is 108 J , which in FLAC3D is represented by a 100 kW source acting for 1000 seconds. A cube of side 5 m is used to represent one octant of the medium, as shown in Figure 5. For the times at which the solution is obtained here (1 day and 3 days), the model is large enough to reduce the effect of the boundaries.

The data file for this problem is shown:

thermomechanical.dat

model new
fish automatic-create off
model large-strain off
;
; Data file for thermomechanical verification of FLAC3D
;
model config thermal
;
; define 5 m by 5 m by 5 m region to represent one octant of
; infinite domain
zone create brick point 0 0,0,0 point 1 5,0,0 point 2 0,5,0 point 3 0,0,5 ...
size 15 15 15 ratio 1.1 1.1 1.1

fish def parameters
global bulk_ = 5.6e8
global shear_ = 4.2e8
global dens_ = 2000.0
global alpha = 5e-3 ; thermal expansion coefficient
global Cp = 500.0 ; specific heat
global cond_ = 1.0 ; conductivity
global Q0 = 1e8 ; 100 kW * 1000 seconds

global diff_ = cond_ / (dens_*Cp) ; diffusivity
end
[parameters]

; Define material properties
zone cmodel assign el
zone prop dens [dens_] bulk=[bulk_] shear=[shear_]

zone thermal cmodel isotropic
zone thermal property expansion [alpha]

zone thermal analytical prop cond [cond_] diffus [diff_]

; mechanical symmetry
zone face apply vel-x 0 range pos-x 0
zone face apply vel-y 0 range pos-y 0
zone face apply vel-z 0 range pos-z 0

; fixed far boundary (only used for second case)
zone face apply vel-x 0 range pos-x 5
zone face apply vel-y 0 range pos-y 5
zone face apply vel-z 0 range pos-z 5

; define heat source in bottom left corner of octant
; (100 kW for 1000 seconds)
zone thermal analytical source point 0 0 0 stre 1e5 t-start 0
zone thermal analytical source point 0 0 0 stre -1e5 t-start 1000

; results after 1 day
zone thermal analytical solve time=8.64e4

model thermal active off
model solve

fish def get_flac3d(stab_sxx, stab_syy, stab_szz, stab_dis)
;
; Put stresses and displacements into tables
;
zone.field.method.name = 'average'
loop local i (1,20)
local dist = i*0.25
zone.field.name = 'stress-xx'
table(stab_sxx,dist) = zone.field.get(0.0,0.0,dist)
zone.field.name = 'stress-yy'
table(stab_syy,dist) = zone.field.get(0.0,0.0,dist)
zone.field.name = 'stress-zz'
table(stab_szz,dist) = zone.field.get(0.0,0.0,dist)
zone.field.name = 'disp-z'
table(stab_dis,dist) = zone.field.get(0.0,0.0,dist)
end_loop
end
[get_flac3d('SXX FLAC3D','SYY FLAC3D','SZZ FLAC3D','Disp FLAC3D')]

; fill tables with analytical solution
program call "thermomechanical-solution.fis"
[analytical('SXX Analytical','SYY Analytical','SZZ Analytical','Disp Analytical',8.64e4)]

model save "tmtest-1day"

table 'SXX FLAC3D' clear
table 'SYY FLAC3D' clear
table 'SZZ FLAC3D' clear
table 'Disp FLAC3D' clear
table 'SXX Analytical' clear
table 'SYY Analytical' clear
table 'SZZ Analytical' clear
table 'Disp Analytical' clear

; results after 3 days
zone thermal analytical solve time=2.592e5
model solve

[get_flac3d('SXX FLAC3D','SYY FLAC3D','SZZ FLAC3D','Disp FLAC3D')]
[analytical('SXX Analytical','SYY Analytical','SZZ Analytical','Disp Analytical',2.592e5)]

model save "tmtest-3days"


The zone plot from FLAC3D is shown in Figure 5. The discretization is finer near the source, where stress and displacement gradients are greatest.

Figure 6 shows the comparison between $$\sigma_{zz}$$ calculated by FLAC3D, and the analytical solution near the $$z$$-axis ($$x = y$$ = 0). Figure Figure #thm-form-crrf-f3d and Figure 7 show similar plots for $$\sigma_{xx}$$ and $$\sigma_{yy}$$, while Figure 8 shows the displacements. Figure 5: FLAC3D model for instantaneous point source.  Figure 6: Radial stress ($$S_{zz}$$) on $$z$$-axis ($$x = y$$ = 0) after 1 and 3 days.  Figure 7: Tangential stresses $$S_{xx}$$ (=S_{yy}) along $$z$$-axis ($$x = y$$ = 0) after 1 and 3 days.  Figure 8: $$z$$-displacement along $$z$$-axis after 1 and 3 days.

## An Example Problem – Waste Repository Drift Model

A model of a drift in a hypothetical nuclear waste repository is created. A heat source in the floor of the drift represents a waste canister. The properties of the host rock are:

 thermal conductivity 2.2 W/m/°C density 2300 kg/m3 specific heat 900 J/kg/°C linear thermal-expansion coefficient 10 × 10−6/°C bulk modulus 8.3 GPa shear modulus 6.3 GPa

The initial stress is hydrostatic at 5 MPa, and the heat source is assumed to have only one component with a decay constant of 4 × 10−10 s−1 and an initial strength of 30 kW.

For this example, it is assumed that there are other heat sources (spaced 15 m apart along the length of the drift), and that there are other drifts (spaced 20 m apart), also with heat sources in the floor. The mechanical boundary conditions incorporate this symmetry, but the thermal symmetry is simulated by grids of regularly spaced heat sources (Figure 9). All but one member of each grid are located outside the FLAC3D grid.

The data file to solve this problem is given below, and Figure 10 shows a close-up FLAC3D plot of the model with temperatures and displacements. Figure 9: FLAC3D model of repository drift with grids of thermal sources. Figure 10: Temperatures and displacements on drift center line.

drift.dat

model new
;
; Repository Model
;
model config thermal
model large-strain on

block thermal analytical on
;
; set up tunnel
;
zone create radial-tunnel point 0 0,0,0  point 1 10,0,0 ...
point 2 0,15,0 point 3 0,0,11.5 ...
point 8 3,0,0  point 9 0,0,1.5 ...
point 10 3,15,0 point 11 0,15,1.5 ...
size 10 15 ...
ratio 1 1 1 1.2

zone reflect norm 0 0 -1

; subdivide for mesh refinement

zone densify range pos-x 0 3 pos-y 6 9 pos-z -5.5 -1.5
zone attach by-face
;
; Confine all boundaries to represent symmetry
;
zone face apply vel-x 0 range pos-x 0
zone face apply vel-x 0 range pos-x 10
zone face apply vel-y 0 range pos-y 0
zone face apply vel-y 0 range pos-y 15
zone face apply vel-z 0 range pos-z -11.5
zone face apply vel-z 0 range pos-z 11.5
;
zone cmodel assign el
zone prop bulk=8.333e9, shear=6.25e9, dens=2300
;
zone thermal cmodel assign iso
zone thermal prop exp=10e-6
zone thermal analytical prop cond=2.2 diffus=1.1e-6
;
zone ini stress xx -5e6 yy -5e6 zz -5e6
;
; Histories
;
zone his dis-z pos 7.5,0.4,10.1
;
; Initial equilibrium
;
model solve
model save "repeqm"
;
; After 1 year
;

;
; Define sources to represent a large area around
; modeled drift
;
zone thermal analytical source grid point-1 100,157.5,-2.5 point-2 100,-142.5,-2.5 ...
point-3 -100,-142.5,-2.5 n-12=21 n-23=21 stre=30e3
zone thermal analytical source grid point-1 100,157.5,-3.5 point-2 100,-142.5,-3.5 ...
point-3 -100,-142.5,-3.5 n-12=21 n-23=21 stre=30e3
zone thermal analytical source grid point-1 100,157.5,-4.5 point-2 100,-142.5,-4.5 ...
point-3 -100,-142.5,-4.5 n-12=21 n-23=21 stre=30e3
zone thermal analytical source grid point-1 100,157.5,-5.5 point-2 100,-142.5,-5.5 ...
point-3 -100,-142.5,-5.5 n-12=21 n-23=21 stre=30e3

zone thermal analytical source components comp-num=1 frac=1 decay= 4e-10

; solve for one year
zone thermal analytical solve time [365*24*60*60]
model thermal active off
model solve

model save "repyr1"