Apostila 2 - método de Runge-Kutta - 2015

Prévia do material em texto

1 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
 
UNIVERSIDADE DE RIBEIRÃO PRETO 
CURSO DE ENGENHARIA QUÍMICA 
 
CÁLCULO DE REATORES 2 
 
 
APOSTILA 2 
 
MÉTODO DE RUNGE-KUTTA APLICADO A REATORES 
QUÍMICOS 
 
 
 
 
 
 
 
 
Prof. Dr. Murilo Daniel de Mello Innocentini 
Curso de Engenharia Química 
Universidade de Ribeirão Preto – UNAERP 
Currículo Lattes: http://lattes.cnpq.br/5681181471077426 
 
 
 
 
 
 
RIBEIRÃO PRETO - SP 
AGOSTO - 2015 
 
 2 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
1. Introdução 
 
 Existem diversos tipos de sistemas em Engenharia Química que podem ser 
representados por problemas de valor inicial. Um reator em batelada pode ser descrito a partir 
das concentrações em um instante de tempo (tipicamente t = 0) e das equações de balanço 
material e de energia. Um problema típico consiste em determinar o tempo necessário para se 
obter uma dada concentração de produto. A variável independente tempo (t) aparece 
freqüentemente em problemas deste tipo. Outras variáveis podem ser utilizadas; um exemplo 
é a determinação do perfil de temperatura T = f(z) ao longo de um trocador de calor, que pode 
ser feita a partir da temperatura de entrada (em z = 0) e das equações de transferência de calor. 
Neste caso, a variável independente é a posição ao longo do eixo do trocador de calor. 
Um exemplo extremamente simples de um problema de valor inicial é 
 
 
)y,t(f
dt
dy

 (1) 
Com condição inicial: t = 0  y = yo 
 
 Em aplicações práticas, muitas vezes necessitamos de mais de uma variável para 
descrever o sistema. Muitos sistemas em Engenharia Química podem ser descritos por 
sistemas de equações diferenciais do tipo: 
 
)t,...z,y,x(f
dt
dx

 (2) 
 
)t,...z,y,x(g
dt
dy

 (3)
 
)t,...z,y,x(w
dt
dz

 (4) 
Com condição inicial: t = 0  x = xo, y = yo, z = zo, ... 
 
 Observe que, em sistemas com condição inicial em t = to diferente de zero, basta 
efetuar a mudança de variável  = t-to para obter a condição inicial em  = 0. 
 Dentre os métodos para a resolução numérica de uma integral, o método de Runge-
Kutta tem a vantagem de ser relativamente simples e de dar resultados mais precisos do que 
outros métodos. Embora as demonstrações não sejam realizadas aqui, sugere-se a leitura de 
material didático específico para a compreensão do método (Moderna Introdução às 
Equações Diferenciais, Richard Bronson, Coleção Schaum, McGraw-Hill, 1977, pg. 285). 
 
 3 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
 
Figura 1. representação gráfica do método de Runge-Kutta de 4ª ordem. Em cada passo a 
derivada é calculada 4 vezes: uma vez no ponto inicial, duas vezes em pontos intermediários e 
uma vez em um ponto (estimado) no final. Dessas derivadas o valor da função final (mostrado 
na figura como um ponto preenchido) é calculado. 
 
2. Método de Runge-Kutta de 4ª ordem para 1 equação diferencial 
 
Considere que a função derivada possa ser representada genericamente por: 
 
 
)y;x(f
dx
dy

 (5)
 
sendo x a variável independente. 
 
A solução numérica para essa equação é dada por: 
 
 
 4321n1n kk2k2k
6
1
yy 
 (6) 
sendo: 
 
)y;x(hfk nn1 
 (7) 
 
)
2
k
y;
2
h
x(hfk 1nn2 
 (8) 
 
)
2
k
y;
2
h
x(hfk 2nn3 
 (9) 
 
)ky;hx(hfk 3nn4 
 (10) 
E h é o intervalo (passo) na variável independente x: 
 
hxx n1n 
 (11) 
 
 4 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
As equações (6) a (11) são aplicadas quantas vezes necessárias para a obtenção da resposta. 
 
Exemplo 1: considere a seguinte função:
2xy
dx
dy

, com condição inicial: x = 0  y = 1. 
Determine graficamente a solução y  x no intervalo x entre 0 e 1. Compare com a solução 
analítica. Use passo h = 0,2. 
 
Resolução: 
 
Neste caso: 
2xy)y;x(f 
, com xo = 0 e yo = 1. Pelas equações (7) a (11): 
 
Para xo = 0 e yo = 1: 
 
0)1)(0)(2,0()1;0(hf)y;x(hfk 2oo1 
 
02,0)1)(1,0)(2,0()1;1,0(hf)
2
0
1;
2
2,0
0(hf)
2
k
y;
2
h
x(hfk 21oo2 
 
020402,0)01,1)(1,0)(2,0()01,1;1,0(hf)
2
02,0
1;
2
2,0
0(hf)
2
k
y;
2
h
x(hfk 22oo3 
 
042,0)0204,1)(2,0)(2,0()0204,1;2,0(hf)0204,01;2,00(hf)ky;hx(hfk 23oo4 
 
 4321o1o kk2k2k
6
1
yy 
 
 042,0)0204,0(2)02,0(20
6
1
1y1 
 
 
020,1y1 
 
hxx o1o 
 
2,00x1 
 
2,0x1 
 
Repetindo os cálculos: 
 
xo= 0 x1= 0.2 x3= 0.4 x4= 0.6 x5= 0.8
yo= 1.000 y1= 1.020 y3= 1.087 y4= 1.220 y5= 1.471
k1= 0.000 k1= 0.042 k1= 0.095 k1= 0.178 k1= 0.346
k2= 0.020 k2= 0.065 k2= 0.129 k2= 0.240 k2= 0.486
k3= 0.020 k3= 0.067 k3= 0.133 k3= 0.251 k3= 0.529
k4= 0.042 k4= 0.095 k4= 0.178 k4= 0.346 k4= 0.799
x1= 0.2 x2= 0.4 x4= 0.6 x5= 0.8 x6= 1
y1= 1.020 y2= 1.087 y4= 1.220 y5= 1.471 y6= 2.000
 
 
 5 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
A solução analítica para esse problema é dada por: 
 
2xy
dx
dy

 
xdx
y
dy
2


 
x
x
y
y
2
oo
xdx
y
dy
 
2
x
2
x
y
1
y
1
2
o
2
o

 
2
x
2
x
y
1
1
y
2
o
2
o


 
2
0
2
x
1
1
1
y
22


 
2x2
2
y


 
Comparando as soluções numérica e analítica: 
x y num y ana Erro (%)
0.0 1.000 1.000 0.0000
0.2 1.020 1.020 0.0400
0.4 1.087 1.087 -0.0040
0.6 1.220 1.220 -0.0400
0.8 1.471 1.471 -0.0280
1.0 2.000 2.000 0.0000
 
0.00
0.50
1.00
1.50
2.00
2.50
0.0 0.2 0.4 0.6 0.8 1.0
x
y
Numérico
Analítico
 
3. Método de Runge-Kutta de 4ª ordem para sistema de 2 equações diferenciais 
 
Considere o seguinte sistema de equações diferenciais: 
 
 
)z;y;x(g
dx
dz
)z;y;x(f
dx
dy


 (12) 
 
onde x é a variável independente. 
 
 6 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
A solução numérica para esse sistema é dada por: 
 
 4321n1n kk2k2k
6
1
yy 
 (13) 
 
 4321n1n ll2l2l
6
1
zz 
 (14) 
sendo: 
 
)z;y;x(hfk nnn1 
 (15) 
 
)z;y;x(hgl nnn1 
 (16) 
 
)
2
l
z;
2
k
y;
2
h
x(hfk 1n
1
nn2 
 (17) 
 
)
2
l
z;
2
k
y;
2
h
x(hgl 1n
1
nn2 
 (18) 
 
)
2
l
z;
2
k
y;
2
h
x(hfk 2n
2
nn3 
 (19) 
 
)
2
l
z;
2
k
y;
2
h
x(hgl 2n
2
nn3 
 (20) 
 
)lz;ky;hx(hfk 3n3nn4 
 (21) 
 
)lz;ky;hx(hgl 3n3nn4 
 (22) 
 
E h é o intervalo (passo) na variável independente x: 
 
 
hxx n1n 
 (23) 
 
As equações (13) a (23) são aplicadas quantas vezes necessárias para a obtenção da resposta. 
 
Exemplo 2: considere o seguinte sistema de equações diferenciais: 
 
xz2
dx
dy

 
 
2x6
dx
dz

 
com condição inicial: x = 0  y = 0 e z = 2. Determine a solução y  x e z  x no intervalo x 
entre 0 e 1. Compare com a solução analítica. Use passo h = 0,1. 
 
 
 
 
 
 7 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
Resolução: 
Neste caso: 
xz2)z;y;x(f 
 e 
2x6)z;y;x(g 
com xo = 0; yo = 0; zo = 2. 
 
Pelas equações (13) a (23): 
 
Para xo = 0, yo = 0 e zo = 2: 
 
0)2)(0)(1,0()2;0;0(hf)z;y;x(hfkooo1 
 
 
0)0(6)1,0()2;0;0(hg)z;y;x(hgl 2ooo1 
 
 
02,0)2)(05,0(2)1,0()2;0;05,0(hf)
2
0
2;
2
0
0;
2
1,0
0(hf)
2
l
z;
2
k
y;
2
h
x(hfk 1o
1
oo2 
 
0015,0)05,0(6)1,0()2;0;05,0(hg)
2
0
2;
2
0
0;
2
1,0
0(hg)
2
l
z;
2
k
y;
2
h
x(hgl 21o
1
oo2 
 
)00075,2;01,0;05,0(hf)
2
0015,0
2;
2
02,0
0;
2
1,0
0(hf)
2
l
z;
2
k
y;
2
h
x(hfk 2o
2
oo3 
 
 
02,0)00075,2)(05,0(2)1,0(k3 
 
 
)00075,2;01,0;05,0(hg)
2
0015,0
2;
2
02,0
0;
2
1,0
0(hg)
2
l
z;
2
k
y;
2
h
x(hgl 2o
2
oo3 
 
 
0015,0)05,0(6)1,0(l 23 
 
 
)0015,2;02,0;1,0(hf)0015,02;02,00;1,00(hf)lz;ky;hx(hfk 3o3oo4 
 
04,0)0015,2)(1,0(2)1,0(k 4 
 
 
)0015,2;02,0;1,0(hg)0015,02;02,00;1,00(hg)lz;ky;hx(hgl 3o3oo4 
 
006,0)1,0(6)1,0(l 24 
 
 
 
 8 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
 4321o1o kk2k2k
6
1
yy 

 04,0)02,0(2)02,0(20
6
1
0y1 
  
02,0y1 
 
 4321o1o ll2l2l
6
1
zz 

 006,0)0015,0(2)0015,0(20
6
1
2z1 
 
002,2z1 
 
hxx o1o 
 
1,00x1 
 
1,0x1 
 
 
Repetindo os cálculos: 
 
xo= 0 x1= 0.1 x2= 0.2 x3= 0.3 x4= 0.4 x5= 0.5 x6= 0.6 x7= 0.7 x8= 0.8 x9= 0.9
yo= 0.000 y1= 0.020 y2= 0.080 y3= 0.182 y4= 0.327 y5= 0.523 y6= 0.779 y7= 1.108 y8= 1.531 y9= 2.075
zo= 2.000 z1= 2.002 z2= 2.016 z3= 2.054 z4= 2.128 z5= 2.25 z6= 2.432 z7= 2.686 z8= 3.024 z9= 3.458
k1= 0.0000 k1= 0.0400 k1= 0.0806 k1= 0.1232 k1= 0.1702 k1= 0.2250 k1= 0.2918 k1= 0.3760 k1= 0.4838 k1= 0.6224
l1= 0.0000 l1= 0.0060 l1= 0.0240 l1= 0.0540 l1= 0.0960 l1= 0.1500 l1= 0.2160 l1= 0.2940 l1= 0.3840 l1= 0.4860
k2= 0.0200 k2= 0.0602 k2= 0.1014 k2= 0.1457 k2= 0.1958 k2= 0.2558 k2= 0.3302 k2= 0.4250 k2= 0.5467 k2= 0.7032
l2= 0.0015 l2= 0.0135 l2= 0.0375 l2= 0.0735 l2= 0.1215 l2= 0.1815 l2= 0.2535 l2= 0.3375 l2= 0.4335 l2= 0.5415
k3= 0.0200 k3= 0.0603 k3= 0.1017 k3= 0.1464 k3= 0.1970 k3= 0.2575 k3= 0.3326 k3= 0.4282 k3= 0.5509 k3= 0.7085
l3= 0.0015 l3= 0.0135 l3= 0.0375 l3= 0.0735 l3= 0.1215 l3= 0.1815 l3= 0.2535 l3= 0.3375 l3= 0.4335 l3= 0.5415
k4= 0.0400 k4= 0.0804 k4= 0.1221 k4= 0.1673 k4= 0.2189 k4= 0.2809 k4= 0.3582 k4= 0.4568 k4= 0.5833 k4= 0.7458
l4= 0.0060 l4= 0.0240 l4= 0.0540 l4= 0.0960 l4= 0.1500 l4= 0.2160 l4= 0.2940 l4= 0.3840 l4= 0.4860 l4= 0.6000
x1= 0.1 x2= 0.2 x3= 0.3 x4= 0.4 x5= 0.5 x6= 0.6 x7= 0.7 x8= 0.8 x9= 0.9 x10= 1.0
y1= 0.020 y2= 0.080 y3= 0.182 y4= 0.327 y5= 0.523 y6= 0.779 y7= 1.108 y8= 1.531 y9= 2.075 y10= 2.773
z1= 2.002 z2= 2.016 z3= 2.054 z4= 2.128 z5= 2.250 z6= 2.432 z7= 2.686 z8= 3.024 z9= 3.458 z10= 4.000
 
 
A solução analítica para esse sistema é dada por: 
2x6
dx
dz

 
dxx6dz 2

 
x
x
2
z
z oo
dxx6dz
 









3
x
3
x
6zz
3
o
3
o
 









3
0
3
x
62z
33  
2x2z 3 
 
 
xz2
dx
dy


)2x2(x2
dx
dy 3 

x4x4
dx
dy 4 

 dxx4x4dy 4 
 
 
  
x
x
4
y
y oo
dxx4x4dy
 
 
 
 9 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 









2
x
5
x
2
x
5
x
4yy
2
o
5
o
25
o










2
0
5
0
2
x
5
x
40y
2525 









2
x
5
x
4y
25 
 
 
25 x2x8,0y 
 
 
Comparando as soluções numérica e analítica: 
 
x y num z num y ana z ana Erro y (%)Erro z (%)
0.00 0.0000 2.0000 0.0000 2.0000 0.000 0
0.10 0.0200 2.002 0.0200 2.002 0.015 0
0.20 0.0802 2.016 0.0803 2.016 0.060 2E-14
0.30 0.1817 2.054 0.1819 2.054 0.130 2E-14
0.40 0.3275 2.128 0.3282 2.128 0.222 4E-14
0.50 0.5233 2.250 0.5250 2.250 0.331 0
0.60 0.7787 2.432 0.7822 2.432 0.455 4E-14
0.70 1.1079 2.686 1.1145 2.686 0.584 3E-14
0.80 1.5311 3.024 1.5421 3.024 0.714 3E-14
0.90 2.0749 3.458 2.0924 3.458 0.837 3E-14
1.00 2.7735 4.000 2.8000 4.000 0.948 0
 
0.0
0.5
1.0
1.5
2.0
2.5
3.0
0.0 0.2 0.4 0.6 0.8 1.0
x
y
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
z
analítico
numérico
 
 
 
 
 
 
 
 10 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
4. Método de Runge-Kutta de 4ª ordem para sistema de 3 equações diferenciais 
 
Considere o seguinte sistema de equações diferenciais: 
 
)w;z;y;x(j
dx
dw
)w;z;y;x(g
dx
dz
)w;z;y;x(f
dx
dy



 (24)
 
onde x é a variável independente. 
 
A solução numérica para esse sistema é dada por: 
 
 
 4321n1n kk2k2k
6
1
yy 
 (25) 
 
 4321n1n ll2l2l
6
1
zz 
 (26) 
 
 4321n1n mm2m2m
6
1
ww 
 (27) 
onde: 
 
)w;z;y;x(hfk nnnn1 
 (28) 
 
)w;z;y;x(hgl nnnn1 
 (29) 
 
)w;z;y;x(hjm nnnn1 
 (30) 
 
)
2
m
w;
2
l
z;
2
k
y;
2
h
x(hfk 1n
1
n
1
nn2 
 (31) 
 
)
2
m
w;
2
l
z;
2
k
y;
2
h
x(hgl 1n
1
n
1
nn2 
 (32) 
 
)
2
m
w;
2
l
z;
2
k
y;
2
h
x(hjm 1n
1
n
1
nn2 
 (33) 
 
)
2
m
w;
2
l
z;
2
k
y;
2
h
x(hfk 2n
2
n
2
nn3 
 (34) 
 
)
2
m
w;
2
l
z;
2
k
y;
2
h
x(hgl 2n
2
n
2
nn3 
 (35) 
 
)
2
m
w;
2
l
z;
2
k
y;
2
h
x(hjm 2n
2
n
2
nn3 
 (36) 
 
)mw;lz;ky;hx(hfk 3n3n3nn4 
 (37) 
 
)mw;lz;ky;hx(hgl 3n3n3nn4 
 (38) 
 
)mw;lz;ky;hx(hjm 3n3n3nn4 
 (39) 
E h é o intervalo (passo) na variável independente x: 
 
 11 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
 
 
hxx n1n 
 (40) 
 
As equações (24) a (40) são aplicadas quantas vezes necessárias para a obtenção da resposta. 
 
 
Exemplo 3: A reação irreversível em fase líquida elementar 2 A  B deve ser realizada em 
um reator descontínuo de volume útil 20 litros. Se o reator for carregado com uma solução 3 
M do reagente A puro, então qual será a concentração após 20 minutos? Compare as 
soluções numérica e analítica. Considere que a constante de velocidade para essa reação é k2 
= 0,1 L/mol.min. 
 
Resolução: 
 
Para essa situação: 
 
equação do reator: 
 
)r(
dt
dC
A
A 
 (1.1) 
equação da reação: 
 
2
A2A Ck)r( 
 (1.2) 
 
de (1.1) e (1.2): 
 
2
A2
A Ck
dt
dC

  
2
A2
A Ck
dt
dC

 (1.3) 
 
com condição inicial: t = 0  CA = CAo = 3 M 
 
 
 
 
 12 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
Método numérico: 
 
fazendo x = t e y = CA e substituindo o valor de k2 temos em (1.3): 
 
 
2y1,0
dx
dy

 com condição inicial: x = 0  y = 3 
 
Aplicando o método de Runge-Kutta de 4ª ordem (para 1 equação diferencial) com passo h = 
0,2 min, temos como exemplo as 5 primeiras iterações: 
 
xo= 0 x1= 0.2 x3= 0.4 x4= 0.6 x5= 0.8
yo= 3.000 y1= 2.830 y3= 2.679 y4= 2.542 y5= 2.419
k1= -0.180 k1= -0.160 k1= -0.143 k1= -0.129 k1= -0.117
k2= -0.169 k2= -0.151 k2= -0.136 k2= -0.123 k2= -0.111
k3= -0.170 k3= -0.152 k3= -0.136 k3= -0.123 k3= -0.112
k4= -0.160 k4= -0.143 k4= -0.129 k4= -0.117 k4= -0.107
x1= 0.2 x2= 0.4 x4= 0.6 x5= 0.8x6= 1.0
y1= 2.830 y2= 2.679 y4= 2.542 y5= 2.419 y6= 2.308
 
 
O valor de CA para t = 20 min é, de acordo com o método numérico: CA = 0,429 M. 
 
Método analítico: 
 
Pela equação (1.3): 
 
2
A2
A Ck
dt
dC

 
dtk
C
dC
22
A
A 
 
 
t
0
2
C
C
2
A
A dtk
C
dCA
Ao
 
 
 
tk
C
1
C
1
2
AoA

 
Ao
2
A
C
1
tk
1
C


  
3
1
t1,0
1
CA


 
 
 
 
 13 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
Para t = 20 min: 
 
 
3
1
)20(1,0
1
CA


  
M429,0CA 
 
 
A figura a seguir apresenta a concordância excepcional entre os perfis de concentração 
em função do tempo para os métodos de resolução analítico e numérico. 
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
0 5 10 15 20 25 30
Tempo (min)
Co
nc
en
tra
çã
o, 
CA
 (M
) Numérico
Analítico
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 14 
Cálculo de Reatores 2 –UNAERP - Prof. Murilo D.M. Innocentini 
Lista de exercícios 
 
1) A reação irreversível de segunda ordem 2 A  B + C é realizada isotermicamente em um 
reator batelada (BSTR). 20 litros de uma solução contendo 40 moles de reagente A puro são 
alimentados no reator. Sabendo-se que a constante de velocidade da reação é k2 = 0,5 
L/mol.min, determine por método numérico a concentração do reagente A no reator após 15 
minutos. Compare a solução com aquela obtida analiticamente. 
 
2) Deseja-se realizar a reação irreversível elementar A + B  C em um reator semi-batelada. 
Uma solução de 10 litros do reagente A em concentração 2 molar deverá ser inicialmente 
introduzida no reator. Na seqüência, o reagente B será alimentado, de modo contínuo e lento 
no reator, em uma vazão de 100 mL/min, correspondendo a uma vazão molar de B de 1 
mol/min durante 20 minutos. Após esse período de alimentação, o reator é mantido sob 
agitação durante outros 20 minutos. Determine o perfil de concentração de A e de B no reator 
durante os 40 minutos de reação. Compare esse perfil com aquele caso todos os reagentes 
fossem alimentados de uma vez só no reator (operação em batelada pura). Considere que a 
constante de velocidade é dada por: k2 = 0,3 L/mol.min. Use passo para resolução numérica h 
= 0,5 min. 
 
3) A produção de brometo de metila (CH3Br) ocorre pela reação irreversível elementar em 
fase líquida: CNBr + CH3NH2  CH3Br + NCNH2, que é realizada em um reator semi-
batelada. Uma solução aquosa de metil amina (CH3NH2) em concentração de 0,025 mol/L 
deve ser alimentada em vazão de 0,05 L/s em uma solução aquosa de CNBr contido em um 
reator de vidro. O volume inicial do fluido no recipiente é de 5 L com CNBr em concentração 
de 0,05 mol/L. A constante de velocidade da reação é k2 = 2,2 L/mol.s. Pede-se: 
 
a) Obtenha os perfis de número de moles de reagentes e produtos (NA, NB, NC e ND) no reator 
em função do tempo (400 segundos). 
b) Obtenha os perfis de concentração de reagentes e produtos (CA, CB, CC e CD) no reator em 
função do tempo (400 segundos). 
c) Obtenha o perfil de grau de conversão (xA) do reagente CNBr no reator em função do 
tempo (400 segundos). 
d) Obtenha o perfil de taxa de reação (-rA) no reator em função do tempo (400 segundos). 
 
Utilize o método de Runge-Kutta com passo h = 5 s. Pode ser utilizada programação 
FORTRAN, BASIC, C++, PASCAL ou planilha Excel. Apresente todos os gráficos 
separadamente em planilha Excel.