transferencia de massa

Apostila de Cálculo de Reatores II

Lindoval Domiciano Fernandes

25 de janeiro de 2008

ii

Sumário

1 Introdução 11.1 Modelos de Adsorção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.1 Modelo de Langmuir . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.2 Modelo BET . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Cinética de Reações Catalíticas Heterogêneas . . . . . . . . . . . . . . . . 41.2.1 Cinética de Langmuir-Hinshelwood-Hougen-Watson . . . . . . . . . 41.2.2 Cinética de Eley-Rideal . . . . . . . . . . . . . . . . . . . . . . . . 16

2 Resistência à Transferência de Massa Externa 212.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.2 Resistência Externa à Transferência de Massa . . . . . . . . . . . . . . . . 21

2.2.1 Reação de Primeira Ordem . . . . . . . . . . . . . . . . . . . . . . 222.2.2 Reação de Segunda Ordem . . . . . . . . . . . . . . . . . . . . . . 23

3 Resistência à Transferência de Massa Interna 253.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2 Sistemas Isotérmicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.2.1 Resolução Analítica da Difusão em Ausência de Resistência Externa 273.2.2 Resolução Analítica Difusão com Resistência Externa . . . . . . . . 313.2.3 Resolução por Diferenças Finitas . . . . . . . . . . . . . . . . . . . 363.2.4 Resolução por Colocação Ortogonal . . . . . . . . . . . . . . . . . . 39

3.3 Sistemas Não-isotérmicos . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.3.1 Resolução por Diferenças Finitas . . . . . . . . . . . . . . . . . . . 453.3.2 Resolução por Colocação Ortogonal . . . . . . . . . . . . . . . . . . 48

4 Reatores de Leito Fixo 534.1 Modelo Unidimensional Básico . . . . . . . . . . . . . . . . . . . . . . . . 534.2 Modelo Pseudo-Homogêneo Bidimensional . . . . . . . . . . . . . . . . . . 56

4.2.1 Resolução por Colocação Ortogonal . . . . . . . . . . . . . . . . . . 594.2.2 Resolução por Diferenças Finitas . . . . . . . . . . . . . . . . . . . 624.2.3 Resolução por Colocação Ortogonal . . . . . . . . . . . . . . . . . . 694.2.4 Resolução por Diferenças Finitas . . . . . . . . . . . . . . . . . . . 73

iv SUMÁRIO

5 Reatores Não Ideais 815.1 Distribuição do Tempo de Residência . . . . . . . . . . . . . . . . . . . . . 81

5.1.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 815.1.2 Distribuição do Tempo de Residência . . . . . . . . . . . . . . . . . 815.1.3 Medida da Distribuição do Tempo de Residência . . . . . . . . . . 815.1.4 Características da Distribuição do Tempo de Residência . . . . . . 835.1.5 Função Distribuição Tempo de Residência Normalizada . . . . . . 865.1.6 Distribuição Interna do Tempo de Residência . . . . . . . . . . . . 86

5.2 A Distribuição do Tempo de Residência para Reatores Ideais . . . . . . . 875.2.1 Distribuição do Tempo de Residência para Reatores Batelada e

Plug Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875.2.2 CSTR Simples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 885.2.3 Reator de Fluxo Laminar . . . . . . . . . . . . . . . . . . . . . . . 895.2.4 Função Distribuição do Tempo de Residência para Séries de PFR/CSTR

Ideais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 905.3 Modelagem de um Reator com a Função Distribuição do Tempo de Resi-

dência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 915.4 Modelos com Zero Parâmetros . . . . . . . . . . . . . . . . . . . . . . . . . 91

5.4.1 Modelo de Segregação . . . . . . . . . . . . . . . . . . . . . . . . . 925.4.2 Modelo de Máxima Antecipação da Mistura . . . . . . . . . . . . . 93

5.5 Modelos com um Parâmetro . . . . . . . . . . . . . . . . . . . . . . . . . . 945.5.1 Modelo de Tanques em Série . . . . . . . . . . . . . . . . . . . . . 945.5.2 Modelo de Dispersão . . . . . . . . . . . . . . . . . . . . . . . . . . 96

5.6 Modelos com 2 Parâmetros . . . . . . . . . . . . . . . . . . . . . . . . . . 995.6.1 CSTR Real Modelado com um Volume de Troca . . . . . . . . . . 99

A Cálculo 103A.1 Integrais Inde�nidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103A.2 Equações Diferenciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

A.2.1 Integração Numérica . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Lista de Figuras

1.1 Isoterma de Langmuir . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Isoterma BET, para c=100. . . . . . . . . . . . . . . . . . . . . . . . . . . 3

3.1 Per�l de concentração para modulo de Thiele igual a 0,2, 1,0, 2,0, 10 e 20.Geometria plana. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.2 Fator de efetividade em função do modulo de Thiele. Geometria plana. . . 293.3 Per�l de concentração para modulo de Thiele igual a 0,2, 1,0, 2,0, 10 e 20.

Geometria cilíndrica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.4 Fator de efetividade em função do modulo de Thiele. Geometria cilíndrica. 303.5 Per�l de concentração para modulo de Thiele igual a 0,2, 1,0, 2,0, 10 e 20.

Geometria esférica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.6 Fator de efetividade em função do modulo de Thiele. Geometria esférica. . 323.7 Per�l de concentração para modulo de Thiele igual a 0,2, 1,0, 2,0, 10 e 20,

com número de Biot igual a 1,0. Geometria plana. . . . . . . . . . . . . . 333.8 Fator de efetividade em função do modulo de Thiele para número de Biot

igual a 1,0. Geometria plana. . . . . . . . . . . . . . . . . . . . . . . . . . 333.9 Per�l de concentração para modulo de Thiele igual a 0,2, 1,0, 2,0, 10 e 20.

Geometria cilíndrica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.10 Fator de efetividade em função do modulo de Thiele. Geometria cilíndrica. 343.11 Per�l de concentração para modulo de Thiele igual a 0,2, 1,0, 2,0, 10 e 20,

para número de Biot igual a 1,0. Geometria esférica. . . . . . . . . . . . . 363.12 Fator de efetividade em função do modulo de Thiele, para número de Biot

igual a 1,0. Geometria esférica. . . . . . . . . . . . . . . . . . . . . . . . . 363.13 Partícula catalítica não isotérmica . . . . . . . . . . . . . . . . . . . . . . 52

4.1 Per�l de concentração para uma pressão parcial de alimentação igual a0,01 bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.2 Per�l de temperatura para uma pressão parcial de alimentação igual a 0,01bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.3 Per�l de concentração para uma pressão parcial de alimentação igual a0,015 bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.4 Per�l de temperatura para uma pressão parcial de alimentação igual a0,015 bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55


vi LISTA DE FIGURAS




4.9 Per�l de concentração para pressão parcial de reagente na alimentaçãoigual 0,005 bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.10 Per�l de Temperatura para pressão parcial de reagente na alimentaçãoigual a 0,005 bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.11 Per�l de concentração para pressão parcial de reagente na alimentaçãoigual a 0,008 bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63


4.13 Per�l de concentração para pressão parcial de reagente na alimentaçãoigual a 0,009 bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63


4.15 Per�l de conversão de B para uma pressão parcial de reagente na alimen-tação igual a 0,003 bar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4.16 Per�l de conversão de C para uma pressão parcial de reagente na alimen-tação igual a 0,003 bar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4.17 Per�l de temperatura para uma pressão parcial de reagente na alimentaçãoigual a 0,003 bar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4.18 Per�l de conversão de B para uma pressão parcial de reagente na alimen-tação igual a 0,0045 bar . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4.19 Per�l de conversão de C para uma pressão parcial de reagente na alimen-tação igual a 0,0045 bar . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74


4.21 Per�l de conversão de B para uma pressão parcial de reagente na alimen-tação igual a 0,005 bar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

4.22 Per�l de conversão de C para uma pressão parcial de reagente na alimen-tação igual a 0,005 bar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75


5.1 Função distribuição do tempo de residência para um reator PFR ideal(linha cheia) e para um reator PFR real (linha pontilhada) . . . . . . . . . 84

5.2 Função distribuição do tempo de residência para um reator CSTR ideal . 845.3 Função distribuição do tempo de residência para um reator de �uxo lami-

nar ideal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 845.4 Distribuição do tempo de residência para um conjunto de n reatores CSTR

ideais em série. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.5 Modelo de um reator com um volume de troca . . . . . . . . . . . . . . . 100

Lista de Tabelas

3.1 Área super�cial de um sólido não poroso. . . . . . . . . . . . . . . . . . . 253.2 Área super�cial de um sólido poroso. . . . . . . . . . . . . . . . . . . . . . 26

viii LISTA DE TABELAS

Capítulo 1

Introdução

1.1 Modelos de Adsorção1.1.1 Modelo de Langmuir

Langmuir foi o primeiro a desenvolver um modelo para a descrição da adsorção quí-mica em 19181. Este modelo se baseia nas seguintes hipóteses:

1. As entidades adsorvidas estão ligadas à superfície em pontos �xos denominadoscentros ativos, ou seja, as espécies adsorvidas não têm mobilidade.

2. Cada centro ativo só pode acomodar uma entidade adsorvida.

3. A adsorção é energeticamente idêntica em todos os sítios e independente da presençaou ausência de espécies adsorvidas em sua vizinhança.

Seja o equilíbrio de uma dada espécie química entre a fase gasosa e a superfície:

(velocidade de adsorção) = (velocidade de dessorção)

ka exp(−Ea

RT

)PΘv = kd exp

(−Ed

RT

)Θ (1.1)

Como:

Θv + Θ = 1 (1.2)

ka exp(−Ea

RT

)P (1−Θ) = kd exp

(−Ed

RT

)Θ (1.3)

Θ1−Θ

=ka

kdexp

(−Ea − Ed

RT

)P = KP (1.4)

∆Hads = Ea − Ed (1.5)1Langmuir, J. Am. Chem. Soc., 40:1361, (1918)

2 Introdução

0,0 0,2 0,4 0,6 0,8 1,0 0,0

0,2

0,4

0,6

0,8

1,0

Qua

ntid

ade

adso

rvid

a

Pressão parcial

Figura 1.1: Isoterma de Langmuir

Θ =KP

1 + KP(1.6)

A equação acima representa a Isoterma de Langmuir, que está ilustrada na Figura 1.1.Das hipóteses levantadas para o desenvolvimento da Isoterma de Langmuir a 3a é a queapresenta maiores erros. Freqüentemente se observa que o calor de adsorção diminui amedida em que o grau de cobertura aumenta, indicando uma heterogeneidade dos sítios.Os sítios mais fortes, que apresentam um calor de adsorção mais alto, são os primeirosa serem ocupados. Assim, os sítios restantes apresentarão um calor de adsorção maisbaixo.

1.1.2 Modelo BETEstendendo o modelo de Langmuir para adsorção em múltiplas camadas, Brunauer,

Emmet e Teller2 chegaram a equação denominada BET. Para seu desenvolvimento osautores se basearam nas seguintes hipóteses:

1. Em cada camada, a velocidade de adsorção é igual a velocidade de dessorção.

2. O calor de adsorção a partir da segunda camada é constante e igual ao calor latentede condensação.

3. Quando a pressão parcial do adsorbato se iguala a sua pressão de saturação, o vaporse condensa como um líquido ordinário e o número de camadas se torna in�nito.

2Brunauer, Emmet e Teller, J. Am. Chem. Soc., 60:309, (1938).

1.1 Modelos de Adsorção 3

0,0 0,2 0,4 0,6 0,8 1,0 0

2

4

6

8

10

Qua

ntid

ade

adso

rvid

a

Pressão parcial

Figura 1.2: Isoterma BET, para c=100.

(na

nam

)=

c (P/P0)(1− P/P0) (1− P/P0 + cP/P0)

(1.7)

(P/P0)na (1− P/P0)

=1

namc

+(

c− 1na

mc

)(P/P0) (1.8)

onde:P: pressão parcial do adsorbato.P0: pressão de saturação do adsorbato.na: quantidade adsorvida.na

m: quantidade adsorvida para formar uma monocamada.c: parâmetro dependente da diferença entre o calor de adsorção para formar a primeira

camada e o calor latente de condensação.

c = exp{

E1 −El

RT

}(1.9)

El: Entalpia de adsorção a partir da segunda camada, corresponde ao calor latentede condensação.

E1: Entalpia de adsorção para formação da primeira camada.O grá�co da quantidade adsorvida em função da pressão relativa, quando o parâmetro

c for alto (>2), nos fornece uma isoterma do tipo II, como representado na Figura 1.2.

4 Introdução

1.2 Cinética de Reações Catalíticas Heterogêneas1.2.1 Cinética de Langmuir-Hinshelwood-Hougen-Watson

Neste modelo, assume-se que uma das etapas do processo é a reação química entre asespécies adsorvidas. Assim, a velocidade global da reação depende da concentração destasespécies. Para relacionar as concentrações destas com a pressão parcial dos componentesde reação na fase gasosa utiliza-se o modelo de Langmuir.

Desta forma, três etapas são necessárias para o desenvolvimento do modelo cinético:

• Adsorção dos reagentes.

• Reação super�cial.

• Dessorção dos produtos.

No modelo de Langmuir-Hinshelwood é assumido que umas destas etapas é a deter-minante do processo.

Etapa determinante: É aquela que consome a maior parte da força motriz (poten-cial químico) disponível. É a etapa intrinsecamente mais lenta. Na reação em estadoestacionário todas três etapas terão a mesma velocidade.

Adsorção QuímicaAdsorção e Dessorção Simples Cada espécie ocupa um centro ativo, ou seja, aadsorção ocorre sem dissociação ou associação.

A + v kaA−−⇀↽−−kdA

Av

rads = kaAPAΘv (1.10)

rdes = kdAΘA (1.11)O somatório dos sítios livres e ocupados é igual a unidade:

ΘA + Θv = 1 (1.12)

Θv = 1−ΘA (1.13)No equilíbrio rads = rdes, então temos:

kaAPA (1−ΘA) = kdAΘA (1.14)

ΘA =KAPA

1 + KAPA(1.15)

1.2 Cinética de Reações Catalíticas Heterogêneas 5

Adsorção de Vários Componentes Cada espécie ocupa um centro ativo, assimtemos:

Θv = 1−N∑

i=1

Θi (1.16)

Para cada uma das espécies, temos no equilíbrio:

kaiPi

1−

N∑

j=1

Θj

= kdiΘi (1.17)

Θi = KiPi

1−

N∑

j=1

Θj

(1.18)

Fazendo-se o somatório:

N∑

i=1

Θi =

[N∑

i=1

KiPi

]1−

N∑

j=1

Θj

(1.19)

1−N∑

i=1

Θi = 1−[

N∑

i=1

KiPi

]1−

N∑

j=1

Θj

(1.20)

1−N∑

i=1

Θi +

[N∑

i=1

KiPi

]1−

N∑

j=1

Θj

(1.21)

1−N∑

j=1

Θj +

[N∑

i=1

KjPj

]1−

N∑

j=1

Θj

(1.22)

1−N∑

j=1

Θj =1

1 +∑N

j=1 KjPj

(1.23)

Θi =KiPi

1 +∑N

j=1 KjPj

(1.24)

Adsorção Dissociativa Cada espécie na fase gasosa ocupa dois centros ativos.

A2 + 2 v kaA−−⇀↽−−kdA

2Av

No equilíbrio rads = rdes, então temos:

kaAPA (1−ΘA)2 = kdAΘ2A (1.25)

KAPA (1−ΘA)2 = Θ2A (1.26)

6 Introdução

ΘA =√

KAPA

1 +√

KAPA(1.27)

Adsorção Associativa Cada duas espécies na fase gasosa ocupa um centro ativo.

2A + v kaA−−⇀↽−−kdA

A2v

No equilíbrio rads = rdes, então temos:

kaAP 2A (1−ΘA) = kdAΘA (1.28)

KAP 2A (1−ΘA) = ΘA (1.29)

ΘA =KAP 2

A

1 + KAP 2A

(1.30)

Reação Irreversível ElementarReação global:

A(g) −−→ B(g)

Mecanismo:

A(g) + v kaA−−⇀↽−−kdA

Av 1

Av k−→ Bv 2

Bv kdB−−−⇀↽−−−kaB

B(g) 3

A velocidade de adsorção de A é dada por:

r1 = kaA(PAΘv −ΘA/KA) (1.31)

A velocidade da reação química super�cial é dada por:

r2 = kΘA (1.32)

A velocidade dessorção de B é dada por:

−r3 = kaB(PB −ΘB/KB) (1.33)

No estado estacionário as velocidades das etapas (1), (2) e (3) serão iguais. A velo-cidade global da reação é dada por:

(−rA) = r2 = kΘA (1.34)Para as etapas (1) e (3) utilizamos suas relações de equilíbrio:


ΘA = KAPAΘv (1.35)ΘB = KBPBΘv (1.36)

O somatório das frações de sítios ocupados por A e B e de sítios livres é igual aunidade:

ΘA + ΘB + Θv = 1 (1.37)

Substituindo 1.35 e 1.36 em 1.37

KAPAΘv + KBPBΘv + Θv = 1 (1.38)

Explicitando Θv:

Θv =1

1 + KAPA + KBPB(1.39)

(−rA) = kΘA = kKAPAΘv =kKAPA

1 + KAPA + KBPB(1.40)

Reação Irreversível Elementar com AssociaçãoReação global:

2A(g) −−→ B(g)

Mecanismo:


Av 1

Av k−→ Bv 2


B(g) 3


r1 = kaA(P 2AΘv −ΘA/KA) (1.41)

A velocidade da reação química super�cial é:

r2 = kΘA (1.42)

A velocidade dessorção de B é:

−r3 = kaB(PBΘv −ΘB/KB) (1.43)


8 Introdução

(−rA) = r2 = kΘA (1.44)Para as etapas (1) e (3) utilizamos suas relações de equilíbrio:

ΘA = KAP 2AΘv (1.45)

ΘB = KBPBΘv (1.46)


ΘA + ΘB + Θv = 1 (1.47)Substituindo 1.45 e 1.46 em 1.47

KAP 2AΘv + KBPBΘv + Θv = 1 (1.48)

Explicitando Θv:

Θv =1

1 + KAP 2A + KBPB

(1.49)

(−rA) = kΘA = kKAPAΘv =kKAP 2

A

1 + KAP 2A + KBPB

(1.50)

Reação Irreversível Elementar com DissociaçãoReação global:

A(g) −−→ B(g)

Mecanismo:

A2(g) + 2 v kaA−−⇀↽−−kdA

2Av 1

Av k−→ Bv 2


B(g) 3


r1 = kaA(PAΘv −Θ2A/KA) (1.51)


r2 = kΘA (1.52)A velocidade dessorção de B é:




(−rA) = r2 = kΘA (1.54)

Para as etapas (1) e (3) utilizamos suas relações de equilíbrio:

ΘA =√

KAPAΘv (1.55)ΘB = KBPBΘv (1.56)


ΘA + ΘB + Θv = 1 (1.57)

Substituindo 1.55 e 1.56 em 1.57√


Explicitando Θv:

Θv =1

1 +√

KAPA + KBPB(1.59)

(−rA) = kΘA = k√

KAPAΘv =k√

KAPA

1 +√

KAPA + KBPB(1.60)

Reação Reversível ElementarReação global:

A(g)−−⇀↽−− B(g)

Mecanismo:


Av 1

Avk+1−−⇀↽−−k−1

Rv 2


B(g) 3




r2 = k(ΘA −ΘB/Krs) (1.62)

10 Introdução

A velocidade dessorção de B é:


No estado estacionário as velocidades das etapas (1), (2) e (3) serão iguais. A veloci-dade global da reação será dada pela velocidade da etapa determinante do processo. Aetapa determinante será aquela que irá consumir a maior parte da força motriz disponívelpara o processo (potencial químico).

A reação química super�cial é a etapa determinante: A velocidade global doprocesso é dada por:

(−rA) = r2 = k(ΘA −ΘB/Krs) (1.64)


ΘA = KAPAΘv (1.65)ΘB = KBPBΘv (1.66)


ΘA + ΘB + Θv = 1 (1.67)



Explicitando Θv:

Θv =1

1 + KAPA + KBPB(1.69)

(−rA) = k(ΘA −ΘB/Krs) = k

(KAPAΘv − KBPBΘv

Krs

)(1.70)

(−rA) =k (KAPA −KBPB/Krs)

1 + KAPA + KBPB(1.71)

A constante de equilíbrio global da reação é dada por:

K =(

PB

PA

)

eq

=(

ΘBKAΘv

ΘAKBΘv

)

eq

(1.72)

K =KA

KB

(ΘB

ΘA

)

eq

=KAKrs

KB(1.73)

Então:


(−rA) =k (KAPA −KAPB/K)1 + KAPA + KBPB

(1.74)

(−rA) =kKA (PA − PB/K)1 + KAPA + KBPB

(1.75)

Adsorção de A é a etapa determinante: A velocidade global do processo é dadapor:

(−rA) = r1 = kaA(PAΘv −ΘA/KA) (1.76)


ΘA = ΘB/Krs (1.77)ΘB = KBPBΘv (1.78)


ΘA + ΘB + Θv = 1 (1.79)


KBPBΘv

Krs+ KBPBΘv + Θv = 1 (1.80)

Explicitando Θv:

Θv =1

(1 + KBPB/Krs + KBPB)(1.81)

(−rA) = kaA (PAΘv −ΘA/KA) = kaA

(PA − KBPB

KAKrs

)Θv (1.82)

(−rA) =kaA

(PA − KBPB

KAKrs

)

(1 + KBPB/Krs + KBPB)(1.83)

A constante de equilíbrio global da reação é dada pela Equação 1.73. Então:

(−rA) =kaA (PA − PB/K)

(1 + KAPB/K + KBPB)(1.84)

(−rA) =kaA (PA − PB/K)

[1 + (KA/K + KB) PB](1.85)

12 Introdução

A dessorção de B é a etapa determinante: A velocidade global do processo é dadapor:

(−rA) = r3 = −kaB(PBΘv −ΘB/KB) (1.86)


ΘB = KrsΘA (1.87)ΘA = KAPAΘv (1.88)


ΘA + ΘB + Θv = 1 (1.89)


KAPAΘv + KrsKAPAΘv + Θv = 1 (1.90)

Explicitando Θv:

Θv =1

1 + KAPA + KrsKAPA(1.91)

(−rA) = kaB(ΘB/KB − PBΘv) = kaB (KrsKAPAΘv − PBΘv) (1.92)

(−rA) =kaB (KrsKAPA − PB)1 + KAPA + KrsKAPA

(1.93)


(−rA) =kaB(KKBPA − PB)

1 + KAPA + KKBPA(1.94)

(−rA) =kaBK (KBPA − PB/K)[1 + (KA + KKB) PA]

(1.95)

Reação Reversível Bimolecular

Reação global:

A(g) + B(g)−−⇀↽−− R(g)

Mecanismo:



Av 1

B(g) + v kaB−−−⇀↽−−−kdB

Av 2

Av + Bvk+1−−⇀↽−−k−1

Rv + v 3

Rv kdR−−⇀↽−−kaR

R(g) + v 4



A velocidade de adsorção de B é dada por:

r2 = kaB(PBΘv −ΘB/KB) (1.97)


r3 = k(ΘAΘB −ΘRΘv/Krs) (1.98)

A velocidade dessorção de R é:

−r4 = kaR(PRΘv −ΘR/KR) (1.99)

No estado estacionário as velocidades das etapas (1), (2), (3) e (4) serão iguais. Avelocidade global da reação será dada pela velocidade da etapa determinante do pro-cesso. A etapa determinante será aquela que irá consumir a maior parte da força motrizdisponível para o processo (potencial químico).


(−rA) = r2 = k (ΘAΘB −ΘRΘv/Krs) (1.100)

Para as etapas (1), (2) e (4) utilizamos suas relações de equilíbrio:

ΘA = KAPAΘv (1.101)ΘB = KBPBΘv (1.102)ΘR = KRPRΘv (1.103)

O somatório das frações de sítios ocupados por A, B e R e de sítios livres é igual aunidade:

ΘA + ΘB + ΘR + Θv = 1 (1.104)

Substituindo 1.101, 1.102 e 1.103 em 1.104

14 Introdução

KAPAΘv + KBPBΘv + KRPRΘv + Θv = 1 (1.105)

Explicitando Θv:

Θv =1

1 + KAPA + KBPB + KRPR(1.106)

(−rA) = k (ΘAΘB −ΘR/Krs) = k

(KAPAKBPBΘ2

v −KRPRΘv

Krs

)(1.107)

(−rA) =k (KAKBPAPB −KRPR/Krs)(1 + KAPA + KBPB + KRPR)2

(1.108)


K =(

PR

PAPB

)

eq

=(

ΘRKAKBΘ2v

ΘAΘbKRΘv

)

eq

(1.109)

K =KAKB

KR

(ΘRΘv

ΘAΘB

)

eq

=KAKBKrs

KR(1.110)

Então:

(−rA) =k (KAKBPAPB −KAKBPR/K)(1 + KAPA + KBPB + KRPR)2

(1.111)

(−rA) =kKAKB (PAPB − PR/K)

(1 + KAPA + KBPB + KRPR)2(1.112)

Adsorção de A (ou de B) é a etapa determinante: A velocidade global do processoé dada por:

(−rA) = r1 = kaA(PAΘv −ΘA/KA) (1.113)

Para as etapas (2), (3) e (4) utilizamos suas relações de equilíbrio:

ΘA =ΘRΘv

KrsΘB(1.114)

ΘB = KBPBΘv (1.115)ΘR = KRPRΘv (1.116)


ΘA + ΘB + ΘR + Θv = 1 (1.117)

Substituindo 1.114, 1.115 e 1.116 em 1.117


ΘRΘv

KrsΘB+ KBPBΘv + KRPRΘv + Θv = 1 (1.118)

KRPRΘv

KrsKBPB+ KBPBΘv + KRPRΘv + Θv = 1 (1.119)

Explicitando Θv:

Θv =1(

1 + KRPRKrsKBPB

+ KBPB + KRPR

) (1.120)

(−rA) = kaA(PAΘv −ΘA/KA) =kaA

(PA − KRPR

KrsKBPBKA

)(1 + KRPR

KrsKBPB+ KBPB + KRPR

) (1.121)


(−rA) =kaA

(PA − PR

K

)(1 + KAPR

KPB+ KBPB + KRPR

) (1.122)

A dessorção de R é a etapa determinante: A velocidade global do processo é dadapor:

(−rA) = r4 = −kaR(PRΘv −ΘR/KR) (1.123)Para as etapas (1), (2) e (3) utilizamos suas relações de equilíbrio:

ΘR = KrsΘAΘB (1.124)ΘA = KAPAΘv (1.125)ΘB = KBPBΘv (1.126)


ΘA + ΘB + ΘR + Θv = 1 (1.127)Substituindo 1.125, 1.126 e 1.124 em 1.127

KAPAΘv + KAPAΘv + KrsKAPAKBPBΘv + Θv = 1 (1.128)Explicitando Θv:

Θv =1

1 + KAPA + KrsKAKBPAPB(1.129)

(−rA) = kaR(ΘR/KR − PRΘv) = kaR(KrsKAKBPAPBΘv − PRΘv) (1.130)

16 Introdução

(−rA) =kaB(KrsKAKBPAPB − PR)

1 + KAPA + KrsKAKBPAPB(1.131)


(−rA) =kaB(KKRPAPB − PR)

1 + KAPA + KKRPAPB(1.132)

(−rA) =kaBKKR

(PAPB − PR

KKR

)


1.2.2 Cinética de Eley-RidealNeste modelo assume-se que a reação ocorre entre uma espécie adsorvida e outra na

fase gasosa. Para correlacionar a concentração da espécie adsorvida com sua pressãoparcial na fase gasosa utiliza-se o modelo de Langmuir.

Reação Reversível BimolecularReação global:

A(g) + B(g)−−⇀↽−− R(g)

Mecanismo:


Av 1

Av + B(g)

k+1−−⇀↽−−k−1

Rv 2

Rv kdR−−⇀↽−−kaR

R(g) 3




r2 = k(ΘAPB −ΘR/Krs) (1.135)

A velocidade dessorção de R é:

−r3 = kaR(PRΘv −ΘR/KR) (1.136)

No estado estacionário as velocidades das etapas (1), (2) e (3) serão iguais. A veloci-dade global da reação será dada pela velocidade da etapa determinante do processo. Aetapa determinante será aquela que irá consumir a maior parte da força motriz disponívelpara o processo (potencial químico).



(−rA) = r2 = k(ΘAPB −ΘR/Krs) (1.137)


ΘA = KAPAΘv (1.138)ΘR = KRPRΘv (1.139)

O somatório das frações de sítios ocupados por A e R e de sítios livres é igual aunidade:

ΘA + ΘR + Θv = 1 (1.140)


KAPAΘv + KRPRΘv + Θv = 1 (1.141)

Explicitando Θv:

Θv =1

1 + KAPA + KRPR(1.142)

(−rA) = k(ΘAPB −ΘR/Krs) = k

(KAPAPBΘv − KRPRΘv

Krs

)(1.143)

(−rA) =k (KAPAPB −KRPR/Krs)

1 + KAPA + KRPR(1.144)


K =(

PR

PAPB

)

eq

=(

ΘRKAΘv

ΘAKRΘvPB

)

eq

(1.145)

K =KA

KB

(ΘR

ΘAPB

)

eq

=KAKrs

KR(1.146)

Então:

(−rA) =k (KAPAPB −KAPR/K)

1 + KAPA + KRPR=

kKA (PAPB − PR/K)1 + KAPA + KRPR

(1.147)

18 Introdução

Adsorção de A é a etapa determinante: A velocidade global do processo é dadapor:

(−rA) = r1 = kaA(PAΘv −ΘA/KA) (1.148)Para as etapas (2) e (3) utilizamos suas relações de equilíbrio:

ΘA =ΘR

KrsPB(1.149)

ΘR = KRPRΘv (1.150)


ΘA + ΘR + Θv = 1 (1.151)Substituindo 1.149 e 1.150 em 1.151

KRPRΘv

KrsPB+ KRPRΘv + Θv = 1 (1.152)

Explicitando Θv:

Θv =1(

1 + KRPRKrsPB

+ KRPR

) (1.153)

(−rA) = kaA (PAΘv −ΘA/KA) = kaA

(PA − KRPR

KAKrsPB

)Θv (1.154)

(−rA) =kaA

(PA − KRPR

KAKrsPB

)(1 + KRPR

KrsPB+ KRPR

) (1.155)


(−rA) =kaA

(PA − PR

KPB

)(1 + KAPR

KPB+ KRPR

) (1.156)

A dessorção de B é a etapa determinante: A velocidade global do processo é dadapor:

(−rA) = r3 = −kaR(PRΘv −ΘR/KR) (1.157)Para as etapas (1) e (2) utilizamos suas relações de equilíbrio:

ΘR = KrsΘAPB (1.158)ΘA = KAPAΘv (1.159)


O somatório das frações de sítios ocupados por A e R e de sítios livres é igual aunidade:

ΘA + ΘR + Θv = 1 (1.160)Substituindo 1.158 e 1.159 em 1.160

KAPAΘv + KrsKAPAPBΘv + Θv = 1 (1.161)Explicitando Θv:

Θv =1

1 + KAPA + KrsKAPAPB(1.162)

(−rA) = r3 = kaR(ΘR/KR − PRΘv) = kaR(KrsKAPAPBΘv/KR − PRΘv) (1.163)

(−rA) =kaR(KrsKAPAPB/KR − PR)1 + KAPA + KrsKAPAPB

(1.164)


(−rA) =kaR(KPAPB − PR)


(−rA) =kaRK(PAPB − PR/K)


20 Introdução

Capítulo 2

Resistência à Transferência deMassa Externa

2.1 IntroduçãoPara obter as equações para a velocidade de reação do capítulo anterior os efeitos

de transferência de calor e de massa foram desprezados, ou seja, assumiu-se que estaetapa era su�cientemente rápida e, assim, a concentração dos componentes de reação ea temperatura próximo aos centros ativos do catalisador eram iguais àquelas da correnteprincipal do �uido. Dois efeitos de transporte devem ser levados em conta: o transportede calor e massa para a superfície da partícula e o transporte de calor e massa no interiorda partícula. Neste capítulo estudaremos os efeitos de transferência de massa e de calorexternos.

2.2 Resistência Externa à Transferência de MassaSeja uma partícula catalítica e não porosa em uma corrente de um determinado rea-

gente A. A velocidade de escoamento do �uido nas vizinhanças da partícula irá variar coma posição ao redor desta, sendo que na superfície esta velocidade será nula. Assim, desen-volverá um gradiente acentuado de velocidade próximo a partícula. Para distâncias dapartícula relativamente curtas a velocidade será praticamente constante e o escoamentodeve ser turbulento. A região onde estão concentrados os gradientes de velocidade deescoamento do �uido denomina-se camada limite hidrodinâmica. Devido ao escoamentoturbulento na corrente principal do �uido, a concentração de reagente nesta região tendea ser homogênea. Próximo à superfície, devido à reação química que está consumindoeste reagente, haverá também um gradiente pronunciado de concentração deste. Esta re-gião onde estão concentrados os gradientes de concentração denomina-se camada limitede concentração.

Uma forma simpli�cada de se modelar a transferência de massa externa é assumirque esta região próxima a superfície da partícula (camada limite) está estagnada e temespessura δ. Assim, para sistemas diluídos ou que apresentem contra-difusão equimolar,temos:

22 Resistência à Transferência de Massa Externa

NAz =DAB

δ(CA0 − CAs) (2.1)

onde CAs é a concentração de A na superfície da partícula de catalisador e CA0 é aconcentração de A na corrente principal do �uido.

A espessura da camada limite (δ) geralmente não é conhecida, mas varia em funçãoda �uidodinâmica do meio, ou seja, depende da velocidade de escoamento, da viscosidadee da densidade do �uido, etc. Esta equação pode ser escrita de forma mais genérica como:

NAz = kc (CA0 − CAs) (2.2)

onde kc é o coe�ciente de transferência de massa, que pode ser calculado por correla-ções. A forma usual é expressar estas correlação como uma função de números adimen-sionais. Assim, por exemplo, Frössling propos a seguinte correlação para a transferênciade massa para uma esfera:

Sh = 2 + 0, 6Re1/2Sc1/3 (2.3)

onde

Sh =kcdp

DAB

Sc =µ

ρDAB

Re =v∞ρdp

µ

(2.4)

2.2.1 Reação de Primeira OrdemSeja a reação de isomerização:

Ak // B

se passando na superfície de uma esfera sólida de um catalisador. A velocidade dereação super�cial por unidade de área de catalisador será dada por:

−r′′As = ksCAs (2.5)

O �uxo de massa para a superfície será dado por:

NAr = kc (CA − CAs) (2.6)

No estado estacionário teremos a igualdade entre as velocidades de reação super�ciale de transferência de massa:

NAr = −r′′As (2.7)

2.2 Resistência Externa à Transferência de Massa 23

kc (CA − CAs) = ksCAs (2.8)

Explicitando CAs:

CAs =kcCA

ks + kc(2.9)

E a velocidade de reação global �ca:

NAr = −r′′As = ksCAs =

kcksCA

ks + kc=

CA1ks

+ 1kc

= k′CA (2.10)

1k′

=1kc

+1ks

(2.11)

Assim, a velocidade de reação global depende da soma de duas resistências: a detransferência de massa externa e de reação reação química super�cial. A maior resistênciairá controlar o processo global. Assim, quando a velocidade de reação química super�cialfor rápida em comparação com a velocidade de transferência de massa, ou seja, ks >> kc,a velocidade global da reação será:

−r′′As =

CA1ks

+ 1kc

' kcCA (2.12)

CAs =kcCA

kc + ks' 0 (2.13)

Neste caso, a resistência à transferência de massa é maior que a resistência a reaçãoquímica super�cial. Dizemos que este sistema é controlado pela transferência de massaexterna.

Quando a velocidade de reação química super�cial for lenta em comparação com avelocidade de transferência de massa, ou seja, kc >> ks

−r′′As =

CA1ks

+ 1kc

' ksCA (2.14)

CAs =kcCA

kc + ks' CA (2.15)

Neste caso, a resistência à transferência de massa é menor que a resistência a reaçãoquímica super�cial. Dizemos que este sistema é controlado pela reação super�cial.

2.2.2 Reação de Segunda OrdemSeja a reação de isomerização:

Ak // B

se passando na superfície de uma esfera sólida de um catalisador. A velocidade dereação super�cial por unidade de área de catalisador será dada por:

24 Resistência à Transferência de Massa Externa

−r′′As = ksC

2As (2.16)

O �uxo de massa para a superfície será dado por:

NAr = kc (CA − CAs) (2.17)No estado estacionário teremos a igualdade entre as velocidades de reação super�cial

e de transferência de massa:

NAr = −r′′As (2.18)

kc (CA − CAs) = ksC2As (2.19)

ksC2As + kcCAs − kcCA = 0 (2.20)

Explicitando CAs:

CAs =−kc ±

√k2

c + 4kskcCA

2ks(2.21)

Capítulo 3

Resistência à Transferência deMassa Interna

3.1 Introdução

Uma reação catalisada por um sólido ocorre na interface sólido/gás. Assim, a velo-cidade da reação depende da área super�cial do catalisador. A forma usual de se obteruma elevada área super�cial é com o uso de materiais porosos. Para ser efetivo em ca-tálise heterogênea, um sólido deve apresentar áreas super�ciais da ordem de 100 m2/gou superiores. Para conseguir esta área utilizando sólidos não porosos deveríamos traba-lhar com partículas extremadamente pequenas, que seriam muito dí�cies de lidar. Porexemplo, na Tabela 3.1 abaixo apresentamos a área super�cial obtida com um sólidonão poroso de densidade igual a 1 g/cm3 formado por partículas esféricas, enquanto naTabela 3.2 apresentamos a área obtida em um sólido preenchido por poros cilíndricos.Comparando-se estas tabelas, pode-se notar que a presença de porosidade é muito maisefetiva para gerar sólidos apresentando áreas super�ciais elevadas.

Devido à presença destes poros, teremos a difusão dos reagentes até o interior dosólido. Dependendo da velocidade de difusão em relação à velocidade de reação quí-mica, podemos ter um gradiente de concentração mais ou menos acentuado ao longo dapartícula de catalisador.

Tabela 3.1: Área super�cial de um sólido não poroso.Diâmetro da partícula (cm) Área super�cial (m2g−1)

1× 100 6× 10−4

1× 10−3 6× 10−1

1× 10−4 6× 100

1× 10−5 6× 101

1× 10−6 6× 102

26 Resistência à Transferência de Massa Interna

Tabela 3.2: Área super�cial de um sólido poroso.Raio do poro (cm) Área super�cial (m2g−1)

2× 100 1× 10−4

2× 10−2 1× 10−2

2× 10−4 1× 100

2× 10−5 1× 101

2× 10−7 1× 103

3.2 Sistemas IsotérmicosA difusão e reação química em um sólido poroso pode ser descrita pela equação da

continuidade molar:

∂CA

∂t+−→∇ · −→NA = R

′′′A (3.1)

Para um sistema em estado estacionário:−→∇ · −→NA = R

′′′A (3.2)

Para um sistema diluído ou que apresente contra-difusão equimolar:−→NA = −DAe

−→∇ · CA (3.3)Onde DAe é a difusividade efetiva, que leva em conta que podem ocorrer outros

mecanismos difusivos no interior de um poro, além da difusividade molecular, como adifusão de Knudsen. Além disto, a difusividade efetiva também leva em conta que nemtodo o volume do sólido está disponível para a difusão, somente sua fração de espaçosvazios, e que os poros do sólido não são retos e podem apresentar seções variáveis.

Para sistemas unidirecionais:

NAr = −DAedCA

dr(3.4)

E, então, teremos, para uma cinética de ordem m:

1r′s

ddr′

[r′sNA,r

]= −kvC

mA (3.5)

onde s depende da geometria, sendo igual a 0 para geometria plana, 1 para cilíndricae 2 para esférica.

DAe

r′sd

dr′

[r′s dCA

dr′

]= kvC

mA (3.6)

Condições de contorno:

r′= R CA = CA0

r′= 0

dCA

dr= 0

(3.7)

3.2 Sistemas Isotérmicos 27

Adimensionalizando:

y =CA

CA0

r =r′

R

(3.8)

1rs

ddr

(rs dy

dr

)= φ2ym (3.9)


dy

dr

∣∣∣∣r=0

= 0

dy

dr

∣∣∣∣r=1

= Bim(1− y)

onde:

φ2 =kvC

m−1A0 R2

DAe

Bim =kcR

DAe

3.2.1 Resolução Analítica da Difusão em Ausência de Resistência Ex-terna

A equação diferencial apresentada no ítem anterior somente apresentará solução ana-lítica para ordem de reação igual a 1.

1rs

ddr

(rs dy

dr

)= φ2y (3.10)


dy

dr

∣∣∣∣r=0

= 0

y|r=1 = 1

(3.11)

Geometria PlanaNeste caso:

d2y

dr2− φ2y = 0 (3.12)

Que tem como solução:


y = A cosh(φr) + Bsenh (φr) (3.13)

Aplicando as condições de contorno:

dy

dr= Aφsenh (φr) + Bφ cosh(φr) (3.14)

dy

dr

∣∣∣∣r=0

= Aφsenh (0) + Bφ cosh(0) = B = 0 (3.15)

Assim:

y = A cosh(φr) (3.16)

Como:

y|r=1 = A cosh(φ) = 1 (3.17)

Então:

A =1

cosh(φ)(3.18)

y =cosh(φr)cosh(φ)

(3.19)

O fator de efetividade será:

η =

∫ +R−R kvCAd(Ar

′)

∫ +R−R kvCA0d(Ar′)

=2kvA

∫ +R0 CAdr

′

2KvACA0

∫ +R0 dr′

(3.20)

Adimensionalizando:

η =2kvACA0R

∫ 10 ydr

2KvACA0

∫ 10 dr

(3.21)

η =∫ 1

0ydr =

∫ 1

0

cosh(φr)cosh(φ)

dr (3.22)

η =senh (φ)φ cosh(φ)

=tanh(φ)

φ(3.23)

Nas Figuras 3.1 e 3.2 são mostrados os per�s de concentração em função da posiçãona partícula para vários módulos de Thiele diferentes e o fator de efetividade em funçãodo modulo de Thiele.


0,0 0,2 0,4 0,6 0,8 1,0 0,0

0,2

0,4

0,6

0,8

1,0

Con

cent

raçã

o

Raio

Figura 3.1: Per�l de concentração paramodulo de Thiele igual a 0,2, 1,0, 2,0, 10e 20. Geometria plana.

1 10 100

0,0

0,2

0,4

0,6

0,8

1,0

Fato

r de

efe

tivid

ade

Modulo de Thiele

Figura 3.2: Fator de efetividade em fun-ção do modulo de Thiele. Geometriaplana.

Geometria CilíndricaNeste caso:

1r

ddr

(rdy

dr

)− φ2y = 0 (3.24)

r2 d2y

dr2+ r

dy

dr− φ2yr2 = 0 (3.25)


y = AI0(φr) + BK0(φr) (3.26)


y|r=0 = �nito (3.27)

Como K0(0) = −∞, B = 0 e:

y = AI0(φr) (3.28)

y|r=1 = AI0(φ) = 1 (3.29)

A =1

I0(φ)(3.30)

y =I0(φr)I0(φ)

(3.31)



0,0 0,2 0,4 0,6 0,8 1,0 0,0

0,2

0,4

0,6

0,8

1,0 C

once

ntra

ção

Raio

Figura 3.3: Per�l de concentração paramodulo de Thiele igual a 0,2, 1,0, 2,0, 10e 20. Geometria cilíndrica.

1 10 100 0,0

0,2

0,4

0,6

0,8

1,0

Fato

r de

efe

tivid

ade

Modulo de Thiele

Figura 3.4: Fator de efetividade em fun-ção do modulo de Thiele. Geometria ci-líndrica.

η =

∫ R0 kvCAd(πr

′2L∫ R0 kvCA0d(πr′2L

=2πkvL

∫ R0 CAr

′dr

′

2πkvLCA0

∫ R0 r′dr′

(3.32)

Adimensionalizando:

η =2πkvLCA0R

2∫ 10 yrdr

2πkvLCA0R2∫ 10 rdr

=2πkvLCA0R

2∫ 10 yrdr

2πkvLCA0R2 12

(3.33)

η = 2∫ 1

0yrdr = 2

∫ 1

0

I0(φr)I0(φ)

rdr (3.34)

η =2I1(φ)φI0(φ)

(3.35)


Geometria EsféricaNeste caso:

1r2

ddr

(r2 dy

dr

)− φ2y = 0 (3.36)

d2y

dr2+

2r

dy

dr− φ2y = 0 (3.37)


y =A cosh(φr) + Bsenh (φr)

r(3.38)



y|r=0 = �nito (3.39)

Como cosh(0)0

= ∞, A = 0 e:

y =Bsenh (φr)

r(3.40)

y|r=1 = Bsenh (φ) (3.41)

B =1

senh (φ)(3.42)

y =senh (φr)rsenh (φ)

(3.43)


η =

∫ R0 kvCAd

(43πr

′3)

∫ R0 kvCA0d

(43πr′3

) =4πkv

∫ R0 CAr

′2dr′

4 πkvCA0

∫ R0 r′2dr′

(3.44)

Adimensionalizando:

η =4πkvCA0R

3∫ 10 yr2dr

4πkvCA0R3∫ 10 r2dr

=4πkvCA0R

3∫ 10 yr2dr

4πkvCA0R3 13

(3.45)

η = 3∫ 1

0yr2dr = 3

∫ 1

0

senh (φr)senh (φ)

rdr (3.46)

η =3 [φ coth (φ)− 1]

φ2(3.47)


3.2.2 Resolução Analítica Difusão com Resistência Externa

1rs

ddr

(rs dy

dr

)= φ2y (3.48)


dy

dr

∣∣∣∣r=0

= 0

dy

dr

∣∣∣∣r=1

= Bim (1− y)

(3.49)


0,0 0,2 0,4 0,6 0,8 1,0 0,0

0,2

0,4

0,6

0,8

1,0 C

once

ntra

ção

Raio

Figura 3.5: Per�l de concentração paramodulo de Thiele igual a 0,2, 1,0, 2,0, 10e 20. Geometria esférica.

1 10 100

0,0

0,2

0,4

0,6

0,8

1,0

Fato

r de

efe

tivid

ade

Modulo de Thiele

Figura 3.6: Fator de efetividade em fun-ção do modulo de Thiele. Geometria es-férica.

Geometria PlanaNeste caso:

d2y

dr2− φ2y = 0 (3.50)


y = A cosh(φr) + Bsenh (φr) (3.51)Aplicando as condições de contorno:

dy

dr= Aφsenh (φr) + Bφ cosh(φr) (3.52)

dy

dr

∣∣∣∣r=0

= Aφsenh (0) + Bφ cosh(0) = B = 0 (3.53)

Assim:

y = A cosh(φr) (3.54)Como:

dy

dr

∣∣∣∣r=1

= Aφsenh (φ) = Bim (1−A cosh(φr)) (3.55)

Então:

A =Bim

φsenh (φ) + Bim cosh(φ)(3.56)

y =Bim cosh(φr)

φsenh (φ) + Bim cosh(φ)(3.57)


0,0 0,2 0,4 0,6 0,8 1,0 0,0

0,2

0,4

0,6

0,8

1,0

Con

cent

raçã

o

Raio

Figura 3.7: Per�l de concentração paramodulo de Thiele igual a 0,2, 1,0, 2,0,10 e 20, com número de Biot igual a 1,0.Geometria plana.

1E-3 0,01 0,1 1 10 100 0,0

0,2

0,4

0,6

0,8

1,0

Fato

r de

efe

tivid

ade

Modulo de Thiele

Figura 3.8: Fator de efetividade em fun-ção do modulo de Thiele para número deBiot igual a 1,0. Geometria plana.

O fator de efetividade global será:

ηg =∫ 1

0ydr =

∫ 1

0

Bim cosh(φr)φsenh (φ) + Bim cosh(φ)

dr (3.58)

ηg =Bimsenh (φ)

φ [φsenh (φ) + Bim cosh(φ)](3.59)

1ηg

=φ

tanh(φ)+

φ2

Bim=

1η

+φ2

Bim(3.60)


Geometria CilíndricaNeste caso:

1r

ddr

(rdy

dr

)− φ2y = 0 (3.61)

r2 d2y

dr2+ r

dy

dr− φ2yr2 = 0 (3.62)


y = AI0(φr) + BK0(φr) (3.63)Aplicando as condições de contorno:

y|r=0 = �nito (3.64)


0,0 0,2 0,4 0,6 0,8 1,0 0,0

0,2

0,4

0,6

0,8

1,0 C

once

ntra

ção

Raio

Figura 3.9: Per�l de concentração paramodulo de Thiele igual a 0,2, 1,0, 2,0, 10e 20. Geometria cilíndrica.

0,01 0,1 1 10 100 0,0

0,2

0,4

0,6

0,8

1,0

Fato

r de

efe

tivid

ade

Modulo de Thiele

Figura 3.10: Fator de efetividade em fun-ção do modulo de Thiele. Geometria ci-líndrica.

Como K0(0) = −∞, B = 0 e:

y = AI0(φr) (3.65)

dy

dr

∣∣∣∣r=1

= AφI1(φ) = Bim [1−AI0(φ)] (3.66)

A =Bim

φI1(φ) + BimI0(φ)(3.67)

y =BimI0(φr)

φI1(φ) + BimI0(φ)(3.68)


ηg = 2∫ 1

0yrdr =

∫ 1

0

BimI0(φr)φI1(φ) + BimI0(φ)

rdr (3.69)

ηg =2BimI1(φ)

φ2I1(φ) + BimφI0(φ)(3.70)

1ηg

=φ2

2Bim+

φI0(φ)2I1(φ)

=1η

+φ2

2Bim(3.71)



Geometria EsféricaNeste caso:

1r2

ddr

(r2 dy

dr

)− φ2y = 0 (3.72)

d2y

dr2+

2r

dy

dr− φ2y = 0 (3.73)


y =A cosh(φr) + Bsenh (φr)

r(3.74)


y|r=0 = �nito (3.75)

Como cosh(0)0

= ∞, A = 0 e:

y =Bsenh (φr)

r(3.76)

dy

dr

∣∣∣∣r=1

= B [φ cosh(φ)− senh (φ)] = Bim [1−Bsenh (φ)] (3.77)

B =Bim

φ cosh(φ) + (Bim − 1)senh (φ)(3.78)

y =Bimsenh (φr)

r [φ cosh(φ) + (Bim − 1)senh (φ)](3.79)


ηg = 3∫ 1

0yr2dr =

∫ 1

0

Bimsenh (φr)φ cosh(φ) + (Bim − 1)senh (φ)

rdr (3.80)

ηg =3Bim [φ cosh(φ)− senh (φ)]

φ2 [φ cosh(φ) + (Bim − 1) senh (φ)](3.81)

ηg =3Bim [φ− tanh(φ)]

φ2 [φ + (Bim − 1) tanh(φ)](3.82)

1ηg

=φ2 [φ− tanh(φ)]

3Bim [φ− tanh(φ)]+

Bimφ2 tanh(φ)3Bim [φ− tanh(φ)]

(3.83)

1ηg

=φ2

3Bim+

φ2

3 [φ coth(φ)− 1]=

1η

+φ2

3Bim(3.84)



0,0 0,2 0,4 0,6 0,8 1,0 0,0

0,2

0,4

0,6

0,8

1,0 C

once

ntra

ção

Raio

Figura 3.11: Per�l de concentração paramodulo de Thiele igual a 0,2, 1,0, 2,0, 10e 20, para número de Biot igual a 1,0.Geometria esférica.

0,01 0,1 1 10 100 0,0

0,2

0,4

0,6

0,8

1,0

Fato

r de

efe

tivid

ade

Modulo de Thiele

Figura 3.12: Fator de efetividade em fun-ção do modulo de Thiele, para número deBiot igual a 1,0. Geometria esférica.

3.2.3 Resolução por Diferenças FinitasFazendo:

dy

dr=

ui+1 − ui−1

2h(3.85)

d2y

dr2=

ui+1 − 2ui + ui−1

h2(3.86)

h =1

N + 1(3.87)

ri = ih (3.88)Em r = 0 −→ i = 0

s

r

dy

dr= lim

r→0

(s

r

dy

dx

)= s

d2y

dr2(3.89)

Então, em r = 0, temos:

d2y

dr2+

s

r

dy

dr= (s + 1)

d2y

dr2= φ2ym (3.90)

(s + 1)u1 − 2u0 + u−1

h2− φ2um

0 = 0 (3.91)

u1 − 2u0 + u−1 − h2

s + 1φ2um

0 = 0 (3.92)

Mas, pela primeira condição de contorno:

dy

dr

∣∣∣∣r=0

=u+1 − u−1

2h= 0 (3.93)


u−1 = u+1 (3.94)E assim:

2 (u1 − u0)− h2

s + 1φ2um

0 = 0 (3.95)

Para i = 1 · · ·Nui+1 − ui + ui−1

h2+

s

ih

ui+1 − ui−1

2h− φ2um

i = 0 (3.96)

(2i + s)ui+1 − 4iui + (2i− s)ui−1 − 2ih2φ2umi = 0 (3.97)

Para i = N + 1 −→ r = 1, 0, pela segunda condição de contorno (Equação 3.2).

uN+2 − uN

2h= Bim(1− uN+1) (3.98)

uN+2 = uN + 2hBim(1− uN+1) (3.99)Assim:

(2i + s)uN+2 − 4iuN+1 + (2i− s)uN − 2ih2φ2umN+1 = 0 (3.100)

(2i + s)[uN + 2hBim(1− uN+1)]− 4iuN+1 + (2i− s)uN − 2ih2φ2umN+1 = 0 (3.101)

4iuN + 4ihBim − 4ihBimuN+1 + 2shBim − 2shBimuN+1 − 4iuN+1 − 2ih2φ2umN+1 = 0(3.102)

2hBim(2i + s)− (4i + 4ihBim + 2shBim)uN+1 + 4iuN − 2ih2φ2umN+1 = 0 (3.103)

Assim, temos um sistema de (N+2) equações e (N+2) incógnitas (yi). Para resoluçãoaplicamos o método de Newton-Rapson:

uk+1 = uk + ∆uk (3.104)

Onde:

H×∆uk = −f (3.105)

f =

2 (u1 − u0)− h2

s+1φ2um0

(2i + s)ui+1 − 4iui + (2i− s)ui−1 − 2ih2φ2umi...

2hBim(2i + s)− (4i + 4ihBim + 2shBim)uN+1 + 4iuN − 2ih2φ2umN+1


Hi,j =∂fi

∂uj(3.106)

H =

∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣

B0 C0 0 0 0 0 . . . 0 0 0A1 B1 C1 0 0 0 . . . 0 0 00 A2 B2 C2 0 0 . . . 0 0 00 0 A3 B3 C3 0 . . . 0 0 0... ... ... ... ... ... . . .

... ... ...0 0 0 0 0 0 . . . AN BN CN

0 0 0 0 0 0 . . . 0 AN+1 BN+1

∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣Onde:

Ci =∂fi

∂ui+1

Bi =∂fi

∂ui

Ai =∂fi

∂ui−1

(3.107)

B0 =∂f0

∂u0= −2− h2

s+1φ2mum−10

C0 =∂f0

∂u1= 2

Ci = (2i + s)

Bi = −4i− 2ih2φ2mum−1i

Ai = (2i− s)

AN+1 =∂fN+1

∂uN= 4i

BN+1 =∂fN+1

∂uN+1= −(4i + 4hiBim + 2hsBim)− 2ih2mum−1

N+1

(3.108)

Então temos uma matriz tridiagonal. Um sistema da forma:

b1 c1 0 0 0 0 . . . 0 0 0a2 b2 c2 0 0 0 . . . 0 0 00 a3 b3 c3 0 0 . . . 0 0 00 0 a4 b4 c4 0 . . . 0 0 0... ... ... ... ... ... . . .

... ... ...0 0 0 0 0 0 . . . an−1 bn−1 cn−1

0 0 0 0 0 0 . . . 0 an bn

x1

x2

x3

x4...

xn−1

xn

=

r1

r2

r3

r4...

rn−1

rn


Tem como solução:

α1 =c1

b1(3.109)

γ1 =r1

b1(3.110)

αi =ci

bi − aiαi−1i = 2 . . . n (3.111)

γi =ri − aiγi−1

bi − aiαi−1i = 2 . . . n (3.112)

xn = γn (3.113)xi = γi − αixi+1 i = n− 1 . . . 1 (3.114)

3.2.4 Resolução por Colocação OrtogonalSubstituindo:

u = r2

A equação diferencial do balanço de massa �ca:

4ud2y

du2+ 2(s + 1)

dy

dx= φ2ym (3.115)

dy

du

∣∣∣∣u=1

=Bim

2(1− y) (3.116)

Com a primeira condição de contorno sendo satisfeita automaticamente.A solução pode ser aproximada por um polinômio de grau N. Seja ui os pontos de

colocação e yi os valores da função nos pontos de colocação, o polinômio que passa pelospontos (yi, ui) pode ser escrito como o polinômio interpolador de Lagrange:

yN =N+1∑

j=1

lj(u)yj (3.117)

Onde:

lj(u) =pN+1(u)

(u− uj)p(1)N+1(uj)

(3.118)

pN+1(u) =N+1∏

j=1

(u− uj) (3.119)

As derivadas primeira e segunda da solução aproximada serão dadas por:


dyN

du=

N+1∑

j=1

l(1)j (u)yj (3.120)

d2yN

du2=

N+1∑

j=1

l(2)j (u)yj (3.121)

As derivadas primeira e segunda da solução aproximada nos pontos de colocação serãodadas por:

dyN

du

∣∣∣∣u=ui

=N+1∑

j=1

l(1)j (ui)yj =

N+1∑

j=1

Ai,jyj (3.122)

d2yN

du2

∣∣∣∣u=ui

=N+1∑

j=1

l(2)j (ui)yj =

N+1∑

j=1

Bi,jyj (3.123)

Fazendo como que a equação diferencial seja válida nos pontos de colocação temos:

4ui

N+1∑

j=1

Bi,jyj + 2(s + 1)N+1∑

j=1

Ai,jyj = φ2ym (3.124)

N+1∑

j=1

[4uiBi,j + 2(s + 1)Ai,j ] yj − φ2ym = 0 (3.125)

Fazendo:Ci,j = 4uiBi,j + 2(s + 1)Ai,j (3.126)

Então, temos:

N+1∑

j=1

Ci,jyj − φ2ym = 0 (3.127)

N∑

j=1

Ci,jyj + Ci,N+1yN+1 − φ2ym = 0 (3.128)

Mas, pela primeira condição de contorno (Equação 3.116):

N+1∑

j=1

AN+1,jyj =Bim

2(1− yN+1) (3.129)

N∑

j=1

AN+1,jyj + AN+1,N+1yN+1 =Bim

2− Bim

2yN+1 (3.130)


yN+1 =Bim/2−∑N

j=1 AN+1,jyj

Bim/2 + AN+1,N+1(3.131)

Assim:

N∑

j=1

Ci,jyj + Ci,N+1

Bim/2−∑Nj=1 AN+1,jyj

Bim/2 + AN+1,N+1− φ2ym = 0 (3.132)

N∑

j=1

[Ci,j − Ci,N+1AN+1,j

Bim/2 + AN+1,N+1

]yj +

(Bim/2)Ci,N+1

Bim/2 + AN+1,N+1− φ2ym = 0 (3.133)

Fazendo:

Di,j = Ci,j − Ci,N+1AN+1,j

Bim/2 + AN+1,N+1(3.134)

vi =(Bim/2)Ci,N+1

Bim/2 + AN+1,N+1(3.135)

N∑

j=1

Di,jyj + vi − φ2ym = 0 (3.136)

Assim, temos um sistema de N equações e N incognitas (yi). Para resolução aplicamoso método de Newton-Rapson:

uk+1 = uk + ∆uk (3.137)

Onde:

H×∆uk = −f (3.138)

fi =N∑

j=1

Di,jyj + vi − φ2ym (3.139)

Hi,j =∂fi

∂yj= Di,j − φ2 ∂ym

i

∂yiδi,j = Di,j −mφ2ym−1

i δi,j (3.140)

O fator de efetividade pode ser calculado como:

η =∫ 1

0ymdu

s+12 (3.141)

η =s + 1

2

∫ 1

0ymu

s−12 du (3.142)


Escolhendo os pontos de colocação como as raízes do polinômio de Jacob P(α,β)N (u),

onde α = 0 ou 1, 0 e β = s−12 esta integral pode ser calculada por quadratura utilizando-se

os pontos yi determinados anteriormente.Com α = 0, temos o método dos momentos e a integral pode ser calculada por

quadratura de Gauss-Jacob:

η =N∑

j=1

wjym (3.143)

sendo wj os pesos da quadratura de Gauss.Com α = 1, 0, temos o método de Galerkin e a integral pode ser calculada por

quadratura de Radau:

η =N∑

j=1

wjym (3.144)

sendo wj os pesos da quadratura de Radau.

3.3 Sistemas Não-isotérmicosOs balanços de massa e de energia para uma partícula catalítica em condições não-

isotérmicas fornece as seguintes equações:

1r′s

ddr′

(Der

′s dCA

dr′

)= rA

(CA, T

′) (3.145)

− 1r′s

ddr′

(λer

′s dT′

dr′

)= (−∆H) rA

(CA, T

′) (3.146)

onde:

(rA) = k0v exp

[−Ea

RT ′

]Cm

A (3.147)


r = 0

dCA

dr′= 0

dT′

dr′= 0

(3.148)

r = R

DedC

dr′= kc (CA0 − CA)

λedT

′

dr′= hf (T0 − T )

(3.149)

Adimensionalizando:

3.3 Sistemas Não-isotérmicos 43

y =CA

CA0(3.150)

r =r′

R(3.151)

T =T′

T0(3.152)

Os balanços de massa e de energia para uma partícula catalítica em condições não-isotérmicas na forma adimensional então �cam:

1rs

ddr

(rs dy

dr

)= φ2 exp

[γ

(1− 1

T

)]ym (3.153)

− 1rs

ddr

(rs dT

dr

)= βφ2 exp

[γ

(1− 1

T

)]ym (3.154)


r = 0

dy

dr= 0

dT

dr= 0

(3.155)

r = 1

dy

dr= Bim (1− y)

dT

dr= Bi (1− T )

(3.156)

onde:

β =(−∆H) DeCA0

λeT0

γ =Ea

RT0

φ = R

√k0

v exp (γ) Cm−1A0

De

Bi =hfR

λe

Bim =kcR

De

(3.157)

Multiplicando a Equação 3.153 por −β e somando a Equação 3.154, temos:


1rs

ddr

[rs

(dy

dr+ β

dT

dr

)]= 0 (3.158)

ddr

[rs

(β

dy

dr+

dT

dr

)]= 0 (3.159)

Integrando:

rs

(β

dy

dr+

dT

dr

)= A (3.160)

βdy

dr+

dT

dr=

A

rs(3.161)

Como, pela condição de contorno (Equação 3.155):

dT

dr

∣∣∣∣r=0

=dy

dr

∣∣∣∣r=0

= 0 (3.162)

A

rs

∣∣∣∣r=0

= 0 (3.163)

Isto implica que A = 0, então:

βdy

dr+

dT

dr= 0 (3.164)

Mas em r = 1 pela condição de contorno (Equação 3.155):

β {Bim [1− y (1)]}+ {Bi [1− T (1)]} = 0 (3.165)

T (1) + βy (1) = 1 + βy (1) + βBimBi

[1− y (1)] (3.166)

Integrando a Equação 3.164 de r=r até 1, temos:

T (r) + βy(r) = T (1) + βy(1) (3.167)

Substituindo a Equação 3.166:

T (r) + βy(r) = 1 + βy (1) + βBimBi

[1− y (1)] (3.168)

T (r) = 1 + β [y (1)− y (r)] + βBimBi

[1− y (1)] (3.169)

Então, a Equação de balanço de massa (Equação 3.153) �ca:

1rs

ddr

(rs dy

dr

)= φ2< [y(1), y(r)] (3.170)

Com:


< [y(1), y(r)]

= exp

γ

1− 1

1 + β [y (1)− y (r)] + βBimBi

[1− y (1)]

ym (3.171)

d2y

dr2+

s

r

dy

dr= φ2< [y(1), y(r)] (3.172)


dy

dr

∣∣∣∣r=0

= 0 (3.173)

dy

dr

∣∣∣∣r=1

= Bim(1− y) (3.174)

3.3.1 Resolução por Diferenças FinitasFazendo:

dy

dr=

ui+1 − ui−1

2h(3.175)

d2y

dr2=

ui+1 − 2ui + ui−1

h2(3.176)

h =1

N + 1(3.177)

ri = ih (3.178)Em r = 0 −→ i = 0

s

r

dy

dr= lim

r→0

(s

r

dy

dx

)= s

d2y

dr2(3.179)

Então, em r = 0, temos:

d2y

dr2+

s

r

dy

dr= (s + 1)

d2y

dr2= φ2< [uN+1, u0] (3.180)

(s + 1)u1 − 2u0 + u−1

h2− φ2< [uN+1, u0] = 0 (3.181)

u1 − 2u0 + u−1 − h2

s + 1φ2< [uN+1, u0] = 0 (3.182)

Mas, pela primeira condição de contorno:

dy

dr

∣∣∣∣r=0

=u+1 − u−1

2h= 0 (3.183)

u−1 = u+1 (3.184)


E assim:

2 (u1 − u0)− h2

s + 1φ2< [uN+1, u0] = 0 (3.185)

Para i = 1 · · ·Nui+1 − ui + ui−1

h2+

s

ih

ui+1 − ui−1

2h− φ2< [uN+1, ui] = 0 (3.186)

(2i + s)ui+1 − 4iui + (2i− s)ui−1 − 2ih2φ2< [uN+1, ui] = 0 (3.187)Para i = N + 1 −→ r = 1, 0, pela segunda condição de contorno (Equação 3.173)

uN+2 − uN

2h= Bim(1− uN+1) (3.188)

uN+2 = uN + 2hBim(1− uN+1) (3.189)Assim:

(2i + s)uN+2 − 4iuN+1 + (2i− s)uN − 2ih2φ2< [uN+1] = 0 (3.190)

(2i + s)[uN + 2hBim(1− uN+1)]− 4iuN+1 + (2i− s)uN − 2ih2φ2< [uN+1] = 0(3.191)

4iuN + 4ihBim − 4ihBimuN+1 + 2shBim − 2shBimuN+1 − 4iuN+1 − 2ih2φ2< [uN+1] = 0(3.192)

2hBim(2i + s)− (4i + 4ihBim + 2shBim)uN+1 + 4iuN − 2ih2φ2< [uN+1] = 0(3.193)

Assim, temos um sistema de (N+2) equações e (N+2) incognitas (yi). Para resoluçãoaplicamos o método de Newton-Rapson:

uk+1 = uk + ∆uk (3.194)

Onde:

H×∆uk = −f (3.195)

f =

2 (u1 − u0)− h2

s+1φ2< [uN+1, u0](2i + s)ui+1 − 4iui + (2i− s)ui−1 − 2ih2φ2< [uN+1, ui]...2hBim(2i + s)− (4i + 4ihBim + 2shBim)uN+1 + 4iuN

−2ih2φ2< [uN+1]


Hi,j =∂fi

∂uj(3.196)

H =

∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣

B0 C0 0 0 0 0 . . . 0 0 0A1 B1 C1 0 0 0 . . . 0 0 00 A2 B2 C2 0 0 . . . 0 0 00 0 A3 B3 C3 0 . . . 0 0 0... ... ... ... ... ... . . .

... ... ...0 0 0 0 0 0 . . . AN BN CN

0 0 0 0 0 0 . . . 0 AN+1 BN+1

∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣Onde:

Ci =∂fi

∂ui+1

Bi =∂fi

∂ui

Ai =∂fi

∂ui−1

(3.197)

B0 =∂f0

∂u0= −2− h2

s+1φ2 ∂<∂u0

B0 = −2− h2

s+1φ2mum−10 exp

{γ

[1− 1

1 + β (uN+1 − u0) + β BimBi (1− uN+1)

]}×

{m

u0− γβ[

1 + β (uN+1 − u0) + β BimBi (1− uN+1)

]2

}

C0 =∂f0

∂u1= 2

(3.198)

Ci = (2i + s)

Bi = −4i− 2ih2φ2 ∂<∂ui

Bi = −4i− 2ih2φ2umi exp

{γ

[1− 1

1 + β (uN+1 − ui) + β BimBi (1− uN+1)

]}×

{m

ui− γβ[

1 + β (uN+1 − ui) + β BimBi (1− uN+1)

]2

}

Ai = (2i− s)(3.199)


AN+1 =∂fN+1

∂uN= 4i

BN+1 =∂fN+1

∂uN+1= −(4i + 4hiBim + 2hsBim)− 2ih2φ2 ∂<

∂uN+1

BN+1 = −(4i + 4hiBim + 2hsBim)− 2ih2φ2umN+1 exp

{γ

[1− 1

1 + β BimBi (1− uN+1)

]}×

{m

uN+1− γβ[

1 + β BimBi (1− uN+1)

]2

}

(3.200)

3.3.2 Resolução por Colocação OrtogonalSubstituindo:

u = r2

4ud2y

du2+ 2(s + 1)

dy

dx= φ2< [y(1), y(r)] (3.201)

(3.202)dy

du

∣∣∣∣u=1

=Bim

2(1− y) (3.203)

A solução pode ser aproximada por um polinômio de grau N. Seja ui os pontos decolocação e yi os valores da função nos pontos de colocação, o polinômio que passa pelospontos (yi, ui) pode ser escrito como o polinômio interpolador de Lagrange:

yN =N+1∑

j=1

lj(u)yj (3.204)

Onde:

lj(u) =pN+1(u)

(u− uj)p(1)N+1(uj)

(3.205)

pN+1(u) =N+1∏

j=1

(u− uj) (3.206)


dyN

du=

N+1∑

j=1

l(1)j (u)yj (3.207)

d2yN

du2=

N+1∑

j=1

l(2)j (u)yj (3.208)



dyN

du

∣∣∣∣u=ui

=N+1∑

j=1

l(1)j (ui)yj =

N+1∑

j=1

Ai,jyj (3.209)

d2yN

du2

∣∣∣∣u=ui

=N+1∑

j=1

l(2)j (ui)yj =

N+1∑

j=1

Bi,jyj (3.210)

Fazendo como que a equação diferencial seja válida nos pontos de colocação temos:

4ui

N+1∑

j=1

Bi,jyj + 2(s + 1)N+1∑

j=1

Ai,jyj = φ2< [yN+1, yi] (3.211)

N+1∑

j=1

[4uiBi,j + 2(s + 1)Ai,j ] yj − φ2< [yN+1, yi] = 0 (3.212)

Fazendo:

Ci,j = 4uiBi,j + 2(s + 1)Ai,j (3.213)

N+1∑

j=1

Ci,jyj − φ2< [yN+1, yi] = 0 (3.214)

N∑

j=1

Ci,jyj + Ci,N+1yN+1 − φ2< [yN+1, yi] = 0 (3.215)

Mas, pela condição de contorno (Equação 3.203):

N+1∑

j=1

AN+1,jyj =Bim

2(1− yN+1) (3.216)

N∑

j=1

AN+1,jyj + AN+1,N+1yN+1 =Bim

2− Bim

2yN+1 (3.217)

yN+1 =Bim/2−∑N

j=1 AN+1,jyj

Bim/2 + AN+1,N+1(3.218)

Assim:

N∑

j=1

Ci,jyj + Ci,N+1

Bim/2−∑Nj=1 AN+1,jyj

Bim/2 + AN+1,N+1− φ2< [yN+1, yi] = 0 (3.219)


N∑

j=1

[Ci,j − Ci,N+1AN+1,j

Bim/2 + AN+1,N+1

]yj +

Bim/2Ci,N+1

Bim/2 + AN+1,N+1− φ2< [yN+1, yi] = 0 (3.220)

Fazendo:


Bim/2 + AN+1,N+1(3.221)

vi =Bim

2 Ci,N+1

Bim/2 + AN+1,N+1(3.222)

N∑

j=1

Di,jyj + vi − φ2< [yN+1, yi] = 0 (3.223)

Assim, temos um sistema de N equações e N incognitas (yi). Para resolução aplicamoso método de Newton-Rapson:

uk+1 = uk + ∆uk (3.224)

Onde:

H×∆uk = −f (3.225)

fi =N∑

j=1

Di,jyj + vi − φ2< [yN+1, yi] (3.226)

Hi,j =∂fi

∂yj= Di,j − φ2 ∂< [yN+1, yi]

∂yjδi,j (3.227)

Hi,j = Di,j − φ2ymi exp

{γ

[1− 1

1 + β (yN+1 − yi) + β BimBi (1− yN+1)

]}×

{m

yi− γβ[

1 + β (yN+1 − yi) + β BimBi (1− yN+1)

]2

}δi,j

(3.228)

O fator de efetividade pode ser calculado como:

η =∫ 1

0< [y(1), y(u)] du

s+12 (3.229)

η =s + 1

2

∫ 1

0< [y(1), y(u)]u

s−12 du (3.230)


Escolhendo os pontos de colocação como as raízes do polinômio de Jacob P(α,β)N (u),

onde α = 0 ou 1, 0 e β = s−12 esta integral pode ser calculada por quadratura utilizando-se

os pontos yi determinados anteriormente.Com α = 0, temos o método dos momentos e a integral pode ser calculada por

quadratura de Gauss-Jacob:

η =N+1∑

j=1

wj< [yN+1, yj ] (3.231)

Com α = 1, 0, temos o método de Galerkin e a integral pode ser calculada porquadratura de Radau:

η =N+1∑

j=1

wj< [yN+1, yj ] (3.232)

Considerando os seguintes valores para os parâmetros: partícula esférica (s = 2),β = 0, 6, γ = 20, φ = 0, 5, Bim = ∞ e Bi = ∞ e resolvendo o sistema com um únicoponto de colocação, temos:

P(1, 1

2)

N (u) =73u− 1 (3.233)

u1 =37

(3.234)

Fazendo com que a equação seja válida em u1, chegamos a equação não-linear:

10, 5(1− y1) = 0, 25y1 exp[20

(1− 1

1 + 0, 6(1− y1)

)](3.235)

Plotando g(y1) = 10, 5(1−y1) versus h(y1) = 0, 25y1 exp[20

(1− 1

1 + 0, 6(1− y1)

)],

Figura 3.13, podemos ver que este sistema apresenta 3 soluções:

y1 = (0, 0255452; 0, 7419215; 0, 9655210) (3.236)


0,0 0,2 0,4 0,6 0,8 1,0 0

4

8

12

16

20

24

28

32

36

g(y)

e h

(y)

y

Figura 3.13: Partícula catalítica não isotérmica

Capítulo 4

Reatores de Leito Fixo

4.1 Modelo Unidimensional BásicoNeste modelo assume-se que gradientes de temperatura e de concentração ocorrem

somente na direção axial. O único mecanismo de transporte operando nesta direção é o�uxo total, que é considerado empistonado. Assim, o balanço de massa e de energia parauma reação simples realizada em um tubo cilíndrico nos dá:

− ddz

(usCA) = raρb (4.1)

usρgcpdT

dz= (−∆H)rAρb − 4

U

dt(T − TR) (4.2)

−dPt

dz= f

ρgu2s

dp(4.3)


z = 0

CA = CA0

T = T0

Pt = Pt0

(4.4)

onde:us: velocidade super�cialρg: densidade do gásρb: densidade do leitorA: velocidade de formação de Acp: capacidade calorí�caTR temperatura da parede do reatorU: coe�ciente de transferência de calor global, dado por

1U

=1αi

+d

λ

Ab

Am+

1αw

Ab

Au(4.5)

αi: coe�ciente de transferência de calor internaαw: coe�ciente de transferência de calor externaAb: área de transferência de calor do lado do leito

54 Reatores de Leito Fixo

λ: condutividade térmica da paredeAu: área de transferência de calor do lado internoAm: média logarítmica de Ab e Au

f : fator de atrito que, segundo Ergun, é dado por:

f =1− ε

ε3

[1, 75 +

150 (1− ε)Re

](4.6)

dp: diâmetro da partícula de catalisadorPt: pressão total

Exemplo Reação de oxidação de o-xileno para produzir o anidrido ftálico:Dados:L = 30 mdt = 2,54 cmT = 335◦-415◦Cdp = 3 mmρp = 1300 kg/m3

FA0 = 4684 kg/m2hρg = 1293 kg/m3

∆H = -1285409 kJ/kgmolcp = 0,992 kJ/kgKU = 0,096 kJ/m2sTR = 352◦CrA = −kP 0

BPA

onde:ln(k) = 19, 837− 13636

T(4.7)

Os balanços de massa e energia, então, �cam:

usdPA

dz+

MmPtρb

ρgkP 0

BPA = 0 (4.8)

usρgcpdT

dz− (−∆H)ρbkP 0

BPA +4U

dt(T − TR) = 0 (4.9)

Condições de contorno:z = 0

{PA = P 0 = 0, 015barT = T0 = TR = 625K

(4.10)

Nas Figuras 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8 são apresentados os per�s de temperaturae concentração em um reator tubular variando-se a concentração de reagente na alimen-tação. Estes per�s foram obtidos pela integração numérica do sistema de equações acimautilizando o método de Range-Kutta. Nestas �guras observa-se o aparecimento de umponto quente (do inglês hot spot) no reator. A intensidade deste ponto quente dependefortemente das condições de alimentação, sendo que alterações minimas desta podemlevar ao descontrole do reator (comparar Figuras 4.6 e 4.8).

4.1 Modelo Unidimensional Básico 55

0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,000

0,001

0,002

0,003

0,004

0,005

0,006

0,007

0,008

0,009

0,010

Con

cent

raçã

o

Comprimento

Figura 4.1: Per�l de concentração parauma pressão parcial de alimentação iguala 0,01 bar.

0,0 0,5 1,0 1,5 2,0 2,5 3,0

626

628

630

632

634

636

638

Tem

pera

tura

Comprimento

Figura 4.2: Per�l de temperatura parauma pressão parcial de alimentação iguala 0,01 bar.

0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,000

0,002

0,004

0,006

0,008

0,010

0,012

0,014

Con

cent

raçã

o

Comprimento


0,0 0,5 1,0 1,5 2,0 2,5 3,0

626

628

630

632

634

636

638

640

642

644

646

648

650

Tem

pera

tura

Comprimento



0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,000

0,002

0,004

0,006

0,008

0,010

0,012

0,014

0,016

Con

cent

raçã

o

Comprimento


0,0 0,5 1,0 1,5 2,0 2,5 3,0 625

630

635

640

645

650

655

660

665

Tem

pera

tura

Comprimento


0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,000

0,002

0,004

0,006

0,008

0,010

0,012

0,014

0,016

0,018

Con

cent

raçã

o

Comprimento


0,0 0,5 1,0 1,5 2,0 2,5 3,0

650

700

750

800

850

900

950

1000

1050

Tem

pera

tura

Comprimento


4.2 Modelo Pseudo-Homogêneo BidimensionalNeste modelo assume-se que gradientes de temperatura e de concentração ocorrem

nas direções axial e radial. O único mecanismo de transporte operando na direção axialé o �uxo total, que é considerado empistonado. Assim, o balanço de massa e de energiapara uma reação simples realizada em um tubo cilíndrico nos dá:

εbDeR

(∂2C

∂r2+

1r

∂C

∂r

)− us

∂C

∂z− ρbrA = 0 (4.11)

λeR

(∂2T

∂r2+

1r

∂T

∂r

)− usρgcp

∂T

∂z− ρb(−∆H)rA = 0 (4.12)


z = 0{

C = C0

T = T0(4.13)

4.2 Modelo Pseudo-Homogêneo Bidimensional 57

r = 0

∂C

∂r= 0

∂T

∂r= 0

(4.14)

r = rt

∂C

∂r= 0

∂T

∂r= − αw

λeR(TR − TW )

(4.15)

Exemplo Reação de oxidação de o-xileno para produzir o anidrido ftálico:

Ak // B

Dados:L = 30 mdt = 2,54 cmdp = 3 mmρp = 1300 kg/m3

FA0 = 4684 kg/m2hρg = 1293 kg/m3

εb = 0,4∆H = 1, 285× 106 kJ/kgmolλer = 0, 78× 10−3 kJ/msKDer = 1, 08× 10−3 cm2/sαw = 0,156 kJ/m2sKP0

B = 0,211 barcp = 0,992 kJ/kgKTR = Tw = 352◦CrA = −kP 0

BPA = kPB0PA0(1− x)onde:

ln(k) = 19, 837− 13636T

(4.16)

O balanço de massa e de energia deste sistema �ca:

εbDeR

(∂2CA

∂r′2+

1r′

∂CA

∂r′

)− us

∂CA

∂z′− ρbrA = 0 (4.17)

λeR

(∂2T

′

∂r′2+

1r′

∂T′

∂r′

)− usρgcp

∂T′

∂z′− ρb(−∆H)rA = 0 (4.18)


z′= 0

{CA = CA0

T′= T0

(4.19)


r′= 0

∂CA

∂r′= 0

∂T′

∂r′= 0

(4.20)

r′= rt

∂CA

∂r′= 0

∂T′

∂r′= − αw

λeR

(T′R − T

′W

) (4.21)

Fazendo:

z =z′

rt(4.22)

r =r′

rt(4.23)

x =CA0 − CA

CA0(4.24)

∂x

∂z= α1

(∂2x

∂r2+

1r

∂x

∂r

)− β1rA (4.25)

∂T

∂z= α2

(∂2T

∂r2+

1r

∂T

∂r

)+ β2rA (4.26)

α1 =εbDer

usrt(4.27)

α2 =λer

ρguscprt(4.28)

β1 =ρbrtMm

ρgusyA0(4.29)

β2 =ρb(−∆H1)rt

ρgusT0cP(4.30)


z = 0{

x = 0T = 1, 0

(4.31)

r = 0

∂x

∂r= 0

∂T

∂r= 0

(4.32)


r = 1, 0

∂x

∂r= 0

∂T

∂r= Bi (TW − T )

(4.33)

Bi =αwR

λeR(4.34)

4.2.1 Resolução por Colocação OrtogonalFazendo u = r2, temos:

∂x

∂z= 4α1

(u

∂2x

∂u2+

∂x

∂u

)− β1rA (4.35)

∂T

∂z= 4α2

(∂2T

∂u2+

∂T

∂u

)+ β2rA (4.36)


z = 0{

x = 0T = 1, 0

(4.37)

u = 1, 0

∂x

∂u= 0

∂T

∂u=

Bi

2(TW − T )

(4.38)

A solução pode ser aproximada por um polinômio de grau N. Seja ui os pontosde colocação e xi(z) e Ti(z) os valores da conversão e da temperatura nos pontos decolocação, respectivamente. Assim, o polinômio que passa pelos pontos [xi(z), Ti(z), ui]pode ser escrito como o polinômio interpolador de Lagrange:

xN (u, z) =N+1∑

j=1

lj (u) xj (z) (4.39)

TN (u, z) =N+1∑

j=1

lj (u) Tj (z) (4.40)

Onde:

lj(u) =pN+1(u)

(u− uj)p(1)N+1(uj)

(4.41)

pN+1(u) =N+1∏

j=1

(u− uj) (4.42)



∂xN (z)∂u

=N+1∑

j=1

l(1)j (u)xj (z) (4.43)

∂2xN (z)∂u2

=N+1∑

j=1

l(2)j (u)xj (z) (4.44)

∂TN (z)∂u

=N+1∑

j=1

l(1)j (u)Tj (z) (4.45)

∂2TN (z)∂u2

=N+1∑

j=1

l(2)j (u)Tj (z) (4.46)


∂xN (z)∂u

∣∣∣∣u=ui

=N+1∑

j=1

l(1)j (ui)xj (z) =

N+1∑

j=1

Ai,jxj (z) (4.47)

∂2xN (z)∂u2

∣∣∣∣u=ui

=N+1∑

j=1

l(2)j (ui)xj (z) =

N+1∑

j=1

Bi,jxj (z) (4.48)

∂TN (z)∂u

∣∣∣∣u=ui

=N+1∑

j=1

l(1)j (ui)Tj (z) =

N+1∑

j=1

Ai,jTj (z) (4.49)

∂2TN (z)∂u2

∣∣∣∣u=ui

=N+1∑

j=1

l(2)j (ui)Tj (z) =

N+1∑

j=1

Bi,jTj (z) (4.50)

Fazendo como que as equações diferenciais sejam válidas nos pontos de colocaçãotemos:

dxi (dz)dz

= 4α1

N+1∑

j=1

uiBi,jxj (z) +N+1∑

j=1

Ai,jxj (z)

− β1rB (4.51)

dTi (z)dz

= 4α2

N+1∑

j=1

uiBi,jTj (z) +N+1∑

j=1

Ai,jTj (z)

+ β2rA (4.52)

Ci,j = 4 (uiBi,j + Ai,j) (4.53)

dxi (z)dz

= α1

N+1∑

j=1

Ci,jxj (z)

− β1rA (4.54)


dTi (z)dz

= α2

N+1∑

j=1

Ci,jTj (z)

+ β2rA (4.55)

dxi (z)dz

= α1

N∑

j=1

Ci,jxj (z) + Ci,N+1xN+1 (z)

− β1rA (4.56)

dTi (z)dz

= α2

N∑

j=1

Ci,jTj (z) + Ci,N+1TN+1 (z)

+ β2rA (4.57)

Mas, pela condição de contorno:

∂x

∂u

∣∣∣∣u=1

=N+1∑

j=1

AN+1,jxj (z) = 0 (4.58)

N∑

j=1

AN+1,jxj (z) + AN+1,N+1xN+1 (z) = 0 (4.59)

xN+1 (z) = −∑N

j=1 AN+1,jxj (z)AN+1,N+1

(4.60)

∂T

∂u

∣∣∣∣u=1

=N+1∑

j=1

AN+1,jTj (z) = Bi/2 [1− TN+1 (z)] (4.61)

N∑

j=1

AN+1,jTj (z) + AN+1,N+1TN+1 (z) = Bi/2 [1− TN+1 (z)] (4.62)

TN+1 (z) =Bi/2−∑N

j=1 AN+1,jTj (z)Bi/2 + AN+1,N+1

(4.63)

dxi (z)dz

= α1

N∑

j=1

Ci,jxj (z)− Ci,N+1∑N

j=1 AN+1,jxj (z)AN+1,N+1

− β1rA (4.64)

dTi (z)dz

= α2

N∑

j=1

Ci,jTj (z) +Ci,N+1

[Bi/2−∑N

j=1 AN+1,j

]Tj (z)

Bi/2 + AN+1,N+1

+ β2rA (4.65)


AN+1,N+1(4.66)

Ei,j = Ci,j − Ci,N+1AN+1,j

Bi/2 + AN+1,N+1(4.67)


vi =Ci,jBi/2

Bi/2 + AN+1,N+1(4.68)

dxi (z)dz

= α1

N∑

j=1

Di,jxj (z)

− β1rA (4.69)

dTi (z)dz

= α2

N∑

j=1

Ei,jTj (z) + vi

+ β2rA (4.70)

Assim, a integração deste sistema de equações ao longo do eixo z, resulta no per�lde temperatura e de concentração. Nas Figuras 4.9, 4.10, 4.11, 4.11, 4.13 e 4.14 estãorepresentados os per�s de concentração de A e de temperatura no reator.

4.2.2 Resolução por Diferenças FinitasAs equações do balanço de massa e de energia:

∂x (z)∂z

= α1

[∂2x (z)

∂r2+

1r

∂x (z)∂r

]− β1rA (4.71)

∂T (z)∂z

= α2

[∂2T (z)

∂r2+

1r

∂T (z)∂r

]+ β2rA (4.72)

Com as condições de contorno:

z = 0{

x = 0T = 1, 0

(4.73)

r = 0

∂x

∂r= 0

∂T

∂r= 0

(4.74)

r = 1, 0

∂x

∂r= 0

∂T

∂r= Bi (TW − T )

(4.75)

Podem ser resolvidas por diferenças �nitas. Assim, fazendo:

dx

dr=

xi+1 − xi−1

2h(4.76)

d2x

dr2=

xi+1 − 2xi + xi−1

h2(4.77)

dT

dr=

Ti+1 − Ti−1

2h(4.78)


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

Con

cent

raçã

o

Comprimento

Figura 4.9: Per�l de concentração parapressão parcial de reagente na alimenta-ção igual 0,005 bar.

0,0 0,5 1,0 1,5 2,0 2,5 3,0 1,00

1,01

1,02

1,03

1,04

1,05

1,06

Tem

pera

tura

Comprimento

Figura 4.10: Per�l de Temperatura parapressão parcial de reagente na alimenta-ção igual a 0,005 bar.

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

Con

cent

raçã

o

Comprimento

Figura 4.11: Per�l de concentração parapressão parcial de reagente na alimenta-ção igual a 0,008 bar.

0,0 0,5 1,0 1,5 2,0 2,5 3,0 1,00

1,05

1,10

1,15

1,20

1,25

1,30

1,35

Tem

pera

tura

Comprimento


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

Con

cent

raçã

o

Comprimento

Figura 4.13: Per�l de concentração parapressão parcial de reagente na alimenta-ção igual a 0,009 bar.

0,0 0,5 1,0 1,5 2,0 2,5 3,0 1,0

1,1

1,2

1,3

1,4

1,5

1,6

Tem

pera

tura

Comprimento



d2T

dr2=

Ti+1 − 2Ti + Ti−1

h2(4.79)

h =1

N + 1(4.80)

ri = ih (4.81)

Temos:

dxi (z)dz

= α1

[xi+1 (z)− 2xi (z) + xi−1 (z)

h2+

1ih

xi+1 (z)− xi−1 (z)2h

]− β1rA (4.82)

dTi (z)dz

= α2

[Ti+1 (z)− 2Ti (z) + Ti−1 (z)

h2+

1ih

Ti+1 (z)− Ti−1 (z)2h

]+ β2rA (4.83)

dxi (z)dz

= α1

{[2i + 1]xi+1 (z)− 4ixi (z) + [2i− 1]xi−1 (z)

2ih2

}− β1rA (4.84)

dTi (z)dz

= α2

{[2i + 1]Ti+1 (z)− 4iTi (z) + [2i− 1]Ti−1 (z)

2ih2

}+ β2rA (4.85)

Em r = 0 −→ i = 0

limr→0

1r

∂x

∂r=

∂2x

∂r2(4.86)

limr→0

1r

∂T

∂r=

∂2T

∂r2(4.87)

dx0 (z)dz

= 2α1∂2x

∂r2− β1rA = 2α1

x1 (z)− 2x0 (z) + x−1 (z)h2

− β1rA (4.88)

dT0 (z)dz

= 2α2∂2T

∂r2+ β2rA = 2α1

T1 (z)− 2T0 (z) + T−1 (z)h2

− β2rA (4.89)


dx (z)dr

=x1 (z)− x−1 (z)

2h= 0 (4.90)

dT (z)dr

=T1 (z)− T−1 (z)

2h= 0 (4.91)

Então:


x−1 (z) = x+1 (z) (4.92)

T−1 (z) = T+1 (z) (4.93)

dx0 (z)dz

=4α1 [x1 (z)− x0 (z)]

h2− β1rA (4.94)

dT0 (z)dz

=4α1 [T1 (z)− T0 (z)]

h2+ β2rA (4.95)

Em r = 1 −→ i = N + 1, pela condição de contorno (Equação 4.74):

dx

dr=

xN+2 (z)− xN (z)2h

= 0 (4.96)

dT

dr=

TN+2 (z)− TN (z)2h

= Bi [Tw − TN+1 (z)] (4.97)

xN+2 (z) = xN (z) (4.98)

TN+2 (z) = TN (z) + 2hBiTw − 2hBiTN+1 (z) (4.99)

Substituindo:

dxN+1 (z)dz

= α1

{[2i + 1]xN (z)− 4ixN+1 (z) + [2i− 1]xN (z)

2ih2

}− β1rA (4.100)

dTN+1 (z)dz

=

α2

{[2i + 1] [TN (z) + 2hBiTw − 2hBiTN+1 (z)]− 4iTi (z) + [2i− 1]Ti−1 (z)

2ih2

}+ β2rA

(4.101)

dxN+1 (z)dz

=2α1 [xN (z)− xN+1 (z)]

h2− β1rA (4.102)

dTN+1 (z)dz

=

α24iTN (z) [2i + 1] 2hBiTw − {[2i + 1] 2hBi + 4i}TN+1 (z)

2ih2+ β2rA

(4.103)

Assim:


dx (z)dz

=

4α1 [x1 (z)− x0 (z)]h2

− β1rA

α1

{[2i + 1]xi+1 (z)− 4ixi (z) + [2i− 1]xi−1 (z)

2ih2

}− β1rA

...2α1 [xN (z)− xN+1 (z)]

h2− β1rA

dT (z)dz

=

4α1 [T1 (z)− T0 (z)]h2

+ β2rA

α2

{[2i + 1] Ti+1 (z)− 4iTi (z) + [2i− 1]Ti−1 (z)

2ih2

}+ β2rA

...α2

4iTN (z) [2i + 1] 2hBiTw − {[2i + 1] 2hBi + 4i}TN+1 (z)2ih2

+ β2rA

Exemplo Seja a reação de oxidação catalítica do o-xileno sobre o pentóxido de vanádio:

Ak3

ÂÂ@@@

@@@@k1 // B

k2 // C

C

Dados:L = 30 mdt = 2,54 cmdp = 3 mmρp = 1300 kg/m3

FA0 = 4684 kg/m2hρg = 1293 kg/m3

εb = 0,4∆H1 = 1, 285× 106 kJ/kgmol∆H3 = 4, 564× 106 kJ/kgmolλer = 0, 78× 10−3 kJ/msKDer = 1, 08× 10−3 cm2/sαw = 0,156 kJ/m2sKP0

B = 0,211 barcp = 0,992 kJ/kgKTR = Tw = 352◦CAs equações da taxa são dadas por:

rA = (k1 + k3)yA0yO2(1− xA) (4.104)

rB = yA0yO2 [k1(1− xA)− k2xB] (4.105)


rC = yA0yO2 [k2xB + k3(1− xA)] (4.106)onde:

ln(k1) =−13588

T ′ + 19, 837 (4.107)

ln(k2) =−15803

T ′ + 20, 86 (4.108)

ln(k3) =−14394

T ′ + 18, 97 (4.109)

O balanço de massa e de energia deste sistema �ca:

εbDeR

(∂2CB

∂r′2+

1r′

∂CB

∂r′

)− us

∂CB

∂z′+ ρbrB = 0 (4.110)

εbDeR

(∂2CC

∂r′2+

1r′

∂CC

∂r′

)− us

∂CC

∂z′+ ρbrC = 0 (4.111)

λeR

(∂2T

′

∂r′2+

1r′

∂T ‘

∂r′

)− usρgcp

∂T′

∂z′+ ρb(−∆H1)rB + ρb(−∆H3)rC = 0 (4.112)


z′= 0

CA = CA0

CB = 0CC = 0T′= T0

(4.113)

r′= 0

∂CB

∂r′= 0

∂CC

∂r′= 0

∂T′

∂r′= 0

(4.114)

r′= rt

∂CB

∂r′= 0

∂CC

∂r′= 0

∂T′

∂r′= − αw

λeR

(T′R − T

′W

)

(4.115)

Fazendo:


z =z′

rt(4.116)

r =r′

rt(4.117)

xA =CA0 − CA

CA0(4.118)

xB =CB

CA0(4.119)

xC =CC

CA0(4.120)

∂xB

∂z= α1

(∂2xB

∂r2+

1r

∂xB

∂r

)+ β1rB (4.121)

∂xC

∂z= α1

(∂2xC

∂r2+

1r

∂xC

∂r

)+ β1rC = 0 (4.122)

∂T

∂z= α2

(∂2T

∂r2+

1r

∂T

∂r

)+ β2rB + β3rC = 0 (4.123)

α1 =εbDer

usrt(4.124)

α2 =λer

ρguscprt(4.125)

β1 =ρbrtMm

ρgusyA0(4.126)


ρgusT0cP(4.127)


ρgusT0cP(4.128)


z = 0

xB = 0xC = 0T = 1, 0

(4.129)

r = 0

∂xB

∂r= 0

∂xC

∂r= 0

∂T

∂r= 0

(4.130)


r = 1, 0

∂xB

∂r= 0

∂xC

∂r= 0

∂T

∂r= Bi (TW − T )

(4.131)

Bi =αwR

λeR(4.132)

4.2.3 Resolução por Colocação OrtogonalFazendo u = r2, temos:

∂xB

∂z= 4α1

(u

∂2xB

∂u2+

∂xB

∂u

)+ β1rB (4.133)

∂xC

∂z= 4α1

(u

∂2xC

∂u2+

∂xC

∂u

)+ β1rC (4.134)

∂T

∂z= 4α2

(∂2T

∂u2+

∂T

∂u

)+ β2rB + β3rC (4.135)


z = 0

xB = 0xC = 0T = 1, 0

(4.136)

u = 1, 0

∂xB

∂u= 0

∂xC

∂u= 0

∂T

∂u=

Bi

2(TW − T )

(4.137)

A solução pode ser aproximada por um polinômio de grau N. Seja ui os pontos decolocação e xBi(z), xCi(z) e Ti(z) os valores da conversão para B e C e da temperaturanos pontos de colocação, respectivamente. Assim, os polinômios que passam pelos pontos[xBi(z), xCi(z), Ti(z), ui] podem ser escritos como polinômios interpoladores de Lagrange:

xBN (u, z) =N+1∑

j=1

lj (u) xj (z) (4.138)

xCN (u, z) =N+1∑

j=1

lj (u) xj (z) (4.139)


TN (u, z) =N+1∑

j=1

lj (u) Tj (z) (4.140)

Onde:

lj(u) =pN+1(u)

(u− uj)p(1)N+1(uj)

(4.141)

pN+1(u) =N+1∏

j=1

(u− uj) (4.142)


dxBN (z)du

=N+1∑

j=1

l(1)j (u)xBj (z) (4.143)

d2xBN (z)du2

=N+1∑

j=1

l(2)j (u)xBj (z) (4.144)

dxCN (z)du

=N+1∑

j=1

l(1)j (u)xCj (z) (4.145)

d2xCN (z)du2

=N+1∑

j=1

l(2)j (u)xCj (z) (4.146)

dTN (z)du

=N+1∑

j=1

l(1)j (u)Tj (z) (4.147)

d2TN (z)du2

=N+1∑

j=1

l(2)j (u)Tj (z) (4.148)


dxBN (z)du

∣∣∣∣u=ui

=N+1∑

j=1

l(1)j (ui)xBj (z) =

N+1∑

j=1

Ai,jxBj (z) (4.149)

d2xBN (z)du2

∣∣∣∣u=ui

=N+1∑

j=1

l(2)j (ui)xBj (z) =

N+1∑

j=1

Bi,jxBj (z) (4.150)


dxCN (z)du

∣∣∣∣u=ui

=N+1∑

j=1

l(1)j (ui)xCj (z) =

N+1∑

j=1

Ai,jxCj (z) (4.151)

d2xCN (z)du2

∣∣∣∣u=ui

=N+1∑

j=1

l(2)j (ui)xCj (z) =

N+1∑

j=1

Bi,jxCj (z) (4.152)

dTN (z)du

∣∣∣∣u=ui

=N+1∑

j=1

l(1)j (ui)Tj (z) =

N+1∑

j=1

Ai,jTj (z) (4.153)

d2TN (z)du2

∣∣∣∣u=ui

=N+1∑

j=1

l(2)j (ui)Tj (z) =

N+1∑

j=1

Bi,jTj (z) (4.154)

Fazendo como que as equações diferenciais sejam válidas nos pontos de colocaçãotemos:

dxBi (z)dz

= 4α1

N+1∑

j=1

uiBi,jxBj (z) +N+1∑

j=1

Ai,jxBj (z)

+ β1rB (4.155)

dxCi (z)dz

= 4α1

N+1∑

j=1

uiBi,jxCj (z) +N+1∑

j=1

Ai,jxCj (z)

+ β1rC (4.156)

dTi (z)dz

= 4α2

N+1∑

j=1

uiBi,jTj (z) +N+1∑

j=1

Ai,jTj (z)

+ β2rB + β3rC (4.157)

Ci,j = 4 (uiBi,j + Ai,j) (4.158)

dxBi (z)dz

= α1

N+1∑

j=1

Ci,jxBj (z)

+ β1rB (4.159)

dxCi (z)dz

= α1

N+1∑

j=1

Ci,jxCj (z)

+ β1rC (4.160)

dTi (z)dz

= α2

N+1∑

j=1

Ci,jTj (z)

+ β2rB + β3rC (4.161)

dxBi (z)dz

= α1

N∑

j=1

Ci,jxBj (z) + Ci,N+1xBN+1 (z)

+ β1rB (4.162)


dxCi (z)dz

= α1

N∑

j=1

Ci,jxCj (z) + Ci,N+1xCN+1 (z)

+ β1rC (4.163)

dTi (z)dz

= α2

N∑

j=1

Ci,jTj (z) + Ci,N+1TN+1 (z)

+ β2rB + β3rC (4.164)

Mas, pela condição de contorno:

∂xB

∂u

∣∣∣∣u=1

=N+1∑

j=1

AN+1,jxBj (z) = 0 (4.165)

N∑

j=1

AN+1,jxBj (z) + AN+1,N+1xBN+1 (z) = 0 (4.166)

xBN+1 (z) = −∑N

j=1 AN+1,jxBj (z)AN+1,N+1

(4.167)

∂xC

∂u

∣∣∣∣u=1

=N+1∑

j=1

AN+1,jxCj (z) = 0 (4.168)

N∑

j=1

AN+1,jxCj (z) + AN+1,N+1xCN+1 (z) = 0 (4.169)

xCN+1 (z) = −∑N

j=1 AN+1,jxCj (z)AN+1,N+1

(4.170)

∂T

∂u

∣∣∣∣u=1

=N+1∑

j=1

AN+1,jTj (z) =Bi

2[1− TN+1 (z)] (4.171)

N∑

j=1

AN+1,jTBj (z) + AN+1,N+1TN+1 (z) =Bi

2[1− TN+1 (z)] (4.172)

TN+1 (z) =Bi/2−∑N

j=1 AN+1,jTj (z)Bi/2 + AN+1,N+1

(4.173)

dxBi (z)dz

= α1

N∑

j=1

Ci,jxBj (z)− Ci,N+1∑N

j=1 AN+1,jxBj (z)AN+1,N+1

+ β1rB (4.174)

dxCi (z)dz

= α1

N∑

j=1

Ci,jxCj (z)− Ci,N+1∑N

j=1 AN+1,jxCj (z)AN+1,N+1

+ β1rC (4.175)


dTi (z)dz

= α2

N∑

j=1

Ci,jTj (z) +Ci,N+1

[Bi/2−∑N

j=1 AN+1,j

]Tj (z)

Bi/2 + AN+1,N+1

+ β2rB + β3rC

(4.176)


AN+1,N+1(4.177)

Ei,j = Ci,j − Ci,N+1AN+1,j

Bi/2 + AN+1,N+1(4.178)

vi =Ci,j(Bi/2

Bi/2 + AN+1,N+1(4.179)

dxBi (z)dz

= α1

N∑

j=1

Di,jxBj (z)

+ β1rB (4.180)

dxCi (z)dz

= α1

N∑

j=1

Di,jxCj (z)

+ β1rC (4.181)

dTi (z)dz

= α2

N∑

j=1

Ei,jTj (z) + vi

+ β2rB + β3rC (4.182)

Assim, a integração deste sistema de equações ao longo do eixo z, resulta no per�l detemperatura e de concentração. Nas Figuras 4.15, 4.16, 4.17, 4.18, 4.19, 4.20, 4.21, 4.22,4.23 são apresentados os per�s de conversão de B e de C e de temperatura no reator.

4.2.4 Resolução por Diferenças FinitasAs equações do balanço de massa e de energia:

∂xB

∂z= α1

(∂2xB

∂r2+

1r

∂xB

∂r

)+ β1rB (4.183)

∂xC

∂z= α1

(∂2xC

∂r2+

1r

∂xC

∂r

)+ β1rC (4.184)

∂T

∂z= α2

(∂2T

∂r2+

1r

∂T

∂r

)+ β2rB + β3rC (4.185)

Com as condições de contorno:

z = 0

xB = 0xC = 0T = 1, 0

(4.186)


0 2 4 6 8 10 0,0

0,2

0,4

0,6

0,8

Con

cent

raçã

o

Comprimento

Figura 4.15: Per�l de conversão de Bpara uma pressão parcial de reagente naalimentação igual a 0,003 bar

0 2 4 6 8 10 0,00

0,05

0,10

0,15

0,20

0,25

Con

cent

raçã

o

Comprimento

Figura 4.16: Per�l de conversão de Cpara uma pressão parcial de reagente naalimentação igual a 0,003 bar

0 2 4 6 8 10 1,000

1,005

1,010

1,015

1,020

1,025

1,030

1,035

1,040

Tem

pera

tura

Comprimento

Figura 4.17: Per�l de temperatura parauma pressão parcial de reagente na ali-mentação igual a 0,003 bar

0 2 4 6 8 10 0,0

0,2

0,4

0,6

0,8

Con

cent

raçã

o

Comprimento


0 2 4 6 8 10 0,00

0,05

0,10

0,15

0,20

0,25

0,30

0,35

0,40

Con

cent

raçã

o

Comprimento


0 2 4 6 8 10 1,00

1,02

1,04

1,06

1,08

1,10

1,12

Tem

pera

tura

Comprimento



0 2 4 6 8 10 0,0

0,2

0,4

0,6

0,8

Con

cent

raçã

o

Comprimento


0 2 4 6 8 10 0,0

0,2

0,4

0,6

0,8

1,0

Con

cent

raçã

o

Comprimento


0 2 4 6 8 10 1,0

1,2

1,4

1,6

1,8

2,0

2,2

Tem

pera

tura

Comprimento



r = 0

∂xB

∂r= 0

∂xC

∂r= 0

∂T

∂r= 0

(4.187)

r = 1, 0

∂xB

∂r= 0

∂xC

∂r= 0

∂T

∂r= Bi (TW − T )

(4.188)

Podem ser resolvidas por diferenças �nitas. Assim, fazendo:

dxB

dr=

xBi+1 − xBi−1

2h(4.189)

d2xB

dr2=

xBi+1 − 2xBi + xBi−1

h2(4.190)

dxC

dr=

xCi+1 − xCi−1

2h(4.191)

d2xC

dr2=

xCi+1 − 2xCi + xCi−1

h2(4.192)

dT

dr=

Ti+1 − Ti−1

2h(4.193)

d2T

dr2=

Ti+1 − 2Ti + Ti−1

h2(4.194)

h =1

N + 1(4.195)

ri = ih (4.196)Temos:

dxBi (z)dz

= α1

[xBi+1 (z)− 2xBi (z) + xBi−1 (z)

h2+

1ih

xBi+1 (z)− xBi−1 (z)2h

]+ β1rB

(4.197)

dxCi (z)dz

= α1

[xCi+1 (z)− 2xCi (z) + xCi−1 (z)

h2+

1ih

xCi+1 (z)− xCi−1 (z)2h

]+ β1rC

(4.198)


dTi (z)dz

= α2

[Ti+1 (z)− 2Ti (z) + Ti−1 (z)

h2+

1ih

Ti+1 (z)− Ti−1 (z)2h

]+ β2rB + β3rC

(4.199)

dxBi (z)dz

= α1

{[2i + 1] xBi+1 (z)− 4ixBi (z) + [2i− 1]xBi−1 (z)

2ih2

}+ β1rB (4.200)

dxCi (z)dz

= α1

{[2i + 1]xCi+1 (z)− 4ixCi (z) + [2i− 1]xCi−1 (z)

2ih2

}+ β1rC (4.201)

dTi (z)dz

= α2

{[2i + 1]Ti+1 (z)− 4iTi (z) + [2i− 1]Ti−1 (z)

2ih2

}+ β2rB + β3rC (4.202)

Em r = 0 −→ i = 0

limr→0

1r

∂xB

∂r=

∂2xB

∂r2(4.203)

limr→0

1r

∂xC

∂r=

∂2xC

∂r2(4.204)

limr→0

1r

∂T

∂r=

∂2T

∂r2(4.205)

dxB0 (z)dz

= 2α1∂2xB

∂r2+ β1rB = 2α1

xB1 (z)− 2xB0 (z) + xB−1 (z)h2

+ β1rB (4.206)

dxC0 (z)dz

= 2α1∂2xC

∂r2+ β1rC = 2α1

xC1 (z)− 2xC0 (z) + xC−1 (z)h2

+ β1rC (4.207)

dT0 (z)dz

= 2α2∂2T

∂r2+ β2rA = 2α1

T1 (z)− 2T0 (z) + T−1 (z)h2

+ β2rA (4.208)


dxB

dr=

xB1 (z)− xB−1 (z)2h

= 0 (4.209)

dxC

dr=

xC1 (z)− xC−1 (z)2h

= 0 (4.210)

dT

dr=

T1 (z)− T−1 (z)2h

= 0 (4.211)


Então:

xB−1 (z) = xB+1 (z) (4.212)

xC−1 (z) = xC+1 (z) (4.213)

T−1 (z) = T+1 (z) (4.214)

dxB0 (z)dz

=4α1 [xB1 (z)− xB0 (z)]

h2+ β1rB (4.215)

dxC0 (z)dz

=4α1 [xC1 (z)− xC0 (z)]

h2+ β1rC (4.216)

dT0 (z)dz

=4α1 [T1 (z)− T0 (z)]

h2+ β2rB + β3rC (4.217)

Em r = 1 −→ i = N + 1, pela condição de contorno (Equação 4.187):

dxB

dr=

xBN+2 (z)− xBN (z)2h

= 0 (4.218)

dxC

dr=

xCN+2 (z)− xCN (z)2h

= 0 (4.219)

dT

dr=

TN+2 (z)− TN (z)2h

= Bi [Tw − TN+1 (z)] (4.220)

xBN+2 (z) = xBN (z) (4.221)

xCN+2 (z) = xCN (z) (4.222)

TN+2 (z) = TN (z) + 2hBiTw − 2hBiTN+1 (z) (4.223)

Substituindo:

dxBN+1 (z)dz

= α1

{[2i + 1]xBN (z)− 4ixBN+1 (z) + [2i− 1]xBN (z)

2ih2

}+ β1rB (4.224)

dxCN+1 (z)dz

= α1

{[2i + 1] xCN (z)− 4ixCN+1 (z) + [2i− 1] xCN (z)

2ih2

}+ β1rC (4.225)


dTN+1 (z)dz

=

α2

{[2i + 1] [TN (z) + 2hBiTw − 2hBiTN+1 (z)]− 4iTi (z) + [2i− 1]Ti−1 (z)

2ih2

}+ β2rA

(4.226)

dxBN+1 (z)dz

=2α1 [xBN (z)− xBN+1 (z)]

h2+ β1rB (4.227)

dxCN+1 (z)dz

=2α1 [xCN (z)− xCN+1 (z)]

h2+ β1rC (4.228)

dTN+1 (z)dz

=

α24iTN (z) [2i + 1] 2hBiTw − {[2i + 1] 2hBi + 4i}TN+1 (z)

2ih2+ β2rA

(4.229)

Assim:

dxB (z)dz

=

4α1 [xB1 (z)− xB0 (z)]h2

+ β1rB

α1

{[2i + 1]xBi+1 (z)− 4ixBi (z) + [2i− 1]xBi−1 (z)

2ih2

}+ β1rB

...

2α1 [xBN (z)− xBN+1 (z)]h2

+ β1rB

dxC (z)dz

=

4α1 [xC1 (z)− xC0 (z)]h2

+ β1rC

α1

{[2i + 1] xCi+1 (z)− 4ixCi (z) + [2i− 1]xCi−1 (z)

2ih2

}+ β1rC

...2α1 [xCN (z)− xCN+1 (z)]

h2+ β1rC


dT (z)dz

=

4α1 [T1 (z)− T0 (z)]h2

+ β2rB + β3rc

α2

{[2i + 1]Ti+1 (z)− 4iTi (z) + [2i− 1]Ti−1 (z)

2ih2

}+ β2rb + β3rC

...α2

{4iTN (z) [2i + 1] 2hBiTw − [(2i + 1) 2hBi + 4i]TN+1 (z)

2ih2

}+ β2rB + β3rC

(4.230)

Capítulo 5

Reatores Não Ideais

5.1 Distribuição do Tempo de Residência5.1.1 Introdução

Os reatores tratados até o momento tinham sido modelados como reatores ideais:batelada perfeitamente agitada, reator contínuo perfeitamente agitado e reator tubularcom �uxo empistonado (plug �ow). Em reatores reais nem sempre a mistura vai serperfeita ou o escoamento vai ser empistonado. Podem ocorrer pontos de estagnação,caminhos preferenciais, etc. Assim, por exemplo, em reatores de leito �xo podem ocorrercaminhos preferenciais fazendo com que o �uxo se desvie de empistonado (1, 2). Emreatores de lama as bolhas de gás podem ter tamanhos diferentes. As bolhas maioressairão mais rapidamente. Em tanques agitados podem ocorrer regiões estagnadas oucurto-circuitos (by pass).

5.1.2 Distribuição do Tempo de ResidênciaEm um reator de �uxo empistonado ideal (PFR - plug �ow reactor) todos os átomos do

material que deixa o reator permaneceram dentro deste pelo mesmo período de tempo.Analogamente, em um reator batelada ideal todos os átomos dos materiais dentro doreator �caram por um tempo idêntico neste. O tempo que os átomos passaram dentrodo reator denomina-se tempo de residência.

Os dois exemplos anteriores são os únicos casos em que temos todos os átomos como mesmo tempo de residência. Em todos os outros tipos de reatores teremos uma distri-buição de tempos de residência. A distribuição do tempo de residência de um reator éuma característica da mistura que ocorre neste em nível macroscópico.

5.1.3 Medida da Distribuição do Tempo de ResidênciaA distribuição do tempo de residência é medida injetando um material inerte traçador

na entrada do reator em um tempo t = 0 e medindo sua concentração na corrente a�uenteem função do tempo. Este traçador, além de ser inerte e facilmente detectado, deverter propriedades similares àquelas da mistura reacional e ser completamente solúvel namistura. Ele também não pode adsorver nas paredes ou em outras superfícies no reator.

82 Reatores Não Ideais

Normalmente, são utilizados corantes ou materiais radiativos como traçadores. As duasformas mais comuns de injeção são o pulso e o degrau.

Pulso

Uma quantidade de traçador, N0, é rapidamente injetada na alimentação do reator esua concentração na saída é monitorada em função do tempo. A quantidade de traçadordeixando o reator entre o tempo t e t + ∆t será:

∆N = C(t)v∆t (5.1)

Onde:C(t): concentração do traçador na saída do reatorv: vazão volumétrica do e�uente∆N representa a quantidade de traçador que �cou no reator por um tempo entre t e

t + ∆t. Dividindo pela quantidade de traçador adicionada, N0, obtemos:

∆N

N0=

C(t)vN0

∆t (5.2)

que representa a fração do material com um tempo de residência no reator entre t et + ∆t. Para a injeção em pulso temos:

E(t) =C(t)vN0

(5.3)

Assim:

∆N

N0= E(t)∆t (5.4)

A quantidade E(t) é chamada função distribuição do tempo de residência, que des-creve de forma quantitativa quanto tempo cada elemento do �uido levou dentro do reator.A quantidade total de traçador adicionada, N0, pode ser obtida integrando ∆N :

N0 =∫ ∞

0vC(t)dt (5.5)

Assim, podemos de�nir E(t) como:

E(t) =C(t)∫∞

0 vC(t)dt(5.6)

Então:∫ ∞

0E(t)dt = 1, 0 (5.7)

5.1 Distribuição do Tempo de Residência 83

DegrauA concentração de traçador na saída do reator pode ser relacionada com a entrada

pela integral de convolução:

Csai(t) =∫ t

0Cent

(t− t

′)E(t

′)dt

′ (5.8)

Para injeção em degrau temos:

C0(t) = 0 para t < 0 (5.9)

C0(t) = constante para t ≥ 0 (5.10)Então:

Csai(t) = C0

∫ t

0E(t

′)dt

′ (5.11)

Csai(t)C0

=∫ t

0E(t

′)dt

′= F (t) (5.12)

Podemos obter E(t) derivando expressão acima:

E(t) =ddt

[Csai(t)

C0

]= F

′(t) (5.13)

onde F(t) é a função distribuição do tempo de residência acumulativa.

5.1.4 Características da Distribuição do Tempo de ResidênciaDistribuição do Tempo de Residência Típicos para Diferentes Tipos de Reator

Nas Figuras 5.1, 5.2 e 5.3 são apresentadas as funções distribuição do tempo deresidência para os reatores PFR, CSTR e de �uxo laminar ideais.

Tempo de Residência MédioO tempo de residência médio é igual ao primeiro momento da função distribuição do

tempo de residência:

tm =

∫∞0 tE(t)dt∫∞0 E(t)dt

=∫ ∞

0tE(t)dt (5.14)

Supondo que o reator está inicialmente cheio de moléculas A e em t = 0 passou-sea alimentá-lo com moléculas B. Inicialmente o volume do reator era igual ao volumeocupado pelas moléculas A. Após o início da alimentação de B, o volume de moléculasque deixam o reator no intervalo de tempo dt será igual a vdt. Assim, o volume demoléculas A saindo do reator em um intervalo de tempo dt será:


E(t

)

t

Figura 5.1: Função distribuição do tempo de residência paraum reator PFR ideal (linha cheia) e para um reator PFR real(linha pontilhada)

E(t

)

t

Figura 5.2: Função distribuição do tempo de residência paraum reator CSTR ideal

E(t

)

t

Figura 5.3: Função distribuição do tempo de residência paraum reator de �uxo laminar ideal

5.1 Distribuição do Tempo de Residência 85

dV = [1− F (t)] vdt (5.15)

onde:F (t): representa a fração das moléculas e�uentes que estiveram no reator por um

tempo menor que t.1−F (t): representa a fração das moléculas e�uentes que estiveram no reator por um

tempo maior que t, ou seja, para um dado tempo t é a fração das moléculas ainda noreator.

Para calcular o volume total do reator, integramos de 0 a ∞:

V =∫ ∞

0[1− F (t)] vdt = v

∫ ∞

0[1− F (t)] dt (5.16)

Integrando por partes:∫

xdy = xy −∫

ydx (5.17)

∫ ∞

0[1− F (t)] dt = t [1− F (t)]|∞0 −

∫ ∞

0td [1− F (t)] (5.18)

∫ ∞

0[1− F (t)] dt =

∫ ∞

0tdF (t) (5.19)

V = v

∫ ∞

0tdF (t) (5.20)

Como:

dF (t) = E(t)dt (5.21)

V

v= τ =

∫ ∞

0tE(t)dt (5.22)

ou seja:

τ = tm (5.23)

(tempo espacial) = (tempo de residência médio) (5.24)

Assim, o volume de um reator pode ser facilmente determinado pela equação:

V = vtm = v

∫ ∞

0tE(t)dt (5.25)


Outros Momentos da Distribuição do Tempo de Residência

O segundo momento em relação a média representa a variança, ou seja, a largura dadistribuição do tempo de residência.

σ2 =∫ ∞

0(t− tm)2 E(t)dt (5.26)

O terceiro momento em relação a média representa a simetria da distribuição dotempo de residência (skewness).

s3 =1

σ3/2

∫ ∞

0(t− tm)3 E(t)dt (5.27)

5.1.5 Função Distribuição Tempo de Residência NormalizadaFreqüentemente, a função da distribuição do tempo de residência normalizada é uti-

lizada no lugar de E(t). De�nindo o parâmetro Θ como:

Θ =t

τ(5.28)

A função normalizada pode ser de�nida como:

E(Θ) = τE(t) (5.29)

e portanto:∫ ∞

0E(Θ)dΘ = 1, 0 (5.30)

5.1.6 Distribuição Interna do Tempo de ResidênciaA função E(t) representa a distribuição do tempo de residência das moléculas que

estão saindo do reator, pois E(t)∆t representa o material que está saindo do reator eque passou dentro deste um tempo entre t e t + ∆t. Por isto, E(t) é conhecida como afunção distribuição externa do tempo de residência. A distribuição interna do tempo deresidência, I(α), é uma função tal que I(α)∆α é a fração de material que está dentro doreator e que �cou neste por um período de tempo entre α e α + ∆α. I(α) caracterizao tempo que o material esteve (e ainda está) no reator em um determinado tempo emparticular.

A relação entre a distribuição externa do tempo de residência e a interna pode serdemonstrada pela análise de um reator contínuo em estado estacionário que está cheiode moléculas A. No tempo t = 0 inicia-se a injeção de moléculas B para substituir asmoléculas A. O volume de moléculas dentro do reator que esteve dentro deste por umperíodo de tempo entre α e α + ∆α, é:

dV = V [I(α)] dα (5.31)

5.2 A Distribuição do Tempo de Residência para Reatores Ideais 87

Em t = 0, volume de moléculas B que entram no reator vai ser v0dα. O volume demoléculas que entraram no sistema entre t = 0 e t = dα e que continuam no sistema notempo t = α é dado por:

dV = v0dα [1− F (α)] (5.32)Igualando as Equações 5.32 e 5.31:

V [I(α)] dα = v0dα [1− F (α)] (5.33)

I(α) =v0 [1− F (α)]

V=

[1− F (α)]τ

(5.34)

Pela de�nição de F (α) Equação 5.12:

I(α) =1τ

[1−

∫ α

0E(α)dα

](5.35)

Diferenciando:

dI(α) =−E(α)dα

τ(5.36)

Explicitando E(α):

E(α) = − ddα

[τI(α)] (5.37)

5.2 A Distribuição do Tempo de Residência para ReatoresIdeais

5.2.1 Distribuição do Tempo de Residência para Reatores Batelada ePlug Flow

Neste caso, todos os átomos que saem do reator passaram o mesmo tempo no interiordeste. A função distribuição do tempo de residência pode ser representada por um deltade Dirac:

E(t) = δ(t− τ) (5.38)Propriedades de δ(x):

δ(x) = 0 para x 6= 0 (5.39)

δ(x) = ∞ para x = 0 (5.40)∫ +∞

−∞δ(x) = 1 (5.41)

∫ +∞

−∞g(x)δ(x− τ)dx = g(τ) (5.42)


5.2.2 CSTR SimplesEm um reator CSTR ideal a concentração de qualquer substância na saída do reator

é igual a concentração em seu interior.Assim, fazendo-se um balanço de massa de um traçador injetado como um pulso,

temos:

(entrada)− (saída) = (acúmulo) (5.43)

0− vC =ddt

(V C) (5.44)

dC

C= − v

Vdt =

dt

τ(5.45)

ln [C(t)] = a− t

τ(5.46)

C(t) = a exp(− t

τ

)(5.47)

Condição inicial:

C(0) = C0 (5.48)

C(t) = C0 exp(− t

τ

)(5.49)

Pela de�nição da função distribuição do tempo de residência (Equação 5.6):

E(t) =C0 exp (−t/τ)∫∞

0 C0 exp (−t/τ) dt(5.50)

E(t) =exp (−t/τ)

τ(5.51)

Na forma adimensional:

E(Θ) = exp (−Θ) = e−Θ (5.52)

O tempo de residência médio será:

tm =∫ ∞

0tE(t)dt =

∫ ∞

0texp (−t/τ)

τdt = τ (5.53)

O segundo momento em torno da média, �ca:

σ2 =∫ ∞

0(t− τ)2

exp (−t/τ)τ

dt = τ2 (5.54)

5.2 A Distribuição do Tempo de Residência para Reatores Ideais 89

5.2.3 Reator de Fluxo LaminarNeste caso, o per�l de velocidade é parabólico:

U = Umáx

[1−

( r

R

)2]

= 2U

[1−

( r

R

)2]

=2v0

πR2

[1−

( r

R

)2]

(5.55)

O tempo necessário para a passagem de um elemento de �uido em uma posição radialigual a r é:

t(r) =L

U(r)=

πR2L

v0

1

2[1− (r/R)2

] (5.56)

t(r) =V

v0

1

2[1− (r/R)2

] (5.57)

t(r) =τ

2[1− (r/R)2

] (5.58)

A fração de �uido passando entre r e r + dr é:

dv

v0=

U(r)2πrdr

v0=

L

t(r)

(2πrdr

v0

)(5.59)

Diferenciando a Equação 5.58:

dt(r) =τ

2R2

2rdr[1− (r/R)2

]2 =4

τR2

τ/2[1− (r/R)2

]

2

rdr (5.60)

Substituindo t dado pela Equação 5.58

dt(r) =4t2

τR2rdr (5.61)

Explicitando dr

dr =τR2

4t2rdt (5.62)

Combinando as Equações 5.58 e 5.62, temos:

dv

v0=

L

t(r)

(2πr

v0

)τR2

4t2rdt (5.63)

dv

v0=

πR2L

v0

τ

2t3dt =

V

v0

τ

2t3dt =

τ2

2t3dt (5.64)

dv

v0=

τ2

2t3dt (5.65)

O tempo mínimo para que um elemento de volume do �uido passe pelo reator será:


t =L

Umáx=

L

2U

πR2

πR2=

V

2v0=

τ

2(5.66)

Assim, a função distribuição do tempo de residência completa será dada por:

E(t) = 0 para t <τ

2(5.67)(5.68)

E(t) =τ2

2t3para t ≥ τ

2(5.69)

A função distribuição acumulativa do tempo de residência para t ≥ τ/2 é:

F (t) =∫ t

τ2

τ2

2t3dt =

τ2

2

∫ t

τ2

dt

t3= 1− τ2

4t2(5.70)

O tempo de residência médio será:

tm =∫ ∞

τ2

tE(t)dt =∫ ∞

τ2

dt

t2= τ (5.71)

A forma adimensionalizada da função distribuição do tempo de residência é:

E(Θ) = 0 para Θ <12

(5.72)(5.73)

E(Θ) =1

2Θ3para Θ ≥ 1

2(5.74)

5.2.4 Função Distribuição do Tempo de Residência para Séries dePFR/CSTR Ideais

Em alguns tanques agitados existe uma zona altamente agitada próxima ao impelidorque pode ser modelada como um CSTR de mistura perfeita. Dependendo da localizaçãoda entrada e da saída, a mistura reacional pode seguir um caminho tortuoso antes oudepois da zona de agitação perfeita. Este caminho tortuoso pode ser modelado comoum reator PFR ideal. Assim, este tipo de reator pode ser modelado como uma série deCSTR e PFR.

Seja τs o tempo de residência médio no CSTR e τp o tempo de residência médio noPFR. Se um pulso de um traçador for injetado na entrada do sistema, no primeiro reator(CSTR), a concentração na saída do CSTR será dada por:

C(t) = C0 exp(− t

τs

)(5.75)

A saída do segundo reator (PFR) será atrasada em um tempo igual a τp. Assim, afunção distribuição do tempo de residência deste sistema será:

5.3 Modelagem de um Reator com a Função Distribuição do Tempo deResidência 91

E(t) = 0 para t < τp (5.76)(5.77)

E(t) =exp

(− t−τp

τs

)

τspara t ≥ τp (5.78)

Colocando primeiro do PFR e depois o CSTR obtemos o mesmo resultado, ou seja,a distribuição do tempo de residência independe da posição dos reatores.

Para reações de primeira ordem, a conversão obtida será a mesma, não importa aordem dos reatores. No entanto, para reações de segunda ordem, vamos obter conversõesdependentes da ordem com que os reatores estão. Assim, a função distribuição do tempode residência não fornece a descrição completa de um reator real, necessitamos de maisinformações para podermos predizer a conversão nestes casos.

5.3 Modelagem de um Reator com a Função Distribuiçãodo Tempo de Residência

Existem situações em que o �uido no reator não é nem perfeitamente misturado nemempistonado. Assim, surge uma pergunta: como podemos utilizar a distribuição do tempode residência para predizer a conversão de um reator real?

A função distribuição do tempo de residência nos diz quanto tempo cada elementode �uido permaneceu dentro do reator (macro-mistura), mas não nos diz nada a respeitoda transferência de massa entre os elementos de �uido (micro-mistura). Para reações deprimeira ordem temos:

dX

dt= k (1−X) (5.79)

ou seja, a conversão é independente da concentração do reagente. Por isto podemospredizer a conversão de uma reação de primeira ordem em um reator real utilizandounicamente a função distribuição do tempo de residência. Para uma reação de segundaordem temos:

dX

dt= kCA0 (1−X)2 (5.80)

ou seja, a conversão depende da concentração do reagente. Por isto não podemospredizer a conversão de uma reação de segunda ordem em um reator real unicamentecom a função distribuição do tempo de residência. Para tal, devemos conhecer o graude mistura em nível microscópico, ou seja, o grau com que os elementos de �uido queapresentam diferentes tempos de residência transferem massa entre si.

5.4 Modelos com Zero ParâmetrosExistem dois extremos para a mistura em nível microscópico: (1) todas as molécu-

las apresentando o mesmo tempo de residência permanecem juntas a medida em que


atravessam o reator e não se misturam até que saiam deste (segregação completa). (2)moléculas apresentando tempos de residência diferentes são misturadas completamentetão logo elas entram no reator (mistura microscópica completa). Para uma dada dis-tribuição do tempo de residência, estes dois extremos de mistura em nível molecularresultarão nos limites superior e inferior de conversão em um reator não-ideal. Para ciné-tica com ordem de reação acima de um, o modelo de segregação irá predizer a conversãomais alta e o de mistura completa a mais baixa. O inverso ocorrerá para cinéticas comordens de reação inferiores a um.

5.4.1 Modelo de SegregaçãoNeste modelo, assume-se que elementos de �uido com diferentes tempos de residência

não trocaram massa entre si. Desta forma, seria como se juntássemos todas as moléculascom o mesmo tempo de residência em um �glóbulo�. Estes �glóbulos� seriam então mistu-rados na saída. Assim, para determinar a conversão média na saída, devemos fazer umamédia das conversões dos vários �glóbulos� que estão saindo. Seja um reator PFR comvárias saídas laterais, com o �uido saindo por estas de forma a reproduzir a distribuiçãodo tempo de residência.

Conversão médiados �glóbulos�

que �caram de ta t + dt no reator

=

Conversãoalcançada porum �glóbulo�

depois de passarum tempo t no

reator

×

Fração dos�glóbulos� que�caram por umtempo entre t et + dt no reator

(5.81)

dX = X(t)× E(t)dt (5.82)

A conversão média vai ser obtida pela integração da equação acima:

X =∫ ∞

0X(t)× E(t)dt (5.83)

Cada glóbulo vai se comportar como um reator batelada individual. O tempo dereação em cada um destes glóbulos é igual ao tempo em que este �cou no reator. Então,a conversão X(t) pode ser calculada pela expressão do balanço de massa em um reatorbatelada. Assim, para uma reação de primeira ordem temos:

−dNA

dt= −rAV = −kCAV (5.84)

dX

dt= k (1−X) (5.85)

Integrando:

X(t) = 1− exp (−kt) (5.86)

Assim:

5.4 Modelos com Zero Parâmetros 93

X =∫ ∞

0[1− exp (−kt)]× E(t)dt (5.87)

X = 1−∫ ∞

0exp (−kt)×E(t)dt (5.88)

5.4.2 Modelo de Máxima Antecipação da MisturaNo modelo anterior foi analisado um extremo no grau de mistura a nível microscópico,

ou seja, a mistura mínima. Neste caso, vamos analisar o outro extremo.Seja um reator PFR com entradas laterais, com o �uido entrando por estas de tal

maneira que reproduza a distribuição do tempo de residência.Assim que entra no reator o �uido é completamente misturado radialmente, mas não

axialmente. Seja λ o tempo que o �uido leva para se mover da entrada até a saída, ouseja, a expectativa de vida do �uido no reator.

A vazão volumétrica do �uido alimentado lateralmente ao reator no intervalo detempo entre λ e λ + ∆λ é igual a v0E(λ)∆λ. A vazão volumétrica do reator no ponto λé:

v(λ) = v0

∫ ∞

λE(λ)dλ = v0 [1− F (λ)] (5.89)

O volume de �uido com expectativa de vida no reator entre λ e λ + ∆λ é igual a:

|∆v| = v0 [1− F (λ)] ∆λ (5.90)

A taxa de geração da substância A neste volume é:

rA |∆V | = rAv0 [1− F (λ)]∆λ (5.91)

O balanço de massa de A entre λ e λ + ∆λ é dado por:

Quantidadede A queentra noreator emλ + dλ

+

Quantidadede A queentra noreator

lateralmente

−

Quantidadede A que saido reator em

λ

+

Quantidadede A geradapor reaçãoquímica nointerior do

volume entreλ e λ + dλ

= 0 (5.92)

v0 [1− F (λ)]CA|λ+∆λ − v0 [1− F (λ)]CA|λ + v0CA0E(λ)∆λ + rAv0 [1− F (λ)]∆λ = 0(5.93)

Dividindo por v0∆λ e tomando o limite para ∆λ → 0:

d [1− F (λ)]CA

dλ+ CA0E(λ) + rA [1− F (λ)] = 0 (5.94)


dCA

dλ= −rA + (CA − CA0)

E(λ)1− F (λ)

(5.95)

Em termos de conversão:

dX

dλ=

rA

CA0+

E(λ)1− F (λ)

X (5.96)

Condição de contorno:

CA = CA0 para λ →∞ (5.97)

5.5 Modelos com um Parâmetro5.5.1 Modelo de Tanques em Série

Neste modelo temos como parâmetro ajustável o número de reatores CSTR ideaisem série. Este número é ajustado para resultar na mesma distribuição do tempo deresidência do reator real.

Inicialmente vamos analisar o comportamento da saída de três reatores CSTR deigual volume em série para a injeção de um pulso de traçador no primeiro reator.

A fração de traçador deixando o terceiro reator que �cou no sistema entre t e t + ∆té dado por:

E(t)∆t =vC3(t)∆t

N0=

C3(t)∫∞0 C3(t)dt

∆t (5.98)

Então:

E(t) =C3(t)∫∞

0 C3(t)dt(5.99)

Para achar a concentração de traçador na saída do terceiro reator em função do tempocomeçamos pelo balanço de massa no primeiro reator:

V1dC1

dt= −vC1 (5.100)

Integrando:

C1 = C0 exp(− vt

V1

)= C0 exp

(− t

τ1

)(5.101)

Como a vazão volumétrica é constante (v = v0) e os volumes dos reatores são iguais(V1 = V2 = V3 = Vi), os tempos espaciais serão também iguais (τ1 = τ2 = τ3 = τi), ondeτi é o tempo espacial de um reator simples que formam o sistema e não do sistema comoum todo.

No segundo reator temos:

VidC2

dt= vC1 − vC2 (5.102)

5.5 Modelos com um Parâmetro 95

VidC2

dt= vC0 exp

(− t

τi

)− vC2 (5.103)

dC2

dt+

C2

τi=

C0

τiexp

(− t

τi

)(5.104)

Esta equação diferencial é resolvida com a utilização do fator de integração et/τi ecom a condição inicial C2 = 0 em t = 0:

C2 =C0t

τiexp

(−t

τi

)(5.105)

No terceiro reator temos:

VidC3

dt= vC2 − vC3 (5.106)

VidC3

dt= v

C0t

τiexp

(−t

τi

)− vC3 (5.107)

dC3

dt+

C3

τi=

C0t

τiexp

(−t

τi

)(5.108)

Esta equação diferencial é resolvida com a utilização do fator de integração et/τi ecom a condição inicial C2 = 0 em t = 0:

C3 =C0t

2

2τ2i

exp(−t

τi

)(5.109)

Assim:

E(t) =C3(t)∫∞

0 C3(t)dt(5.110)

E(t) =t2

2τ3i

exp(−t

τi

)(5.111)

Generalizando para n reatores CSTR em série, temos:

E(t) =tn−1

(n− 1)!τni

exp(−t

τi

)(5.112)

Como o volume total do sistema é nVi então τi = τ/n:

E(t) =n (nΘ)n−1

(n− 1)!τexp (−nΘ) (5.113)

τE(t) = E(Θ) =n (nΘ)n−1

(n− 1)!exp (−nΘ) (5.114)

O número de tanques em série, n, pode ser calculado pela variança adimensional doexperimento com o traçador:


σ2Θ =

σ2

τ2=

∫ ∞

0(Θ− 1)2 E(Θ)dΘ =

∫ ∞

0

(Θ2 − 2Θ + 1

)E(Θ)dΘ (5.115)

σ2Θ =

∫ ∞

0Θ2E(Θ)dΘ− 1 =

∫ ∞

0Θ2 n (nΘ)n−1

(n− 1)!exp (−nΘ)dΘ− 1 (5.116)

σ2Θ =

1n

(5.117)

n =1

σ2Θ

=τ2

σ2(5.118)

Que representa o número de tanques em série necessários para modelar a distribuiçãodo tempo de residência de um reator real.

Caso a reação apresente cinética de primeira ordem, podemos utilizar a equação:

X = 1− 1(1 + τik)n (5.119)

com:

τi =V

v0n(5.120)

para predizer a conversão. Para outras cinéticas devemos fazer o balanço em cadareator.

5.5.2 Modelo de DispersãoNeste modelo, existe uma dispersão axial do material, o qual é governada por uma

analogia com a lei de Fick, juntamente com o �uxo. Assim, além do �uxo bruto, uAc,cada componente da mistura é transportado através de uma seção transversal a uma taxaigual a −DaAc

dCdz devido à difusão turbulenta e molecular.

A vazão molar do traçador por convecção e por dispersão será:

FT = −DaAcdC

dz+ uAcCT (5.121)

O balanço molar do traçador resulta em:

−∂FT

∂z= Ac

∂CT

∂t(5.122)

Da∂2CT

∂z2− ∂

∂z(uCT ) =

∂CT

∂t(5.123)

Fazendo:

5.5 Modelos com um Parâmetro 97

E (� )00 . 511 . 522 . 533 . 54

�0 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 4Figura 5.4: Distribuição do tempo de residência para um conjunto de n reatores CSTRideais em série.

Ψ =CT

CT0(5.124)

(5.125)λ =

z

L(5.126)(5.127)

Θ =ut

L(5.128)

Para injeção em pulso, CT0 é de�nido como o número de moles de traçador injetadono sistema pelo volume do reator.

A equação diferencial então �ca:

Da

uL

∂2Ψ∂λ2

− ∂Ψ∂λ

=∂Ψ∂Θ

(5.129)

onde uLDe

é uma forma do número de Peclet:

Pe =uL

De=

velocidade de transporte por convecçãovelocidade de transporte por dispersão (5.130)

1Pea

∂2Ψ∂λ2

− ∂Ψ∂λ

=∂Ψ∂Θ

(5.131)

Condições de contorno:Existem dois casos a considerar:


• Condições de contorno para o reator fechado: não existe dispersão ou variação deradial de concentração na entrada ou na saída do reator

• Condições de contorno para o reator aberto: existe dispersão na entrada e na saídado reator

Condição de contorno fechado-fechado

Neste caso, temos �uxo empistonado (sem dispersão) um pouco antes da entrada(z = 0−) e um pouco depois da saída (z = L+). Porém, em z = 0 e z = L temosdispersão e reação. Então:

FT

(0−, t

)= FT

(0+, t

)(5.132)

uAcCT

(0−, t

)= −AcDa

∂CT

∂z

∣∣∣∣z=0+

+ uAcCT

(0+, t

)(5.133)

Como:

CT

(0−, t

)= CT0 (5.134)

CT0 = −Da

u

∂CT

∂z

∣∣∣∣z=0+

+ CT

(0+, t

)(5.135)

Na seção de saída a concentração de traçador é continua e não existe gradiente deconcentração.

em: z = L

CT (L−) = CT (L+)

∂CT

∂z= 0

(5.136)

Condições de contorno na forma adimensional:

em: λ = 0 − 1Pea

∂Ψ∂λ

+ Ψ =CT (0−, t)

CT0= 1, 0

em: λ = 1∂Ψ∂λ

= 0

(5.137)

Esta equação diferencial com estas condições de contorno não possui solução analítica,devendo ser resolvida numericamente.

O tempo de residência médio para este sistema será igual ao tempo espacial e avariança adimensional será dada por:

σ2

t2m=

2Pea

− 2Pea

[1− exp(−Pea)] (5.138)

5.6 Modelos com 2 Parâmetros 99

Condições de contorno aberto-abertoQuando um traçador é injetado em um leito empacotado em uma posição de mais

de 2 ou 3 diâmetros de partícula antes da entrada e medido alguma distância depois dasaída, as condições de contorno aberto-aberto se aplicam.

FT

(0−, t

)= FT

(0+, t

)(5.139)

−Da∂CT

∂z

∣∣∣∣z=0−

+ uCT

(0−, t

)= −Da

∂CT

∂z

∣∣∣∣z=0+

+ uCT

(0+, t

)(5.140)

CT

(0−, t

)= CT

(0+, t

)(5.141)

Na saída:

CT

(L−, t

)= CT

(L+, t

)(5.142)

Utilizando estas condições de contorno, a equação diferencial pode ser resolvida ana-liticamente, resultando em:

Ψ(1, 0) =CT (L, t)

CT0=

1√πΘ3/Pe

exp

[− (1−Θ)2

4Θ/Pe

](5.143)

Neste caso, o tempo de residência médio será diferente do tempo espacial:

tm =(

1 +2

Pe

)τ (5.144)

A variança será:

σ2

t2m=

2Pe

+8

Pe2(5.145)

5.6 Modelos com 2 ParâmetrosUm reator real pode ser modelado pela combinação de reatores ideais. Neste item

vamos �xar o número de parâmetros em dois.

5.6.1 CSTR Real Modelado com um Volume de TrocaUm reator que apresente uma região altamente agitada e outra menos agitada pode

ser modelado como dois reatores CSTR que trocam material entre si (Figura 5.5).Seja β a fração da vazão volumétrica total que é trocada entre os reatores CSTR

ideais 1 e 2:

v1 = βv2 (5.146)

Seja α a fração do volume total ocupada pela região altamente agitada:


0

CA1

v0

v0

v1

v1

CA2V1

V2

v

Figura 5.5: Modelo de um reator com um volume de troca

V1 = αV (5.147)

então:

V2 = (1− α) V (5.148)

O balanço de massa no reator 2 �ca:

βv0CA1 − βv0CA2 − (1− α) V kCA2 = 0 (5.149)

CA2 =βCA1

β + (1− α) τk(5.150)

Balanço de massa no reator 1:

v0CA0 − βv0CA2 − v0CA1 − βv0CA1 − kCA1αV = 0 (5.151)

v0CA0 − βv0βCA1

β + (1− α) τk− v0CA1 − βv0CA1 − kCA1αV = 0 (5.152)

CA1 =CA0

1 + β + ατk −{

β2

β+(1−α)τk

} (5.153)

Em termos de conversão:

X = 1− CA1

CA0=

(β + ατk) [β + (1− α) τk]− β2

(1 + β + ατk) [β + (1− α) τk]− β2(5.154)

Para valores altos de τk:

5.6 Modelos com 2 Parâmetros 101

CA1 ' CA0

1 + β + ατk − 0' CA0

ατk(5.155)

Para valores baixos de τk:

CA1 ' CA0

1 + β + 0− β2/β' CA0 (5.156)

Os parâmetros α e β podem ser obtidos através de experimentos com traçadores


Apêndice A

Cálculo

A.1 Integrais Inde�nidas∫

adx = ax (A.1)

∫af(x)dx = a

∫f(x)dx (A.2)

∫(u± v ± w ± · · · )dx =

∫udx±

∫vdx±

∫wdx± · · · (A.3)

∫udv = uv −

∫vdu (A.4)

∫undu =

un+1

n + 1, n 6= −1 (A.5)

∫du

u= ln |u| (A.6)

∫eudu = eu (A.7)

∫audu =

au

ln a(A.8)

∫sen udu = − cosu (A.9)

∫cosudu = sen u (A.10)

∫tanudu = − ln cosu (A.11)

104 Cálculo

∫cotudu = ln sen u (A.12)

∫sec2 udu = tanu (A.13)

∫sen 2udu =

12

(u− sen u cosu) (A.14)

∫cos2 udu =

12

(u + sen u cosu) (A.15)

∫senh udu = coshu (A.16)

∫coshudu = senh u (A.17)

∫tanhudu = ln coshu (A.18)

∫xeaxdx =

eax

a

(x− 1

a

)(A.19)

∫x2eaxdx =

eax

a

(x2 − 2x

a+

2a2

)(A.20)

∫xneaxdx =

eax

a

[xn − nxn−1

a+

n(n− 1)xn−2

a2− · · · (−1)nn!

an

](A.21)

∫xJ0(x)dx = xJ1(x) (A.22)

∫x2J0(x)dx = x2J1(x) + xJ0(x)−

∫J0(x)dx (A.23)

∫J1(x)dx = −J0(x) (A.24)

∫xJ1(x)dx = −xJ0(x) +

∫J0(x)dx (A.25)

∫ ∞

0e−ax2

dx =12

√π

a(A.26)

A.2 Equações Diferenciais 105

A.2 Equações DiferenciaisEquação Diferencial Solução

Separação de variáveisf1(x)g1(y)dx + f2(x)g2(y)dy = 0

∫ f1(x)f2(x)

dx +∫ g2(y)

g1(y)dy = c

Equação linear de primeira ordemdy

dx+ P (x)y = Q(x) ye

RPdx =

∫Qe

RPdxdx + c

Equação de Bernoullidy

dx+ P (x)y = Q(x)yn y1−ne(1−n)

RPdx = (1− n)

∫Qe(1−n)

RPdxdx + c

Equação ExataM(x, y)dx + N(x, y)dy = 0

onde ∂M/∂y = ∂N/∂x∫

M∂x +∫ (

N − ∂

∂y

∫M∂x

)dy = c

Equação Homogêneady

dx= F

(y

x

)ln x =

∫ d( y

x

)

F ( yx)− ( y

x)Equação linear homogênea de 2a ordem

d2y

dx2+ a

dy

dx+ by = 0 Sejam m1 e m2 as raízes de m2 + am + b = 0.

Então existem 3 casos:Caso 1: m1 e m2 reais e distintas

y = c1em1x + c2e

m2x

Caso 2: m1 e m2 reais e iguaisy = c1e

m1x + c2xem1x

Caso 3: m1 = p + qı e m2 = p− qıy = epx [c1 cos(qx) + c2sen (qx)]onde p = −a/2 e q =

√b− a2/4

Equação de Bessel

x2 d2y

dx2+ x

dy

dx+

(λ2x2 − n2

)y = 0 y = c1Jn(λx) + Yn(λx)

Equação de Bessel modi�cada

x2 d2y

dx2+ x

dy

dx− (

λ2x2 + n2)y = 0 y = c1In(λx) + Kn(λx)

Equação de Legendre

(1− x2)d2y

dx2− 2x

dy

dx+ n (n + 1) y = 0 y = c1Pn(λx) + Qn(λx)

A.2.1 Integração NuméricaRegra do trapézio

∫ X1

X0

f(X)dX =h

2[f(X0) + f(X1)] (A.27)

onde:

106 Cálculo

h = X1 −X0 (A.28)

Regra de Simpson (1/3)

∫ X2

X0

f(X)dX =h

3[f(X0) + 4f(X1) + f(X2)] (A.29)

onde:

h =X2 −X0

2(A.30)

Regra de Simpson (3/8)

∫ X3

X0

f(X)dX =3h

8[f(X0) + 3f(X1) + 3f(X2) + f(X3)] (A.31)

onde:

h =X3 −X0

3(A.32)

Regra de Simpson (5 pontos)

∫ X4

X0

f(X)dX =h

3[f(X0) + 4f(X1) + 2f(X2) + 4f(X3) + f(X4)] (A.33)

onde:

h =X4 −X0

4(A.34)

Regra de Simpson (N+1 pontos - N/3 é um inteiro)

∫ XN

X0f(X)dX =

3h

8[f(X0) + 3f(X1) + 3f(X2) + 2f(X3) + 3f(X4) + 3f(X5) + 2f(X6) + · · ·

+3f(XN−1) + f(XN )](A.35)

onde:

h =XN −X0

N(A.36)

A.2 Equações Diferenciais 107

Regra de Simpson (N+1 pontos - N/2 é um inteiro)

∫ XN

X0f(X)dX =

h

3[f(X0) + 4f(X1) + 2f(X2) + 4f(X3) + 2f(X4) + · · ·

+4f(XN−1) + f(XN )](A.37)

onde:

h =XN −X0

N(A.38)

108 Cálculo

Referências Bibliográ�cas

1 FOGLER, H. Elements of Chemical Reaction Engineering. [S.l.]: Prentice Hall, 1992.

2 FROMENT, G.; BISCHOFF, K. Chemical Reactor Analysis and Design. [S.l.]: JohnWiley & Sons, 1990.

Environment

transferencia de massa