EQE038 – Simulação e Otimização de Processos Químicos
PFR: Métodos de Discretização
Argimiro R. Secchi
Programa de Engenharia Química – COPPE/UFRJ
Rio de Janeiro, RJ
EQ/UFRJ
março de 2014
1
Modelo Dinâmico de um Reator Tubular
Considerações:
•
Propriedades físicas constantes (, , cp);
•
Fluido Newtoniano;
•
Simetria angular;
•
vr = v = 0
2
Classificação dos modelos de acordo com os princípios físico-químicos
Modelo Microscópico
Modelo de Gradientes Múltiplos
Modelo de Máximo Gradiente
Modelo de Macroscópico
3
Modelo Microscópico:
(  .v )  0
C i
 vz (r , t )
t
 cp

vz  vz (r , t )
T
t
v
t
Ci
z
 (   C i v )  (   D
  c p v z ( r , t )
T
z
(l )
 C i )  Ri
  c p (   v  T )  (   k
  v   v   [   v  v ]   P  
(l )
(l )
 T )  (
(l )
v   v )  Sr
(l )
(t )
(t )
 v g
2
Modelo de turbulência (exemplo simples):
v  c i  J i
(t )
 D
 c p vT   q
 v v  
(t )
(t )
(t )
Ci
 k T
 
(t )
(t )
v
Computational Fluid
Dynamic (CFD)
4
Método de resolução recomendado
Volumes Finitos:
Consiste na realização de balanços de propriedades em volumes elementares
(volumes finitos), ou de forma equivalente na integração sobre o volume
elementar da equação diferencial na forma conservativa (ou forma divergente,
onde os fluxos aparecem dentro das derivadas).
 Uso de softwares de CFD (dinâmica de fluido computacional)
5
Exemplo simples do uso do método dos volumes finitos:
Equação da reação-difusão transiente em uma partícula catalítica esférica









C
t

1   2 C 
2
r
   C (r , t )
2
r r 
r 
C
r
reação de primeira ordem
 0 ; C (1, t )  1
r 0
C ( r , 0)  C 0
dV  4  r dr
2
e

w
C
e
C 

r 2 dr     r 2

t
r 

w
e

2
r 2 C dr
w
6
Exemplo simples do uso do método dos volumes finitos
e
 C r dr
2
Cp 
w
e

 r dr
2
e
3
r
2
e
r
2
w

 C r dr
valor médio
2
w
w
dC p
dt

e
3
r
3
e
dC p
dt
AW 
 rw
3

 2 C 
2
r


Cp


r  w

C
r

r  re
 AW CW  A p C p  AE C E   2 C p
3 rw2
 re3  rw3   rw
 re2
rw2 
Ap 



3
3
 re  rw    re  rw 
3
CE  C p
C
 re
r

r  rw
C p  CW
 rw
p = 2, ..., N  1
AE 
3 re2
 re3  rw3   re
7
Exemplo simples do uso do método dos volumes finitos
Condições de contorno
C
r
C
0
r
rw
Sistema resultante:
dC
 AC  b

re
Ce  C p
rf

1 Cp
rf
onde A é matriz tridiagonal
dt
Ou um sistema não-linear para reações de ordem  1:
dC
 F(C)
dt
8
Modelo de Gradientes Múltiplos:
D D
C i
t
 cp

vz
t
(t )
D
 vz (r , t )
T
t
k k
(l )
C i
z
(t )
k
(l )

(t )

(l )
 (  D  C i )  Ri
  c p v z ( r , t )
T
z
 (  k  T )    v  S r
  P    v z   g
2
Condições de Contorno:
vz , C i0 , T0
zo na de re aç ão
z= 0
z = 0¯
z= L
z = 0+
z = L¯
z = L+
9
Balanço de massa por componente (removendo a notação de média temporal):
C i
t
 Ci
2
 D z (r , t )
z
2
C i 
C i
1  

 Ri
 r D R (r , t )
  vz (r , t )
r r 
r 
z
v z ( r , t ) C i 0 ( t )  v z ( r , t ) C i (0, r , t )  D z ( r , t )
 Ci
z
 Ci
r
 Ci
r
 C i (0, r , t )
z
( L,r,t)  0
Sem reação na saída
( z , 0, t )  0
Simetria
(z, R ,t )  0
Parede impermeável
C i ( z , r , 0)  C i _ in ( z , r )
Condição inicial
10
Balanço de energia (removendo a notação de média temporal):
 cp
T
t
 T
2
 k z (r , t )
z
2

1  
T 
T
r
k
(
r
,
t
)


c
v
(
r
,
t
)
 H r R A
R
p
z


r r 
r 
z
v z ( r , t ) T 0 ( t )  v z ( r , t ) T (0, r , t ) 
T
z
T
r
k z ( r , t )  T (0, r , t )
 cp
z
(L, r, t)  0
Sem reação na saída
( z , 0, t )  0
Simetria
kR (R,t)
T
r
( z , R , t )  U Tw  T ( z , R , t ) 
T ( z , r , 0)  Tin ( z , r )
Troca de calor pela parede
Condição inicial
11
Balanço de quantidade de movimento (removendo a notação de média temporal):
vz
1   v z 



r

t
z
r r  r 
vz
r
P
(0, t )  0
Simetria
vz ( R , t )  0
Parede imóvel
v z ( r , 0)  v z _ in ( r )
Condição inicial
12
Simplificação adicional:
- coeficientes difusivos efetivos constantes
- velocidade constante
C i
t
 cp
2
 DL
T
t
v z C i 0 ( t )  v z C i (0, r , t )  D L
 Ci
z
 Ci
r
 Ci
r
( L,r,t)  0
 Ci
z

2
 T
C i
D R   C i 
r

v
 Ri
z


r r  r 
z
2
 kL
z
2

 C i (0, r , t )
z
k R   T 
T
r


c
v
 H r RA
p
z


r r  r 
z
v z T 0 ( t )  v z T (0, r , t ) 
T
z
( z , 0, t )  0
T
r
(z, R ,t )  0
C i ( z , r , 0)  C i _ in ( z , r )
kR
k L  T (0, r , t )
 cp
z
(L, r, t)  0
( z , 0, t )  0
T
r
( z , R , t )  U Tw  T ( z , R , t ) 
T ( z , r , 0)  Tin ( z , r )
13
Métodos de resolução de S.E.D.P.
Para simulação dinâmica o mais indicado é o método das linhas:
- discretização espacial (diferenças finitas, volumes finitos, elementos finitos)
- integração temporal
Exemplo: diferenças finitas na direção axial e colocação ortogonal na radial
dC i , j , k
dt
 cp
dT j , k
dt
onde:
 C i , j  1, k  2 C i , j , k  C i , j 1, k 
 C i , j  1, k  C i , j 1, k 
 n 1

 DL 

4
D
u
B

A
C

v

  Ri , j ,k
R   k
k ,m
k ,m 
i , j ,m 
z 
2
z
2z
 m 0





 T j  1, k  2 T j , k  T j 1, k 
 T j  1, k  T j 1, k 
 n 1

 kL 
  4 k R    u k B k , m  Ak , m  T j , m    c p v z 
   H r R A , j ,k
2
z
2z
 m0





n 1
u r
2
lm (u ) 
u  up
u
p0
m
 up
Ak . m 
2
dl m ( u k )
Bk ,m 
du
pm
j = 1,2,...,N
k = 1,2,...,n
(  , )
n 1
y (u )  P
d lm (u k )
n 1
(u ) 
du
l
m 0
m
2
(u ) y m
14
Métodos de resolução de S.E.D.P.
Condições de contorno
 C i ,2 ,k  C i 0 (t ) 
v z C i 0 ( t )  v z C i ,1, k  D L 

2z


v z T0 ( t )  v z T1, k 
C i , N 1, k  C i , N , k  0
T N  1, k  T N , k  0
n 1

n 1

A0 , m C i , j , m  0
m0
A0 , m T j , m  0
m0
n 1
n 1

k L  T 2 , k  T0 ( t ) 


 cp 
2z

k R  An  1, m T j , m  U  T w  T j , n 1 
An  1, m C i , j , m  0
m0
m 0
C i , j , k  C i _ in , j , k
T j , k  Tin , j , k
j = 1,2,...,N
k = 1,2,...,n
Resulta em um sistema de equações algébrico-diferenciais, ou
somente diferenciais se as condições de contorno (lineares)
forem incorporadas nas equações diferenciais
15
Outro Modelo de Gradientes Múltiplos: (ignorando gradientes radiais)
PFR com dispersão axial
R
C i
t
 Ci
2
 DL
z
2
 vz
C i
z
 Ri
v z C i 0 ( t )  v z C i (0, t )  D L
C i
z
Ci ( z, t) 
C
i
( r , z , t ) r dr
0
R
 r dr
 C i (0, t )
0
z
(L, t)  0
C i ( z , 0)  C i _ in ( z )
R
 cp
T
t
 T
2
 kL
z
2
  c p vz
v z T 0 ( t )  v z T (0, t ) 
T
z
T
z

2U
R
 Tw  T    H r
k L  T (0, t )
 cp
RA
T ( z, t) 
 T (r , z, t ) r dr
0
R
 r dr
0
z
(L, t)  0
T ( z , 0)  Tin ( z )
16
Métodos de resolução do S.E.D.P.
Método das linhas:
- discretização espacial (D.F., V.F., colocação ortogonal em elementos finitos)
- integração temporal
Exemplo: diferenças finitas na direção axial
dC i , j
dt
 cp
 C i , j  1  2 C i , j  C i , j 1 
 C i , j  1  C i , j 1 
 DL 
  vz 
  Ri , j
2

z
2

z




 T j  1  2 T j  T j 1
 kL 
2
dt
z

dT j
v z C i 0 ( t )  v z C i ,1
 C i ,2  C i 0 ( t ) 
 DL 

2z


C i , N 1  C i , N  0
C i , j  C i _ in , j

 T j  1  T j 1
   c p vz 
2z


j = 1,2,...,N
 2U
 Tw  T j    H r R A , j

R

v z T0 ( t )  v z T1 
k L  T 2  T0 ( t ) 


 cp 
2z

T N 1  T N  0
T j  Tin , j
17
Modelo de Gradiente Máximo: (ignorando dispersão axial)
PFR sem dispersão axial
C i
t
 vz
C i
z
 Ri
C i (0, t )  C i 0 ( t )
C i ( z , 0)  C i _ in ( z )
 cp
T
t
  c p v z
T
z

2U
R
 Tw  T    H r
RA
T (0, t )  T0 ( t )
T ( z , 0)  Tin ( z )
Resulta em um sistema de equações diferenciais ordinárias ao
aplicar o método das linhas
18
Modelo de Macroscópico:
L
V
dC i
dt
 C i 0 (t ) v z S  C i v z S  Ri V
C i (t ) 
C
i
( z , t ) S dz
0
L
 S dz
0
C i (0 )  C i _ in
L
 cp V
dT
dt
  c p v z S T 0 ( t )   c p v z S T  U At ( T w  T )   H r R A V
T (t ) 
 T ( z , t ) S dz
0
L
 S dz
T (0)  Tin
0
Resulta em um sistema de equações diferenciais ordinárias
19
Esquema da aplicação
do MCO ao EMSO
• Contornos
• Pontos Internos
• Alfa e Beta
Simulador
EMSO
Plugin
Rotina de cálculo
das raízes e matrizes
• Raízes de Jacobi
• Matriz A e B
DD as Plugin (Type=“OCFEM”, Boundary="BOTH”, InternalPoints=5
alfa=1, beta=1)
Plugin: ocfem_emso.dll
20
Reator de leito fixo com dispersão axial
(reação de ordem m)
y


y
x
1  y
2

Pe  x
2
 Da y
m
Condições de contorno:

1 y
P e x
y
x
 1  y (  , 0)
ou
y (  , 0)  1
x0
0
x 1
Condição inicial:
y (0, x )  0
21
Sistema de Equações Diferenciais Parciais
– Método das Linhas com D.F. e Colocação Ortogonal –
Exemplo: adicionar o Plugin
ocfem_emso.dll e executar os flowsheets
dos arquivos FDM_ss.mso, OCM_ss.mso
e OCFEM_ss.mso e comparar os
resultados das discretizações. Fazer o
mesmo para o caso transiente nos
arquivos FDM_din.mso e OCM_din.mso.
22
Comparando os Resultados
1.2000
1.2000
1.0000
1.0000
0.8000
0.8000
0.6000
0.6000
0.4000
0.4000
0.2000
0.2000
0.0000
0.0000
0.2000
0.4000
0.6000
0.8000
1.0000
0.0000
0.0000
0.2000
0.4000
0.6000
0.8000
1.0000
MCO pelo EMSO   1   1
Número de pontos internos: 5
Método das Diferenças Finitas
Número de pontos internos: 6000
y(x=1) = 0.151475 (erro de 0.038%)
y(x=1) = 0.15155 (erro de 0.087%)
y(x=1) = 0.151418 (exato)
23
Estudo de Caso
• Produção de anidrido acético em reator PFR adiabático
– O anidrido acético, em especial, é um dos mais importantes produtos
químicos, e freqüentemente é preparado pela reação de um ácido acético com
um composto chamado ceteno, obtido a partir do aquecimento da acetona a
700-770oC.
– Etapa importante é o craqueamento em fase vapor de acetona para ceteno e
metano:
C H 3C O C H 3  C H 2 C O  C H 4
– A segunda etapa é a reação do ceteno com ácido acético.
C H 2 C O  C H 3C O O H   C H 3C O  O
Ref: G. V. Jeffreys, A Problem in Chemical Engineering Design: The Manufacture of
Acetic Anhydride, 2nd ed. (London: Institution of Chemical Engineers, 1964)
24
Estudo de Caso
• Definição do problema
– A primeira etapa da produção é realizada através de uma reação
na fase vapor de acetona em um reator tubular (PFR) adiabático.
onde A = acetona; B = ceteno e C = metano
A  B C
– A reação é de 1a ordem em relação à cetona na reação de
craqueamento, com constante de Arrhenius dada por:
• k – segundos
• T – Kelvin
34222 

k  exp  34.34 

T


25
Estudo de Caso
• Descrição do Processo
– Geometria do reator
• um reator tubular adiabático contínuo;
• contém um banco de 1000 tubos de 1 in sch. 40 correspondendo uma
seção transversal de 0,557 m2;
• com comprimento total de 2,28 m;
– Condições de operação
• temperatura da carga 762oC (1035 K);
• pressão de operação de 1,6 atm
• vazão de 8000 kg/h (137,9 kmol/h);
– Composição
• Acetona (ACETONE), Ceteno (KETENE) e Metano (METHANE)
• alimentação de acetona pura
– Cinética
•
•
•
•
ocorre uma reação de primeira ordem,
fator pré-exponencial (k0): 8,2 x 1014 s-1
energia de ativação (E/R): 34222 K
calor da reação de -80,77 kJ/mol
26
Estudo de Caso
– Produção de anidrido acético –
Exemplo: executar o FlowSheet do
arquivo PFR_Adiabatico.mso e construir
os perfis de temperatura e composição no
estado estacionário. Mostrar também a
evolução do perfil de temperatura da
condição inicial até o estado estacionário.
Discutir o tipo e a qualidade da
discretização usada neste exemplo.
27
Exercícios
1) Resolver o problema da reação-difusão em uma partícula
catalítica esférica, descrito por:
y
1 y  2 y 
2 1/ 2
 2
r


y


t
r r 
r 
y
r
0
y (t , r )
r 1
1
y (t , r ) t  0  0
r 0
 2
(m ódulo de T hiele)
2) Apresentar os resultados do exemplo do slide 27.
28
Bibliografia
• Himmelblau, D. M. & Bischoff, K. B., "Process Analysis and Simulation - Deterministic Systems", John Wiley &
Sons, 1968.
• Finlayson, B. A., "The Method of Weighted Residuals and Variational Principles with Application in Fluid Mechanics,
Heat and Mass Transfer", Academic Press, 1972.
• Villadsen, J. & Michelsen, M. L., "Solution of Differential Equation Models by Polynomial Approximation", PrenticeHall, 1978.
• Davis, M. E., "Numerical Methods and Modeling for Chemical Engineers", John Wiley & Sons, 1984.
• Denn, M., "Process Modeling", Longman, New York, 1986.
• Luyben, W. L., "Process Modeling, Simulation, and Control for Chemical Engineers", McGraw-Hill, 1990.
• Silebi, C.A. & Schiesser, W.E., “Dynamic Modeling of Transport Process Systems”, Academic Press, Inc., 1992.
• Biscaia Jr., E.C. “Método de Resíduos Ponderados com Aplicação em Simulação de Processos”, XV CNMAC,
1992
• Ogunnaike, B.A. & Ray, W.H., “Process Dynamics, Modeling, and Control”, Oxford Univ. Press, New York, 1994.
• Rice, R.G. & Do, D.D., “Applied Mathematics and Modeling for Chemical Engineers”, John Wiley & Sons, 1995.
• Maliska, C.R. “Transferência de Calor e Mecânica dos Fluidos Computacional”, 1995.
• Bequette, B.W., “Process Dynamics: Modeling, Analysis, and Simulation”, Prentice Hall, 1998.
• Fogler, H.S., “Elementos de Engenharia de Reações Químicas”, Prentice Hall, 1999.
29