Transients in long cylinder with step temperature change
11-22-2016, 03:29 PM (This post was last modified: 11-23-2016 11:05 AM by Ángel Martin.)
Post: #1 Ángel Martin Senior Member Posts: 1,068 Joined: Dec 2013
Transients in long cylinder with step temperature change
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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)