Minicomputers and Large Scale Computations - ACS Publications

1 Microcomputer

Plus Saul'yev M e t h o d Solves

Simultaneous Partial Differential Equations of the Diffusion T y p e w i t h H i g h l y N o n l i n e a r

Downloaded by 185.14.195.179 on November 7, 2016 | http://pubs.acs.org Publication Date: June 1, 1977 | doi: 10.1021/bk-1977-0057.ch001

B o u n d a r y Conditions R. KENNETH WOLFE, DAVID C. COLONY and RONALD D. EATON University of Toledo, Toledo, OH 43606 Important today is the ability to answer rapidly and inexpen sively the complex questions posed by an increasingly complex society. Mathematics has played an important role in s c i e n t i f i c problem solving. Practical solutions today rely heavily on com puterized numerical approaches. This paper extends for use with the Hewlett-Packard 67/97 a numerical method due to Saul'yev (1,2). His method i s very similar to the popular method of Schmidt (3) used in graphical, numerical and computer computations to study transient heat con duction problems. This paper w i l l i l l u s t r a t e the use of a small minicomputer (microprocessor) to apply the Saul'yev approach to a simple case and also to a more complex case. The complex case is that of a hot s o l i d slab bounded on one side by a cooler semi - i n f i n i t e s o l i d and exposed at the hot surface to solar radiation, cloud cover and forced or free convective heat losses to a i r . A Simple Case Consider a solid cylinder with faces fixed at two different tem peratures. The sides of the cylinder are insulated. Temperature and time are then related through the extension of Fourier's law to the parabolic partial differential equation: 2

1) 01

8T 9X2

_

3T 3t Τ = T(x,t) = temperature at a point χ and a time t. χ = distance, in feet t = time in hours α = k/pc = thermal d i f f u s i v i t y k = thermal conductivity, BTU/hrft °F/ft p = density, l b / f t C = heat capacity, BTU/lb °F 3

m

m

1 Lykos; Minicomputers and Large Scale Computations ACS Symposium Series; American Chemical Society: Washington, DC, 1977.

2

MINICOMPUTERS

A N D LARGE

SCALE

COMPUTATIONS

Equation 1) i s the d i f f u s i v i t y equation which applies to heat transfer as well as to the transport of matter. Assume that the cylinder has an i n i t i a l constant temperature of Τχ and that at time zero one face i s instantaneously brought to a temperature T . The time-temperature relationship can then be determined analytc a l l y by any of several methods to be: s

2)

T(x,t) = ( T x - T s J e r f t x / v C T + T

$

erf(z) = error function 0

_ , exp(-z ) , 1 Downloaded by 185.14.195.179 on November 7, 2016 | http://pubs.acs.org Publication Date: June 1, 1977 | doi: 10.1021/bk-1977-0057.ch001

2

r

1_

+

z ' 2z3

νπ

Tj3

1-3-5

22 5 " 23 7

K

Z

Z

, ··'

;

'

Equation 2) i s developed in most heat transfer texts. Two good references which treat this problem are Chapman (4) and Carslaw and Jaeger ( 5 j . But direct application of this analytical treatment i s not often f e a s i b l e . Other methods of solution are required. Approximations of the error function are given in Abramowitz and Stegun (6) which are highly accurate and very fast on computers. Schmidt's Numerical Method Consider, as shown in Figure 1, a cylinder that i s divided into hypothetical elements which are frequently called nodes in heat transfer l i t e r a t u r e . To develop Schmidt's numerical method, an energy balance around an element i i s written: [Heat flow from i - 1 ] + [Heat flow from i + l ] 3)

kA(T _ i

1

- T.)

+

kA(T.

ΔΧ

A Δχ at T^ T.

+1

- T.)

=

P

=

H e a t

accumula

tion i n i

CAAx(T' - T. )

ΔΧ

ΔΪ

= area perpendicular to flow, f t = element length, f t = small increment of time, hours = temperature of element i at time t+ t = temperature of element i at time t . 2

Rearranging equation 3) gives: τ

[(τ^ - V

+

(T

i+1

-

Τ ι

)] =

τ;-τ.

Lykos; Minicomputers and Large Scale Computations ACS Symposium Series; American Chemical Society: Washington, DC, 1977.

1.

W O L F E

E T AL.

3

Microcomputer Plus Saul'yev


ΔΧ

τ

t

2

τ

Δχ

'i-l

3

1

'i+l

isothermal 1 ines

depth , D-

Figure 1.

Imaginary division of slab with finite depth and large surface dimensions in relation to the depth


4

4)

MINICOMPUTERS

a N T

i-i

+

(1-2αΝ)Τ. + aNT.

SCALE

COMPUTATIONS

= Τ! ,

+1

a = thermal

A N DL A R G E

diffusivity,

Ν = At/AX

2

For numerical s t a b i l i t y , the coefficients of a l l temperature variables must be non-negative. This means as far as equation 4) is concerned that Ν must be selected such that aN

8

>

ο ο

s

05

1.

W O L F E

E T A L .

7


Saul'yev method. Other methods, such as Schmidt's, require the use of new registers for this purpose. It can be stated that the agreement between the analytical results and the numerical results is excellent. Since the Saul'yev method is an alternating direction method, an even number of rows i s required for accurate r e s u l t s . Each pair of rows i s called a pass; the total number of passes, P, i s equal to t/(2At). Programming the HP 67/97 to obtain the results in Table I , with printout of intermediate r e s u l t s , required 56 program steps out of 224 available. Two memory registers were used to store aN and to maintain a count of the number of rows computed. Downloaded by 185.14.195.179 on November 7, 2016 | http://pubs.acs.org Publication Date: June 1, 1977 | doi: 10.1021/bk-1977-0057.ch001

A More Complex Problem In this section we examine a more complex case and develop some extended formulas. Figure 2 shows a recent situation where the authors (7_) needed to provide a method for f i e l d engineers to predict time-temperature cooling curves for hot asphaltic surfaces placed during highway construction, in order to support decisions on whether or not to allow paving work during marginal weather conditions. Figure 2 also shows the various factors which influence the heat transfer and the temperature-time history of a hot pavement layer. The mathematical model which i s applicable to this s i t u a tion i s summarized below: Governing Equation for Hot Layer 10)

aT 9^ 2

α ι

_ =

3T

3t

Governing Equation for Cold Base

8U

11)

2

ax "

012

=

7

9U

3x

Surface Energy Balance 12) ^ a j l ^ t l

= - aMH + . 6 5 V ( T ( 0 , t ) - T 8

+ ε σ ( Τ ( 0 , ΐ ) + 460)

a i >

)

4

Hot Layer - Cold Layer Interface 13)

TU,t) = UU,t)

contact condition

14)

k BTij^tl

k MhH

x

=

2

e condition at medium A to medium Β interface

e n e r g y

b a l a n C


8

MINICOMPUTERS

A N D LARGE

SCALE

COMPUTATIONS

H = solar flux Cloud Transmission Factor j= M = .15 for clouds = 1.0 for no clouds Fraction of cloud cover = W (visually estimated) \adiation=

e a A

( i T

+ 4 6 0

>

4


Tn=T„.

\onnection-

h A ( T l

-

air

T o )

i n i t i a l uniform

q

T

s o l a r =aMWH

lΊ . T

x=0

2

temperature of media A y at time zero.* " . T

Medium A HOT

k

5

•T 6

Τβ = i n i t i a l tempera ture of base at time zero.* Medium Β

T

7

T

8

T

9

Τ/\.β

Medium to medium change = T in thib case

=

6

Tio Tn

Step size = bAx = 9ΔΧ

COLD BASE Tl2

t Step size = bAx = 9ΔΧ

I 13

* Temperatures T. are specified at a l l node points at time zero. Figure 2.

Hot shb on cold semi-infinite base with surface radiation, convection, and insolation


1.

W O L F E

E T

9


AL.


I n i t i a l Conditions 15)

T(x,0) = T

16)

U(x,0) = U

0

TO

'V,

+

τ

) Mi

Μ

- s f f o )T, • 0 * ^

) η

b Backward 26

>

TO

0, T(0,t-*0) T(x,t) = T^rUx 1+r When the two media have the same thermophysical properties, r=l and T(0,0) = Ti+Ui = 1/2(T +U ) 1

1

1+1 It i s clear from the foregoing equation that the arithmetic average correctly represents the interface temperature only when the thermophysical properties of the hot and cold media are equal.



12

MINICOMPUTERS

A N D LARGE

SCALE

COMPUTATIONS

Surface Equations The surface of a hot l a i d asphaltic con crete pavement layer is acted upon by several environmental influences. These influences include solar radiation (insolation) and cloud cover along with a i r velocity and turbulence. Solar radiation varies in intensity with time of day and the season of the year, while wind conditions are even more variable. Solar radiation can be measured with a pyrheliometer, but in practice, the equations and tables in the ASHRAE Handbook of Fundamentals (8) can be used to accurately determine the solar f l u x . The effect of cloud cover is perhaps more d i f f i c u l t to e v a l uate since height, thickness, water droplet size and percentage of cloud cover a l l influence the transmission of solar energy. For the present purposes, the cloud cover is assumed either to exist or not to e x i s t . If i t e x i s t s , solar radiation is reduced by 85% in a l l computations. The following diagram depicts the surface element and node construction appropriate to the problem under discussion:

(air)

1 2

^One-half element assigned to Τχ

AX-

ΔΧ

The energy balance at the surface element y i e l d s : - radiation loss + solar gain + gain from T

- connective loss

2

= energy gain/loss in Τ χ element 29)

- ε σ Α ί Τ ^ β Ο ) + aMAH + kA(T -Ti) 4

2

Δχ

- ΗΑ(Τχ-Τ . ) a

i

r

= 1/2ΑΔχρΟ(Τ{-Τχ) "Ix

Rearranging terms gives: 30)

2haNAx

(Το-Τχ) + 2 Ν ( Τ - Τ χ ) 2

+ 2aaNMAxH k =T{

2εσαΝΔχ

(Τχ+460) + Τχ 4


1.

W O L F E

13

Microcomputer Plus SauVyev

E TA L .

Forward and backward interpretations of the above are: Forward 31)

2αΝΔχ k

(hTft + aMH - εσ(Τ +460) ) Ι+

1

+ 2αΝ(Τ -Τχ) + Τι = (1 + ^ Ν Δ χ 2

}

J

{

Backward


32)

2αΝΔχ k

(hT + aMH - ε σ ί Τ ι ^ β Ο ) ) 4

0

+ (i . 2oNA_j! j

Τ

ι

+

2 α Ν Τ

,

=

( 1 + 2 α Ν ) τ

.

In equations 31) and 32) the radiation term, (T +460) , has not been applied in a forward and backward sense, since the solu tion for T'. would otherwise be unduly complicated. Tests show this omission to y i e l d negligible errors. The largest discrepancy occurs in T only during the early minutes after time zero. I+

1

x

Solution of Problem in Figure 2 The above equations provide a means for solving the problem depicted in Figure 2. To further i l l u s t r a t e this problem, the following data are hypothesized: Environmental conditions Solar radiation Η = 200 BTU/hr (obtained from ASHRAE Handbook of Fundamental Tables (8). M = 1 or . 1 5 , assume cloud cover with M=.15 Air velocity = 10 MPH A i r Temp. = 80°F h = convection coefficient α = . 6 5 v = .65(10) · = 4.10 Air temperature = 70°F Surface radiation = εσ(Τ+460) = .95 χ 1.731-10~ (T+460) = 1.644-10" (T+460) 8

8

4

9

Hot Solid Absorptivity a for solar flux = .85 Emissivity for solar radiation = .95 I n i t i a l temperature = 300°F k = 1.5 ρ = 150 C = .25 ι α = 0.04 Δχ= 0 . 5 " = ft.

y

4

4

Cold Solid I n i t i a l temperature = 70°F k = 3 C = .25 ρ = 150 α = 0.08 Δχ from node 5 to node 10 = .5 in Δχ from node 10 to node 13 = 9 χ .4 = 3.6 inches

Elapsed time = 15 minutes At =


14

MINICOMPUTERS

A N D LARGE

SCALE

COMPUTATIONS

Total depth of base = 1.6 + 10.8 = 12.4 inches. Table II gives the results of computations using the HP 67/97 at one minute increments up to 15 minutes. These results have been compared with highly accurate results from an IBM 360 and they agree within 2°F at a l l points.


Program The program to obtain the results in Table II took 223 steps out of an available 224 steps. The authors have made several hundred computations using the HP67 or HP97. To these authors, the results are extremely satisfactory. The computations are convenient and the method is generally superior to other approaches. The programs are appended to this paper. Current purchase price of an HP67 computer is $400 and that of an HP97 i s $750. No hardware other than one of the foregoing was needed to perform the complex heat transfer calculations which have been described. Monthly maintenance cost of these instruments can be considered negigible. Conclusion Examples have been presented which demonstrate the usefulness of a small, handheld computer for performing numerical solutions of simultaneous partial d i f f e r e n t i a l equations of the diffusion type. Some mathematical development, or extension, of a standard numerical solution method was required to adapt the method to a small computer. But the results obtained compare very closely to those yielded by an IBM 360 computer; and the use of a small computer makes possible rational decisions on the s i t e in real time by a construction project engineer. It has not been possible, hitherto, to support "go - no go" decisions at a paving s i t e with such detailed analysis of environmental data. Acknowledgements The authors wish to thank Mr. Leon Talbert and Mr. W i l l i s Gibboney of the Department of Transportation, State of Ohio, for their assistance and encouragement. The research reported in this paper was supported in part by the Department of Transportation, State of Ohio, and the U.S. Department of Transportation, Federal Highway Administration.



Minutes

in

Time 275

274

263

254

244

235

227

219

256

247

238

230

222

214

207

200

193

188

183

177

172

3

4

5

6

7

8

9

10

11

12

13

14

15

124 131 135 137 139

158 157 159 159

208

155

207 201 195 190 186 182

219 212 205 199 194 189 184

212

205

199

193

187

182

178

212

226

165

168

171

174

177

180

147

138

138

139

150 149

139

140

140

153 152

154

140

157

188

219

234 184

140

157

192

226 140

158

196

235

253 243

201

245

216

111

156

263

256

228

271

286

283

264

70 71

72 76

93

101 103 104 106 107

113 114 115 116 117

126

127

127

127

127

100

112 126

70 70

72 73

96 98

70

72 95

70

70 70

71

70

70

71 72

91

89

98

71

87 111

71

84 125

96

109

70

70

82

124

93

107

122

70

70 80

91

104

120

70

70

70

70 70

70

70

70

14.4

70

70

70

10.8"

70

77

88

101

117

75

73

70

70

80

7.2"

3.6"

84

97

91

108 113

84

99

76

85

250

288

2min

293

296

275

70

70

70

165

300

300

300

1min

300

300

0

3.2"

2.8"

1.2"

0.8"

0.0"

Interface 1.6" 2.0" 2.4"

Inches from Surface

Temperature Profile of Problem in Figure 2 as a Function of Time

t\x

0.4"

Table II.


^

y?

RCLl RCL5

16-34

23 36 16 26 35 36 36

86 86 4£ 45 81 85 -35 36 15

X.

RCLE

-35

t

T

ISZI STOi RCLl RCL3 T

used

Σαχο^Ν 2^ι (node # + 10) of solid-solid inter face «1^2

+

α

h = 1.35 OR 0.8 h = 0.65 v' NH (as entered) Ν = cloud cover factor = 1.0 or 0.15

36 82 -24

RCL2

τ

-35 -55 -£'4 26 4c 35 45 36 0£

86 -24 1£ 26 46 35 45 36 81 36 03 -24

RCL4 RCLZ

r

λ

S Tût RCLO RCLE

+

EtiTI 1 •

f

2ajAt

54 36 84 36 02 -24 -35 35 15 3£ .3 36 le -55 36 15 -21 il -55 -24

ST0L

35 14

PZS RTH $L6L£

16-51 24 21 06 -62

MH

«NODE. int %)de int

αΝ

Ν I

'air

ft

2

ύ£

2

α ι

2ΡΔΧ

Ρ = # of passes

26 4£

K

RC-L3

'

Hr

35 45 3t *5 3c 53 -35 26 •it 35 45 36" 04

Λ

RCLI RC^4

_t

35 4t

STOI R0L9 ST 01 ROLE

17


E T AL.

" used

S 8

MH *used


R so

f~

MINICOMPUTERS

18

A N DLARGE

SCALE

COMPUTATIONS

Program Listing


STEP

KEY TNTRY 113 114 115 lit 11? 118 119 126 121 122 123 124 125 126 127 12S 129 138 131 132 133 134 135 136 137 138 139 148 141 142 143 144 145 146 147

'ύί

1

_

A a 0 5

COMMENTS

STEP

êc*

S τ*-

5 3 F27 Ri

x

è£ ëz

STOi RTH *LBLÎ 1

35 45 24

STO: RCL* STCi RCLI x;i RCLB ISZl

35 46 36 il

570 ί

RCLI RCL6

KZI

It 2 163

λ=νν

164 165 166 167 168

£T0? tuBLD 1 8 STOI

1 88 189 198 91 192 193 194 195 S6

6i -45 16-33 22 ëi ëi -34

QTûS 1 4 CHS iTOS *LBLS RCLD ISZl STOi ISZl *L8L9 RCLC STOi ISZl RCLI 2 4

8

36 46 16-4* 36 12 16 26 4c 35 45 36 46 36 86

λ=/?

149 158 151 152 153 154 155 156 15? 158 159 168 161

n

TEMPERATURES DISTRIBUTED FROM RIO to R23

ëi

e

h = 0.65V OR IF F2 ON « h = 0.53v°-

71 172 173 174 175 176 17? 178 1 79 188 181 182 183 184 1 85 186 18?

0,8

-62 es ëZ 16 23 CZ -41 -31

197 198 99 2 2 81 2 82 2 83 2 264 2 85 66 2 67 2 88 2 89 2 18 2 2 12 2 13

ee

16-41 22 45 ZI 3c 36 14 16 26 46 35 16 26 21 36 35 16 Zô 36

KEY ENTRY 169 178

-£Z 6 5 EHT'

148

Date input input" by pa G Q used

KEY CODE

45 46 CS 1Ζ 45 46 46 ëZ v4

14 15 2 16 2 1?

RCLi 1 ii FSE RCLi PRTk $LBL8 ISZl tLBLi 5 CHS. XZI GTOi $LBLE

i6~4l 36 4t ëi ëë ~4Z 16 5. 36 45

DSZI R1 *LBL7 CLh RCLi * DSZI RCLI

£70? R; RCL6 1 1

R7h

Display node number

Zi ëë it 26 46 21 Ci èl ë>5 16-4* 22 45 21 15 ëi

1

STOI

RCL ί RCL6 STOI Ri RCLi

COMMENTS

KEY CODE

35 36 36 Z5

46 45 66 46 -ci 36 45 -55

COMPUTE AVERAGE TEMPERATURE OF LAYER 1 USING TRAPEZOIDAL RULE

-24 It 25 46 16-31 21 ë7 -51 36 45 16 25 46 36 4t 81 ëi -45 16-42 22 C7 -3i 36 06 ëi ëi -45 -24 24

R-ί

ib-

RTri

B

22 U; •a 14 ë* ëë 46

Compute

220

PRINT TERMPERATURES FROM RIO to R23 LAE ELS

c

Allocate D Print/ nut

0

e

1 2

b

c

used h=c v0.8

2

3

4

7

8

9

1

6

d

—uspri

SET STATUS

FLAGS Ε

useo*

3

FLAGS ON OFF O D D 1 • • 2 • • 3 • •

TRIG DEG GRAD RAD

DISP • • •


FIX SCI ENG

• •

α

WOLFE

E T AL.


APPENDIX Β

User Instructions SAUL'YEV PROGRAM fr-


STEP

1 2 3

Start

•

•

INSTRUCTIONS

Put in data with Saul'yev input/output program Read Saul'yev program Output from Saul'yev input/output program

INPUT DATA/UNITS

OUTPUT DATA/UNITS

KEYS

, :

11 11

!

1 1

"J [ """""]

L "Ί 1 -

1

LID (Z~1IZZ1 i : ! 1 II 1 1 1 1 C~3 '• 1 • CZI i 1 i 1 1 1 • CZI ι :i ι ι :ι ι ι ir I : π ι

π

α

IZUtZZ]

• ι •

11 11

;

11 11

["

]ι

; ι 1 ι

:

i!

;

:

11

• Γ

: u d 11

cri 1 i 1

ι "J ( II 1 ί

I I

• 1

1 ι 1 1 1


MINICOMPUTERS

20

STEP

KEY ENTRY 891


062 005 004 005 006 007 008 005 010 011 012 013 014 015 01 6 917 818 015 020 021 022 023 024 025 026 027 028 025 030 031 032 033 034 035 036 Θ37 03S 035 040 041 042 043 044 045 046 047 048 045 050 051 052 053 054 055 056

*L£LH

ÎSBE 1

STxi DSZJ. RCi. ι ISZI ISZI RCL i * RCL: X

USZ:

ST+i 1 RCL: +

ISZI RCL6 RCLI (iBS

X=V? &70C 2 i x=y:

C7CD 2 3 RCLI

X=r'7 C7Ci 1 1 CHS

Α=ϊ'? CT02 CTOtt *L6L1

SFÊ 2 CHS S7CI CTCa *LfiL2 CSEE CF6 1

ST-e

RCLC X=07

23 ;2 »o ti 36 t'. 35--35 16 25 56 16 26 it 26 56

START Cale. Τ ι MAIN PROGRAM

45 4o 4C4c 4c

45 -55 36 Ci -35 4c it 35--55 45 Ci 36 Ci -55 35-24 45 it 26 46 36 66 36 4 6 16 3i 16--33 22 13 32 Ci 16--33 22 14 '02 ti«i 36 46 i6--33 22 Ci Si 01 -22 16 -33 22 02 22 16 11 21 Ci 16 21 Ct? 02 02 -22 35 46 22 16 i i Zi 62 23 12 16 22 30 Ci 35 -45 CC 36 ce 16 -43

Next Τ T. ? interface* Step size change?

Tl ? 3

Surface T ? x

Prepare for backwarq pass, set flag 0.

Set RI to negative values. END OF BACKWARD PASd Reduce # of passes left by one.

J12

2

Ti»T,

BP* 5

057 058 055 060 061 062 063 064 065 066 067 068 065 070 071 Θ72 073 074 075 076 077 078 075 080 081 082 083 084 085 086 087 088 085 050 051 052 053 054 055 056 057 058 055 100 101 102 103 104 105 106 107 106 105 110 111

RTN 67 OH «LÊL3 CHS STOi

isz;

16 22

CTCa *LEL8

i 0

STCI RCLi RCL?

it

ISZI RCLi

4 6 0 + 4 yx

j 6 4

4 EEA

5 CHS Λ

-

RCLS 8

5

Λ

RC15 A

RCLi •

RCLi 2 Λ

ISZI it RCL i DSZI it F0? 16 CTOfc 22 RCi.;

-

21

λ

Node

24 ii 63 46 0C -

Minicomputers and Large Scale Computations - ACS Publications

Recommend Documents