HP Forums
Transients in long cylinder with step temperature change - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP-41C Software Library (/forum-11.html)
+--- Thread: Transients in long cylinder with step temperature change (/thread-7287.html)

Transients in long cylinder with step temperature change - Ángel Martin - 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!)

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