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:
The infinite thermal boundary is automatically incorporated.
The calculation procedure is extremely quick, compared with finite element or other numerical methods.
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.
The mechanical boundary conditions can be applied correctly, because the standard mechanical logic in FLAC3D is used for the mechanical part of the calculation.
Inhomogeneous and anisotropic mechanical properties can be used.
There are two limitations:
The material is assumed to be thermally homogeneous and isotropic with constant properties.
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)
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
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))
where erfc (\(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)
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 configure 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 property diffusivity=3.2e-6 conductivity=2.4
zone thermal analytical source point 0,0,0 strength 30000 time-start 0
fish define get_temp(stab)
;
; Put temperatures into input table
;
loop foreach local gp gp.list
local radius = math.mag(gp.pos(gp))
if radius >= 2.0
table(stab,radius) = gp.temp(gp)
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, erfc(\(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 configure 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 define 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 components component-number 1 fraction 0.2 component-number 2 fraction 0.8 range group 's1'
zone thermal analytical source components component-number 1 fraction 0.4 component-number 2 fraction 0.6 range group 's2'
zone thermal analytical property diffusivity=3.2e-6 conductivity=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 configure 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 property diffusivity=3.2e-6 conductivity=2.4
zone thermal analytical source point 1,1,1 strength 30000
zone thermal analytical symmetry-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 configure 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 property diffusivity=3.2e-6 conductivity=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 configure 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 property diffusivity=3.2e-6 conductivity=2.4
zone thermal analytical source point 1,1,1 strength 30000
zone thermal analytical isothermal-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 configure 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 property diffusivity=3.2e-6 conductivity=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 configure 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 isotropic
zone gridpoint initialize temperature 100
zone thermal analytical property diffusivity=3.2e-6 conductivity=2.4
zone thermal analytical source point 0,0,0 strength 3000 time-start 0
fish define get_temp(itab)
;
; Put temperatures into input table
;
loop foreach local gp gp.list
local radius = math.mag(gp.pos(gp))
if radius >= 2.0
table(itab,radius) = gp.temp(gp)
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
fish define add_100(stab)
local tab_ = table.find(stab)
loop local i (1, table.size(tab_))
table.y(tab_,i) += 100.0
end_loop
end
[add_100('2e6 s - Analytical')]
; solve to 2e7 s
zone thermal analytical solve time=2e7
[get_temp('2e7 s - FLAC3D')]
[point_source_line(2e7,'2e7 s - Analytical',3000)]
[add_100('2e7 s - Analytical')]
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 configure thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,10,0 point 3 0,0,2 size 1 5 1
zone thermal analytical property diffusivity=3.2e-6 conductivity=2.4
; make points
fish define 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 configure thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,10,0 point 3 0,0,2 size 1 5 1
zone thermal analytical property diffusivity=3.2e-6 conductivity=2.4
zone thermal analytical source line point-1 0.5 0.5 0.5 point-2 1.5 0.5 0.5 ...
number=3 strength 3000 time-start 1e7
zone thermal analytical source line point-1 0.5 0.5 0.5 point-2 1.5 0.5 0.5 ...
number=3 strength 4000 time-start 2e7
zone thermal analytical source line point-1 0.5 0.5 0.5 point-2 1.5 0.5 0.5 ...
number=3 strength 5000 time-start 0
;
zone thermal analytical source line point-1 0.5 0.5 1.0 point-2 1.5 0.5 1.0 ...
number=3 strength 3000 time-start 1e7
zone thermal analytical source line point-1 0.5 0.5 1.0 point-2 1.5 0.5 1.0 ...
number=3 strength 4000 time-start 2e7
zone thermal analytical source line point-1 0.5 0.5 1.0 point-2 1.5 0.5 1.0 ...
number=3 strength 5000 time-start 0
;
zone thermal analytical source line point-1 0.5 0.5 1.5 point-2 1.5 0.5 1.5 ...
number=3 strength 3000 time-start 1e7
zone thermal analytical source line point-1 0.5 0.5 1.5 point-2 1.5 0.5 1.5 ...
number=3 strength 4000 time-start 2e7
zone thermal analytical source line point-1 0.5 0.5 1.5 point-2 1.5 0.5 1.5 ...
number=3 strength 5000 time-start 0
;
;
zone thermal analytical source line point-1 0.5 1.0 0.5 point-2 0.5 1.0 1.5 ...
number=3 strength 3000 time-start 1e7
zone thermal analytical source line point-1 0.5 1.0 0.5 point-2 0.5 1.0 1.5 ...
number=3 strength 4000 time-start 2e7
zone thermal analytical source line point-1 0.5 1.0 0.5 point-2 0.5 1.0 1.5 ...
number=3 strength 5000 time-start 0
;
zone thermal analytical source line point-1 1.0 1.0 0.5 point-2 1.0 1.0 1.5 ...
number=3 strength 3000 time-start 1e7
zone thermal analytical source line point-1 1.0 1.0 0.5 point-2 1.0 1.0 1.5 ...
number=3 strength 4000 time-start 2e7
zone thermal analytical source line point-1 1.0 1.0 0.5 point-2 1.0 1.0 1.5 ...
number=3 strength 5000 time-start 0
;
zone thermal analytical source line point-1 1.5 1.0 0.5 point-2 1.5 1.0 1.5 ...
number=3 strength 3000 time-start 1e7
zone thermal analytical source line point-1 1.5 1.0 0.5 point-2 1.5 1.0 1.5 ...
number=3 strength 4000 time-start 2e7
zone thermal analytical source line point-1 1.5 1.0 0.5 point-2 1.5 1.0 1.5 ...
number=3 strength 5000 time-start 0
;
zone thermal analytical source line point-1 0.5 1.5 0.5 point-2 0.5 1.5 1.5 ...
number=3 strength 3000 time-start 1e7
zone thermal analytical source line point-1 0.5 1.5 0.5 point-2 0.5 1.5 1.5 ...
number=3 strength 4000 time-start 2e7
zone thermal analytical source line point-1 0.5 1.5 0.5 point-2 0.5 1.5 1.5 ...
number=3 strength 5000 time-start 0
;
zone thermal analytical source line point-1 1.0 1.5 0.5 point-2 1.0 1.5 1.5 ...
number=3 strength 3000 time-start 1e7
zone thermal analytical source line point-1 1.0 1.5 0.5 point-2 1.0 1.5 1.5 ...
number=3 strength 4000 time-start 2e7
zone thermal analytical source line point-1 1.0 1.5 0.5 point-2 1.0 1.5 1.5 ...
number=3 strength 5000 time-start 0
;
zone thermal analytical source line point-1 1.5 1.5 0.5 point-2 1.5 1.5 1.5 ...
number=3 strength 3000 time-start 1e7
zone thermal analytical source line point-1 1.5 1.5 0.5 point-2 1.5 1.5 1.5 ...
number=3 strength 4000 time-start 2e7
zone thermal analytical source line point-1 1.5 1.5 0.5 point-2 1.5 1.5 1.5 ...
number=3 strength 5000 time-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 configure thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,10,0 point 3 0,0,2 size 1 5 1
zone thermal analytical property diffusivity=3.2e-6 conductivity=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) number-12=3 number-23=3 ...
strength 3000 time-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) number-12=3 number-23=3 ...
strength 4000 time-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) number-12=3 number-23=3 ...
strength 5000 time-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) number-12=3 number-23=2 ...
strength 3000 time-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) number-12=3 number-23=2 ...
strength 4000 time-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) number-12=3 number-23=2 ...
strength 5000 time-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) number-12=3 number-23=2 ...
strength 3000 time-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) number-12=3 number-23=2 ...
strength 4000 time-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) number-12=3 number-23=2 ...
strength 5000 time-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) number-12=3 number-23=2 ...
strength 3000 time-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) number-12=3 number-23=2 ...
strength 4000 time-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) number-12=3 number-23=2 ...
strength 5000 time-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 configure thermal
zone create brick point 0 0,0,0 point 1 2,0,0 point 2 0,10,0 point 3 0,0,2 size 1 5 1
zone thermal analytical source point 0 0 0 strength 62.8318e3 time-start 0
zone thermal analytical source components component-number 1 fraction 0.2 decay 1.5 ...
component-number 2 fraction 0.8 decay 1.5
zone thermal analytical property diffusivity=.166667 conductivity 5
zone thermal analytical solve time=6
fish define get_temp(stab)
;
; Put temperatures into input table
;
loop foreach local gp gp.list
local radius = math.mag(gp.pos(gp))
if radius >= 2.0
table(stab,radius) = gp.temp(gp)
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
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):
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 configure 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 define 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 elastic
zone property density [dens_] bulk=[bulk_] shear=[shear_]
zone thermal cmodel isotropic
zone thermal property expansion [alpha]
zone thermal analytical property conductivity [cond_] diffusivity [diff_]
; mechanical symmetry
zone face apply velocity-x 0 range position-x 0
zone face apply velocity-y 0 range position-y 0
zone face apply velocity-z 0 range position-z 0
; fixed far boundary (only used for second case)
zone face apply velocity-x 0 range position-x 5
zone face apply velocity-y 0 range position-y 5
zone face apply velocity-z 0 range position-z 5
; define heat source in bottom left corner of octant
; (100 kW for 1000 seconds)
zone thermal analytical source point 0 0 0 strength 1e5 time-start 0
zone thermal analytical source point 0 0 0 strength -1e5 time-start 1000
; results after 1 day
zone thermal analytical solve time=8.64e4
model thermal active off
model solve
fish define 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 configure thermal
model large-strain 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 normal 0 0 -1
; subdivide for mesh refinement
zone densify range position-x 0 3 position-y 6 9 position-z -5.5 -1.5
zone attach by-face
;
; Confine all boundaries to represent symmetry
;
zone face apply velocity-x 0 range position-x 0
zone face apply velocity-x 0 range position-x 10
zone face apply velocity-y 0 range position-y 0
zone face apply velocity-y 0 range position-y 15
zone face apply velocity-z 0 range position-z -11.5
zone face apply velocity-z 0 range position-z 11.5
;
zone cmodel assign elastic
zone property bulk=8.333e9, shear=6.25e9, density=2300
;
zone thermal cmodel assign isotropic
zone thermal property expansion 10e-6
zone thermal analytical property conductivity=2.2 diffusivity=1.1e-6
;
zone initialize stress xx -5e6 yy -5e6 zz -5e6
;
; Histories
;
zone history displacement-z position 7.5,0.4,10.1
;
; Initial equilibrium
;
model thermal active off
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 number-12=21 number-23=21 strength=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 number-12=21 number-23=21 strength=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 number-12=21 number-23=21 strength=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 number-12=21 number-23=21 strength=30e3
zone thermal analytical source components component-number=1 fraction=1 decay= 4e-10
; solve for one year
zone thermal analytical solve time [365*24*60*60]
model solve
model save "repyr1"
References
Abramowitz, M., and I. A. Stegun. Handbook of Mathematical Functions. New York: Dover Publications Inc. (1972).
Carslaw, H. S., and J. C. Jaeger. Conduction of Heat in Solids, 2nd Ed. Oxford: Clarendon Press (1959).
Crank, J. The Mathematics of Diffusion, 2nd Ed. Oxford: Oxford University Press (1975).
Crouch, S. L., andW. C. McClain. “Interim Report on Development of a Semi-Empirical Numerical Model for Simulating the Deformational Behavior of a High-Level RadioactiveWaste Repository in Bedded Salt,” University of Minnesota Report to Oak Ridge National Laboratory, ORNL/TM-5462 (1978).
Dahlquist, G., and A. Bjorck. Numerical Methods. Englewood Cliffs, New Jersey: Prentice-Hall (1974).
Karlekar, B. V., and R. M. Desmond. Heat Transfer, 2nd Ed. St. Paul: West Publishing Co. (1982).
Nowacki, W. Thermoelasticity. New York: Addison-Wesley (1962).
Perrochet, P., and D. Berod. “Stability of the Standard Crank-Nicolson-Galerkin Scheme Applied to the Diffusion-Convection Equation: Some New Insights,” in Water Resources Research, 29(9), 3291-3297 (September 1993).
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Oct 31, 2024 |
