• Author / Uploaded
• Diego Alejandro

• Categories
• Integral
• Ecuaciones
• Análisis matemático
• Física y matemáticas
• Matemáticas

Citation preview

CÁLCULO NUMÉRICO DE LA EFICIENCIA DE ALETAS CON SECCIÓN TRANSVERSAL NO UNIFORME Introducción Las ecuaciones diferenciales resultantes de los balances diferenciales de calor para aletas con sección transversal no uniforme son del tipo de Bessel y cuyas soluciones resultan ser funciones modificadas de Bessel. La solución simplificada sugerida por Harper y Brown (2) es la que usualmente se presenta en forma gráfica en la mayoría de los textos de transferencia de calor y corresponde a la extensión, por medio de una longitud corregida, de la solución obtenida para una aleta con un extremo aislado (4). Esta solución aproximada puede llevar a errores y para situaciones de gran precisión sería bueno disponer de una forma más elaborada de cálculo. El presente cálculo numérico obvia

estas aproximaciones y

adicionalmente otras, que usualmente se hacen para

efectos de poder obtener soluciones analíticas al problema, tales como: coeficientes de película por convección y radiación y conductividad térmica constantes. Debido a que la propuesta numérica de cálculo trabaja directamente con la ecuación diferencial, los anteriores parámetros pueden ser corregidos a lo largo del proceso de integración. El problema de valores de frontera se resuelve mediante una función objetivo, que es planteada en una condición de frontera y sobre la cual actúa um método mejorado de memoria de interpolación IMM, el cual suministra valores cada vez mas adecuados para que, en un proceso iterativo de integración numérica, se satisfaga la condición de frontera exigida. Adicionalmente y para efectos de comparación, se presenta y desarrolla una alternativa numérica que calcula directamente las funciones de Bessel (5), mediante técnicas de integración bastante precisas como las cuadraturas de Gauss y de Laguerre. Esta estrategia numérica puede considerase como una alternativa intermedia en la solución del problema, pero mucho más exacta que la presentada usualmente en los textos. La alternativa general propuesta será presentada con cierto detalle y es importante anotar que para situaciones de coeficientes de convección, radiación, y conductividad térmica constantes la solución convencional debe ser igual a la que origina la propuesta numérica implementada en este trabajo.

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

1

Eficiencia de las aletas Para una condición de flujo unidimensional estable la ecuación diferencial que se origina en un balance de calor es d  dT  A(x) - (hc /k) T-T  P(x) = 0  dx  dx 

A(x), área de la sección transversal de la aleta en x P(x), perímetro de la aleta en x

1

hc , coeficiente de película k, conductividad térmica T, temperatura en la aleta T , temperatura del fluido x, dirección del cambio en la temperatura

Existe una gran variedad de formas geométricas (4), (7), (13), (14). En la Figura 1 se presentan dos tipos de aletas bastante utilizadas, con sección transversal no uniforme: aleta triangular y aleta anular de espesor uniforme.

hc,T hc,T

t

t r L ra,Ta x To

rb,Tb Figura1. Aletas triangular y anular de espesor uniforme

El caso de la aleta triangular recta, (8), genera la siguiente ecuación diferencial 2 d T + 1 dT - 1 ( 2h L /kt )( T- T ) = 0  c 2 x dx x dx

2

Donde, Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

2

L, longitud de la aleta t, espesor de la aleta en la base T, temperatura del fluido que rodea la aleta hc, coeficiente de película del fluido que rodea la aleta La ecuación inmediatamente anterior es una ecuación diferencial de Bessel, (15), cuya solución se describe a continuación. La forma general de la ecuación diferencial de Bessel es

d2 Y (1-2a) dY (a2 -p2c2 ) + +  ( bcXc-1 )2  Y + =0 2 X dX dX X2

3

y su solución esta dada por

Y= X aZp( bX c) Zp = AJn ( AJn (

) + BI-n ( ) p=n=fraccionario ) + BY+n ( ) p=n=0,1,2,3,4,5,.

AIn ( AIn (

) + BI-n ( ) p=n=fracc.(im.) ) + BK+n ( ) p=n=0,1,2,.(im.)

4

Si se comparan término a término, exponentes y coeficientes, para el caso se obtiene (1-2a)= 1

a= 0

2(c-1)= -1  2h L  b2c 2 = -  c   kt 

c= 1/2

1/2

 -8hcL  b=   = mi  kt  a2 -p2c2= 0

5

b= m = (8hcL/kt)½

p=0

En esta situación las funciones de Bessel son imaginarias y el orden n = 0; por lo tanto, el perfil de temperatura es

1

1

6

(T-T ) = C1I0 [mx 2 ] + C 2K o[mx 2]

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

3

donde, m = (8hcL/kt)½ I0, función de Bessel modificada de primera clase y orden cero K0, función de Bessel modificada de segunda clase y orden cero C1 y C2, constantes de integración; valores A y B de la ecuación 4. De la Figura 1 las condiciones de frontera son: Para x = 0 T = T0 (temperatura en el extremo de la aleta)

(T0 - T )= C1I0[0] + C2K 0[0]

7

De (6) la función K0(0) de la ecuación 7 es igual a  y como T tiene un valor finito, entonces C2= 0. Para x = L T = Ts (temperatura en la base de la aleta)

(T s - T  ) = C1I0 [mL1/2 ]

8

De donde se obtiene que el valor de C1 es

C1 =

(Ts -T ) 1/2 I0[mL ]

9

Al reemplazar la ecuación 9 en la ecuación 6, resulta la siguiente expresión para el perfil de temperatura

1

(T-T ) [ 2] = IO mx 1 (T s-T ) IO [mL2 ]

10

El calor real transferido por la aleta al fluido puede calcularse, valorando la derivada en x=L, así

(

dT m = ) dx x=L 2L1/2

1

I1[mL 2] (T  T ) s  1 I0 [mL 2]

11

Lo que conduce finalmente a la siguiente expresión para el cálculo del calor real entregado por la aleta para una profundidad z Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

4

1

1 I [mL 2] Qreal = z (2hck t)2 (T s  T ) 1 1 I0[mL 2]

12

Donde z es la profundidad de la aleta y Qr el calor real. La eficiencia η de la aleta se define como:

η=

[ calor real transferido por la aleta ] [calor teórico transferdido por la aleta ]

13

El calor teórico es el transferido como si toda la aleta estuviese a la temperatura de la base. Por lo tanto, se tiene que

Q t = 2hc .z.L  ( T s-T )

14

donde, Qt = calor teórico L'= [L2 + (t/2)2]½ Al reemplazar las ecuaciones 12 y 14 en la ecuación 13 se obtiene

1

η=

1

( 2hc k t )2 I1[mL 2] 1 2hcL  I0 [mL2 ]

15

La ecuación diferencial para la aleta anular con espesor uniforme es 2 d T + 1dT - ( 2h /kt)(T-T ) = 0  c 2 r dr dr

16

donde, t, espesor de la aleta r, posición radial

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

5

Para el caso de la ecuación diferencial 16 y al comparar término a término en la ecuación 3, resulta

17

(T-T ) = C1I0[mr] + C2K 0[mr]

donde m = (2hc/kt)½ Los valores de C1 y C2 se determinan a partir de las condiciones de frontera, así: En r = ra T = Ta

18

(T a-T ) = C1I0[mr a] + C2K 0[mr a]

Para r = rb se tiene que

-k (dT/dr) =hc (Tb -T )

19 -k (mC1I1 (mrb ) - mC2K1[mr b] ) = hc ( C1I0[mr b]+ C2K 0[mr b] )

Al solucionar el sistema de ecuaciones 18 y 19 se obtiene

C1 =

(T a-T ) ( I0[mr a] + a2K 0[mr a] )

C2 = a2C1

20 a1 = (hc /k m) a1I0[mr b]+I1[mr b] a2 = K1[mr b]-a1K 0[mr b]

El calor real transferido por la aleta al fluido puede calcularse conductivamente, valorando dT/dr en r= ra, y se obtiene

21

Qr = -2π k.t. ra m(C1I1[mra ]-C2K1[mra ]) (T a-T )

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

6

El calor teórico transferido por la aleta es 2 2 Q t = [ 2πhc (rb -ra ) + 2 hc t rb ] (T a-T )

22

Finalmente, la expresión para la eficiencia de la aleta es

η=

-2π k t r a m (C1I1[mr a]-C2K1[mr a]) [ 2πhc (rb2 -ra2 )+2πhc trb ] (TA -T )

23

Cálculo numérico de las funciones modificadas de Bessel

La función de Bessel modificada de primera clase y orden v entero es

Iv [x] =

1





e

x cos( s )

24

cos(vs) ds

0

El cálculo de Iv[x] se efectúa por medio de una cuadratura Gaussiana de diez y seis puntos (1), (6). Los valores de las abcisas y de las funciones de peso son presentados en la rutina de Datos del programa anexo.

La función de Bessel modificada de segunda clase y orden v entero es 

K v [x] =

1 exp (vs-xcosh (vs))(1+ e-2vs ) ds 2 0

25

La valoración numérica de Kv[x] se realiza mediante una cuadratura de Gauss-Laguerre. Para una mayor precisión se dividió el intervalo de integración en dos. Inicialmente utiliza una cuadratura Gaussiana hasta un valor relativamente alto y luego se cambia a una cuadratura de Laguerre, (1), (5). Cálculo numérico de la eficiencia

A manera de ilustración se presentará la técnica numérica, en una forma más o menos detallada, para el caso de la aleta anular. Lo primero es convertir la ecuación diferencial de segundo orden en dos ecuaciones diferenciales de primer orden, así Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

7

y1= T dy1/dr = dT/dr = y2 d2T/dr2 = dy2/dr Al reemplazar en la ecuación diferencial 16 se obtiene dy2/dr = (-1/r)y2 + (2hc/kt)(y1-T) Las dos ecuaciones diferenciales resultantes son: dy1/dr = y2 dy2/dr = (-1/r)y2 + (2hc/kt)(y1-T) La nomenclatura corresponde a la utilizada en el programa anexo de computación. Mediante un análisis numérico sencillo se observa que para este caso es conveniente empezar el proceso de integración en r = rb, lo cual implica suponer la temperatura (si se comenzara en r = ra, sería necesario suponer el valor de la derivada, lo cual sería un poco más difícil). Adicionalmente en r = rb se debe cumplir que –k dt/dr = hc (Tb-T), lo cual permite calcular el valor de la derivada de la temperatura y así puede iniciarse el proceso de integración, que en términos de las variables definidas anteriormente queda como y2 = -(hc/k)(y1- T) La función que debe satisfacerse se ubica en el otro extremo, r = ra, y es F = Ta - y1calc. Para un valor supuesto de la temperatura en r = rb y mediante una técnica de integración (Runge Kutta), se consige un valor de la temperatura en r = ra, lo que suministra un valor de la función F. Un método de interpolación IMM actúa, a partir de dos valores iniciales, para originar valores mejorados de la temperatura en r = rb y así rápidamente llegar a que la función F sea cero. Debido a la alta velocidad de convergencia del IMM, el número de iteraciones es bastante reducido, con respecto a cálculos convencionales. Para un mayor entendimiento del IMM, como método de interpolación, por favor consultar (3), (10), (12).

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

8

Para la solución numérica de la eficiencia en el caso de la aleta triangular, por conveniencia se modificó la ubicación en el sistema de coordenadas, para evitar problemas numéricos de indeterminación, con respecto al empleado en la solución analítica. El sistema de coordenadas parte en x=0 y corresponde a la base de la aleta, lo que conduce a la siguiente ecuación diferencial

2 1 dT 1 dT (2hc / kt)(T-T ) = 0 2 (L-x) dx (L-x) dx

26

Adicionalmente, se cambió la estrategia numérica, pues el proceso de integración comienza en x=0 y la variable supuesta fué el valor de la eficiencia de la aleta, lo que implícitamente condiciona el valor de la primera derivada y así puede iniciarse el proceso de integración. La función a ser cero se planteó en x=L y corresponde a la condición de frontera conductiva-convectiva. Se efectúa, en forma similar, el cambio de variables anotado anteriormente y mediante la técnica de integración y el método de interpolación IMM se busca la solución del problema.

Resultados y conclusiones

En primer lugar se comparan los resultados de la propuesta numérica con respecto a cálculos convencionales presentandos en los textos y con respecto a la alternativa intermedia que calcula directamente las funciones de Bessel, mediante cuadraturas de Gauss y Gauss-Laguerre. Para esta situación es necesario suponer que el coeficiente de película y la conductividad térmica son constantes. Lo anterior permite validar la propuesta numérica, pues para la situación de coeficientes constantes en la ecuación diferencial los resultados de las dos alternativas deben coincidir. En (8) se presenta el cálculo de la eficiencia de una aleta triangular, con los siguientes datos: hc = 15 Btu/h.ft2.oF L = 4 in k = 15 Btu/h.ft.oF t = 1 in T= 100 oF Ts = 1100 oF

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

9

Los valores obtenidos mediante las dos alternativas numéricas fueron idénticos. El calor transferido por unidad de profundidad es de 5069.60 Btu/h-ft y la eficiencia de la aleta η = 0.5030. Del libro de Pitts y Sissom (11), para una aleta anular con ra = 1/2 in rb = 1 in t = 0.009 in hc = 1.5 Btu/h.ft2.oF k = 93 Btu/h.ft.oF T= 80 oF Ta= 330 oF El libro suministra un valor de la eficiencia de η  0.94 y textualmente anota la dificultad para la obtención de valores exactos de las funciones de Bessel y además la condición de frontera en r = rb se considera como un extremo aislado. Mediante las dos alternativas numéricas el valor de η es 0.9653. Es importante anotar que el número de iteraciones es reducido (tres ó cuatro), para un paso de integración de 0.01. Para el caso inmediatamente anterior de la aleta anular y mediante la propuesta numérica que permite tener en cuenta el cambio en el valor de los coeficientes de película por convección y radiación se obtuvo un valor para η = 0.90. Para el cálculo de hcr se utiliza una expresión empírica para la convección libre en placas horizontales (4), de la forma hcr = 0.29(T-T)0.25 y la ecuación de Stefann Boltzman para valorar el coeficiente por radiación. Para este caso el valor de la emisividad se tomó como 0.9.

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

10

Figura 2 Eficiencia de aleta triangular 1 0,9

Eficiencia,

0,8 0,7 0,6 0,5 0,4 0,3 0,2 0,1 0 0

1

2

3 L(2h/kt)

4

5

0.5

Figura 3 Eficiencia de aleta radial 1 0,9

Eficiencia,

0,8 0,7

(ra/rb)=1

0,6

(ra/rb)=2

0,5

(ra/rb)=3

0,4

(ra/rb)=4

0,3 0,2 0,1 0 0

1

2

3

(ra-rb)[2hc/kt]

4

5

0.5

La propuesta numérica de cálculo resulta ser atractiva, pues levanta las restricciones que usualmente se hacen para la valoración de la eficiencia de aletas y las condiciones de frontera del modelo matemático son las reales del fenómeno físico. Adicionalmente, corrige en el proceso de integración los coeficientes de la respectiva ecuación diferencial, garantizando así una mayor precisión en el cálculo del perfil de temperatura y por ende Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

11

de la eficiencia de la aleta. La técnica numérica propuesta es ágil, exacta y muy versátil. El programa elaborado puede utilizarse como un módulo numérico muy práctico y prueba de ello son las Figura 2 y 3, originadas mediante un cambio iterativo de un parámetro geométrico o de una propiedad térmica. Las anteriores gráficas permiten el cálculo de la eficiencia de aletas triangular recta y anular de espesor uniforme. El listado del programa tiene la propuesta numérica para el caso de coeficientes constantes. Adicionalmente se indican las leves modificaciones que son necesarias en el caso de coeficientes variables, para una aleta anular con espesor uniforme.

Bibliografía

1. Carnahan, B., Luther, H.A. Cálculo numérico. Métodos, aplicaciones. Editorial Rueda, Madrid, 1979. 2. Harper, W.P. Brown, D.R. Mathematical equations for heat conduction in the fins of air-cooled engines, NACA Rep. 158, 1922. 3. Hildebrand, F.B., Introduction to Numerical Analysis. Second Edition. McGraw-Hill. 1974. 4. Holman, J.P. Transferencia de calor. McGraw-Hill book Co. Compañía Editorial Continental S.A. Segunda edición, 1980. 5. Jahnke, E., Emde, Losch, F. Tables of Higher Functions, 6th ed., McGraw-Hill, 1960. 6. Jenson, V.G., Jeffreys, G.V. Métodos matemáticos aplicados en Ingeniería Química. Ed. Alhambra. Madrid. Primera edición,1969. 7. Karlekar, B.V., Desmond, R.M. Transferencia de calor. Nueva Editorial Interamericana S.A., 1985. 8. Kreith, F. Principles of heat transfer. Harper international Edition. Fourth edition, 1986. 9. Ozisik, N.M. Tranferencia de calor. McGraw-Hill Latinoa. S.A. 1985.

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

12

10.Pang, T., Knopf, F.C., Numerical analysis of single variable problems with the use of continued fractions. Computers Chem. Eng. 10, 87, 1986. 11.Pitts, D.R. Sissom, L.E. Heat transfer. McGraw-Hill Latinoam. S.A., 1977. 12.Shacham, M. An improved memory method for the solution of a nonlinear Chem. Eng. Sc. 44, p.1495, 1989. 13.Welty, J.R. Tranferencia de calor aplicada a la Ingeniería. John Wiley & Sons. Editorial Limusa, 1981. 14.Welty, J.R., Wilson, R.E., Wicks, C.E. Fundamentals of momentum, heat and mass transfer. Wiley International Editorial, 1969. 15.Wylie, R.C. Advanced engineering mathematics. McGraw-Hill Book Co., 1975.

SCREEN 12 REM"-------------------------------------------------------------------PRINT " Este programa calcula analítica y numéricamente la eficiencia de " PRINT " aletas con seccion transversal no uniforme, tales como una aleta " PRINT " anular de espesor uniforme y triangular recta. Para el cálculo de " PRINT " las funciones modificadas de Bessel en la solución analítica se " PRINT " utilizaron cuadraturas de Gauss y Gauss-Laguerre. Para el cálculo " PRINT " numérico se crearon funciones en las condiciones de frontera no " PRINT " conocidas. Al utilizar un método de interpolación -IMM- y con la " PRINT " ayuda simultánea de una técnica de integración de Runge Kutta de " PRINT " cuarto orden se busca satisfacer las funciones respectivas." REM"-------------------------------------------------------------------DEFDBL A-Z DEF fnf (x) = EXP(x9 * COS(x)) * COS(n0 * x) DEF fng (x) = (1 + EXP(-s2 * x)) * EXP(n0 * x - x8 * (EXP(x) + EXP(-x))) REM Para el calculo de In y Kn de la funciones de Bessel REM se utiliza una integracion Gaussiana y de Laguerre DIM u(16), w(16), v(16), z(16), ii(5), kk(5), iii(5), kkk(5) DIM y(5), dery(5), phi(5), savy(5), x(25), t(25) DIM aa(25), b(25), c(25) pi = 3.141592654# comenzar: PRINT PRINT PRINT "Opcion=1. Solucion analitica" PRINT "Opcion=2. Solucion numerica" PRINT PRINT

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

13

INPUT "Opcion"; opcion IF opcion = 1 THEN GOTO analitica IF opcion = 2 THEN GOTO numerica IF opcion < 1 OR opcion > 3 THEN GOTO comenzar REM ------------- SOLUCION ANALITICA DE LA EFICIENCIA ---------analitica: GOSUB datos PRINT "Opcion=1. Aleta radial" PRINT "Opcion=2. Aleta triangular" PRINT PRINT INPUT "Opcion"; opcion IF opcion = 1 THEN GOTO radial IF opcion = 2 THEN GOTO triangular IF opcion < 1 OR opcion > 3 THEN GOTO analitica REM ------ SOLUCION ANALITICA ALETA RADIAL-----radial: INPUT "hc="; hc INPUT "ra="; ra INPUT "rb"; rb INPUT "t="; bb INPUT "k="; kk INPUT "Ta="; ta INPUT "Too="; too m1 = (2 * hc / (kk * bb)) ^ .5 n0 = 0 x9 = m1 * ra GOSUB funi ii(0) = q GOSUB funk kk(0) = q n0 = 1 x9 = m1 * rb GOSUB funi ii(1) = q GOSUB funk kk(1) = q n0 = 0 GOSUB funi iii(0) = q GOSUB funk kkk(0) = q n0 = 1 x9 = m1 * ra GOSUB funi iii(1) = q GOSUB funk kkk(1) = q REM ----Evaluacion de constantes--a1 = (hc / (kk * m1)) a2 = (ii(1) + a1 * iii(0)) / (kk(1) - a1 * kkk(0)) tta = ta - too c1 = tta / (ii(0) + a2 * kk(0)) c2 = a2 * c1 REM Calculo de la eficiencia qt = (2 * hc * (ta - too) * 3.14159 * (rb ^ 2 - ra ^ 2)) + (hc * 2 * 3.14159 * rb * bb) * (ta - too)

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

14

qr = -kk * 2 * 3.14159 * ra * bb * m1 * (c1 * iii(1) - c2 * kkk(1)) nn = qr / qt PRINT "Eficiencia="; nn END REM ----SOLUCION ANALITICA ALETA TRIANGULAR---triangular: INPUT "hc="; hc INPUT "L="; l1 INPUT "t="; t INPUT "k="; kk INPUT "Ta="; ta INPUT "Too="; too b1 = 1 ll = ((l1 ^ 2 + (t / 2) ^ 2)) ^ .5 m2 = (8 * hc * l1 / (kk * t)) ^ .5 n0 = 0 x9 = m2 * (l1 ^ .5) GOSUB funi ii(0) = q n0 = 1 GOSUB funi ii(1) = q REM Calculo de la eficiencia qt = 2 * hc * ll * b1 * (ta - too) qr = b1 * ((2 * hc * kk * t) ^ .5) * (ta - too) * ii(1) / ii(0) nn = qr / qt PRINT "Eficiencia="; nn END REM ------------- SOLUCION NUMERICA DE LA EFICIENCIA -----------numerica: PRINT "Opcion=1. Aleta radial" PRINT "Opcion=2. Aleta triangular" PRINT PRINT INPUT "Opcion"; opcion IF opcion = 1 THEN GOTO radia IF opcion = 2 THEN GOTO triangu IF opcion < 1 OR opcion > 3 THEN GOTO numerica REM --------- ALETA RADIAL -----radia: nume = 1 INPUT "hc="; hc INPUT "ra="; ra INPUT "rb"; rb INPUT "t="; bb INPUT "k="; kk INPUT "Ta="; ta INPUT "Too="; too m1 = (2 * hc / (kk * bb)) ^ .5 INPUT "NUMERO DE DIVISIONES="; m% dr = -(rb - ra) / m% x1 = ta x(0) = x1 GOSUB integra

Página Análisis Conceptual, Modelamiento y Cálculo Numérico en Ingeniería Química. HRangel J.

15

t(0) = ta - y(1) PRINT t(0) x1 = too x(1) = x1 GOSUB integra t(1) = ta - y(1) PRINT t(1) REM ---- IMM aleta radial ---imm: aa(0) = x(0) aa(1) = (t(1) - t(0)) / (x(1) - x(0)) x(2) = aa(0) - (t(0) / aa(1)) mm% = 2 calculando: x1 = x(mm%) GOSUB integra t(mm%) = ta - y(1) b(0) = x(mm%) FOR i% = 1 TO mm% - 1 b(i%) = (t(mm%) - t(i% - 1)) / (b(i% - 1) - aa(i% - 1)) NEXT i% c(mm%) = (t(mm%) - t(mm% - 1)) / (b(mm% - 1) - aa(mm% - 1)) aa(mm%) = c(mm%) FOR i% = mm% TO 1 STEP -1 c(i% - 1) = aa(i% - 1) - (t(i% - 1) / c(i%)) NEXT i% x(mm% + 1) = c(0) x1 = x(mm% + 1) GOSUB integra t(mm% + 1) = ta - y(1) erro = ABS(t(mm% + 1)) IF erro