UNIVERSIDADE FEDERAL DO PARÁ

INSTITUTO DE TECNOLOGIA

FACULDADE DE ENGENHARIA QUÍMICA

ANA LUIZA DE BARROS SOUZA CAMPOS

MATRÍCULA: 201807540011

3º RELATÓRIO: MODELAGEM E SIMULAÇÃO

BELÉM

2022
 ANA LUIZA DE BARROS SOUZA CAMPOS

MATRÍCULA: 201807540011

3º RELATÓRIO: MODELAGEM E SIMULAÇÃO

Relatório apresentado como 3ª

avaliação para a disciplina de
modelagem e simulação no curso
de engenharia química da
Universidade Federal Do Pará.

BELÉM

2022
 SUMÁRIO

1. INTRODUÇÃO
2. METODOLOGIA
2.1. Questão 3ª avaliação
2.1.1. Balanço de massa e energia
2.1.2. Código Rotina bvp4c
2.1.3. Código Diferenças Finitas
3. RESULTADOS
3.1. Rotina bvp4c
3.1.1. Alternativa i
3.1.2. Alternativa ii
3.1.3. Alternativa iii
3.2. Método diferenças finitas
3.2.1. Alternativa i
4. CONCLUSÃO
REFERENCIAS
1. INTRODUÇÃO

Os reatores químicos são equipamentos onde acontecem reações em escala

industrial para transformação de matérias primas em produtos, há dois dos tipos de
reatores ideais mais comuns cujo são denominados: batelada (BR) e tubular (PFR)
(FOGLER, 2012). Em reatores ideais do tipo Plug-Flow Reactor - PFR, o fluido
reagente se move através da tubulação como um pistão paralelo ao eixo do tubo,
sendo a composição uma função do volume ou do eixo axial do equipamento,
caracteriza-se por um escoamento ordenado dos elementos, de modo que, não ocorre
difusão ao longo do percurso e a velocidade entre dois elementos do fluido não muda
(SANTOS et al, 2019).

Porém, normalmente um reator real opera com desvios do comportamento ideal

devido a presença de zonas de turbulência, formação de regiões de by-pass, regiões
mortas e caminhos preferenciais no interior do reator, onde tais desvios são
responsáveis por alteração tanto na qualidade do produto quanto na sua conversão.
O modo de alimentação também é um fator que pode interferir no comportamento de
um reator, modificando o escoamento em seu interior (ISHIDA,2017).

Alimentações transientes, a título de exemplo, apresentam potencial para

aumentar a turbulência do escoamento, criando zonas de vorticidade que são
responsáveis pelo aumento do transporte de massa e energia e pela promoção de um
melhor contato entre as moléculas de reagente. É interessante saber se as
modificações geradas pelas diferentes alimentações acarretam diferentes conversões
ao final do processo (ISHIDA,2017).

Uma forma comum de analisar um reator para obter resultados sem a realização
de testes em uma planta real, consiste na utilização da modelagem e simulação, além
de ser uma estratégia de otimização, isto é maximizar a geração de produtos
desejados ou de maior valor agregado, gerando o melhor rendimento através dos
menores custos (STEPHANOPOULOS, 1984). A modelagem e simulação de
processos é uma ferramenta bastante utilizada na engenharia, numa vez que pode
prever o comportamento de sistemas sem a necessidade da montagem de uma
estrutura física de processo, reduzindo os custos do estudo. A fluidodinâmica
computacional realiza a simulação numérica dos processos físicos e químicos
envolvidos no escoamento reativo através da implementação de modelos
matemáticos que sejam representativos para o sistema (SILVA & SOARES, 2011).

Para este tipo de modelagem, pode-se destacar o método de diferenças finitas,

um método de discretizarão de equações diferenciais. Isto significa que ele transforma
uma função contínua em uma representação discreta (pontos), de modo que, quanto
mais elementos forem utilizados, maior será a precisão do método, porém o gasto
computacional também irá aumentar. A estratégia deste método consiste em buscar
equações algébricas que aproximem a solução em cada ponto i, na qual a origem
destas aproximações sempre é uma aproximação em série de Taylor em torno de
cada ponto (FONTANA,2018). Além disso, há rotinas para aplicar nestes tipos de
problemas como a bvp4c, do Matlab que resolve problemas de valor de contorno para
EDOs por colocação (Software Matlab).

O objetivo deste trabalho é analisar a adimensionalização da concentração e da

coordenada de distância, bem como comportamento de parâmetros com base na
modelagem de um reator PFR, utilizando duas metodologias distintas de solução.

2. METODOLOGIA

A metodologia da solução utilizada reúne códigos utilizados para a resolução do

problema aplicados no software MATLAB R2022a. Assim, segue abaixo a descrição
do código, junto com este trabalho foi encaminhado em uma pasta com os arquivos
em formato .m (Arquivos do Matlab) para facilitar a conferência da resolução.

2.1. Questão 3ª avaliação (PFR)

Admitindo-se as adimensionalizações para coordenada l; w = l/L e para a

concentração CA, xA = CA/CAf, e para a temperatura T, y = T/Tf, deduza a equação
para concentração adimensional x e para a temperatura admissional y em função da
coordenada adimensional w para uma reação de primeira ordem n= 1 e segunda
ordem n = 2, onde os parâmetros número de Peclet e número de Damköhler são 𝑃𝑒 =
𝐿.𝑣 𝑘.𝐿.𝐶𝑎𝑓𝑛−1
e 𝐷𝑎 = respectivamente.
𝐷𝑖 𝑣
 Figura 1: Representação de um retor PFR

Fonte: Autora.

i) Considerando um reator isotérmico calcule o perfil de concentração do reator

em função da coordenada adimensional, realize um estudo completo para o
reator para Da = 1,5, v = 3x10-4m/s (velocidade), Di = 10-5m2/s (coeficiente
de difusão) e L = 1m (comprimento do reator), calcule o perfil e a conversão na
saída do reator, para ordens de reação n = 1 e n = 2. Para solução do use o
método das Diferenças Finitas dividindo o domínio em 10 pontos e a rotina
bvp4c do Matlab.
ii) Para uma reação de ordem n =1 assumindo o valor de Da = 3,0, analise o
efeito de PeM variando de 0 a 200. Utilize a rotina bvp4c do Matlab.
iii) Para uma reação de ordem n =1 assumindo o valor de PeM = 20, analise o
efeito de Da variando de 0 a 20. Utilize a rotina bvp4c do Matlab.

É importante ressaltar que o número de Péclet (Pe) Em problemas de

transferência de massa, é o produto do número de Reynolds , que descreve o regime
de fluxo , e o número de Schmidt , usado para caracterizar fluxos de fluidos nos quais
há processos simultâneos de convecção de momento e difusão de massa. No âmbito
da engenharia esse valor numérico é frequentemente muito grande, de modo que
quando se modela certas situações com números de Péclet altos, modelos
computacionais mais simples podem ser adotados. Já os números de Damköhler( Da)
são números adimensionais usados em engenharia química para relacionar a escala
de tempo da reação química ( taxa de reação ) à taxa de fenômenos de transporte
que ocorrem em um sistema (CONNOR, N. .THERMAL-ENGINEERING).

2.1.1. Balanço de massa e energia

Primeiramente foi necessário realizar o balanço de massa e energia em termos

dos parâmetros de Peclet e Damköhler no reator PFR, bem como suas variáveis
realizando as devidas adimensionalizações para a possível resolução da
problemática.

𝑑𝑁𝑖 𝑑2 𝐶𝑖
= −𝐷𝑖 (Equação 1)
𝑑𝑙 𝑑𝑙 2

𝑑𝐻𝑖 𝑑𝑛𝑖 𝑑𝑞
∑𝑖 𝑛𝑖 + ∑𝑖 𝐻𝑖 + 𝐴𝑡 𝑑𝑙 = 𝑄′ (Equação 2)
𝑑𝑙 𝑑𝑙

1 𝑑2 𝑦 𝑑𝑦 −𝛾
⁄𝑦
− 𝑑𝜔 + 𝛽. 𝐷𝑎. 𝑒 . 𝑋𝑎 = −𝑄 (Equação 3)
𝑃𝑒𝐻 𝑑𝜔2

Adimensionalizando o sistema temos que para uma ordem n de reação

(Equação 4) e respectivamente as condições de contorno:

1 𝑑2 𝑥 𝑑𝑥
∙ − 𝑑𝜔 − 𝐷𝑎 ∙ 𝑥 𝑛 = 0 (Equação 4)
𝑃𝑒 𝑑𝜔2

𝑑𝑥
𝜔=1→ =0 (Equação 5)
𝑑𝜔

1 𝑑𝑥
𝜔 = 0 → 𝑃𝑒 ∙ 𝑑𝜔 = 𝑥 − 1 (Equação 6)

2.1.2. Código Rotina bvp4c

 Aplicação i

function dydx = bvpfcn2(x,y) %#ok<*INUSL>

global L Pe Da %#ok<*GVMIS>
dydx = zeros (2,1); %#ok<PREALL>

dydx = [y(2)
(y(2) + Da*y(1))/(1/Pe)];

function res = bcfcn2(ya,yb)

global Pe %#ok<GVMIS>

res = [(1/Pe)*ya(2)-ya(1)+1

yb(2)];
end

function g = guess2(x) %#ok<*INUSD>

g = [0
1];
end

clear all %#ok<CLALL>

clc

global L k Pe Da Di v %#ok<GVMIS>
L = 1; k = 4.5e-4; v=3e-4; Di= 1e-5;
Pe=L*v/Di;
Da= 1.5 ; %#ok<*NOPTS>

deltaL = linspace(0,L,10);
solinit = bvpinit(deltaL, @guess2);
sol = bvp4c(@bvpfcn2, @bcfcn2, solinit);
X = sol.x; T = sol.y;
T = T';
X = X';
plot(X, T(:,1))

legend('n=1;alternativa i')

 Aplicação ii

function dydx = bvpfcn2(x,y) %#ok<*INUSL>

global L Pe Da %#ok<*GVMIS>
dydx = zeros(2,1); %#ok<PREALL>

dydx = [y(2)
(y(2) + Da*y(1))/(1/Pe)];

function res = bcfcn2(ya,yb)

global Pe %#ok<GVMIS>

res = [(1/Pe)*ya(2)-ya(1)+1

yb(2)];
end

function g = guess2(x) %#ok<*INUSD>

g = [0
1];
end

clear all %#ok<CLALL>

clc

global L k Pe Da %#ok<GVMIS>
L = 1; k = 4.5e-4; v=3e-4;
Da=3;

for Pe=0.1:200
deltaL = linspace(0,L,10);
solinit = bvpinit(deltaL, @guess2);
sol = bvp4c(@bvpfcn2, @bcfcn2, solinit);
X = sol.x; T = sol.y;
T = T';
X = X';
end

plot(X, T(:,1))

legend('n=1;alternativa ii')

 Aplicação iii

function dydx = bvpfcn2(x,y) %#ok<*INUSL>

global L Pe Da %#ok<*GVMIS>
dydx = zeros(2,1); %#ok<PREALL>

dydx = [y(2)
(y(2) + Da*y(1))/(1/Pe)];

function res = bcfcn2(ya,yb)

global Pe %#ok<GVMIS>

res = [(1/Pe)*ya(2)-ya(1)+1

yb(2)];
end

function g = guess2(x) %#ok<*INUSD>

g = [0
1];
end

clear all %#ok<CLALL>

clc

global L k Pe Da %#ok<GVMIS>
L = 1; k = 4.5e-4; v=3e-4;
Pe=20 ;

for Da=0.1:20
deltaL = linspace(0,L,10);
solinit = bvpinit(deltaL, @guess2);
sol = bvp4c(@bvpfcn2, @bcfcn2, solinit);
X = sol.x; T = sol.y;
T = T';
X = X';
end
plot(X, T(:,1))

legend('n=1;alternativa iii')

2.1.3. Código Diferenças Finitas

 Alternativa i

n=1

function [ dydt ] = pfrdf( t ,y ) %#ok<INUSL>

global deltax X0 X9 %#ok<GVMIS>

Pe=3;
Da=1.5;
X1 = y(1); X2 = y(2); X3 = y(3); X4 = y(4); X5 = y(5); X6 = y(6);
X7 = y(7); X8 = y(8);
X9=(4*X8-X7)/3;
X0=(4*X1-X2+1)/4;

dydt = zeros(8,1);

dydt(1) = (1/Pe)*(X2 - 2*X1 + X0)/deltax^2 - Da*X1;

dydt(2) = (1/Pe)*(X3 - 2*X2 + X1)/deltax^2 - Da*X2;
dydt(3) = (1/Pe)*(X4 - 2*X3 + X2)/deltax^2 - Da*X3;
dydt(4) = (1/Pe)*(X5 - 2*X4 + X3)/deltax^2 - Da*X4;
dydt(5) = (1/Pe)*(X6 - 2*X5 + X4)/deltax^2 - Da*X5;
dydt(6) = (1/Pe)*(X7 - 2*X6 + X5)/deltax^2 - Da*X6;
dydt(7) = (1/Pe)*(X8 - 2*X7 + X6)/deltax^2 - Da*X7;
dydt(8) = (1/Pe)*(X9 - 2*X8 + X7)/deltax^2 - Da*X8;

end
clear all
clc
global deltax X0 X9 %#ok<GVMIS>
L = 1.0; np = 10; nt = 100;
X10 = 0.3097; X0 = 0.7334;

deltax = L/(np-1);

w = (0:deltax:L);
y0 = 1*ones(8,1);
tspan = linspace(0,100,nt);

[tempo,Y] = ode45(@pfrdf,tspan,y0);

for i = 1:nt
X(i,1) = (4*Y(i,1)-Y(i,2)+1)/4;;
for j=2:np-1
X(i,j) = Y(i,j-1);
end
X(i,np) =(4*Y(i,8)-Y(i,7))/3;;
end

plot(w,X(nt,:))

n=2

function [ dydt ] = pfrdf( t ,y ) %#ok<INUSL>

global deltax X0 X9 %#ok<GVMIS>

Pe=3;
Da=1.5;
X1 = y(1); X2 = y(2); X3 = y(3); X4 = y(4); X5 = y(5); X6 = y(6);
X7 = y(7); X8 = y(8);
%X9= (1/Pe)*((+ - 2*X9 + X8)/deltax^2)/Da;
%X0=(1/Pe)*((X0 - 2*X1 + X2)/deltax^2) + 1 - Da*X0;
X9=(4*X8-X7)/3;
X0=(4*X1-X2+1)/4;

dydt = zeros(8,1);

dydt(1) = (1/Pe)*(X2 - 2*X1 + X0)/deltax^2 - Da*X1^2;

dydt(2) = (1/Pe)*(X3 - 2*X2 + X1)/deltax^2 - Da*X2^2;
dydt(3) = (1/Pe)*(X4 - 2*X3 + X2)/deltax^2 - Da*X3^2;
dydt(4) = (1/Pe)*(X5 - 2*X4 + X3)/deltax^2 - Da*X4^2;
dydt(5) = (1/Pe)*(X6 - 2*X5 + X4)/deltax^2 - Da*X5^2;
dydt(6) = (1/Pe)*(X7 - 2*X6 + X5)/deltax^2 - Da*X6^2;
dydt(7) = (1/Pe)*(X8 - 2*X7 + X6)/deltax^2 - Da*X7^2;
dydt(8) = (1/Pe)*(X9 - 2*X8 + X7)/deltax^2 - Da*X8^2;

end

clear all
clc

global deltax X0 X9 %#ok<GVMIS>

L = 1.0; np = 10; nt = 100;
X10 = 0.3097; X0 = 0.7334;

deltax = L/(np-1);

w = (0:deltax:L);
y0 = 1*ones(8,1);
tspan = linspace(0,100,nt);

[tempo,Y] = ode45(@pfrdf,tspan,y0);

for i = 1:nt
X(i,1) = (4*Y(i,1)-Y(i,2)+1)/4;;
for j=2:np-1
X(i,j) = Y(i,j-1);
end
X(i,np) =(4*Y(i,8)-Y(i,7))/3;;
end

plot (w,X(nt,:))

3. RESULTADOS
3.1. Rotina bvp4c
3.1.1. Alternativa i

Figura 2: Perfil de concentração do reator em função da coordenada adimensional para 1ª ordem.

 Fonte: Autora.

Figura 3: Perfil de concentração do reator em função da coordenada adimensional para 2ª ordem.

Fonte: Autora.

Utilizando a rotina bvp4 do software Matlab, ajustou-se de maneira eficaz o

código para visualização do comportamento da concentração em coordenada
adimensional. O valor de Da fornece uma estimativa rápida do grau de conversão que
pode ser alcançado, de modo que, quando este valor numérico é menor que 0.1, uma
conversão menor que 10% , e quando é maior que 10, uma conversão maior que 90%
é esperada, comprovada para ambos os casos de primeira e segunda ordem.
Teoricamente, para reação de primeira ordem, espera-se uma conversão de
praticamente 100%, já para segunda ordem este valor encontra-se entre 90% e 100%.
Portanto nota-se que o ajuste do modelo deste reator para o perfil de concentração
em coordenada adimensional para primeira e segunda ordem, fornece resultados
plausíveis com a literatura.

3.1.2. Alternativa ii

Para uma reação de ordem n =1 assumindo o valor de Da = 3,0, pode-se

analisar o efeito de PeM variando de 0 a 200, utilizando a rotina bvp4c do Matlab. Para
metodologia abordada neste trabalho, o ajuste do modelo deste reator para o número
de Peclet estipulado, fornece resultados plausíveis com a literatura. Observa-se que
quando esta variável aumenta afeta o escoamento do reator diminuindo a dispersão
no mesmo, assim o Pe de 0 a 200 nota-se que com seu aumento diminui a dispersão
no mesmo. Isto pode ser observado na figura abaixo.

Figura 4:Efeito de PeM variando de 0 a 200 com o valor de Da = 3,0.

Fonte: Autora.

3.1.3. Alternativa iii

Para uma reação de ordem n =1 fixando o valor de PeM = 20, analisou-se o

efeito de Da variando de 0 a 20, com a rotina bvp4c do Matlab gerando o gráfico
abaixo. Nota-se que para a variável desta variável a concentração adimensional chega
entre 60% e 70%, cumprindo aumento para coordenada adimensional de distância até
o seu valor total equivalente a 1.
 Figura 5: Efeito de Da variando de 0 a 20 com o valor de Pe=20.

Fonte: Autora.

3.2. Método das diferenças finitas

3.2.1. Alternativa i

Já a resolução utilizando o método das Diferenças Finitas dividindo o domínio

em 10 pontos, atingiu-se um resultado menos influente do que para a rotina bvp4c,
uma vez que por este método de discretização, mais pontos estudados evidenciaria a
eficácia e precisão deste método de aproximação. Esta observação aplica-se tanto
para a reação caracterizada de primeira ordem quanto para a de segunda ordem.
 Figura 6: Perfil de concentração do reator em função da coordenada adimensional para 1ª ordem .

Fonte: Autora.

Figura 7: Perfil de concentração do reator em função da coordenada adimensional para 2ª ordem

Fonte: Autora.

4. CONCLUSÃO

Conclui-se, por conseguinte, que o comportamento de diversos processos da

Engenharia Química vem sendo estudados utilizando as ferramentas de modelagem
e simulação, devido às possibilidades em analisar condições do processo sem
geração de gastos operacionais de equipamento, além de permitir testes a fim de
otimizar o estudo. O software Matlab pôde solucionar o perfil de concentração
adimensional auxílio do método de diferenças finitas , já que se realizou o balanço
total para domínio de acordo com a problemática em questão, bem como com a
utilização da rotina bvp4c, sendo esta considerada uma resolução mais apropriada .
Além disso, pôde-se analisar a influência do número de Peclet e Damköhler no reator
PFR, uma vez que, a variável Pe crescente sugere maior fração de fluxo, porém se
decrescente sugere maiores níveis de mistura.
 REFERÊNCIAS

CONNOR, N. .THERMAL-ENGINEERING. O que é o Número Péclet – Definição.

Acesso em: < https://www.thermal-engineering.org/pt-br/o-que-e-o-numero-peclet-
definicao/>.
FOGLER, H. S. Elementos de engenharias das reações químicas. 4. Ed. Rio de
Janeiro: LTC, 2012
FONTANA, ÉLITON. Métodos Numéricos em Engenharia Química. Universidade
Federal do Paraná, UFPR-Programa de Pós-Graduação em Engenharia Química
PPGEQ, 2018.
ISHIDA, S. P. Análise Da Eficiência De Um Reator Tubular (PFR) Submetido A
Diferentes Modos De Alimentação. Universidade Tecnológica Federal Do Paraná
Departamento De Engenharia Química ,Ponta Grossa 2017.
SANTOS S. B. F.; DAMASCENO, D. O.; SILVA, L. N.; DE LIMA, J. S.; DE SOUSA,
C. F.; MOITINHO, J. A. D.; DE OLIVEIRA, L. G. Modelagem de um Reator
Químico do Tipo PFR não Isotérmico com Jaqueta como Instrumento Didático
no Ensino de Engenharia Química: Um Estudo de Caso. Rev. Virtual Quim.,
2019, 11 (1), 50-59.
Silva, S. F; Soares, A. A. B. O MÉTODO DAS DIFERENÇAS FINITAS APLICADO À
TEORIA DAS VIGAS. Traços, Belém, v. 13, n. 27, p. 9-23, jun. 2011.
STEPHANOPOULOS, G. Chemical process control: an introduction to theory and
practice. New Jersey: Prentice Hall, 1984;