11-22-2016, 03:29 PM
Transients in long cylinder with step temperature change [ TRT ]
Extension of the TXT program from the ETSII Collection
This program calculates the temperature T(r) in points r of an infinite cylinder of radius R, after experiencing a thermal shock – or sudden change of ambient temperature, from To initial to Tf final.
Similar to the flat plate case, the Biot number is calculated from its constituent factors. The same data entry process is used like in the infinite plate, only now it is cylindrical symmetry instead.
The resulting temperature is expressed as an infinite sum as follows:
T(x,t) = Tf + (T0-Tf) SUM { (2/Mn). f(n, r). exp[-a.t.(Mn/R)^2 ]} ; n = 1,2,...
With: f(n,r) = J1(Mn). J0(Mn.r/R) / [ J1^2 (Mn) + J0^2(Mn)]
And (Mn) are the n roots of the equation defined by: (Mn) J1(Mn) = Bi J0(Mn)
Which, leaving the Biot number alone in the second term, can be expressed as the intersection of the Biot number with the function x.J1(x)/J0(x), shown in the graphic (see WolframAlpha), where the asymptotic boundaries will provide a reasonable criteria for the estimations needed by the root-finding routine (more about this later).
Example.
A very long metal rod of radius R=0.14 has a uniform temperature of 1,000 deg C. It is suddenly immersed into a cooling fluid stream at 50 deg C. Calculate the temperature in its center and outer boundary 15, 30 and 60 minutes after the sudden step temperature change. The physical properties of the material are given below:
alpha = 1.66 E-6 m^2/s
h = 20,000 kcal/H.m^2.C = 23,260.0 W/m^2 K
K = 100 kcal/h.m.C = 116.30 W/m.K
The result temperatures are shown in the table below: (warning: very slow convergence!)
A few remarks regarding the implementation.
By direct inspection of the plot in previous URL, it’s clear that this case is much more demanding on the root-finder algorithm than the previous one. As the Biot number value increases, the intersection with the graphic will occur in zones with a very steep slope, making the identification of the root very tricky – so much so that the FOCAL routine “SLV” is not adequate and misses the roots, even if very fine-tuned search intervals are provided – which is also a difficult affair.
To search for each of the Mn roots, the program uses symmetric intervals centered at the initial value, and distance "one", i.e:
[ n*(Mn)init - 0.5 ; n*(Mn)init + 0.5]
With: (Mn)init = sqr{ (3/2) [ sqr (1+4Bi/3) – 1]
In this version we’ve used FROOT instead, also included in the SandMath - which was already required for the Bessel functions, so no more dependencies are added. The estimation for the initial guesses becomes very important for the successful root identification, and the execution times – which are going to be very long regardless; better crank up your turbo emulator for this one!
Another important remark is that repeating the calculations for different values of (t, r) (analysis time and distance to the cylinder axis) has been expedited dramatically for subsequent runs with longer times than previous executions). In that case there's no need to re-calculate or find additional Mn roots beyond those already identified, as the contribution of the terms to the infinite sum will be smaller due to the larger argument in the inverse exponential function:
f(n, r) . exp[- alpha.t.(Mn/R)^2 ]}
This of course is not so straight-forward as one may think, because the series is alternating the sign of its terms so the contributions are not always in the same direction. The program stores the successive roots found in an X-memory file, to be reused when the analysis is repeated with longer values of cooling time. In this way, if the X-mem value is zero then the corresponding root needs to be sought for.
The program listing is given below. Note that the ALPHA registers are used by the infinite sum routine to calculate the partials and to store the current term. Because the MCODE function JBS also uses the ALPHA registers for scratch, we’ll use the function A<>RG to preserve ALPHA in {R17-R20} while the general term is being calculated.
XROM “?” is a simple data-entry utility functions to save bytes. Be careful if you use arithmetic functions with the value in X – that would alter the expected stack configuration and may be disruptive to the program.
And here's the complete listing - feedback is always welcome.
Extension of the TXT program from the ETSII Collection
This program calculates the temperature T(r) in points r of an infinite cylinder of radius R, after experiencing a thermal shock – or sudden change of ambient temperature, from To initial to Tf final.
Similar to the flat plate case, the Biot number is calculated from its constituent factors. The same data entry process is used like in the infinite plate, only now it is cylindrical symmetry instead.
The resulting temperature is expressed as an infinite sum as follows:
T(x,t) = Tf + (T0-Tf) SUM { (2/Mn). f(n, r). exp[-a.t.(Mn/R)^2 ]} ; n = 1,2,...
With: f(n,r) = J1(Mn). J0(Mn.r/R) / [ J1^2 (Mn) + J0^2(Mn)]
And (Mn) are the n roots of the equation defined by: (Mn) J1(Mn) = Bi J0(Mn)
Which, leaving the Biot number alone in the second term, can be expressed as the intersection of the Biot number with the function x.J1(x)/J0(x), shown in the graphic (see WolframAlpha), where the asymptotic boundaries will provide a reasonable criteria for the estimations needed by the root-finding routine (more about this later).
Example.
A very long metal rod of radius R=0.14 has a uniform temperature of 1,000 deg C. It is suddenly immersed into a cooling fluid stream at 50 deg C. Calculate the temperature in its center and outer boundary 15, 30 and 60 minutes after the sudden step temperature change. The physical properties of the material are given below:
alpha = 1.66 E-6 m^2/s
h = 20,000 kcal/H.m^2.C = 23,260.0 W/m^2 K
K = 100 kcal/h.m.C = 116.30 W/m.K
The result temperatures are shown in the table below: (warning: very slow convergence!)
Code:
Point t = 15 min t = 30 min t = 1 hour
--------------------------------------------------------------------
center (x=0 cm) 945.7185485 704.2922460 343.4690201
Outer edge (x=0.14 m) 102.5288706 80.51769740 63.05841690
A few remarks regarding the implementation.
By direct inspection of the plot in previous URL, it’s clear that this case is much more demanding on the root-finder algorithm than the previous one. As the Biot number value increases, the intersection with the graphic will occur in zones with a very steep slope, making the identification of the root very tricky – so much so that the FOCAL routine “SLV” is not adequate and misses the roots, even if very fine-tuned search intervals are provided – which is also a difficult affair.
To search for each of the Mn roots, the program uses symmetric intervals centered at the initial value, and distance "one", i.e:
[ n*(Mn)init - 0.5 ; n*(Mn)init + 0.5]
With: (Mn)init = sqr{ (3/2) [ sqr (1+4Bi/3) – 1]
In this version we’ve used FROOT instead, also included in the SandMath - which was already required for the Bessel functions, so no more dependencies are added. The estimation for the initial guesses becomes very important for the successful root identification, and the execution times – which are going to be very long regardless; better crank up your turbo emulator for this one!
Another important remark is that repeating the calculations for different values of (t, r) (analysis time and distance to the cylinder axis) has been expedited dramatically for subsequent runs with longer times than previous executions). In that case there's no need to re-calculate or find additional Mn roots beyond those already identified, as the contribution of the terms to the infinite sum will be smaller due to the larger argument in the inverse exponential function:
f(n, r) . exp[- alpha.t.(Mn/R)^2 ]}
This of course is not so straight-forward as one may think, because the series is alternating the sign of its terms so the contributions are not always in the same direction. The program stores the successive roots found in an X-memory file, to be reused when the analysis is repeated with longer values of cooling time. In this way, if the X-mem value is zero then the corresponding root needs to be sought for.
The program listing is given below. Note that the ALPHA registers are used by the infinite sum routine to calculate the partials and to store the current term. Because the MCODE function JBS also uses the ALPHA registers for scratch, we’ll use the function A<>RG to preserve ALPHA in {R17-R20} while the general term is being calculated.
XROM “?” is a simple data-entry utility functions to save bytes. Be careful if you use arithmetic functions with the value in X – that would alter the expected stack configuration and may be disruptive to the program.
Code:
1 LBL "?"
2 RCL IND X
3 "|="
4 ARCL X
5 "|-?"
6 CF 22
7 PROMPT
8 FS?C 22
9 STO IND X
10 END
And here's the complete listing - feedback is always welcome.
Code:
1 LBL "TRT"
2 SIZE?
3 21
4 X>Y?
5 PSIZE
6 E
7 CF 04
8 "SAVE-RT? YN"
9 PMTK in OS/X Module
10 X#Y?
11 GTO 04
12 SF 04
13 "mN" file Name
14 SF 25
15 PURFL purge if exists
16 9 up to nine roots
17 CF 25
18 CRFLD create new file
19 LBL 04
20 "TINI"
21 11
22 XROM "?"
23 "TEND"
24 10
25 XROM "?"
26 ST- 11
27 "H"
28 E
29 XROM "?"
30 STO 14
31 "K"
32 2
33 XROM "?"
34 ST/ 14
35 "RO"
36 3
37 XROM "?"
38 ST* 14 Bi
39 "ALPHA"
40 0
41 XROM "?"
42 LBL C
43 0
44 FS? 04
45 SEEKPT
46 "R"
47 4
48 XROM "?"
49 "TIME"
50 12
51 XROM "?"
52 RCL 14
53 0,75
54 /
55 E
56 +
57 SQRT
58 E
59 -
60 1.5
61 *
62 SQRT
63 STO 13
64 "aJ"
65 ASTO 06
66 XROM "sI" infinite SUM
67 ST+ X
68 RCL 11
69 *
70 RCL 10
71 +
72 "T(X,T)="
73 ARCL X
74 PROMPT
75 GTO C
76 LBL "aJ"
77 A<>RG preserve alpha
78 17 in {R17 - R20}
79 FC? 04 roots in X-Mem?
80 GTO 04 no, need to search
81 GETX yes, get current
82 X#0? valid root?
83 GTO 05 yes!
84 RCLPT no, backtrack pointer
85 E
86 -
87 SEEKPT
88 LBL 04 need to search
89 RCL 18 n
90 RCL 13 delta
91 *
92 RCL X
93 ,5
94 ST- Z
95 +
96 "JN" function to solve
97 FROOT
98 FS? 04 save option?
99 SAVEX yes, oblige
100 LBL 05
101 STO 07 Mn
102 VIEW X
103 E
104 X<>Y
105 JBS J1
106 STO 08
107 0
108 RCL 07 ln
109 JBS Jo
110 X^2 Jo^2
111 RCL 08 J1
112 X^2 J1^2
113 + Jo^2+J1^2
114 STO 09
115 0
116 RCL 07 Mn
117 RCL 03 R0
118 / Mn/R0
119 RCL 04 r
120 *
121 JBS Jo(ln.r.R0)
122 RCL 08 J1
123 *
124 RCL 09 Jo^2+J1^2
125 /
126 RCL 07 Mn
127 /
128 LASTX Mn
129 RCL 03 R0
130 / Mn /R0
131 X^2 (Mn/Ro)^2
132 RCL 12 t
133 *
134 RCL 00 a
135 * a.t(Mn/Ro)^2
136 CHS
137 E^X
138 *
139 A<>RG restore ALPHA
140 17 from {R17-R10}
141 RTN
142 LBL "JN" function to solve
143 E
144 X<>Y
145 JBS
146 X<>Y
147 ST+ X
148 ,
149 X<>Y
150 JBS
151 ST/ Z
152 RDN
153 ST+X
154 *
155 RCL 14 Bi
156 -
157 END