UNIVERSITY OF LONDON: KING'S COLLEGE

The prediction of temperatures in mass concrete

by numerical computation
by Professor A. D. Ross, B.Sc., Ph.D., A.M.I.C.E. and J. W. Bray, B.Sc.(Eng.)

SUMMARY: In an appendix the accuracy obtainable by these

The effects in mass structure$ of the chemical numerical solutions is discussed.
generation of heat in the concrete are briefly reviewed.
Some reference is made to the practical techniques Introduction
designed to reduce steep temperature gradients and THE DESIGN AND CONSTR UCTION of mass concrete
hence to control cracking due to thermal causes. structures are largely influenced by the tendency of
There remains, however, the problem of predicting the concrete to crack and the adoption of necessary
temperatures in a mass of concrete of known thermal preventive measures. Cracking is mainly the result
properties which is placed under specified conditions. of a non-uniform temperature distribution in the
Formal mathematics can give solutions only in the concrete brought about by heat evolved' chemically
simplest cases but the number of influential factors in in the cement during its setting and hardening
practice necessitates the employment of approximate processes: the conditions will be aggravated to some
step-by-step methods of numerical calculation. extent by the action of. drying shrinkage but this is
The method proposed by Schmidt for the uni- confined to the surface layers and in large masses
dimensional flow of heat is outlined and illustrated by a is of only secondary importance. The rate at which
numerical solution for the case of a mass, initially at a heat, once generated, is conducted to the exposed
high uniform temperature, which loses heat from two surfaces depends on the thermal diffusivity, which,
parallel boundaries a finite distance apart. The treat- in the case of concrete, has a relatively low value.
ment of a similar mass in which there is generation of When concrete is laid in small quantities or as a
heat is then considered and graphical solutions for these thin slab the heat evolved is conducted away at a
simple cases are also presented. sufficiently rapid rate to prevent the development of
It is shown that a more elegant computational tech- very high temperatures and no ill effect is produced
nique can be achieved when heat is generated by in the concrete. However, when concrete is placed
lowering the boundary temperature and, further, that in large masses, the rate of heat dissipation is slow,
the numerical working can be reduced by a simple much slower in fact than the initial rate of heat
mod({ication of the tabular method of setting out the generation in normal cements. The high temperature
calculations. The utility of these devices is illustrated gradient developed at the surface produces an
by application to the same example for comparison excessive degree of tensile strain in the surface layers
with the original solution. of the concrete. I n the early stages the tendency
The method is then extended to two and three dimen- for cracks to form is relieved by the high creep
sions for which the appropriate values of a non-dimen- property of the concrete, but later, as the concrete
sional number connecting dijjusivity space and time loses its plasticity, cracking becomes inevitable unless
intervals are established. To illustrate the power of the suitable measures are taken. Quite apart from the
step-by-step calculation, the numerical results of a aesthetic point of view, cracking is very undesirable
simple two-dimensional problem are compared with the as it conduces to a rapid disintegration of the surface
rigorous mathematical solution and show excellent and, if it is not restricted to the surface layers, the
agreement. stability of the structure as a whole is likely to be
It is shown that the method can readily be applied to affected.
more complex cases. A solution is presented for the In practical construction a number of special
case of two lifts of concrete placed at different times on precautions may be taken to avoid steep temperature
an underlying rock base which has a thermal capacity gradients and so to control thermal cracking. The
but no heat generation. more usual devices include the use of low-heat or
9

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 Magazine of Concrete Research: January 1949
modified Portland cement, a restriction on the height boundaries parallel to the direction of heat flow are.
of -lift, a rigid, placing schedule, the use of leaner either insulated or too remote from the region under
mixes in the heart of the structure, the pre-cooling consideration to cause any significant distortion of the
of mixing water and aggregates and accelerated flow lines. In explaining the principles of the method
cooling by internal circulating pipes. These pre- it is most convenient to consider the problem of
cautions can be effective in the control of cracking simple cooling first, leaving the consideration of heat
but large quantities of heat are still evolved and generation until later.
basically the problem is altered only in degree. It is The following are the more important symbols
therefore desirable to be able to predict temperatures used in this work :
under any given set of conditions so th/lt the probable
efficacy of these measures of ~ontrol may be assessed
o denotes temperature in 0 F.
k denotes thermal conductivity in B.Th.V. per
before testing them empirically in construction. ft. per day per 0 F.
An investigation of the problem involves in the first c denotes specific heat.
instance a determination of the temperature dis- p denotes density in lb. per cu. ft.
tribution in the concrete at any time after placing. h2 denotes thermal diffusivity in sq. ft. per day.
This may be estimated from a knowledge of the heat
generation characteristic of the cement, the thermal A rod of the material is considered having a uniform
properties of the concrete, the initial and boundary cross-sectional area A sq. ft. and with its axis parallel
conditions and details of the procedure of construc- to the direction of flow. It is imagined to be divided
tion; but the problem is intrinsically a difficult one. up into a number of equal elements of length l::,.x.
Complications arise from the following causes: Figure 1 shows thr,ee adjacent elemerrts p, q and r.
the rate of ,heat generation varies with time in no The temperatures at the centres of the elements
simple mathematical manner, the concrete is built up Op; Oq and Or are assumed to be approximately equal
in successive lifts which are cast at different times at to the mean temperatures and the temperature
probably different temperatures, the ambient tempera- gradients between the centres are assumed to be
ture is continually varying both diurnally and linear.
seasonally, shuttering is left in place for a few days The rate of flow of heat from r to q,
providing a semi-insulating cover, and artificial
cooling may be employed in the first few months to il = k A (fJr- fJq)
l::,.x .
restrict internal temperatures. In a limited number of The rate of flow of heat from q to p.
cases thi!- problem may be solved by rigorous mathe-
matical methods, but usually it is too complex to be i2=~(fJq-OP)
l::,.x
tackled in this way. Step-by-step methods have been
devised by Carlson 1 and Schmidt 2 for dealing with The rate at which heat is absorbed at q,
i=il -i2
heat flow problems, and these provide approximate
but very powerful methods of solution. kA kA
The two methods are very similar in principle but =-(fJr-(jq)- -«iq-(jp)
l::,.x l::,.x
Carlson's method is less convenient in application
and will not be considered here. kA
= l::,.x«(jp+(jr-2(jq)
Schmidt's method is simple but entails somewhat
heavy arithmetic when applied to practical problems. This rate of heat flow is assumed to remain constant
The authors introduce two modifications which reduce during an interval of time 61.
the volume of numerical working. Again, although
practical problems are, in fact, three-dimensional,
they have been treated as cases oflinear flow. In this
paper it is shown that this step-by-step method can be
extended to two and three dimensions.
Schmidt's method for simple cooling
II
~ r-r}
I
I
I
I , ' II
I I
This method was devised for solving problems in
which heat flows in one direction only, i.e. linear
flow. In practice this corresponds to cases where the

1Carlson, R. E. 'A simple method of computation of

temperatures in concrete structures.' J. Am. Cone.
Inst., Vol. 34 (1937), p.89.
2 Schmidt, E. 'Foppls Festschrift', p.179. Springer,
~ 1>11--.+--1,-1>11 ~.-t-I.~
Berlin, 1924. Figure 1.

10

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 711c p/'cdictioll or telll/l('/"(/Iu/,e ill II/(/I".\' cOl/erct'::.

Then the total heat supplied to q in time L,I is i L,I. 1'0 -[---1 ---
This results in a rise of temperature at q of
L, Oq=_~Qt_
1
: !
cp A L,x 0'8r~--i---- ,-
i !

_k ~ (Op+ Or-20q) 0-6!-

cp (L,X)2 I - True distribution
, !
• Calculated points
h 2 L,I W ; i
c£ 0'4~--"-------,--
= --- (Op+Or-2()q) ::> -
(L,X)2
~
«
c£
where jz2 =k is the thermal diffusivity w
Q.
pc
::E
The temperature at q at the end of the time interval w
~----- ,
will then be
.... 3 4 5
-- j

6 7
Oq'= Oq+L, Oq
DISTANCE FROM BOUNDARY - ft.
jz2L,1
=Oq+(Ex)2 (Op+Or-20q) (I) Figure 2.
If L,t is chosen in relation to L,x so that
h 2L,t _ l
(6X)2 - 2 • •• (2)

then Oq'= Oq+t (Op+ Or-20q)

or Oq'=t (Op+Or) .. (3)
This means that the arithmetic mean of the tem- ~
peratures at two alternate points at any instant gives ~'"
..e
the temperature at the intermediate point after an a.
interval of time equal to (6X)2j2h2.
~L-____L -______~______-L______~i~s~to~n~cc~_
The first step in computation is to decide the number
of elements such that the boundaries pass through ba
the centres of the first and last elements. An equal
number of columns is drawn up on the calculation Figure 3:
sheet and in the first line of each column is written calculate the temperatures for one half of the slab,
the initial temperature of the corresponding element. as indicated in Table I.
The second line representing the temperatures after The final temperature distribution is shown by the
the first interval is obtained by averaging alternate plotted points in Figure 2 together with the true
temperatures in the first line. The temperature distribution calculated by rigorous methods. In this
distribution at the end of the second interval is problem the initial temperature distribution suffers a
obtained by averaging the figures in the second· line ; discontinuity at the boundary where the temperature
and so on. This method is illustrated by the following drops suddenly frpm one to zero. The effect of this
problem. is to produce rather poor results in the first few time
intervals, the extent of the error involved being
PROBLEM I indicated by the equal values obtained for the tem-
A very large slab of concrete, 14 ft. thick, is perature at any given point for successive intervals
initially at unit* temperature throughout. It is of time. These errors, however, can be avoided by
desired to determine the distribution of temperature the use of tables and charts which have been prepared
at the end of four days when the boundaries are main- for the case of a semi-infinite solid. It can be shown
tained at zero* temperature. The diffusivity may be that up to a time of about X2j64 days, where X is
taken as 1 sq. ft. per day, which is a good average the distance between the boundaries in feet, the
value for normal concretes. temperature at any point is very nearly the same
Let 6x= I ft.; then 6t=(6x)2j2jz2 =0.5 day. for a finite as for a semi-infinite body. This means
The temperature distribution is symmetrical about tha t, in a practical case, the temperature up to
the central plane and hence it is only necessary to X2
64 days can be obtained from the solution for the
·These arbitrary values are adopted merely for clarity.
The procedure is unaltered when practical values of the semi-infinite case and, thereafter, the calculation
internal and boundary temperatures are inserted. continued by the numerical method discussed above.
I[

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 Magazine of Concrete Research: January 1949

TABLE 1

Distances from one boundary (feet)

TIme

_(d~~_+_+-i_~)!_~__I_(~)_+I_ (:) -_~_ (~) __~__

0'5 ! 0 '500 (1) 1 (I) 1 (1) 1 (1)
1'0 0 ('500) ·750 (1) 1 (1) 1 (I) 1
1'5 0 ·375 (,750) '875 (1) 1 (1) 1 (1)
2'0 0 ('375) '625 ('875) '938 (1) 1 (1) 1
2'5 0 '313 ('625) '781 ('938) '969 (1) 1 (1)
3'0 0 ('313) ·547 (-781) '875 ('969) '984 (1) ·984
3'5 0 '273 (,547) '711 ('875) '930 ('984) '984 ('984)
4'0 0 I
('273) ·492 ('711) '820 ('930) '957 ('984) '957
I
I I I I I

The brackets In this table are for future reference and should be Ignored in this section of the work.
Graphical treatment
An alternative method given by Schmidt for solving
the simple cooling problem is to obtain the mean of
alternate temperatures graphically. Figure 3 shows
the temperature distribution in a material at a given
instant of time; aoa=Op, bob=Oq and coc=Or. Let
1·0

09

0·8
I / ~ ~ ;Z;
I IA W7
- --::--
' ..... ~

""
a straight line joining a and c cut the ordinate bob at d.
Then bob=t (aoa+coc)
0·7
II 1;/ ~
j/

=t (Op+ Or) 06
Ill. Ifjj
= Oq'
Therefore bod represents the new temperature at q os
IJI/i ~
after an interval of time 6t, equal to (6x)2/1h 2 •
The graphical method presents no difficulty and 0-4
,1!l1!
Figure 4 shows a solution of Problem 1 worked out
in this way.
I, If 3

... 2)
Cooling with internal heat generation 0:0·
:::l

:50·II
!;;:
When heat generation is to be taken into account,

~
it is necessary to determine the adiabatic rise in Q.

temperature which would be produced in the concrete .......

:l
.0
by the heat generated during each time interval. o 234 678910
DISTANCE FROM BOUNDARY- ft.
This is added to the individual temperatures in the
tabulation after averaging. The following problem
illustrates the procedure.
Figure 4.
PROBLEM 2
A very large slab of concrete, 20 ft. thick, is initially
at zero temperature throughout. It is desired to find
the temperature distribution at the end of four days Let .6x=1 ft. and .6t=0·5 days, as before.
if, during that time, the boundaries are maintained Successive increments of adiabatic temperature
at zero temperature and the adiabatic temperature will be 11'5, 10'0, 6'5, 5'0,4'0, 3'5,2'5 and 2'0 degrees.
characteristic is as follows : - The calculations are shown in Table 2.

Time (days) .. .. .. 0 0'5 1'0 1'5 2,0 2'5 3'0 3'5 4'0
-------------- ---- ---- ---- ---- - - - - - - - - - ---
Adiabatic temp. (0 F.) .. 0 11'5 21'5 28'0 33'0 37·0 40'5 43·0 45'0

12

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 The prediction of temperature in mass concrete

TABLE 2

I Distance from one boundary (feet)

Time
(days) 0 1 2 3 4 5 6 7 8
----- ----- ---- I -----
o (a l ) 0 0 0 0 0 0 0 0 0
(b l ) 0 0 0 0 0 0 0 0
(c l ) 11'5 11'5 11'5 11'5 11'5 11'5 11·5 1l·5
0'5 (a2 ) 0 11'5 11'5 11'5 11'5 11'5 11'5 11'5 11'5
(b 2 ) 5'8 11'5 11'5 11·5 11·5 11'5 11'5 11'5
10'0 10·0 10'0 10'0 10'0 10'0 10'0 10'0
1'0 0 15'8 21'5 21'5 21'5 21'5 21·5 21'5 21'5
10'8 18'6 21'5 21'5 21'5 21'5 21'5 21'5
6'5 6'5 6'5 6'5 6'5 6'5 6'5 6'5
1'5 0 17'3 25-1 28'0 28·0 28'0 28·0 28'0 28'0
12·6 22'6 26·6 28'0 28'0 28·0 28'0 28'0
5'0 5·0 5'0 5'0 5'0 5'0 5'0 5'0
2'0 0 I 17'6 27'6 31'6 33'0 33'0 33'0 33'0 33'0
13'8 24'6 30'3 32'3 33'0 33·0 33'0 33'0
4'0 4'0 4'0 4'0 4'0 4'0 4'0 4·0
2'5 0 17'8 28'6 34'3 36'3 37'0 37'0 37'0 37·0
14'3 26'1 32'4 35'7 36'6 37·0 37'0 37·0
3'5 3'5 3'5 3'5 3'5 3'5 3'5 3'5
3'0 0 17'8 29'6 35'9 39'2 40'1 40'5 40'5 40'5
14'8 26·9 34'4 38'0 39'8 40'3 40'5 40'5
2'5 2'5 2'5 2'5 2'5 2'5 2'5 2'5
3'5 0 17'3 29'4 36·9 40'5 42'3 42'8 43·0 43'0
14'7 27'1 34·9 39'6 41'7 42'7 42'9 43'0
2'0 2'0 2'0 2'0 2'0 2'0 2'0 2'0
4'0 0 16'7 29'1 36'9 41'6 43'7 44'7 44'9 45'0
ROW a l represents the initial temperature distribution.
" b l is obtained by averaging the temperatures at alternate stations.
" CI gives the increment of adiabatic temperature to be added to the value in row bl> resulting in
" a 2 which is the temperature distribution at the end of the first time interval and at the beginning of the
next.
As before it is only necessary to consider one half us consider also the problem of the same solid in
of the slab. In fact only the first eight feet have been which no heat is generated but in which the boundary
taken, as the computed temperatures in the middle temperature is varied instead. It can be shown that
five feet of the slab are unaffected by the boundary the heat flow in these two cases is exactly the same,
conditions and rise adiabatically. This is in agree- provided that the variation in boundary temperature
ment with the criterion given above for neglecting in the second case is at all times equal to the variation
one boundary since the time period of four days is in adiabatic temperature in the first case, only nega-
X 2 (20)2 tive in sign. If rp is the rise of adiabatic temperature
less than 64 =6"4 at any given time t, then the boundary temperature
=6'25 days. must be decreased by this amount in time t.
The fact that the heat flow is the same in both cases
Modifications of tabular calculations implies that the djfference of temperature between
As can be judged from Table 2, the tabular calcula- any two points at a given time is the same. This
tions may become rather unwieldy and the Authors provides us with an alternative way of solving heat
have developed two modifications which render the generation problems by reference to the equivalent
working more concise. cases of boundary temperature variation. In most
(1) The kind of problem with which we have to cases this results in a considerable saving of work as
deal is one in which the rate of heat generation at a fewer steps are required in the computation; it is
given instant is the same throughout the solid. Let no longer necessary to determine successive incre-
\3

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 Magazine of Concrete Research.' January 1949

TABLE 3

Distances from boundary (feet)

Time
(days) 0 1 2 3 4 5 6 7 8
---- -----
o (a) 0 0 0 0 0 0 0 0 0
0'5 (b) 11'5 0 0 0 0 0 0 0 0
1'0 (b) 21'5 5·7 0 0 0 0 0 0 0
1'5 28'0 10'7 2'9 0 0 0 0 0 0
2'0 33'0 15'4 5'4 1'4 0 0 0 0 0
2'5 37'0 19'2 8·4 2'7 0'7 0 0 0 0
3'0 40'5 22'7 10'9 4'6 1·3 0·4 0 0 0
3'5 43'0 25'7 13'7 6'1 2'5 0·7 0·2 0 0
4'0 (b) 45'0 28'3 15'9 8'1 3'4 1'3 0·3 0'1 0
4'0 (c) 0 16'7 29'1 36'9 41'6 43'7 44'7 44'9 45'0
I
..
ROW (a) represents the InItIal temperature dIstnbutIOn .
" (b) is obtained by inserting the adiabatic temperature in Column 0 and averaging the tempt:nuures at
alternate stations.
" (c) gives the final temperature distribution obtained by subtracting the tabular values from the full
adiabatic rise.
ments of adiabatic temperature nor to add these of one set are enclosed in brackets; the remainder
increments to the temperatures at the end of every constitute the second set. These sets are independent
time interval. The only extra work entailed is in in that any term of one set is entirely dependent for
obtaining the final temperatures by subtracting the its value on the other terms in the same set and is
adiabatic temperature from the tabulated values in unaffected by, and does not itself affect, any of the
the last line. Instead of varying the boundary tem- figures in the other set. Reference to Figure 2 shows
perature in the negative direction it is more convenient one set to be much more accurate than the other.
to omit the minus signs, in which case the final It is therefore reasonable to use the first set of results
temperatures are obtained by subtraction of the only and ignore the second. In the same way Table 3
tabulated values from the adiabatic temperature. consists of two independent sets of calculations but
Table 3 shows Problem 2 solved in this way and this time they do not consist of the same terms. The
indicates the reduction in the figuring. amount of labour may therefore be halved by working
out only one of these sets and ignoring the other.
(2) If Table 1 is examined it is seen to consist of To make this point clear, Problem 2 is worked out
two equal but independent sets of results. The terms on these lines in Table 4 below; comparison of
TABLE 4

Time Distances from boundary (feet)

(days)
o I 1 2 3 4 5 6 7 8
--1----- - - - .---- ---- -----
0 - . 0 0 0 - 0 -
0'5 Il'5 - 0 - 0 - 0 - 0
1'0 - I 5'7 - 0 - 0 - 0 -
1'5 28·0 - 2'9 - 0 - 0 - 0
2'0 - 15'4 - 1·4 - 0 - 0 -
2'5 37·0 - 8'4 - 0'7 - 0 - 0
3'0 - 2')'7 - 4'6 - 0'4 - 0 -
3'5 43'0 - 13'7 - 2'5 - 0·2 - 0
4'0 45'0 28'3 - 8'1 - 1'3 - 0'1 -
4'0 0 16'7 - 36'9 - 43'7 - 44'9 -
II I I

14

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 Tables 3 and 4 will show the set of values omitted.
It is to be noted that at any given time only half the
full number of results are obtained by this method;
normally these results will be quite sufficient to
represent the temperature distribution, but if inter-
mediate values are required, they may be estimated
fairly accurately by taking the mean of the values in
the preceding and succeeding time intervals for the
given point.
The graphical method of solution for the case of
heat generation is also facilitated by the use of
equivalent boundary temperature variation. Problem
2 is worked out graphically in Figure 5; both modifica-
tions have been adopted and consequently the diagram
is the graphical counterpart of Table 4.
T'W' -1110"
More complex cases Figure 6.
Thus far, only simple cases have been considered. Space intervals of 1 ft., corresponding to time
But the methods described can be employed to give intervals of t day, were employed in the computation
solutions in more complex problems as, for example, and the resulting temperature history for the first
the case of a dam built on an underlying rock stratum seven days at a point 1 ft. below the upper surface
and consisting of a number of lifts of concrete placed of the first lift is shown in Figure 6. The temperature
at appropriate intervals. At any given time the rate distribution after 6 days from first placing is shown
at which heat is· generated by the concrete varies in Figure 7. These diagrams are characteristic in
from lift to lift, the rock base generates no heat but that they show the two maxima occurring at a point
allows diffusion to occur, and, of course, the con- at different times and also indicate the severe tempera-
ditions and configuration are altered with the placing ture gradient at an exposed boundary.
of each successive lift. The results for two lifts only are presented here but
To illustrate the application of these methods, the calculations can readily be extended to include
this problem has been solved using the following further lifts and longer time periods. It has be&n
numerical data: assumed above that the heat flow is linear. This is
Height of lifts: 6 ft. the assumption usually made but in consequence the
Interval between placing of each lift: 4 days. solution can only be a rough approximation to the
Diffusivity of rock and concrete: 1 sq. ft. per day. actual temperature history. If a more accurate know-
Adiabatic temperature rise for the concrete as ledge of temperature distribution is required it is
indicated in Figure 6. necessary to consider heat flow in two dimensions at
DISTANCE fROM BOUNDARY - fl.
least. This is considered in the following section.
?---~---} 6 8 9 10
Extension of Schmidt's method
~~
~7~~ ~
i ~
j--t-- 4o
So far we have dealt solely with problems of linear
heat flow, that is, cases in which the temperature is
10
dependent upon only one geometrical co-ordinate x.
1---
~v 1-- --- t - - t--- -- -
r- I3o In the more general case the temperature will also
vary with the co-ordinates in the y oand z directions.
20 I - ""--
1

---~ Problems are usually called one-, two- or three-

1 II~ I
I dimensional according to whether the heat flow is

-+-
20
dependent on one, two or three co-ordinates respec-
o~
I
r--- -- tively. It will now be shown that the method can be
extended to deal with both two- and three-dimensional
10
problems.
40
'/ I
I
Two-dimensional case
I
--- L_ o In this case heat flow takes place in parallel planes
I and the boundaries parallel to these planes are either
0 -'-----
insulated or are too remote to affect the temperature
distribution in the region of the solid under con-
Figure 5. sideration.
15

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 Magazine of Concrete Research: January J949
40

f 1\
30
::~I: :::::
-
'Onere
r:....
L
:/
II \
\
l.
0

/
V L r--
::-I
'" 0
r-- ... ... IL I t.Lifl. -1~Lltt
l---t
...
gl

.
0:
'"
V \
.'"
~
0
10
V
.......
" 2 0 2 4 6 8
DISTANCE ABOVE & BEL?W FORMATION LEVEL -f~
to 12

TEMPERATURE DISTIlIBUTION AFTER 6 DAYS,

~Ax---o,-tl-- ---+~-
Figure 7
Let heat flow take place in planes parallel to the AX AX

xy-plane and consider a slab of material bounded by

two such planes at a distance z apart. Imagine this Figure 8.
slab to be divided into a number of square elements, and the amount of work Involved somewhat heavier.
the sides being parallel to the x and y axes and of The calculation sheet is divided up into rectangles
length l:::,.x. One such element q is shown in Figure 8 of a convenient size representing the elements in order.
with its four adjacent elements m, n, p and r; let their The initial temperatures are written at the top of each
mean temperatures at a given instant be ()q, ()m, ()n, rectangle and the subsequent temperatures obtained
()p and ()r.
by averaging are written down in order underneath.
Using assumptions and reasoning similar to that The modifications to the method in the linear case
in the linear case, the resultant rate of heat flow into q, may be applied equally well nere; in fact, the saving
i=i1 - i2 +ia - i4 of work secured by their use is likely to be much
=
l:::,.x l
f
k (zl:::,.x) «()r- ()q)-( ()q- ()p) +( ()n- ()p)-( ()q- ()m)}
greater.
PROBLEM 3
Find the temperature distribution in a sq,uare
=k z «()p+ ()r+ ()m+ ()n-4()q) section block of concrete of 20 ft. side when cooled
After an interval of time l:::,.t, the rise of tempera- for eight days from an inttial uniform unit tempera-
ture at q, ture with the boundaries maintained at zero tempera-
il:::,.l ture. ,
6 ()q cpZ(l:::,.X)2 Let 6x=2 ft. and with h2 =1 ft2 per day as before
then l:::,.t=t 22=1'0 day.
There are four lines of symmetry as shown in
= k l:::,.1 «()p+ ()r+ ()m+ ()n--4()q) Figure 9 and hence temperatures need only be
Cp(l:::,.X)2
calculated for one-eighth of the whole section, i.e.
h2 6t the part shown hatched in the figure. The calculations
=(6X)2 «(Jp+(Jr+(Jm+(Jn-4()q)
are given in Table S.
= N «()p+ ()r+ ()m+ ()n-4()fj) In this Table each space represents a station spaced
2 at 2 ft. from adjacent stations in both directions and
h l:::,.t a non- d"ImenSIOnaI num ber.
where N= (6X)2' the temperatures in successive time intervals are
given in successive lines. The zero figure in row (I)
The temperature at q at the end of time interval 6t and column (f) are the boulldary temperatures and the
will therefore be unit figures at the head of spaces (a)(l), (c)(l),
()q' = ()q+ l:::,. ()q
(e)(l), (b)(2), etc., represent the initial temperatures.
=N«()p+ ()r+ ()m+ ()n)+(1-4N)()q (4) In subsequent lines the figures give the result of
Then, if N=t, averaging the four surrounding temperatures, for
()q' =1, «()p + ()r + ()m + ()n) (5) example, 0'75 in line 2 of (b)(2) is the result of
The computational procedure in this case is averaging the temperatures in (b)(l), (b)(3), (a)(2)
obviously more complicated than in the linear case, and (c)(2); that is, the average of 0, 1, 1, and I.
16

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 The prediction 0/ temperature in mass concrete

TABLE 5

(a) (b) (c)

---- ---- - - - -------
(d) (e) I
([)

0 0 0 o 0 o (1)
---- - - - - - - - ---- - - - ~--

1'00 - 1'00 - 1'00

- 0·75 - 0'75 - x
0'62 - 0'62 - 0'37
- 0'55 - 0'47 -
0'49 - 0·47 - 0'23 0 ( 2)
- 0'44 - 0'34 -
0'41 - 0'38 - 0'17
- 0'38 - 0'27 -
0'36 - 0'31 - 0'13
1 - - - - - - - - - - - - ---- - - - ---
- 1'00 - 1'00
1·00 - 1'00 - Figure 9.
- 0·94 - 0'87
Again, for example, the figure 0'82 in line 7 of (c)(4)
0'87 - 0'86 - is obtained by averaging 0'92 in (b)(4), 0.92 in (c)(5),
- 0'82 - 0'66 0 (3)
0'77 - 0.72 - 0'72 in (c)(3) and 0·72 in (d)(4), the figures in (c)(5)
- 0'71 - 0'53 and (d)(4) being obtained from symmetry.
Isotherms obtained from the final results shown in
0·69 - 0·61 - Figure 8, are indistinguishable from those obtained
- 0'63 - 0'44
- - - - - - ----- --- by more rigorous methods *.
1'00 - 1'00 Three-dimensional case
- 1'00 - This is the most general case and no restrictions
1'00 - 1'00 on the boundaries are implied. For the purposes of
- 0·98 - computation it is imagined that the solid is divided
0·96 - 0·92 0 (4) into cubes of side DX, the sides being parallel to the
- 0·92 - xy, yz and zx planes. This time, each element q has
0'90 - 0·82 six adjacent elements, labelled p, r, m, n, v and w
- 0'85 - as in Figure 10. By the same kind of argument used
0'83 - 0'73 above it can be shown that
1 - - - - - - - - - - - - - - - - --- ()q' =N«()p+ ()r+ ()m+ (In+ (Jv+ ()w) +(l-6N) ()q, (6)
- 1'00
1·00 - and if N=~
6'
- 1'00
1.00 - ()ql=~«()p+()r+(Jm~(Jn+(Jv+()w) .......... (7)
- 0'99 0 ( 5)
It appears that the best method of setting out the
0'99 - calculations is to have a number of sheets drawn up
- 0'95
into rectangles as in the previous case, each sheet
0'95 -
representing a different xy plane.
- 0'90
In conclusion it may be said that the method is
- - - - - - - - - - - - - - - - - --- very flexible ana additional factors such as a fluctuat-
1-00
ing boundary temperature and a non-uniform initial
- distribution of temperature in the concrete can be
1'00
incorporated in the computations without serious
-
(6) difficulty. Nevertheless, even with the modifications
1'00 0
- suggested, the arithmetic is likely to be somewhat
0'99 arduous particularly if two- or three-dimensional flow
- is considered. For this reason, the development of an
0.95 *See 'Carslaw and Jaeger. .. Conduction of Heat in
Solids." Figure 17, p.IS3.

17

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 Magazine of Ccmcrele Research: January 1949

y directions in the manner of Figure 10, the distance of

all points from q being !::,x. Further, let OqH and Oq'

~.
represent the temperatures at q, !::,t sees before, and
!::,t sees after, a given instant. Then equation (8)
becomes, when expressed in finite difference form,

z Op+ Or-20q Om+ On-20q Ov+ Ow-2l}q

(!::,X)2 + (!::'X)2 +
(!::'X)2
I Oq' - Oq'
h2 2!::,t
Equations (1), (4) and .(6) on the other hand give
Op+Or-20q _ I Oq'-Oq
p (one-dimensional)
(!::,X)2 - h2 ------z:"t
Op+ Or+ Om+ On-40q Oq'- Oq
w (!::,X)2 h2------;s:t
(two-dimensional)
and
Op+ Or+ Om+ On+ Ov+ Ow-60q I Oq'- Oq
Figure 10.
(!::,X)2 - = iii ----;s:r-
(three-dimensional)
electrical analogue is now in hand. It is too early as Comparing these with the finite difference equation
yet to report this work but it is believed that the we see that there is a difference in the expression for
device will give rapid solutions to a wide variety of
thermal problems in mass concrete. o Q: moreover from the purely mathematical point
?il
APPENDIX 1 of view the expression given by the Schmidt type of
equation is much the poorer approximation. On this
Accuracy of numerical calculations argument it would seem preferable to calculate tem-
It is instructive to compare the general equations, j::eratures on the basis of the finite difference equation.
(1), (4) and (6) above with those obtained by the use This, however, leads to a more complex expression
of finite differences. Stated mathematically, all for Oq', resulting in a less practicable method of
simple cooling problems amount to the solution of computation. Nevertheless the foregoing remarks
the differential equation concerning the approximation involved in Schmidt's
02 0 02 0 02 0 1 00 method necessitates an enquiry into the order of
- + "Cly2
ox 2
~-- + -
Cl 2
= ~ -
h 2 ot
........ (8)
accuracy to be expected from calculations using this
Z
for given initial and boundary conditions. method.
Let q represent any point with an array of points For this purpose it is essential to confine our inveSll-
q, r, m, n, v and w disposed around it in the x, y, and z gation in the first place to the simple and yet quite
TAB LE 6
Distance from boundary
Time
0 !::,x 2!:::..x 3!:::..x 4!:::..x 5!:::..x 6!:::..x .. . m!::,x
- - - - - - - - - - - - - - - ----- -----
0 0 (1 ) 1 (1) 1 (1) 1
I!:::..t 0 '5 (1) 1 (1) 1 (1)
2!:::..t 0 ('5) '75 (1) 1 (1) 1
3!:::..1 0 '375 ('75) '875 (1) 1 (1)
46t 0 ('375) '625 ('875) '9375 (1) 1
5!:::..t 0 '3125 ('625) '7813 ('9375) ·9688 (1)
6!:::..1 0 ('3125) '5469 ('7813) '8750 ('9688) '9844

n!:::..t

18

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 fundamental case 'Of the CQQling 'Of a semi-infinite 1-0 ~--'-~~'--~~~~=-~I

SQlid frQm a given uniform temperature, the bQun- ~r- I'

0-91---+--+--+-~+-,..q -~-t----j--+--+----j
daries being held at zero temperature. To keep the ,Y'
reasQning as general as PQssible, no specific time 'Or 08 / i
space intervals are chQsen and all work is carried 'Out
'On a non-dimensional basis. The result 'Of the applica-
0'71-_+-~t-il~- / . I f ,~-~
tion of Schmidt's method is shQwn in Table 6, the
06 -
temperature values being given as fractiQns 'Of the
'Original temperature' which will be referred tQ as 1- O-SI---+--It--~I-
Om... Six space and six time intervals have teen
chQsen. 0-41- ----1-/'-+-+
FQr comparison with the true result shQwn in
Figure 11, non-dimensional plotting is used, values 'Of 03 - Lt
! -- - --L+ --1
iii
0/0 max. being plQtted on a base of xl2yh 2 t, where
0-2 / :t--+
x is the distance measured from the bQundary and I
is the time from start 'Of cOQling.
NQW if m =the number of space intervals from the
bQundary
0y o 05
-r-t-i--t-
10() 15 2-0
I
2-5
, m
n =the number 'Of time intervals from the 2Jh1i or An
start of cooling
Figure 11.
then x =m fix and t=n fit
In the case 'Of heat generatiQn, it is assumed in
x m fix Schmidt's methQd that the heat is suddenly eVQlved
··2Yh'2t 2yh2nfit at the end 'Of each time interval. If this were actually
the case then by superpQsitiQn the same remarks as
1_ ~ ~ tQ accuracy would apply here aISQ. Since heat
2 vi h;;;;t y n generatiQn is in fact a cQntinuQus prQcess it is not
PQssible tQ use the previQus results tQ draw any
1 - m cQnclusiQn as tQ the error in calculatiQns invQlving
- y2 ----,-
2 yn heat generation. There appears tQ be nQ a priori
m reaSQn fQr including the heat generation term after
the prQcess of averaging; it might just as well be
y2n introduced befQre and in fact it might be sUPPQsed
This is then the equivalent base fQr the plQtting 'Of that the most accurate results would be 'Obtained by
the tabulated results. Thus, for example, we find including half the heat generated before and half
that ~. after the averaging process. The calculation of par-
o= max
0'5 when m = and n = 2, i.e.
ticular cases shQWS clearly, hQwever, that the first
m 1 methQd gives the better result. TQ 'Obtain some idea
when ---== = 0'5. of the errQr invQlved, the case 'Of a semi-infinite solid
y2n y2X2
FQr ease of reference the bracketed values in is taken, where the rate of heat generation (R) is
Table 6, which corresPQnd tQ 'Odd values 'Of the sum constant. Again, we adopt a nQn-dimensional
(m+n), are plQtted as circles in Figure 11, whereas apprQach to the problem, plotting values of
the remaining values are represented by dots. It is ~ )on a base Qf_xco=
( Rt ('Or _m )
clear that the fQrmer set of results show little diver- 2yh 2 t Y2n.
gence from the true curve, the error not exceeding In wQrking out Table 7, R has been taken for
abQut 4 per cent. at any point, and getting much cQnvenience equal tQ (II fi 1). The apprQximate values
smaller as n increases. together with the true curve are plotted in Figure 12.
It can be shQwn that in the case of a parallel slab It will be seen that, with the exceptiQn of the two
'Of finite width the temperature values given by unit values cQrresPQnding to m = I, n = 1 and m = 2,
Schmidt's method are compQsed 'Of a number 'Of n=2, the error invQlved nQwhere exceeds 4 per cent.
values given by this methQd fQr the semi-infinite Here we have only considered the case of linear
case which we have just considered. Hence, in this flow when heat is eVQlved at a constant rate corres-
case alsQ, if we always take the figures cQrresponding ponding tQ a linear adiabatic characteristic. Nothing
tQ 'Odd values 'Of (m+n), the maximum error will definite is inferred as tQ the accuracy tQ be expected
never exceed 4 per cent. and will prQbably be much in the general case but it seems prQbable that, provid-
less. ing the time intervals chQsen (with the exception 'Of

19

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.
 TABLE 7

Distance from boundary

T

__ '~'__ II~= :==1= ~X T 2~X T 3~X =_4~X == 5~X =6~X _ - .•. m/ix

IL,t I 0 I 1 1 1 1 1 1
2L,t I 0 1'5 2 2 2 2 2
3L,t I, 0 I 2'0 2'75 3 3 3 3
4L,t II 0 i 2'375 3'5 3'875 4 4 4
5L,t II 0 I 2'75 4'125 4·75 4'9375 5 5
6L,t I 0 3'0625 4'75 5'531 5'875 5'9687 6
I

I
the first two time intervals) are sufficiently small for an overall saving of work.
the heat generation characteristic to be well repre- Furthermore, in the two- and three-dimensional
sented, the error is nowhere likely to exceed 4 per cent. cases, the overall dimensions in the three directiol\s
are often very different and it is useful fo adopt
APPENDIX 2 correspondingly different values of the elemental
Other values of N lengths in these directions, say, L,x, L,y and L,z.
It will be noticed that the method as so far pre- For a given value of L,t, this will imply different
sented depends on fixed relationships between space values of N in the three directions, say, Nx, Ny and Nz
and time increments L,x and L,', the actual relation- Then the general equations for Oq' become
ships depending upon whether the problem is one-,
Oq'=Nx(6p+Or)+(l-2Nx) 6q
two- or three-dimensional; this means that, given L, t
Oq'=Nx(6p+ ()r)+Ny(Om+ On)+[1-2(Nx+Ny)] {}q
then L,x must follow, or alternatively, given L,x
Oq' = Nx(Op + OrHNy«()m+ ()n)+Nz«()v+ ()w)
then L,t must follow. ]n particular cases this may be
+il-2(Nx+Ny+Nz)] ()q.'
a source of inconvenience, the number of either space
or time intervals required being excessive. The method If it is considered desirable, Nx, Ny and Nz may be
may be made more elastic in this respect by choosing chosen so that Oq is eliminated from the calculations.
N (h2L,t)
= Cf-,ix)2
.
to gIve . bl e space an d tIme
sUlta .

intervals, ignoring the specific values required by the

Schmidt method. The general equations for Oq'
given above, are:-
(I) Oq' =N (Op+ Or) + (l--2N) Oq 0·7 e----+---.'f--+.-+---+
(one-dimensional)
(2) Oq'=N(Op+Or+Om-! On) + (l-4N) Oq
•
(two-dimensional) Rt. o.s ~

(3) Oq'=N(f)p+ Or+ Om+ On+ Or+ Ow)+(l-6N)Oq

(three-dimensional)
Now, for example, if we make N=O'I, then these
equations give
-- --I

(I) Oq'=O'1 (Op+Or);t-0'8 Oq

(2) Oq' =0,1 (Op + Or+ Om + On) + 0·6 Oq
(3) Oq'=O'1 (Op+ Or+ Om i- On+ Ov+ Ow)+0'4 Oq.
Obviously this modification complicates the com- • m
putional procedure and increases the amQunt of work 2J1:f"t or r.
involved in the calculations: but, depending on the
requirements of any given problem, there may be Figure 12.

20

Downloaded by [ Mr Michel Di Tommaso] on [16/08/18]. Copyright © ICE Publishing, all rights reserved.