Reference / Paper · 1969

Partial Differential Equations for Heat Conduction Analysis of Frozen Layer Shifting

Technical Information Series No. 11 from Hitachi demonstrates how the heat-conduction partial differential equation governing the shifting of a frozen soil layer is converted to a set of finite-difference equations and solved on a Hitachi analog-hybrid computer. The document derives boundary conditions for a one-dimensional frozen/unfrozen soil column, scales the equations to machine variables, and presents temperature-transition curves (time and depth as parameters) obtained from the simulation. A full block diagram for the HITAC 505 High-Precision Analog Computer implementation is included, showing 7 integrators, 19 summing amplifiers, 12 sign-changers, 7 potentiometers, 2 dead-zone units, and 15 diodes.

Manufacturer: Hitachi
System: HITAC 505 High-Precision Analog Computer
Year: 1969
Type: Reference / Paper
Language: English
Learning track: specific applications
Pages: 8

HITAC 505 High-Precision Analog Computer
Hitachi
heat conduction
partial differential equations
finite difference methods
frozen layer / permafrost

← Back to the Reference Library

HITACHI Analog-eHybrid Computer Technical Information Series No.11 Partial Differential Equations for Heat Conduction Analysis of Frozen Layer Shifting. 1969 Hitachi, Ltd. Partial Differential Equations for Heat Conduction Analysis of Frozen Layer Shifting In analyzing partial differential equations for heat conduction on the analog computer, the boundary face between two different substances has been usually assumed not to shift. In the present example, however, since temperature distribution in the process of freezing is expressed as a function of time, the boundary face between frozen and non-frozen layer is unintentionally displaced with time, presenting a very interesting way of analysis by the analog computer. 1. Physical Conditions In the process of freezing, in which the frozen layer proceeds with time, shifting of the boundary face between frozen and non-frozen layers occurs, and the position of boundary face and temperature distribution are sought for as a function of time. In this case, heat conduction is assumed to occur unidirectionally. éti a0. kK Ote _ ge Ke O2ti 1 0 xX? a2 te fe) Initial conditions For where 6=0; At (o<x<) (1) (@< x) (2) Froze Heat flux. Hi ti =te =Q; (3) 0 (4) <t Fig. 1 , 6 Hi @ pe - ow tart; j at at py 98 “Oy a nr t; =te -Oo te =C, thermal conductivity of ice thermal conductivity of water specific heat of ice specific heat of water density of ice density of water initial temperature atmospheric temperature thermal conduction coefficient of air heat of fusion heat of fusion time temperature of frozen laver temperature of non-frozen laver distance from the surface of frozen layer to the boundary face thermal diffusibility of ice thermal diffusibility of water heat density distance from the surface of frozen layer n tayer Non-frozen layer te ix, & 1.05 Keal/m hr °C 0.344 Keal/m hr°C 0.475 Keal/kg °C 0, 889 Keal/kg °C 980 kg/m? 1,020 kg/m? +3,9°C -20°C 5.3 Keal/m? hr °C 65.2 Keal/kg 80.0 Keal/kg hr °C °C Keal/m? m The freezing conditions are assumed as follows: the medium is in contact with the atmosphere of temperature -20°C and infinitely large heat capacity, and convection in the water phase is neglected. 2. Conversion In the equation of heat conduction at art (9) a6 K@x2 . a putting Kap (10) at azt PCp a6 =A 9 x2 (11) Let us consider freezing process in unit area, Heat contained in unit volume is Q=OCpt (12) Putting (12) into (11), 2 t ao or) (13) Assume that the medium to be frozen, with thickness L, is divided into n | parts, and the average temperature of each partis ti, tis, tis ccc -tin If each part has volume dV (= dx x1 m7’), heat Qn contained in dv is Qn=f Ope tin dv DAtin (14) if the part is pure water and Qn=A, Cp, tin dv- 65.2Cor 80.0 )x dv, (15) if the part is of pure ice. --80. 0x 1000dV In order to avcid comolication in cealeulation, oT sa. Tadody’ volume incrés Suet freesiv gg ta ce zlectes, rT ; Qn and heat con ae : - Ff fs regarded to pe 7. From Eqs. (l4) a0i 010, 00 ts nl tte: . against Q as shown in Fis. 2 P. Cp, dV If Eq. (13) is converted to a difference equation, Fig. 2 The relationship dQi _ A ai), +t), 16) between Qn and (4t)n Since |t|< 207, so|(At)| max<20 x L0 5. (At)i Taking a favorable scale factor, the computer variable is set to [ 40 It is easily understood that the maximum of | Qi] does not exceed 10‘ Kcal/dV if the width of the partition range is set to 0.1 m. Hence, the equation after scaling becomes as in (17). 1 dQi A - & na = othe {CR iti, (Adi (Atdicd a (17) In order to obtain a sopresentation compatible with the machine units of the com- puter, the graph of Fig. 2 is changed to that of Fig. 3. EAL i/40 The boundary condition is, at x = 0, at — ta27tia (Shoko BE D . ae Slope=p, Gp, dV xqqg7 0 948 Putting (18) into (5), " yj, $422 t 1 La cta-t,) (19) — —0.652 N—0. 800 dx r Qi A0* or htat A, ty “dx = 2 rE fax +h 20) @ 10* Ar _ Slope=jop¢ ayo 64 Accordingly, " A ih bs aan x A, xti Ss Ay hta 4 Ay “ax 4 CAt)2 Fic.3 Th lati hi 40a, “dxth) (A, “ax*t ) 40 ig. e relationship between Qn and ; (At)n in accordance with machine Ati ~5775 , 105 (At )e (21) unit representation (dV = 0.1 m*) 40 . 32 16 40) Operation of Analog Computer Although the accuracy is much improved by dividing the range as finely as possible in handling such difference equations as mentioned in the above, in the present experi- ment, a water depth of 0.8 m is equally divided into 8 steps at 0.1m intervals, ap- proximating the total range with 7 difference equations and an algebraic equation. There are two ways of setting boundary conditions at x = Lm, t = const. at ax = 0 and The latter means that a perfectly insulating wall is placed at x = L, which may be a good approximation if the wall is regarded as constituting a part of the refrigerating chamber. However, the former has more practical significance as an approximation for infinite water depth. Anyway, the more the range is divided, the less error between different ways of approximation. In the present report, attention was focused on how the solution waveform is affected by the way of setting boundary conditions Equations. (At), — 3.775 10.5 (Ate 22 yo 32 7 716 * 40 (22) Do de 40 (hte 5g Ate CA! (23) 10° d4 betcdx 3 $0 - 40 $0 ° 1 dQ _ 40 Cts 5 tds os ets (24) 10* d¢6 10% dx)° 4 Ty bo _1 dQ _ 40 ! Cat is yet 4 Cate (25) 16 dé 10°(dx)?) 49 ~ 40 10 1 dQs _ 40 (hts 5 Cats , Cats | ivf dg. 10%dx) | ra 2 (26) d 40 At 5 (At) (At) vriy ~ cree [at ~ +99 |} (27 1 dQ _ 40 Cats 5 (At , CAtdo | itt} de tis | 40 or rr <r (28) i od 40 At At (At) 7 i = sate et a i (29) -3- at The last equation should be replaced for the boundary condition (Fx) 1 dQs __ 40 CAt)o _CAt)s 10° dg 10%d x)? 40 10 CAt)o 1 t 2 7 a0 =-3 Ag x bie Fy X 0:344x3.9= 0.0336 (41) i = (Qi) (The function form is shown in Fig. 3.) =9 with x=L (29') (30) (31) The results obtained by analyzing these equations are illustrated in Figs. 4 — 7. It should be noted that considerable differences occur depending upon the boundary conditions. As is clear from the fact that the boundary conditions affect the solution less at the vicinity of the surface, it may be supposed that the more the range is divided, the less error due to difference in boundary conditions. +4°C @t =3.9C at x =0.8m @Heat of fusion 65.2Kcal kg 0.7m =) ol +47 200hr 9 gm 250hr 150hr wi 0.5m 0.4m Temperature transition with sentn as cori a @ Stig at x =0.8m ax @) Heat of fusion : 65. 2Kcal /kg — 20°C) Fig. 5 Temperature transition with depth as parameter ~4- @ +t =3.9°C at x =0.8m @ Heat of fusion : 80.0Keal/kg +4t 0.7m 0 50hr 100hr \_ 10h 0.6m 200hr ————— 250hr 0.5m ~° 0.4m 0.3m —10 0.2m 0.1m -15 @Depth om — 200 Fig.6 Temperature transition with depth as parameter at t D gy 0 at x =0.8m +4°C ~ @ Heat of fusion : 80.0Kcal.kg SS btm 0 50hr 100hr\ 9 Gm | 150hr 200hr 250hr 0.5m ~ 0.4m 0.3m ~10> 0.2m 0.1m @ Depth 154 Om —20°C 5 Fig.7 Temperature transition with depth as parameter 1 Heat of fusion : &80Kcal ke. 7. at x =0.8m 250hr —20°C Fig.8 Temperature distribution with time as parameter -5- at pb— ax ° at x =0.8m @e----- t =3.9C at x =0.8m 0.7 F @ Heat of fusion65.2kcal/kg, eT °° c= QO ms Ez a : @ Heat of fusion&0. Okeal kg 0.3} ta 0.2 0.1 ( 6 0 50hr 100hr 150hr 200hr Fig.9 Displacement of boundary face 4. Discussion (1) (3) In the present experiment, the point of 0.8 m depth was always held either ata temperature of 3.9°C or at temperature gradient 0, for computing temperature changes at points of 0~0.7 m depth. These data, however, are not restricted to the 0.8 m boundary, depending upon time scale keeping operation. For instance, temperature transition at 1 m depth obtained by setting 0.8 m depth at 3,9°C can be read as that »t 1 m depth obtained by keeping 8 m depth at 3.9°C, if the time axis is compressed 1/10 times. Since the data »presentec here are obtained by approximation with division-in-8, the solution waveforms taxe a stepwise transition, differing considerably from those expected on — foorystieal phenomena. This is inevitable so long as approximation is mace us Denne, anc in order to obtain smooth curves, approximatio: s aan easing the number of divisions, However, when finer civisi- - tne error may be augmentec cor x5 setting potentiometers. It seems ces oS ivis as lar ge as the computer capacity permits, irrespe tiv. of teratiinal errors, when the partial differential equations are analvzec. The displacement of the boundary be. = and water is plotted as a function of time in Fig.9, in which the curves present some disaccord depending upon the boundary conditions approximating infinite deoth, while they agree fairly well with each other in time shorter than 30 =rs. This is also ascribed to the problem of division number stated in (2): the larger the number of divisions, the less error due to boundary condition setting. In this analysis, displacement of boundary layers within the same section is not taken into consideration. This computation was executed in response to a request from the Engineering Faculty, Kanazawa University. The model employed was HITACHI 505 High Precision Analog Computer. If potentiometers are saved as much as possible, the necessary composition is as follows: Integrating amplifier 7 Summing amplifier 9 Sign changer 16 Potentiometer 20 Dead zone unit 7 Diode 16 <b 0.180 40 —REF 40 10*(dx)?2 __ 80 Q, 2 04 (Able 10M dee 10 40° Lo, 40 10° (dx)? _ Ct), 40 0*idx)?s 80 _ Odds 40 (At) 10* (dx)? A 5 “a0 OU 40 (ans. lOrtaxB 40 Ys 40 80 LAL, — (Ate (Abe 10* (dx)? 40 1 80 Gb, 10*(dx)?2 40° VV 1 40 ADs 10% dxl?p ~~ 40 At); 40 80 (At) 10*(dx)?8 -~\ 0 C- 40 (At), 10*(dx)#2 ~ 40 PCpzto 10° —REF B: time scale factor( =0.4) Fig. 10 Block Diagram a: time scale factor (=0. 4)