ModelagemeSimulacao-Rocha-2019

Prévia do material em texto

NATAL 
2019 
 
 
 
UNIVERSIDADE FEDERAL DO RIO GRANDE DO NORTE 
DEPARTAMENTO DE ENGENHARIA QUÍMICA 
 
 
 
 
ARTHUR LUCAS DO NASCIMENTO ROCHA 
 
 
 
 
 
 
 
MODELAGEM E SIMULAÇÃO DE UMA COLUNA DE 
DESTILAÇÃO MULTICOMPONENTE OPERANDO 
EM REGIME ESTACIONÁRIO 
 
 
 
 
NATAL 
2019 
 
 
 
ARTHUR LUCAS DO NASCIMENTO ROCHA 
 
 
 
 
 
 
 
 
MODELAGEM E SIMULAÇÃO DE UMA COLUNA DE 
DESTILAÇÃO MULTICOMPONENTE OPERANDO 
EM REGIME ESTACIONÁRIO 
 
 
 
 
 
Trabalho de conclusão de curso apresentado à 
Universidade do Rio Grande do Norte, como requisito para 
o recebimento do Bacharel em Engenharia química. 
 
Orientador: André Luis Lopes Moriyama 
 
 
 
NATAL 
2019 
ARTHUR LUCAS DO NASCIMENTO ROCHA 
 
MODELAGEM E SIMULAÇÃO DE UMA COLUNA DE DESTILAÇÃO 
MULTICOMPONENTE OPERANDO EM REGIME ESTACIONÁRIO 
 
Trabalho de conclusão de curso aprovado como requisito parcial à obtenção do 
título de Bacharel, Curso de Engenharia química, Centro de tecnologia, Universidade 
Federal do Rio Grande do Norte. Universidade Federal do Rio Grande do Norte, pela 
seguinte banca examinadora: 
___________________________________ 
Prof. Dr. André Luis Lopes Moriyama 
Orientador – Departamento de Engenharia Química – UFRN 
___________________________________ 
MSc. Maitê Medeiros de Santana e Silva 
Departamento de Engenharia Química – UFRN 
___________________________________ 
Eng. Químico Gustavo Henrique Louredo 
Departamento de Engenharia Química – UFRN 
 
 
 
Natal, 20 de agosto de 2019 
 
 
Sumário 
Lista de figuras ................................................................................................................. 5 
Lista de tabelas ................................................................................................................. 6 
Lista de quadros ................................................................................................................ 6 
 Resumo ...................................................................................................................... 7 
 Abstract ..................................................................................................................... 8 
 Introdução.................................................................................................................. 9 
 Revisão bibliográfica............................................................................................... 11 
4.1 Modelo matemático .............................................................................. 12 
4.1.1 Equações constitutivas ..................................................................... 12 
4.1.2 Equações MESH .............................................................................. 14 
4.2 Método matemático para a resolução do problema .............................. 16 
4.2.1 Método de Newton-Raphson ........................................................... 16 
 Metodologia ............................................................................................................ 20 
5.1 Detalhes no modelo .............................................................................. 20 
5.1.1 Modelo NRTL e Parâmetros de interação ....................................... 21 
5.1.2 Hipóteses .......................................................................................... 22 
5.1.3 Entalpia de residual e entalpia de excesso ....................................... 23 
5.2 Organização do programa em Fortran .................................................. 23 
5.2.1 Disposição das equações e variáveis ............................................... 23 
5.2.2 Variáveis do problema ..................................................................... 23 
5.2.3 Fluxograma do programa ................................................................. 24 
5.2.4 Inicialização ..................................................................................... 25 
5.3 Simulação do processo utilizando ProSimPlus ..................................... 27 
5.4 Estudo paramétrico ............................................................................... 28 
 Resultados ............................................................................................................... 29 
6.1 Resultados com as condições padrão .................................................... 29 
 
6.2 Validação do programa utilizando o ProSimPlus ................................. 32 
6.3 Variação da potência do condensador................................................... 33 
6.4 Variação na composição da alimentação (Alimentação 2) ................... 35 
6.5 Abordagem termodinâmica da destilação ............................................. 36 
6.5.1 Ponto azeotrópico ............................................................................ 36 
6.5.2 Zonas de destilação e natureza dos pontos singulares ..................... 38 
6.6 Estudos complementares ....................................................................... 39 
6.6.1 Retiradas intermediárias .................................................................. 39 
6.6.2 Alimentações simultâneas ............................................................... 40 
6.7 MRLS01 vs MRLS21 ........................................................................... 42 
 Conclusão ................................................................................................................ 44 
 Referências bibliográficas ....................................................................................... 45 
Propriedades termodinâmicas dos componentes puros ................................... 46 
 
 
 
 
 
 
 
5 
 
 
Lista de figuras 
Figura 1 Coluna de destilação. ......................................................................................... 9 
Figura 2 Estágio teórico da coluna de destilação. .......................................................... 11 
Figura 3 Matriz tridiagonal por blocos completa. .......................................................... 18 
Figura 4 Identificação dos blocos na matriz total. .......................................................... 19 
Figura 5 Fluxograma do código do programa. ............................................................... 24 
Figura 6 Flowsheet do ProSimPlus com a coluna de destilação estudada. .................... 28 
Figura 7 Perfil das frações molares de acetona ao longo da coluna. .............................. 30 
Figura 8 Diagrama ternário das frações por prato. ......................................................... 31 
Figura 9 Perfil de temperatura ao longo da coluna. ........................................................ 31 
Figura 10 Gráfico dos perfis de temperatura ao longo da coluna pelo programa e pelo 
simulador. ....................................................................................................................... 32 
Figura 11 Comparação entre os perfis de fração molar na fase líquida calculada pelo 
ProSimPlus e pelo programa. ......................................................................................... 32 
Figura 12 Influência da potência no condensador sobre o destilado. ............................. 33 
Figura 13Influência da potência no condensador sobre as vazões de destilado e resíduo.
 ........................................................................................................................................ 34 
Figura 14 Influência da potência no condensador sobre a temperatura no condensador.
 ........................................................................................................................................ 34 
Figura 15 Influência da potência no condensador sobre o resíduo................................. 35 
Figura 16 Diagrama ternário contendo os pontos singulares. ........................................ 37 
Figura 17 Diagramaternário com as frações no líquido e no vapor por prato ............... 38 
Figura 18 Diagrama ternário para a coluna de destilação com duas retiradas 
intermediárias. ................................................................................................................ 40 
Figura 19 Diagrama ternário com as frações líquidas e vapor por prato para o caso com 
duas alimentações. .......................................................................................................... 41 
Figura 20 Tempo de compilação segundo o número de pratos teóricos por método. .... 42 
Figura 21 Diferença entre os tempos de compilação segundo o número de pratos teóricos.
 ........................................................................................................................................ 43 
 
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319709
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319710
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319711
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319712
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319713
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319714
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319715
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319716
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319717
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319718
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319718
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319719
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319719
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319720
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319721
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319721
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319722
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319722
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319723
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319724
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319728
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319729
file:///C:/Users/arthu/Desktop/Lucas/ENSGTI/Modélisation%20d'opérations%20unitaires%202/TCC/TCC%20-%20Arthur%20Lucas%20do%20Nascimento%20Rocha.docx%23_Toc17319729
6 
 
 
Lista de tabelas 
Tabela 1 Parâmetros de interação binária τij ∙ RT. ......................................................... 21 
Tabela 2 Parâmetros de interação binária αij. ................................................................ 22 
Tabela 3 Condições operatórias da coluna de destilação. .............................................. 29 
Tabela 4 Resultados com as condições padrão. .............................................................. 29 
Tabela 5 Resultados do programa com a variação da composição de alimentação 
(alimentação 2). .............................................................................................................. 36 
Tabela 6 Resultados obtidos com a adição de retiradas intermediárias. ........................ 39 
Tabela 7 Resultados obtidos com duas alimentações simultâneas. ................................ 40 
 
Lista de quadros 
Quadro 1 Descrição das variáveis. ................................................................................. 15 
Quadro 2 Disposição das variáveis e equações em seus respectivos vetores. ................ 23 
 
 
 
 
7 
 
 
 Resumo 
Este trabalho de conclusão de curso se propõe a modelar e simular uma coluna de 
destilação multicomponente em regime estacionário a partir do desenvolvimento de um 
programa computacional escrito em linguagem de programação Fortran. A destilação é 
uma das operações unitárias (separação) mais utilizadas na indústria, daí o interesse em 
abordá-la. Esse processo de separação baseia-se no equilíbrio líquido-vapor de 
componentes em mistura e tem a diferença de volatilidade relativa como força motriz para 
separar os componentes da mistura o que ocorre pela formação de duas fases (líquido e 
vapor) com concentrações diferentes (fase vapor enriquecida do componente mais 
volátil). A validação dos resultados obtidos foi realizada utilizando uma ferramenta de 
simulação fiável e consolidada no mercado (ProSimPlus), o erro relativo percentual dos 
resultados foi da ordem de 0,01%. O trabalho também visou conceber uma ferramenta 
computacional capaz de adaptar-se a diferentes casos reais, o que foi consolidado a partir 
de estudos paramétricos e interpretação física dos resultados. O presente trabalho 
permitiu, acima de tudo, colocar em pratica os conhecimentos em várias disciplinas 
diferentes cursadas ao longo da graduação em engenharia química. 
 
Palavras-chave: Modelagem, Simulação, Fortran, Equilíbrio líquido-vapor, Coluna de 
destilação. 
 
8 
 
 
 Abstract 
This work proposes to model and simulate a stationary multi-component 
distillation column by developing a computational program written in Fortran 
programming language. Distillation is one of the most used unitary operations 
(separation) in the industry, hence the interest in addressing it. This separation process is 
based on the liquid-vapor equilibrium of the components in the mixture and the difference 
between relative volatility is the driving force to separate the components of the mixture 
which occurs by the formation of two phases (liquid and vapor) with different 
concentrations (vapor phase enriched by the most volatile component). The validation of 
the results obtained was performed using a reliable and consolidated simulation tool 
(ProSimPlus), the relative percentage error of the results was of the order of 0.01%. The 
work also aimed to design a computationaltool capable of adapting to different real cases, 
which was consolidated using parametric studies and physical interpretation of the results. 
Above all, the present work allowed to put into practice the knowledge in several different 
disciplines studied throughout the graduation in chemical engineering. 
 
Keywords: Modeling, Simulation, Fortran, Liquid-vapor equilibrium, Distillation 
column. 
 
 
9 
 
 
 Introdução 
Quando é necessário separar uma mistura de componentes, esteja ela em estado 
líquido, vapor ou líquido-vapor, a destilação (ilustrada pela Erro! Fonte de referência 
não encontrada.) é a operação unitária mais utilizada devido sua facilidade de operação 
e o vasto conhecimento existente adquirido nessa área, muito desse conhecimento obtido 
devido sua utilização pela indústria petroquímica que junto com outros setores industriais 
fazem uso dessa operação (ex.: Indústria de alimentos, cosméticos, sucroalcooleira etc.). 
A base teórica para este processo de separação vem da termodinâmica de equilíbrio de 
fases e da transferência de matéria. 
De forma simples, a destilação consiste em fornecer calor à solução até que esta 
atinja uma temperatura na qual a mistura se torna bifásica, o que caracteriza o equilíbrio 
líquido-vapor (ELV). Essa temperatura é denominada temperatura de ebulição e depende 
da composição e da pressão do sistema. Nessas condições, o vapor formado é mais rico 
no componente mais volátil e é essa diferença de concentração que permite a separação 
dos constituintes. 
Os sistemas que são mais facilmente separáveis são sistemas nos quais os 
componentes tem uma grande diferença na volatilidade relativa (αAB). Por outro lado, 
quando essa variável possui valores similares para as substâncias da mistura, a separação 
se torna mais difícil e, então, requer-se um processo mais elaborado (destilação 
fracionada). Essa dificuldade aumenta mais ainda para o caso dos azeótropos, uma 
mistura que apresenta um comportamento de um componente puro, assim, a concentração 
Figura 1 Coluna de destilação. 
 
Fonte: Portal laboratórios virtuais de processos químicos. 
10 
 
 
de todos os componentes do azeótropo é igual nas fases vapor e liquida, isso faz com que 
seja necessário a utilização de técnicas de separação como destilação azeotrópica ou 
destilação com duas colunas a pressões diferentes. 
Nesse trabalho, uma coluna de destilação a pressão constante e regime 
estacionário será modelada matematicamente e esse modelo será simulado a partir de um 
programa escrito em Fortran. Além disso, serão realizados também estudos para avaliar 
a influência das condições operatórias sobre as variáveis ao longo da coluna, temperatura 
e composição. 
Os principais objetivos deste trabalho são modelar a coluna de destilação 
multicomponente, desenvolver um programa computacional capaz de resolver esse 
modelo, validar o programa utilizando uma ferramenta de simulação consolidada, 
delimitar as zonas de destilação e estudar alguns parâmetros operacionais. 
11 
 
 
 Revisão bibliográfica 
O presente trabalho tem como foco a modelagem e simulação (usando Fortran) de 
uma coluna de destilação à pressão constante, em regime contínuo e estacionário. Para 
isso, é necessário escrever o modelo que descreve o fenômeno assim como implementar 
o método matemático para a resolução 
O modelo matemático em questão será construído à partir das equações do tipo 
MESH, ou de conservação (M-Balanço de massa, E-Equilíbrio termodinâmico, S- 
Somatório das frações, mássicas ou molares e H – Balanço de energia), junto à essas 
equações são adicionadas também equações do tipo constitutivas (relações entre variáveis 
termodinâmicas), como por exemplo, cálculo de entalpia liquida e vapor, lei de pressão 
de saturação, modelo de energia de Gibbs de excesso (Determinação dos coeficientes de 
atividade e entalpia de excesso) ou equações de estado (Cálculo de fugacidade). 
O modelo será escrito de forma genérica, ou seja, não se limita as características 
de um sistema fixo, logo, parâmetros como, entradas, saídas, perdas e ganhas térmicas, 
número de pratos (estágios) teóricos e número de componentes poderão ser alterados. A 
F a seguir representa um prato intermediário da colona de destilação: 
A coluna de destilação em questão possui “Np” estágios teóricos, a numeração 
desse tipo de coluna, por convenção, começa do topo com o prato 1 que representa o 
condensador e termina no prato Np que é o refervedor. 
Figura 2 Estágio teórico da coluna de destilação. 
 
Fonte: Própria. 
 
12 
 
 
4.1 Modelo matemático 
4.1.1 Equações constitutivas 
4.1.1.1 Lei de Antoine 
A equação constitutiva escolhida para calcular a pressão de vapor saturante, essa 
é uma correlação empírica, onde suas constantes são obtidas experimentalmente e 
dependem da substância química. Essa relação matemática está apresentada abaixo: 
Os valores das constantes A, B, C, D e E podem ser encontrados em bases de 
dados em livros ou em bases de dados de simuladores. 
4.1.1.2 Modelo de coeficiente de atividade e fugacidade 
Para descrever o equilíbrio líquido-vapor, existem dois tipos de abordagem. A 
primeira é a abordagem Φ-Φ, ou abordagem homogênea onde o Φ é o coeficiente de 
fugacidade variável que representa a variação entre o comportamento real de um 
componente e seu comportamento no caso ideal. A outra abordagem chama-se abordagem 
heterogênea ou γ – Φ, onde o γ é o coeficiente de atividade e também está relacionado ao 
desvio do caso ideal, porém para a fase liquida (SMITH, 2007). 
Em termos gerais o equilíbrio de fases é caracterizado principalmente pelo 
equilíbrio entre os potenciais químicos(µ) e, uma vez que este é proporcional a fugacidade 
dos componentes químicos, esse equilíbrio pode ser descrito como a igualdade entre as 
fugacidades de uma dada substância nas duas fases: 
 𝜇𝑖
𝑓𝑎𝑠𝑒1(𝑇, 𝑃, 𝑧𝑓𝑎𝑠𝑒1) = 𝜇𝑖
𝑓𝑎𝑠𝑒2(𝑇, 𝑃, 𝑧𝑓𝑎𝑠𝑒2) (Eq. 2) 
O índice i designa um constituinte genérico, enquanto T, P e z são a temperatura, 
pressão e composição respectivamente. Visto que: 
 𝑑𝜇𝑖 = 𝑅𝑇𝑑𝑙𝑛(𝑓𝑖) (Eq. 3) 
Sendo 𝑓𝑖, a fugacidade do componente i, então: 
 𝑓𝑖
𝑓𝑎𝑠𝑒1
= 𝑓𝑖
𝑓𝑎𝑠𝑒2
 (Eq. 4) 
A partir da definição de coeficiente de fugacidade descrito abaixo, é possível 
alterar a equação acima para obter uma equação com variáveis conhecidas (SMITH, 
2007): 
 
Φ𝑖
𝑓𝑎𝑠𝑒
=
𝑓𝑖
𝑓𝑎𝑠𝑒
𝑍𝑖𝑃
 (Eq. 5) 
 
𝐿𝑛(𝑃𝑠𝑎𝑡(𝑇)) = 𝐴 +
𝐵
𝑇
+ 𝐶 ∙ 𝑙𝑛(𝑇) + 𝐷𝑇𝐸 (Eq. 1) 
13 
 
 
A variável Zi indica composição do componente i e pode será substituída por y ou 
x, para representar as composições nas fases vapor e liquida. 
Utilizando uma abordagem Φ-Φ, a equação de equilíbrio, Eq. 4, pode ser reescrita 
da seguinte forma: 
 𝑦𝑖Φ𝑖
𝑣𝑎𝑝𝑜𝑟P = 𝑥𝑖Φ𝑖
𝑙𝑖𝑞𝑢𝑖𝑑𝑜P (Eq. 6) 
A definição de coeficiente de atividade calcula a fugacidade do componente i na 
fase liquida de acordo com a expressão abaixo: 
 
γ𝑖 =
𝑓𝑖
𝑙𝑖𝑞𝑢𝑖𝑑𝑜
𝑥𝑖𝑓𝑖
0𝑙𝑖𝑞𝑢𝑖𝑑𝑜
 (Eq. 7) 
A variável 𝑓𝑖
0𝑙𝑖𝑞𝑢𝑖𝑑𝑜
 é a fugacidade padrão do componente i na fase liquida e não 
leva em consideração a presença de outras substâncias na mistura. Aplicando a Eq. 7 e a 
Eq. 5, obtém-se a equação de equilíbrio líquido-vapor (ELV) segundo uma abordagem γ 
– Φ (SMITH, 2007): 
 𝑦𝑖Φ𝑖
𝑉P = 𝑥𝑖𝛾𝑖
𝐿𝑓𝑖
0𝐿 (Eq. 8) 
Nota-se que, na Eq. 8, os sobrescritos líquido e vapor foram trocados pelas letras 
L e V, para simplificar a notação. 
Ambas abordagens podem ser utilizadas para representar o equilíbrio líquido 
vapor, porém dependendo do sistema uma abordagem é preferível em relação à outra 
(MARIAS, 2019). Por exemplo, o coeficiente de fugacidade é calculado por meio das 
equações de estado (Peng-Robinson, RK, SRK etc.) e na maioria dos casos estas 
descrevem bem o comportamento da fase vapor, entretanto, as equações de estado não 
funcionam muito bem para descrever o comportamentoda fase líquida (com exceção de 
líquidos não polares). Logo, os modelos de energia livre de Gibbs de excesso (modelos 
de coeficiente de atividade) são frequentemente utilizados para as fases líquidas não 
ideais. 
4.1.1.3 Cálculo de entalpia 
Das variáveis de estado que estão presentes no modelo as entalpias líquidas e 
vapor também se encaixam na categoria de equações constitutivas. 
Para calcular a entalpia de um composto químico, faz-se necessário estabelecer 
condições de referência, ou seja, temperatura, pressão, estado físico e entalpia de 
referência. Essas variáveis de referência podem ser arbitrárias, porém é essencial mantê-
las para todos os cálculos de entalpia do sistema. Além disso, algumas hipóteses e 
14 
 
 
modelos também interveem no cálculo das entalpias de misturas. Nesse contexto, é 
necessário escolher um tipo de abordagem para o equilíbrio (heterogênea ou homogênea) 
a fim de calcular as entalpias residuais ou de excesso. 
No presente trabalho, o estado de referência é: 
• Temperatura (Tref) = 298 K; 
• Pressão = 1 atm; 
• Estado físico = Vapor; 
• Entalpia de referência (href) = 0 j/mol. 
Visando essas condições e não aplicando os modelos de entalpia residual, a 
entalpia do vapor é dada pelo modelo: 
 
𝐻𝑣 = ∑ 𝑦𝑖[𝐶𝑝𝑖,𝑣(𝑇 − 𝑇𝑟𝑒𝑓) ] + 𝐻
∗
𝑛𝑐
𝑖=1
 (Eq. 11) 
Onde nc é o número de substâncias, yi é a fração do componente i na fase vapor, 
Cpiv é a capacidade térmica a pressão constante do componente i puro em estado vapor, 
T é a temperatura do sistema, H* é a entalpia residual e Hv é a entalpia total do vapor. 
A expressão para calcular a entalpia líquida é apresentado abaixo: 
𝐻𝑙 = ∑ 𝑥𝑖[𝐶𝑝𝑖,𝐿(𝑇 − 𝑇𝑒𝑏𝑖) − ∆𝐻𝑣𝑎𝑝𝑖 + 𝐶𝑝𝑖,𝑣( 𝑇𝑒𝑏𝑖 − 𝑇𝑟𝑒𝑓)]
𝑛𝑐
𝑖=1
+ 𝐻𝑒𝑥 (Eq. 12) 
No caso da entalpia líquida, xi é a fração de i na fase liquida, Tebi é a temperatura 
de ebulição de i, ∆𝐻𝑣𝑎𝑝𝑖 é o calor de vaporização de i e Hex é a entalpia de excesso. 
O cálculo de entalpias encerra a lista de equações constitutivas e, a seguir, as 
equações do tipo MESH. 
4.1.2 Equações MESH 
As equações do tipo MESH, balanços de massa total e parciais, equações de 
equilíbrio, a soma das frações na fase líquida e na fase vapor e o balanço de energia são 
aplicados por prato teórico que recebe o índice j e o índice i é mantido para designar os 
componentes (FELDER, 2005),(LUYBEN, 1996). 
De maneira geral, o conjunto dessas equações é apresentado abaixo para todos os 
pratos, porém no encerramento desta seção, alguns pratos são destacados. 
4.1.2.1 Balanço de massa total 
 𝐿𝑗−1 − (𝐿𝑗 + 𝐷𝐿𝑗) − (𝑉𝑗 + 𝐷𝑉𝑗) + 𝑉𝑗+1 + 𝐹𝑗 = 0 (Eq. 13) 
15 
 
 
4.1.2.2 Balanço de massa parcial 
𝐿𝑗−1𝑥𝑖,𝑗−1 − (𝐿𝑗 + 𝐷𝐿𝑗)𝑥𝑖,𝑗 − (𝑉𝑗 + 𝐷𝑉𝑗)𝑦𝑖,𝑗 + 𝑉𝑗+1𝑦𝑖,𝑗 + 𝐹𝑗𝑍𝑖,𝑗 = 0 (Eq. 14) 
4.1.2.3 Equilíbrio 
 𝑦𝑖,𝑗 − 𝑚𝐾𝑖,𝑗𝑥𝑖,𝑗 = 0 (Eq. 15) 
4.1.2.4 Soma (diferença das somas) 
 
∑ 𝑥𝑖,𝑗
𝒏
𝒊=𝟏
− ∑ 𝑦𝑖,𝑗
𝒏
𝒊=𝟏
= 0 (Eq. 16) 
A diferença das somas das frações foi preferida às somas das frações iguais a 1, devido a 
estabilidade numérica dessa equação em relação às outras duas. 
4.1.2.5 Balanço de energia 
𝐿𝑗−1𝐻𝐿𝑗−1 − (𝐿𝑗 + 𝐷𝐿𝑗)𝐻𝐿𝑗 − (𝑉𝑗 + 𝐷𝑉𝑗)𝐻𝑉𝑗 + 𝑉𝑗+1𝐻𝑉𝑗+1 + 𝐹𝑗𝐻𝐹𝑗
+ 𝑄𝑗 = 0 
(Eq. 17) 
Quadro 1 Descrição das variáveis. 
Lj Vazão de líquido saindo do prato j 
Vj Vazão de vapor saindo do prato j 
xi,j Fração do componente i na fase líquida e no prato j 
yi,j Fração do componente i na fase vapor e no prato j 
DLj Saída de líquido intermediária no prato j 
DVj Saída de vapor intermediária no prato j 
Fj Alimentação no prato j 
Zi,j Fração do componente i na alimentação do prato j 
γi,j Coeficiente de atividade do componente i no prato j 
Pi,j
sat Pressaão de saturação do componente i no prato j 
HLj Entalpia da fase liquida do prato j 
HVj Entalpia da fase vapor do prato j 
HFj Entalpia da alimentação do prato j 
Qj Calor perdido ou recebido no prato j 
mKi,j Modelo de equilíbrio aplicado para o componente i no prato j 
Fonte: Própria. 
Nesse conjunto de equações, dois pratos teóricos merecem uma atenção maior, 
eles são o prato equivalente ao refervedor (j = número de pratos (Np)) e o prato 
correspondente ao condensador (j = 1). Para esses dois casos, é raro encontrar 
16 
 
 
alimentações, assim como não existem pratos acima do condensador nem pratos abaixo 
do refervedor, logo Vj+1 e Lj-1 não existem. Qj, nesses casos, representam o calor 
adicionado no refervedor e o calor retirado do condensador. 
Da maneira como o modelo foi escrito, é possível representar colunas de 
destilação com um número genérico de pratos, componentes, entradas e saídas, assim 
como representar perdas ou ganhos de calor ao longo da coluna, o que faz com que o 
modelo seja bem completo podendo, potencialmente, representar operações de destilação 
reais. 
4.2 Método matemático para a resolução do problema 
Uma vez que o modelo foi construído, é necessário recorrer mais uma vez à 
matemática para resolvê-lo. Pela apresentação do modelo, pode-se deduzir erroneamente 
que o problema se trata de um sistema de equações lineares. Contudo, observando melhor 
cada termo das equações, é possível destacar a não linearidade. 
O sistema de equações pode ser escrito na forma: 
 𝐴𝑥 = 0 (Eq. 18) 
Onde A é a matriz completa dos coeficientes das variáveis organizadas no vetor 
x. 
As partes não lineares das equações vêm justamente das equações constitutivas 
utilizados nos cálculos de algumas variáveis. Tratando-se de um sistema de equações não 
lineares, o método escolhido para solucionar numericamente o modelo foi o método de 
Newton-Raphson. 
4.2.1 Método de Newton-Raphson 
O método de Newton-Raphson foi escolhido, pois ele apresenta uma boa e rápida 
convergência, entretanto, faz-se necessário realizar uma boa inicialização para as 
variáveis que serão calculas, por boa inicialização, entende-se dar valores iniciais 
coerentes às variáveis que apresentem sentido físico. Assim, acelera-se a resolução do 
problema visto que os valores iniciais não estão muito afastados do resultado. 
O método de Newton-Raphson é um método iterativo que permite encontrar o 
zero de uma função F(x) a partir do desenvolvimento de Taylor de ordem 1 (Eq. 19) dessa 
função (BURDEN, 2010). 
 𝑇(𝑥) = 𝐹(𝑏) + 𝐹′(𝑏)(𝑥 − 𝑏) (Eq. 19) 
17 
 
 
Sendo T(x) a expansão de Taylor de ordem 1 de F em torno de b, encontrar o zero 
da função F é a mesmo que encontrar o zero de T. Substituindo b por xk e levando em 
conta que função f pertence ao Rn, pode-se obter a equação abaixo. 
 0 = 𝐹(𝑥𝑘) + 𝐽𝑎𝑐(𝑥𝑘)(𝑥𝑘+1 − 𝑥𝑘) (Eq. 20) 
Onde F(x) é a função objetivo, xk é o valor de x na iteração k, xk+1 é o valor de x 
para a próxima iteração e Jac(xk) é a matriz jacobiana da função F(x) avaliada na iteração 
k. 
Então resolver o modelo matemático consiste em transformar o sistema não-linear 
em um sistema linear, onde cada iteração é um passo em direção a solução (caso a solução 
convirja). Como na programação o zero “não existe”, faz-se necessário criar um critério 
de tolerância (ε), um número real próximo de zero para a realidade do problema. O 
processo iterativo para quando o “passo”, calculado a partir da norma do vetor [xk+1 - xk], 
é inferior ao critério de tolerância. 
O sistema linear a ser resolvido é o sistema abaixo: 
 𝐽𝑎𝑐(𝑥𝑘)(𝑥𝑘+1 − 𝑥𝑘) = −𝐹(𝑥𝑘) (Eq. 21) 
Durante esse trabalho, os sistemas lineares obtidos em cada iteração serão 
resolvidos por 2 métodos diferentes. Para entender a diferença entre esses métodos, é 
necessário explicar as dimensões do sistema. O vetor da função objetivo possui 
(2*Np*Nc + 3*Np), Nc equações de balanço de material parcial, Nc equações de 
equilíbrio, uma equação de balanço de massa total, um balanço de energia e a diferença 
entre os somatórios de y e x, tudo isso multiplicado pelo número de pratos teóricos Np. 
Consequentemente, a matrizjacobiana trata-se de uma matriz quadrada de tamanho 
(2*Np*Nc + 3*Np)X(2*Np*Nc + 3*Np), logo, dependendo do número de pratos (Np) e 
de componentes (Nc) o esforço de cálculo computacional necessário para resolver o 
sistema é enorme. 
Outro fator que deve ser levado em conta é que nessa matriz existem várias 
entradas nulas e que mesmo que elas não façam diferença no cálculo, o computador 
realiza mesmo assim as operações com esses valores. Pensando nisso, percebeu-se que 
existe um padrão na configuração da matriz e esse padra está representado na F. 
 
18 
 
 
Nesse sentido, um dos métodos (MRLS01) de resolução de sistemas lineares vai 
resolver o sistema utilizando todas as entradas da matriz (incluindo as nulas), enquanto 
que o outro (MRLS01), vai resolver o sistema apenas utilizando os blocos da matriz. 
Esses blocos estão identificados na F. 
 
Figura 3 Matriz tridiagonal por blocos completa. 
 
Fonte: MARIAS, F. (2019). Material de suporte da disciplina : Modélisation des opérations unitaires 
2. 
 
19 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 4 Identificação dos blocos na matriz total. 
 
Fonte: Alterado de: MARIAS, F. (2019). Material de suporte da disciplina : Modélisation des 
opérations unitaires 2 
20 
 
 
 Metodologia 
A partir dos conhecimentos adquiridos durante a formação, principalmente, nas 
áreas de balanços de massa e energia, termodinâmica, operações unitárias e programação, 
a metodologia deste trabalho consiste na elaboração do modelo que descreve o fenômeno, 
escolha do método de resolução deste modelo, programação do modelo em linguagem 
Fortran, simulação do fenômeno com o auxílio de um software (ProSimPlus), 
comparação e análise dos resultados e estudo paramétrico. 
Primeiramente, o modelo que será utilizado foi demonstrado e justificado na seção 
3, em seguida, o mesmo foi feito para o método de resolução, porém faz-se necessário 
ainda especificar alguns aspectos do modelo. 
Do ponto de vista de resolução do modelo, foi criado um programa em Fortran do 
qual o fluxograma será mostrado e explicado nessa seção. Este código resolverá um 
problema “padrão” que também será melhor abordado a seguir. 
A etapa de simulação utilizando um software já consolidado no mercado (francês 
e europeu) é de extrema importância para o trabalho, pois, a partir dessa simulação o 
programa desenvolvido elaborado pode ser validado. 
O estudo paramétrico tem como objetivo representar casos reais que possam 
ocorrer no meio industrial como uma modificação no tipo de mistura a ser destilada, 
variações nas potências no condensador e no refervedor, retiradas em pontos 
intermediários da coluna, adição de mais de uma alimentação etc. 
No âmbito deste trabalho, todos as etapas serão realizadas utilizando uma mistura 
ternária de acetona, querosene e clorofórmio. Esta mistura foi escolhida, pois trata-se de 
um sistema completo do ponto de vista termodinâmico, pois no diagrama ternário desse 
sistema formado com as curvas de resíduo é possível destacar vários pontos singulares e 
existe um azeótropo nessa mistura. 
5.1 Detalhes no modelo 
Durante a revisão bibliográfica, foi-se discutido a respeito dos modelos de 
coeficiente de atividade, coeficiente de fugacidade e tipos de abordagem, contudo, tudo 
foi feito de forma geral. Nesta subseção, as escolhas e justificativas relacionas ao modelo 
serão abrangidas. 
Antes de tudo, vale destacar que a abordagem heterogênea foi escolhida pra 
caracterizar o equilíbrio, pois o sistema escolhido contém substâncias polares o que 
dificulta a representação dele a partir das equações de estado. 
21 
 
 
Retomando a equação 15 e a equação 8, podemos escrever um modelo para a 
constante de equilíbrio K. 
 
𝐾𝑖,𝑗 =
𝛾𝑖
𝐿𝑓𝑖
0𝐿
Φ𝑖
𝑉𝑃
 (Eq. 22) 
O modelo de coeficiente de atividade escolhido foi o modelo NRTL, pois este é 
um dos modelos mais difundidos de energia livre de Gibbs de excesso e também existe 
uma vasta base de dados para ele. 
5.1.1 Modelo NRTL e Parâmetros de interação 
O modelo NRTL é um modelo que calcula a composição local dos componentes 
baseado nos parâmetros de interação binária, estes por sua vez são obtidos de forma 
empírica. 
Segundo o modelo NRTL o coeficiente de atividade de um determinado 
componente i é calculado a partir da expressão abaixo: 
 
ln(𝛾𝑖) =
∑ 𝑥𝑗𝜏𝑗𝑖𝐺𝑗𝑖
𝑛𝑐
𝑗=1
∑ 𝑥𝑘𝐺𝑘𝑖
𝑛𝑐
𝑘=1
+ ∑
𝑥𝑗𝐺𝑖𝑗
∑ 𝑥𝑘𝐺𝑘𝑗
𝑛𝑐
𝑘=1
𝑛𝑐
𝑗=1
(𝜏𝑖𝑗 −
∑ 𝑥𝑚𝜏𝑚𝑗𝐺𝑚𝑗
𝑛𝑐
𝑚=1
∑ 𝑥𝑘𝐺𝑘𝑗
𝑛𝑐
𝑘=1
) (Eq. 23) 
 
Onde n é o número de constituintes da mistura, i e j designam os componentes na 
mistura, 𝜏 , αij e G são parâmetros de interação binária, sendo G calculado pela equação 
10. 
 𝐺𝑖𝑗 = exp (−𝛼𝑖𝑗𝜏𝑖𝑗) (Eq. 24) 
Por definição, Gii = 1, os parâmetros de interação binária foram obtidos da base 
de dado do ProSimPlus. 
Tabela 1 Parâmetros de interação binária τij ∙ RT. 
𝜏𝑖𝑗 ∙ 𝑅𝑇 Acetona Benzeno Clorofórmio 
Acetona 0 -193.34 -643.277 
Benzeno 569.931 0 0 
Clorofórmio 229.457 0 0 
Fonte: Base de dados do simulador ProSimPlus. 
 
 
 
 
22 
 
 
Tabela 2 Parâmetros de interação binária αij. 
𝛼𝑖𝑗 Acetona Benzeno Clorofórmio 
Acetona 0 0.3007 0.3043 
Benzeno 0.3007 0 0 
Clorofórmio 0.3043 0 0 
Fonte: Base de dados do simulador ProSimPlus. 
5.1.2 Hipóteses 
Agora que, o tipo de abordagem foi escolhido, a constante de equilíbrio pôde ser 
representada e o modelo NRTL foi exibido, cabe-se apresentar algumas hipóteses do 
modelo: 
H1: Gás ideal; 
H2: 𝑓𝑖
0𝐿 = 𝑃𝑖
𝑠𝑎𝑡. 
A primeira hipótese permite simplificar a representação termodinâmica da fase 
vapor, uma vez que, o coeficiente de fugacidade Φ𝑖
𝑉 para um gás ideal é igual a 1. Esta 
hipótese é coerente, pois em baixas pressões (condições operatórias) os gases tem um 
comportamento próximo do comportamento de gás perfeito. 
Para o cálculo da fugacidade padrão do componente puro, existem dois modelos, 
o modelo padrão e o modelo com correção de Poynting descritos abaixo (SMITH, 2007): 
 𝑓𝑖
0𝐿 = 𝑃𝑖
𝑠𝑎𝑡Φ𝑖
𝑉 (Eq. 25) 
 
𝑓𝑖
0𝐿 = 𝑃𝑖
𝑠𝑎𝑡Φ𝑖
0𝑉exp (
𝑣𝑖
0𝐿
𝑅𝑇
(𝑃 − 𝑃𝑖
𝑠𝑎𝑡)) (Eq. 26) 
Onde 𝑣𝑖
0𝐿 é o volume molar padrão da espécie i pura e R é a constante universal 
dos gases. A expressão exp (
𝑣𝑖
0𝐿
𝑅𝑇
(𝑃 − 𝑃𝑖
𝑠𝑎𝑡)) é conhecida como fator de correção de 
Poynting. 
Considerando a primeira hipótese (H1) e que a pressão do sistema é baixa o fator 
de Poynting tem valor aproximado igual a 1. Logo a hipótese (H2) é justificável. 
Aplicando H1 e H2 na expressão da constante de equilíbrio, a equação se resume 
a: 
 
𝐾𝑖,𝑗 =
𝛾𝑖
𝐿𝑃𝑖
𝑠𝑎𝑡
𝑃
 (Eq. 27) 
Para calcular 𝑃𝑖
𝑠𝑎𝑡, a lei de Antoine apresentada na seção 4.1.1.1 é utilizada com 
os coeficientes A,B,C,D e E da base de dados do ProSimPlus. 
23 
 
 
5.1.3 Entalpia de residual e entalpia de excesso 
Utilizando a hipótese H1, pode-se afirmar que o termo residual (H*) da Eq. 12 é 
igual a zero. Contudo o termo de excesso (Hex) da Eq.12 é calculado pela expressão abaixo 
(SMITH, 2007): 
 
𝐻𝑒𝑥 = −𝑅𝑇
2 ∑
𝜕ln (𝛾𝑖)
𝜕𝑇
𝑛𝑐
𝑖=1
 (Eq. 28) 
 
5.2 Organização do programa em Fortran 
5.2.1 Disposição das equações e variáveis 
Para que o sistema possa ser resolvido pelo método escolhido, é preciso que não 
haja elementos nulos na diagonal principal e, para isso, disposição das equações e das 
variáveis na forma vetorial deve ser escolhida de maneira a respeitar esse critério. 
Pensando nisso, a configuração apresentada na Quadro 2 foi sugerida. 
Quadro 2 Disposição das variáveis e equações em seus respectivos vetores. 
Variáveis Equações 
Vj Balanço de massa total 
yi,j Balanço de massa parcial de i no prato j 
Tj Balanço de energia do prato j 
xi,j Equação de equilíbrio do prato j 
Lj Diferença das somas de x e y 
Fonte: Própria. 
5.2.2 Variáveis do problema 
Observando o conjunto de equações apresentado no modelo matemáticoe a 
descrição das variáveis, é possível destacar facilmente as variáveis do problema, frações 
molares no líquido, frações molares no vapor, vazões de líquido e vapor e temperatura. 
Esse é o conjunto de variáveis para um prato teórico, logo esse número multiplica-se 
quando de acordo com as dimensões da coluna. Então, faz-se necessário estruturar um 
sistema de índices para encontrar qualquer uma dessas variáveis em qualquer estágio 
teórico. O sistema de índices criado é apresentado a seguir: 
𝐼𝑉 = 1 
𝐼𝐶𝑉 = 𝐼𝑉 + 1 
𝐼𝑇 = 𝐼𝐶𝑉 + 𝑁𝑐 
24 
 
 
𝐼𝐶𝐿 = 𝐼𝑇 + 1 
𝐼𝐿 = 𝐼𝐶𝐿 + 𝑁𝑐 
Com IV, ICV, IT, ICL et IL é possível localizar as variáveis do problema e, se 
deseja-se escolher um prato teórico j qualquer, basta adicionar um fator multiplicativo (j-
1)(2Nc + 3). 
5.2.3 Fluxograma do programa 
Para apresentar de forma visual as etapas da programação o fluxograma abaixo 
(F) foi desenhado e alguns processos serão melhor descritos em detalhe. A partir dele, é 
possível ter uma noção dos diferentes tipos de sub-rotinas e módulos utilizados. 
O primeiro passo na lógica do progama é a definição dos parâmetros globais dos 
problemas. Esses parâmetros são declarados em módulos. Em seguida, ocorre a leitura de 
um arquivo de entrada contendo as informações relativas a esses parâmetros. A paritr daí 
Figura 5 Fluxograma do código do programa. 
 
Fonte: Própria. 
 
25 
 
 
a inicialização pode ser feita. Existem também subrotinas que se encarregam do cálculo 
de coefientes de atividade, entalpia e pressão de saturação. Outro passo paralelo à 
inicialização, é a definição da função objetivo e das variáveis. Fornecendo as informações 
de inicialização e função objetivo, é possível resolver o problema com o auxílio de outras 
subrotinas como MRLS01 e MRLS21, gerando, então, os resultados que correspondem 
ao perfil das variáveis ao longo da coluna de destilação. 
5.2.4 Inicialização 
Como abordado na seção 4.2.1, a inicialização é um aspecto fundamental para a 
convergência do método de Newton-Raphson. Portanto, uma atenção especial deve ser 
dada à essa etapa. 
5.2.4.1 Equilíbrio na temperatura de bolha 
Para determinação das condições iniciais no condensador, faz-se necessário 
calcular a fase vapor em equilíbrio com uma mistura liquida contendo as frações molares 
da alimentação, ou seja, o ponto de bolha da solução de alimentação (pressão constante). 
De forma simples, soma das equações de equilíbrio por componente (Eq. 15) são 
somadas originando a equação abaixo. 
 
1 − ∑
𝑥𝑖𝛾𝑖𝑃𝑖
𝑠𝑎𝑡
𝑃
= 0
𝑛
𝑖=1
 (Eq. 29) 
Essa equação é caracterizada como uma equação não-linear com uma só incógnita, 
a temperatura que está implícita nos modelos de coeficiente de atividade e pressão de 
saturação. A solução dessa equação é a temperatura de bolha (Tbolha), temperatura 
mínima na qual o sistema começa a tornar-se bifásico. As frações molares da fase 
vapor(yibolha) são, diretamente, calculadas pelas equações de equilíbrio por componente. 
5.2.4.2 Equilíbrio na temperatura de orvalho 
Para determinação das condições iniciais no refervedor, faz-se necessário calcular 
a fase líquida em equilíbrio com uma mistura vapor contendo as frações molares da 
alimentação, ou seja, o ponto de orvalho da solução de alimentação (pressão constante). 
Diferentemente do caso anterior, no caso do ponto de orvalho as frações molares 
da fase líquida são desconhecidas, o que torna o problema levemente mais complicado 
que o caso anterior. Ao invés de encontrar a solução para a soma das equações de 
equilíbrio, esta será resolvida em simultâneo com as equações de equilíbrio por 
componente, assim o sistema a resolver é: 
26 
 
 
 
1 − ∑
𝑥𝑖𝛾𝑖𝑃𝑖
𝑠𝑎𝑡
𝑃
= 0
𝑛
𝑖=1
 (Eq. 30) 
 𝑥𝑖𝛾𝑖𝑃𝑖
𝑠𝑎𝑡 − 𝑦𝑖𝑃 = 0 (Eq. 31) 
Avaliando os graus de liberdade nota-se que as variáveis são (Nc)frações molares 
da fase líquida (xiorvalho) e a temperatura de orvalho (Torvalho), em contra partida, tem-
se (Nc)equações de equilíbrio e uma soma das equações de equilíbrio (considerando o 
somatório do yi). Tendo o número de variáveis igual ao número de equações linearmente 
independentes o sistema é possível de resolver. 
Solucionando o sistema de equações não lineares, encontra-se as condições no 
refervedor. 
5.2.4.3 Condições intermediárias 
As variáveis, temperatura e frações molares no vapor e no líquido, foram obtidas 
nas duas últimas subseções, sintetizando: 
Esses valores tem sentido físico, pois no condensador encontra-se a temperatura 
mais baixa da coluna e no refervedor, a mais alta. Para os pratos intermediários (j), a 
inicialização dessas variáveis é feita por interpolação linear, ou seja: 
Até o momento, as temperaturas e fações molares inicializadas, contudo, ainda 
restam as vazões de líquido e de vapor. 
5.2.4.4 Balanços de massa e de energia iniciais 
Para inicializar as variáveis que faltam é preciso, resolver algumas equações de 
balanço. Primeiramente, a partir do balanço de massa global e um balanço parcial é 
possível calcular as vazões de destilado e resíduo: 
Balanço de massa 
total: 
𝐹 = 𝐷 + 𝑊 
(Eq. 33) 
Balanço de massa 
parcial: 
𝐹𝑥𝐹 = 𝐷𝑥𝐷 + 𝑊𝑥𝑤 
(Eq. 34) 
Nas equações acima, D e W são, respectivamente, o destilado e o resíduo da 
destilação enquanto que xd e xw são as frações de um componente arbitrário nessas duas 
𝑇1 = 𝑇𝑏𝑜𝑙ℎ𝑎 𝑥𝑖,1 = 𝑧𝑖,𝑎𝑙𝑖𝑚𝑒𝑛𝑡𝑎çã𝑜 𝑦𝑖,1 = 𝑦𝑖,𝑏𝑜𝑙ℎ𝑎 
𝑇𝑁𝑝 = 𝑇𝑜𝑟𝑣𝑎𝑙ℎ𝑜 𝑥𝑖,1 = 𝑥𝑖,𝑜𝑟𝑣𝑎𝑙ℎ𝑜 𝑦𝑖,1 = 𝑧𝑖,𝑎𝑙𝑖𝑚𝑒𝑛𝑡𝑎çã𝑜 
 𝑇𝑗 − 𝑇1
𝑇𝑁𝑝−𝑇1
=
𝑗 − 1
𝑁𝑝 − 1
 (Eq. 32) 
27 
 
 
correntes de matéria. Determinando os valores de D e W, é possível calcular os fluxos de 
matéria ao longo da coluna. 
Na notação do programa, D equivale a V1, o vapor que sai do condensador e W 
corresponde a vazão de líquido do último prato (refervedor), então, LNp. Para o caso geral, 
o estado físico da alimentação pode ser líquido, vapor ou líquido-vapor. Porém, para 
simplificar o equacionamento, num primeiro momento, foi considerado que a alimentação 
é líquida e no ponto de bolha. Dessa maneira, a vazão de líquido acima do prato de 
alimentação tem uma vazão igual a vazão de líquido do condensador(L1) e, abaixo do 
estágio de alimentação, essa vazão é igual a L1+F. Enquanto, não há modificação na vazão 
de vapor que é constante ao longo da coluna e igual à V1+L1. 
Para calcular L1, utiliza-se o balanço de energia no condensador: 
 
𝐿1 = 
𝑄𝑐𝑜𝑛𝑑𝑒𝑛𝑠𝑒𝑢𝑟
∆𝐻𝑣𝑎𝑝𝑜𝑟𝑖𝑠𝑎𝑡𝑖𝑜𝑛
 (Eq. 35) 
5.3 Simulação do processo utilizando ProSimPlus 
A ProSim é uma empresa de desenvolvimento de softwares e prestadora de 
serviços na simulação de processos, é líder no mercado de softwares de engenharia 
química na Europa. Trata-se de uma empresa consolidada e renomada na área. 
Tendo isso em mente, é de grande valia validar os resultados obtidos com o 
programa utilizando os resultados de um dos simuladores da ProSim, o ProSimPlus. 
28 
 
 
Visando esse objetivo o problema “padrão” foi também executado no ProSimPlus, 
guardando todas as condições iguais (modelo termodinâmico, hipóteses, número de 
pratos, condições de alimentação, potência no condensador e no refervedor etc). O 
flowsheet da simulação (F) é apresentado abaixo. 
A área de trabalho do simulador com a operação unitária selecionada (exibida na 
F), no caso, uma coluna de destilação com condensador parcial. Os dados fornecidos (no 
caso desse trabalho) para esse modelo de coluna no ProSimPlus são a quantidade de 
pratos, a posição de alimentação e as potências no refervedor e condensador. 
5.4 Estudo paramétrico 
Referente ao estudo paramétrico, o objetivo é de representar casos reais de 
operação de colunas de destilação como variações nas potências do condensador e 
refervedor, flutuações na composição da alimentação, alimentações em paralelo, 
obtenção de produtos intermediários e estudo da diferença de tempo de cálculo entreos 
métodos MRLS01 e MRLS21. 
 
Figura 6 Flowsheet do ProSimPlus com a coluna de destilação estudada. 
 
Fonte: Própria. 
 
29 
 
 
 Resultados 
Para testar o programa, foram escolhidas inicialmente algumas condições 
“padrão” como quantidade de componentes, componentes, número de pratos teóricos, 
potências no refervedor e no condensador etc. Essas condições foram retiradas de um 
exercício proposto em (MARIAS, 2019) e são apresentadas na Tabela 3. 
Tabela 3 Condições operatórias da coluna de destilação. 
Número de pratos teóricos 30 
Potência no refervedor 19000 cal/s 
Potência no condensador 15000 cal/s 
Prato de alimentação 14 
Pressão de funcionamento 1 atm 
Temperatura de alimentação Temperatura de bolha da mistura 
Vazão molar parcial: 1-Acetona 0.6 mol/s 
Vazão molar parcial: 2-Benzeno 0.3 mol/s 
Vazão molar parcial: 3-Clorofórmio 0.1 mol/s 
Fonte: Própria. 
Tendo essas condições como dados de entrada assim como os dados 
termodinâmicos para cada componente é possível obter os primeiros resultados de 
simulação. 
6.1 Resultados com as condições padrão 
Os primeiros resultados para obtidos com as condições padrão são apresentados a 
seguir na Tabela 4, num primeiro momento, será dado destaque para as variáveis de 
resposta nos pratos teóricos 1 e 30 (condensador e refervedor). Posteriormente, será 
apresentado um perfil de algumas dessas variáveis ao longo da coluna. 
Tabela 4 Resultados com as condições padrão. 
Prato(j) 1 30 
Vj(mol/s) 0.5575 2.6042 
y1,j 0.9937 0.1909 
y2,j 0.0063 0.5283 
y3,j 2.29E-05 0.2807 
Tj(K) 329.386 344.564 
x1,j 0.9903 0.104 
x2,j 0.0097 0.6701 
30 
 
 
x3,j 5.20E-05 0.226 
Lj(mol/s) 2.1161 0.4425 
Fonte: Própria 
Nota-se que, na saída do condensador, a concentração de acetona é a mais 
importante chegando à ordem de 99,37% molar. Isso é totalmente coerente às 
características do componente que é o mais volátil dos três. Outro destaque é a 
temperatura de equilíbrio no condensador, 329,38 K, que é muito próxima à temperatura 
de ebulição da acetona pura (329,44 K). 
Na mesma linha do pensamento, ao observar o refervedor (prato 30), pode-se 
observar que o benzeno, componente menos volátil, é a substância predominante, porém 
há quantidades não desprezíveis de clorofórmio e acetona. A temperatura de equilíbrio 
do refervedor é maior do que as temperaturas de ebulição da acetona e do clorofórmio, 
mas o fato destes componentes estarem em mistura permite que eles permaneçam na fase 
aquosa. 
A Fe a F mostram o perfil das composições líquidas e vapor ao longo da coluna. 
 
 
 
Figura 7 Perfil das frações molares de acetona ao longo da coluna. 
 
Fonte: Própria. 
 
0
0,2
0,4
0,6
0,8
1
1 6 11 16 21 26F
ra
çõ
es
 m
o
la
re
s 
em
 a
ce
to
n
a
Nº de Pratos
Perfil da fração molar de acetona ao longo da 
coluna
Vapor
Líquido
31 
 
 
Observando o diagrama (Figura 8) e o gráfico (Figura7), é possível reforçar o fato 
de que a acetona é concentra ao subir da posição de alimentação até o condensador (curvas 
vão da esquerda para direita). 
Na F abaixo, o perfil de temperatura é mostrado ao longo da coluna: 
À medida que os componentes voláteis vão tendo sua fração reduzida na mistura 
(aumento no número de pratos), a temperatura aumenta. Do ponto de vista físico, o 
Figura 8 Diagrama ternário das frações por prato. 
 
Fonte: Própria 
 
Figura 9 Perfil de temperatura ao longo da coluna. 
 
Fonte: Própria 
 
328
330
332
334
336
338
340
342
344
346
1 6 11 16 21 26
Te
m
p
er
at
u
ra
 (
K
)
Nº de Pratos
Perfil de temperatura ao longo da coluna
32 
 
 
refervedor é a fonte de calor da coluna de destilação então, é natural que nesse prato a 
temperatura seja a mais elevada. 
6.2 Validação do programa utilizando o ProSimPlus 
Guardando todas as condições para o a simulação e o programa, as F e a F exibem 
a comparação entre os resultados obtidos no simular (linha) e os resultados do algoritmo 
desenvolvido (pontos). 
 
Figura 10 Gráfico dos perfis de temperatura ao longo da coluna pelo programa e pelo 
simulador. 
 
Fonte: Própria. 
* 
320
325
330
335
340
345
350
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29
Te
m
p
er
at
u
ra
 (
K
)
Nº de Pratos
Perfil de temperatura ao longo da coluna
ProSimPlus
Programa
Figura 11 Comparação entre os perfis de fração molar na fase líquida calculada pelo ProSimPlus e 
pelo programa. 
 
Fonte: Própria 
 
0
0,2
0,4
0,6
0,8
1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
Fr
aç
õ
es
 m
o
la
re
s 
n
a 
fa
se
 lí
q
u
id
a
Nº de Pratos
Perfil da fração molar na fase líquida ao longo da 
coluna
Acetona - ProSimPlus Benzeno - ProSimPlus Clorofórmio - ProSimPlus
Acetona - Prorgrama Benzeno - Programa Clorofórmio - Programa
33 
 
 
É perceptível que os erros entre a fração calculada pelo programa e a calculada 
pelo ProSimPlus é mínima da ordem de menos de 0,1%. Tendo isso em vista é possível 
validar que o modelo matemático e a resolução executada pelo programa desenvolvido 
são validados por uma ferramenta profissional de simulação. 
6.3 Variação da potência do condensador 
A variação dos parâmetros operacionais é uma das vantagens da programação, 
uma vez que é possível de estudar a influência desses parâmetros sobre as variáveis e, 
com base nisso, é possível escolher ou encontrar os valores ótimos para responder a 
diferentes demandas. 
Na F abaixo o comportamento da fração de acetona na fase vapor de acordo com 
a potência retirada do condensador é exibido. 
No gráfico, nota-se que o vapor é enriquecido em acetona conforme o calor retirado do 
condensador aumenta. Isso pode ser explicado pelo fato de que, quando se aumenta o 
refluxo e, consequentemente, a transferência de matéria é favorecida. Outro ponto 
importante é que ao atingir 15 kcal/s a fração de acetona na fase vapor não aumenta mais, 
indicando assim que esse é um ponto ótimo para a recuperação de acetona, pois se gastaria 
maior quantidade de utilidade fria (Fluido refrigerante) a 20 kcal/s para obter uma mesma 
recuperação de acetona. 
Figura 12 Influência da potência no condensador sobre o destilado. 
 
Fonte: Própria 
 
 
0
0,2
0,4
0,6
0,8
1
11 12 13 14 15 16 17 18 19
Fr
aç
ão
 m
o
la
r 
d
e 
ac
et
o
n
a 
n
o
 v
ap
o
r
Potência (kcal/s)
Influência da potência no condesador sobre o 
destilado
34 
 
 
Outras variáveis que foram observadas com a variação de potência do 
condensador foram a temperatura e as vazões molares. Exibidas nas figuras 13 e 14 a 
seguir. 
Como seria esperado, quanto mais concentrado em acetona no condensador, mais 
a temperatura se aproxima da sua temperatura de ebulição (329,44 K). Esse tipo de estudo 
permite auxiliar no dimensionamento desse condensador, dando pistas do tipo de fluido 
refrigerante que pode ser usado ou até mesmo do tipo de trocador. 
Figura 14 Influência da potência no condensador sobre a temperatura no condensador. 
 
Fonte: Própria. 
 
328
330
332
334
336
338
340
11 13 15 17 19
Te
m
p
er
at
u
ra
(K
)
Potência (kcal/s)
Influência da potência no condesador 
sobre a temperatura no condesador
Figura 13Influência da potência no condensador sobre as vazões de destilado e resíduo. 
 
Fonte: Própria. 
 
0
0,2
0,4
0,6
0,8
1
1,2
11 13 15 17 19
V
az
ão
 m
o
la
r(
m
o
l/
s)
Potência (kcal/s)
Influência da potência no condesador sobre as 
vazões de destilado e resíduo
Destilado
Resíduo
35 
 
 
A F acima mostra que a vazão de destilado diminui com o aumento da potência, 
esse acontecimento tinha sido previsto na discussão sobre a pureza do destilado. 
Obviamente, pelo balanço de massa a saída do resíduo é compensada por essa diminuição 
na quantidade de destilado (balanço de massa sem reação química). 
Ao concentrar-se no resíduo, pode-se perceber que algo bastante interessante 
também acontece, esse comportamento é mostrado pela F.Pelo gráfico acima, nota-se que a recuperação de benzeno é favorizada por 
potências mais baixas no condensador, isso ocorre, pois o benzeno é o componente mais 
pesado e à medida que a acetona e o clorofórmio saem na região superior da coluna e 
concentra-se na região mais quente da coluna. 
6.4 Variação na composição da alimentação (Alimentação 2) 
Uma variação bastante comum na indústria química é a variação da composição 
de uma alimentação, pois como existe uma linha de processo flutuação nas etapas 
anteriores a destilação podem fazer com a composição que chegue na coluna não seja a 
mesma. Nesse contexto, foi feito um teste de variação de composição para avaliar como 
o modelo reage a esse tipo de diferença cujos resultados são apresentados na Tabela 5. 
A alimentação nomeada “alimentação 2” possui a seguinte vazão molar parcial: 
0,1 mol/s de acetona, 0,2 mol/s de benzeno e 0,7 mol/s de clorofórmio. 
 
 
Figura 15 Influência da potência no condensador sobre o resíduo. 
 
Fonte: Própria 
 
0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
1
10 12 14 16 18 20F
ra
çã
o
 d
e 
b
en
ze
n
o
 n
a 
fa
se
 lí
q
u
id
a
Potência no condensador (kcal/s)
Influência da potência no condesador sobre o 
resíduo
36 
 
 
Tabela 5 Resultados do programa com a variação da composição de alimentação (alimentação 2). 
Prato(j) 1 30 
Vj(mol/s) 0.5506 2.6493 
y1,j 0.0663 0.1838 
y2,j 0.0158 0.3214 
y3,j 0.9179 0.4949 
Tj(K) 336.3 341.65 
x1,j 0.1026 0.1413 
x2,j 0.0244 0.4256 
x3,j 0.873 0.4331 
Lj(mol/s) 2.1258 0.4494 
Fonte: Própria. 
No caso dessa nova composição, é visível que na cabeça da coluna a saída não é 
mais tão pura em acetona, isso ocorre, pois, acetona e clorofórmio têm pontos de ebulição 
muito próximos e a fração molar de clorofórmio é muito importante na alimentação. A 
partir desse resultado, pode-se entender que uma alimentação muito rica em clorofórmio, 
possibilita a obtenção dele como produto majoritário no destilado. O mesmo raciocínio 
pode ser aplicado para a acetona. 
Olhando para o benzeno, verifica-se mais uma vez o enriquecimento de benzeno 
no resíduo, embora a baixa quantidade dele na alimentação. 
6.5 Abordagem termodinâmica da destilação 
Os pontos singulares do sistema são os pontos no diagrama ternário que 
correspondem, aos componentes puros e às misturas azeotrópicas se houveram. Esses 
pontos podem ser classificados como nós estáveis, instáveis ou ponto de cela. 
6.5.1 Ponto azeotrópico 
Para encontrar o ponto azeotrópico e também para testar a generalidade do 
programa elaborado. O código foi compilado utilizando uma alimentação contendo 
somente acetona e clorofórmio. O resultado dessa simulação indicou um ponto 
azeotrópico cuja composição é aproximadamente 0,65 em clorofórmio e 0,35 em acetona, 
a temperatura do azeótropo de 338,2 K. Na literatura, o azeótropo está situado em 0,66 
para a fração em clorofórmio e temperatura de 337,5 K (ROSSI, 2013). 
37 
 
 
O diagrama ternário da Erro! Fonte de referência não encontrada. contém os 
pontos singulares representados, componentes puros e a mistura azeotrópica. 
 
Figura 16 Diagrama ternário contendo os pontos singulares. 
 
Fonte: Própria. 
 
38 
 
 
6.5.2 Zonas de destilação e natureza dos pontos singulares 
Para conseguir determinar a natureza dos pontos singulares e caracterizar as zonas 
de destilação, foram utilizadas três alimentações em pontos diferentes do diagrama e as 
curvas de destilação foram construídas com as frações molares no líquido e no vapor por 
prato (F). 
Figura 17 Diagrama ternário com as frações no líquido e no vapor por prato 
 
Fonte: Própria. 
Os pontos verde e ciano representam, respectivamente, o clorofórmio e a acetona 
puros que por sua vez podem ser identificados como nós instáveis, pois estes têm os 
pontos mais de ebulição mais baixos nas regiões onde estão (329,4 K para a acetona e 
334,3 K para o clorofórmio). Tendo isso em vista, pode-se concluir que estas duas 
substâncias não podem ser obtidas no resíduo. Contudo, é notável que se a composição 
da alimentação está à esquerda da curva que separa as zonas de destilação, o clorofórmio 
é obtido no destilado e, se a composição está à direita, a acetona é obtida na cabeça da 
coluna (LUKACS, 2009). 
O ponto amarelo escuro corresponde ao benzeno puro, este por sua vez pode ser 
identificado como um nó estável, pois, dentre os pontos singulares, o benzeno puro é o 
que tem o ponto de ebulição mais alto (353,2 K). Então, ele sempre será o componente 
prioritário de ser obtido no resíduo, não importa a zona de destilação. 
O azeótropo por sua vez tem uma temperatura de ebulição de 337,5 K, ou seja, ele 
está no meio termo entra os nós instáveis e o nó estável. O que acontece é que esse ponto 
39 
 
 
é denominado ponto de cela, o que caracteriza esse ponto é que as curvas de resíduo 
tendem a se dirigir em direção a ele, porém elas se distanciam para se redirecionar ao nó 
estável. 
A linha azul no meio do diagrama delimita as diferentes zonas de destilação que 
podem ser obtidas no sistema. A existência dessas zonas preconiza a busca pelo estudo 
de sistema uma vez que os produtos podem ser diferentes dependendo da zona, mas 
possuem as mesmas características dentro de uma mesma zona. Por exemplo, na zona da 
esquerda, o clorofórmio será sempre o produto principal no destilado enquanto que, na 
zona da direita, a acetona é obtida na cabeça da coluna, contudo nos dois casos o resíduo 
apresenta a tendência de se enriquecer em benzeno. 
6.6 Estudos complementares 
A fim de explorar um pouco mais as funcionalidades do código, diferentes estudos 
complementares foram feitos. O primeiro estudo avalia o impacto de se ter saídas 
intermediárias, no segundo, o objetivo é de avaliar o efeito de ter-se duas alimentações 
com composições diferentes em simultâneo. 
6.6.1 Retiradas intermediárias 
Para testar se o programa pode ser usado em um caso um pouco mais complexo, 
foram adicionadas duas correntes de saída intermediárias na coluna de destilação, uma 
corrente líquida e outra corrente, vapor. O prato intermediário escolhido foi o prato 10 e 
a vazão de saída de cada uma das retiradas é de 0,1 mol/s. Os resultados obtidos são 
mostrados na Tabela 6 abaixo, destacando o condensador e refervedor. 
Tabela 6 Resultados obtidos com a adição de retiradas intermediárias. 
Prato(j) 1-Padrão 1-Dl et Dv 30-Padrão 30-Dl et Dv 
Vj(mol/s) 0.5575 0.45489 2.6042 2.59317 
y1,j 0.9937 0.9867 0.1909 6.50E-02 
y2,j 0.0063 1.20E-02 0.5283 0.6398 
y3,j 2.29E-05 7.26E-04 0.2807 0.3037 
Tj(K) 329.386 329.511 344.564 347.097 
x1,j 0.9903 0.9799 0.104 3.10E-02 
x2,j 0.0097 1.90E-02 0.6701 0.7588 
x3,j 5.20E-05 1.65E-03 0.226 0.210542 
Lj(mol/s) 2.1161 2.1143 0.4425 0.345109 
Fonte: Própria. 
40 
 
 
Pode-se notar que o código conseguiu calcular o modelo com as retiradas 
intermediárias e, do ponto de vista físico, conclui-se que existe uma redução na 
quantidade de acetona recuperada que passa de, aproximadamente, 0,554 mol/s para 
0,444mol/s, em contrapartida, houve um aumento da ordem de 10% no grau de pureza do 
benzeno que passa de 67,01% a 75,88% molar. 
O diagrama ternário abaixo (Figura 18) mostra como as retiradas afetam a forma 
da curva. 
Figura 18 Diagrama ternário para a coluna de destilação com duas retiradas intermediárias. 
 
Fonte: Própria. 
6.6.2 Alimentações simultâneas 
No caso desse teste, uma segunda alimentação de 1 mol/s foi introduzida (no prato 
5) em paralelo com a alimentação “padrão” (no prato 14). A composição dessa corrente 
é a mesma da alimentação 2. O objetivo é de verificar que mais uma vez o programa seja 
capaz de descrever o problema. 
Os resultados são apresentados a seguir na Tabela 7, mais uma vez, com foco no 
condensador e no refervedor. 
Tabela 7 Resultados obtidos com duas alimentações simultâneas. 
Prato(j) 1 - 2xAlimentações 30- 2xAlimentações 
Vj(mol/s) 0.52942.6564 
y1,j 0.3942 0.4285 
y2,j 0.0749 0.2425 
41 
 
 
y3,j 0.5309 0.3291 
Tj(K) 338.5540 339.3670 
x1,j 0.3620 0.3341 
x2,j 0.0910 0.3130 
x3,j 0.5470 0.3529 
Lj(mol/s) 2.1290 1.4707 
Fonte: Própria. 
Primeiro de tudo, é importante perceber que o balanço de massa total foi 
respeitado, aumentando, assim as vazões molares do destilado e do resíduo. O 
componente majoritário no destilado foi o clorofórmio, mas com pouca diferença em 
relação à acetona, o que é interessante é que o benzeno não é mais o produto majoritário 
no pé da coluna, devido a sua baixa quantidade nas alimentações. 
A figura 19 abaixo apresenta as composições molares ao longo da coluna de 
destilação. 
Figura 19 Diagrama ternário com as frações líquidas e vapor por prato para o caso com duas 
alimentações. 
 
Fonte: Própria 
Nitidamente, é possível perceber a presença das duas alimentações pela mudança 
no perfil das curvas. Outro destaque é que mesmo o benzeno não estando em grande 
quantidade, há um aumento na sua fração no pé da coluna, pois ele é um nó estável. 
Simular este tipo de problema tem um interesse do ponto de vista industrial, pois 
pode-se adequar a posição de entrada da corrente de alimentação de forma a encontrar o 
42 
 
 
melhor prato de alimentação, o prato em que sejam causados o mínimo de perturbações 
ao funcionamento ideal. 
6.7 MRLS01 vs MRLS21 
Por fim, um estudo foi feito para avaliar o desempenho dos dois métodos, porém 
o desafio é criar parâmetros objetivos para fazer esse tipo de avaliação. Pensando nisso, 
a comparação entre os dois métodos será feita a partir do tempo necessário para executar 
uma mesma tarefa. 
Nesse contexto, será observado o tempo de compilação em relação ao número de 
pratos teóricos com um número fixo de componentes (F). 
No gráfico acima percebe-se uma diferença nos tempos de simulação que não são 
tão significativas, porém, ao observar, a evolução dessa diferença com o número de pratos 
teóricos (F), é possível verificar que quando o sistema se torna muito grande essa 
diferença de tempos pode ser significante. 
Figura 20 Tempo de compilação segundo o número de pratos teóricos por método. 
 
Fonte: Própria. 
 
0
5
10
15
20
25
30
35
0 10 20 30 40 50 60 70
Te
m
p
o
(s
)
Nº de pratos
Tempo de simulação vs número de pratos 
teóricos
MRLS01
MRLS21
43 
 
 
Pode-se perceber que a evolução da diferença segue uma tendência com 
características exponenciais, o que reforça a ideia de que, o desempenho do método 
MRLS21 em relação ao MRLS01 aumenta conforme o sistema aumenta. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 21 Diferença entre os tempos de compilação segundo o número de pratos teóricos. 
 
Fonte: Própria 
 
0
1
2
3
4
5
6
0 10 20 30 40 50 60 70
Te
m
p
o
(s
)
Nº de pratos
Diferença entre os tempos de simulção
44 
 
 
 Conclusão 
O desenvolvimento permitiu, em primeiro lugar, de colocar em prática inúmeros 
conhecimentos adquiridos ao longo da formação em Engenharia Química. Áreas como 
balanços materiais e balanços de energia, termodinâmica de equilíbrio, métodos 
matemáticos, programação computacional (Fortran), simulação utilizando uma 
ferramenta consolidada no ramo de processos químico (ProSimPlus) e também a 
modelagem matemática. 
No geral, pode-se concluir que os objetivos foram alcançados, foi construído um 
modelo, este foi resolvido pelo programa desenvolvido durante o trabalho e, em 
comparação por simulação utilizando ProSimPlus, os parâmetros do modelo puderam ser 
testados. 
O programa desenvolvido permite de estudar alguns aspectos da termodinâmica 
de equilíbrio, estudo dos pontos singulares, zonas de destilação, cálculo de pontos 
azeotrópicos etc. Isso é de grande valor para aprendizagem desses fenômenos, inclusive, 
é possível realizar um estudo para validar e adaptar o código para um contexto industrial. 
45 
 
 
 Referências bibliográficas 
 
FELDER, R.M. Princípios elementares dos processos químicos. 3. Ed.: LTC, 2005. 
LUKACS, T. S. Etude de la distillation réactive dans une colonne avec un bac 
intermédiaire avec des réactions consécutives. 2009. 249 p. Tese de doutorado 
(Doutorado) - Institut National Polytechnique de Touloues e Université des Sciences 
Techniques et Economiques de Budapest, Toulouse, 2009. 
LUYBEN, W.L. Process modeling, simulation, and control for chemical engineers. 2 
.Ed.: McGraw-Hill, 1996. 
MARIAS, F. (2019). Material de suporte da disciplina : Modélisation des opérations 
unitaires 2. Pau, France. 
RENEAUME, J.-M. (2017). Material de suporte da disciplina : Distillation. Pau, France. 
SCHWARTZENTRUBER, J. (s.d.). Pôle Numérique Pour La Pédagogie - École des 
Mines d'Albi. Disponível em: https://nte.mines-
albi.fr/ThermoIntro/fr/co/uc_Distillation.html; Acesso em ago. 2019. 
SMITH, J. Introdução à Termodinâmica da Engenharia Química. 7. Ed: LTC, 2007. 
ROSSI, A.S. Cálculo de mapas de curvas residuais aplicando modelos de equilíbrio e 
correção por eficiência. 2013. 129 f. Dissertação (Mestrado em engenharia química) - 
Programa de pós-graduação em engenharia química da Universidade Federal de 
Uberlândia, Uberlândia, 2013. 
BURDEN, R.L. Análise numérica. 9. Ed : Brooks/Cole Cengage learning, 2010. 
Chemstations. (s.d.). Disponível em: https://www.chemstations.com/; Acesso em ago. 
2019. 
Portal laboratórios virtuais de processos químicos. Disponível em: 
http://labvirtual.eq.uc.pt/siteJoomla/index.php?option=com_content&task=view&id=22
3&Ite; Acesso em ago. 2019. 
 
 
 
 
 
 
 
46 
 
 
Anexos 
Propriedades termodinâmicas dos componentes puros 
Propriedades / Componente Acetona Benzeno Clorofórmio 
Temperatura de ebulição 1 atm (K) 329.44 353.24 334.33 
ΔH de vaporisation a Teb e 1 atm 
(J/kmol) 
2.96e7 3.08e7 2.95e7 
Cp – vapor (J/kmol/K) 8.05e4 9.60e4 6.89e4 
Cp – líquido (J/kmol/K) 1.34e5 1.47e5 1.17e5 
Coeficientes da lei de Antoine 
A 69.006 83.107 146.43 
B -5599.6 -6486.2 -7792.3 
C -7.0985 -9.2194 -20.614 
D 6.224e-6 6.984e-6 0.024578 
E 2 2 1 
Fonte: Base de dados do simulador ProSimPlus.