PROGRAMA EQ-ANP
Processamento, Gestão e Meio Ambiente na Indústria
do Petróleo e Gás Natural
DESENVOLVIMENTO DE UM
PROCEDIMENTO PARA SIMULAÇÃO
DINÂMICA DE TAMBORES DE “FLASH”
MULTIFÁSICOS USANDO COMPUTAÇÃO
ALGÉBRICA
Fernanda Mota Gonçalves
Tese de Mestrado
Orientadores
Ofélia de Queiroz Fernandes Araújo, Ph.D.
Marcelo Castier, Ph.D.
Julho de 2001
DESENVOLVIMENTO DE UM PROCEDIMENTO PARA
SIMULAÇÃO DINÂMICA DE TAMBORES DE “FLASH”
MULTIFÁSICOS USANDO COMPUTAÇÃO ALGÉBRICA
Fernanda Mota Gonçalves
Tese submetida ao Corpo Docente do Curso de Pós-Graduação em Tecnologia
de Processos Químicos e Bioquímicos da Escola de Química da Universidade
Federal do Rio de Janeiro, como parte dos requisitos necessários para a
obtenção do grau de Mestre em Ciências.
Orientada por:
________________________________________
Ofélia de Queiroz Fernandes Araújo, Ph.D.
(orientadora – presidente da banca)
________________________________________
Marcelo Castier, Ph.D.
(orientador)
Aprovada por:
________________________________________
Evaristo Chalbaud Biscaia Junior, D.Sc.
________________________________________
Marcio José Estillac de Mello Cardoso, Ph.D.
________________________________________
Maurício Bezerra de Souza Junior, D.Sc.
Rio de Janeiro, RJ - Brasil
Julho de 2001
i
Gonçalves, Fernanda Mota.
Desenvolvimento de um Procedimento para Simulação Dinâmica de
Tambores de “Flash” Multifásicos usando Computação Algébrica/ Fernanda
Mota Gonçalves. Rio de Janeiro: UFRJ/EQ, 2001.
xv, 114 p.; il.
(Dissertação) – Universidade Federal do Rio de Janeiro, Escola de Química,
2001. Orientador(es): Ofélia Queiroz Fernandes Araújo e Marcelo Castier
1. Simulação Dinâmica
2. Equilíbrio de Fases
3. Processos de Separação.
I. EQ/UFRJ
II. Título (série)
ii
Dedico este trabalho a todos os que me incentivaram
a concluir mais esta etapa da minha vida,
em especial àquela que foi muito mais que orientadora,
foi uma grande amiga, Ofélia.
iii
O importante é estar pronto para, a
qualquer momento, sacrificar o que somos
pelo que poderíamos vir a ser.
(Charles Du Bois)
iv
AGRADECIMENTOS
Gostaria de conseguir expressar minha gratidão a todos que
contribuíram para o encerramento desta etapa da minha vida. O interesse pelo
andamento da tese, as palavras de incentivo e a certeza de que eu
conseguiria foram muito importantes.
Inicialmente agradeço aos meus orientadores, Marcelo e Ofélia, que
cumpriram muito além do papel de mestres, me ensinando, incentivando e
sempre dispostos a ajudar.
Agradeço aos meus pais Fernando e Marly, responsáveis pela pessoa
que sou hoje, às minhas irmãs Daniela e Adriana e também às minhas
“branquinhas” Cindy e Megg.
Agradeço ainda a todos os amigos que estiveram presentes mesmo que
distantes. Em especial, aos amigos André (Dedé) e Maria Isabel (Bebel) pelos
telefonemas diários sempre com palavras amigas e de incentivo. Foram muito
mais que alguns minutos de papo, tenham certeza disso.
Agradeço ao Robson, que agüentou com muita paciência e carinho
todos os dias de mau humor e cansaço. Nos momentos de indecisão, me
incentivou, nos momentos de tristeza fez tudo o que podia para me animar.
Sem você não teria conseguido.
Agradeço à galera da piscina pelos excelentes momentos que temos
passado juntos. Em especial agradeço à Dani, minha "sobrinha", pelo carinho.
Agradeço aos professores e funcionários da Escola de Química em
especial aos amigos que fiz na época da graduação e que continuam a me dar
muito carinho: Márcia, Chico, Zé e Aquino, Jane e Maria. Agradeço ainda à
Rosellee e Marlene responsáveis pelos serviços da pós-graduação, apesar de
nunca terem se limitado a tal função. Os momentos que passei com todos
vocês são inesquecíveis.
Agradeço a todos do Laboratório H2CIN (I-221), em especial ao Prof.
José Luiz Medeiros, pela infraestrutura que me foi concedida para desenvolver
este trabalho.
v
Agradeço ainda à Agência Nacional de Petróleo (ANP) pelo apoio
financeiro através do PRH-13.
Enfim, gostaria de agradecer àquele que é responsável pela vida, pelo
amor, pelas alegrias, enfim, é o responsável por tudo. Obrigado Senhor DEUS
pelos inúmeros detalhes visíveis e invisíveis da minha existência.
vi
Resumo da Tese de Mestrado apresentada ao Curso de Pós-Graduação em
Tecnologia de Processos Químicos e Bioquímicos da Escola de Química/UFRJ
como parte dos requisitos necessários para obtenção do grau de Mestre em
Ciências, com ênfase na área de Petróleo e Gás Natural.
DESENVOLVIMENTO DE UM PROCEDIMENTO PARA SIMULAÇÃO
DINÂMICA DE TAMBORES DE “FLASH” MULTIFÁSICOS
USANDO COMPUTAÇÃO ALGÉBRICA
Fernanda Mota Gonçalves
Julho, 2001
Orientadores: Profa. Ofélia de Queiroz Fernandes Araújo, Ph.D.
Prof. Marcelo Castier, Ph.D.
Neste trabalho, um procedimento de simulação dinâmica de tambores de
"flash" foi desenvolvido. Foi considerada a possibilidade de desaparecimento
ou surgimento de fases em equilíbrio além da realização do teste de
estabilidade global de fases, utilizando-se propriedades termodinâmicas
rigorosamente calculadas. O modelo consiste em um conjunto de equações
diferenciais e algébricas (EAD's). Os balanços de massa e energia foram
escritos como equações diferenciais para a energia interna e número de moles
de cada componente. As equações algébricas no sistema constituem relações
de equilíbrio de fases e conservação do número de moles e volume nas fases
presentes no sistema. As propriedades físicas necessárias para o modelo
foram calculadas usando-se a equação de estado de Peng-Robinson.
A implementação computacional do modelo do equipamento e das
propriedades termodinâmicas necessárias foram realizadas utilizando
computação algébrica (Mathematica®, Wolfram Inc.) que permite a geração
automática de códigos em linguagem FORTRAN. Para a simulação, um
ambiente com processamento numérico e interface gráfica foi empregado
(MATLAB®, The MathWorks Inc.). A inteface entre código FORTRAN e o
MATLAB® foi realizada através DLL's (dynamic link libraries) compiladas.
Com o desenvolvimento de tais procedimentos computacionais, a
simulação dinâmica de tambores de “flash” foi realizada em diversas situações,
incluindo o caso de tanques de GLP (Gás Liquefeito de Petróleo), produto de
interesse na indústria de Petróleo e Gás Natural. A abordagem desenvolvida
deverá ser estendida para a modelagem de equipamentos de separação mais
complexos, tais como colunas de destilação.
vii
Abstract of a Thesis presented to Curso de Pós-Graduação em Tecnologia de
Processos Químicos e Bioquímicos - EQ/UFRJ as partial fulfillment of the
requirements for the degree of Master of Science with emphasis on Petroleum
and Natural Gas.
DEVELOPMENT OF A PROCEDURE FOR THE DYNAMIC
SIMULATION OF MULTIPHASE FLASH DRUMS
USING COMPUTER ALGEBRA
Fernanda Mota Gonçalves
July, 2001
Supervisors: Ofélia de Queiroz Fernandes Araújo, Ph.D.
Marcelo Castier, Ph.D.
In this work, a procedure was developed for the dynamic simulation of
flash drums. The possibility of the disappearance of phases, and also of their
appearance, was taken into account by means of the global phase stability test,
and with rigorously calculated thermodynamic properties. The model consists of
a set of differential and algebraic equations (DAEs). The energy and mass
balances were written as differential equations for the internal energy and the
number of moles of each species. The algebraic equations in the system were
the phase equilibrium relationships and mole number and volume conservation
among the phases present in the flash drum. The physical properties required
by the model were calculated using the Peng-Robinson equation of state.
The computational implementation of the drum model and of the
thermodynamic properties were done with the aid of computer algebra
(Mathematica®, Wolfram Inc.), which allowed automatic generation of
FORTRAN codes. For simulation, an environment with numerical and graphical
features (MATLAB®, The MathWorks Inc.) was employed. FORTRAN codes
were interfaced through DLL (dynamic link libraries) compiled versions, callable
from within the MATLAB® programming language.
With the developed computational procedure, the simulated dynamic
behavior of flash drums in several situations was obtained, including the case of
LPG (Liquefied Petroleum Gas) tanks, of interest to the Oil and Gas Industry.
The developed approach is to be extended to model the dynamic behavior of
more complex pieces of separation equipment, such as distillation columns.
viii
ÍNDICE
Capítulo 1 – Introdução
1
Capítulo 2 – Revisão Bibliográfica
4
2.1 – Simulação Dinâmica de Processos
4
2.2 – Simulação de “Flashes”
6
2.3 – Propriedades Termodinâmicas
9
2.4 – Resolução de Equações Algébrico-Diferenciais
10
2.4.1 – O Código ode15s
12
2.4.2 – Resolução de EAD’s no SIMULINK®
14
2.5 – Da Computação Algébrica à Simulação Dinâmica Rigorosa de
Processos: O Problema de “Flash” Multifásico
Capítulo 3 – Formulação
16
19
3.1 – “Flash” UVN
19
3.2 – Implementação Computacional
25
3.2.1 – Integração de códigos FORTRAN ao Ambiente
MATLAB®
28
3.2.2 – Códigos em Ambiente MATLAB®
30
3.2.3 – Códigos em Ambiente SIMULINK®
30
3.3 – Cálculo de Propriedades Termodinâmicas
33
3.3.1 – Propriedades Termodinâmicas
36
3.3.2 – Teste de Estabilidade Global de Fases
37
3.3.3 – Equação de Estado
38
3.4 – “Flash” TVN
39
3.5 – “Flash” TPN
40
Capítulo 4 –Resultados
41
4.1 – Sistema Exemplo 1: BENZENO-ETILBENZENO
41
4.2 – Sistema Exemplo 2: METANO
45
4.3 - Controle para o Sistema Exemplo 1
46
4.4 – Sistema Exemplo 3: GLP
53
ix
Capítulo 5 – Conclusão
65
Referências Bibliográficas
68
Apêndice A1 – Montagem das Equações Algébricas de Equilíbrio
de Fases
Apêndice A2 – Cálculo das Propriedades Termodinâmicas
72
80
A2.1 – Entrada dos dados para geração da rotina flsuvn
80
A2.2 – Rotina flsuvn
82
A2.3 – Entrada dos dados para geração da rotina penrob
89
A2.4 – Rotina penrob
90
x
ÍNDICE DE FIGURAS
Figura 2.1 Diagrama de Blocos do Simulador SIMULINK®
15
Figura 3.1 Algoritmo para Resolução Numérica
27
Figura 3.2 Estrutura MEX FORTRAN
29
Figura 3.3 Modelo em Ambiente SIMULINK®
32
Figura 3.4 “Flash” com Controle de Temperatura
34
Figura 4.1 Sistema Exemplo 1 Submetido a Padrão Temporal
Descrito na Tabela 4.2: Telas de Especificação
e Acompanhamento SIMULINK®
Figura 4.2 Sistema Exemplo 1 sem Controle: Séries Temporais
42
43
Figura 4.3 Sistema Exemplo 1 Submetido a Padrão Temporal
Descrito na Tabela 4.3
47
Figura 4.4 Simulação de Injeção de Gás Metano Puro no
Tambor de “Flash” Submetido a Padrão Temporal
Descrito na Tabela 4.5
48
Figura 4.5 Simulação de Injeção de Gás Metano Puro no
Tambor de “Flash” Submetido a Padrão Temporal
Descrito na Tabela 4.6
Figura 4.6 Sistema Exemplo 1 sob Controle de T e P
49
51
Figura 4.7 Sistema Exemplo 1 sob Controle de T e P:
Séries Temporais
52
Figura 4.8 Adaptação de Ganho do Controlador sob Evento
de Surgimento de Fase
54
Figura 4.9 Desempenho da Estratégia de Controle com
Adaptação de Ganho de Nova Fase Líquida
Figura 4.10 Sistema Exemplo 3: Flash
55
56
Figura 4.11 Tela de Especificação para o Sistema Exemplo 3
em Ambiente SIMULINK®
58
Figura 4.12 Sistema Exemplo 1: Série Temporais
60
Figura 4.13 Sistema Exemplo 3: Temperatura e Carga
Térmica do Sistema
61
xi
Figura 4.14 Sistema Exemplo 3: Resposta da Pressão a
modificações no “Set-Point”
Figura 4.15 Resposta do Controlador de Nível ao Surgimento da Fase
62
63
xii
ÍNDICE DE TABELAS
Tabela 3.1 Ajuste de Passo Máximo para Auxílio de Convergênca
das Equações Agébricas
31
Tabela 4.1 Condições Iniciais do Sistema Exemplo 1
41
Tabela 4.2 Condições de Entrada 1 do Sistema Exemplo 1
44
Tabela 4.3 Condições de Entrada 2 do Sistema Exemplo 1
44
Tabela 4.4 Condições Iniciais do Sistema Exemplo 2
45
Tabela 4.5 Condições de Entrada 1 do Sistema Exemplo 2
45
Tabela 4.6 Condições de Entrada 2 do Sistema Exemplo 2
46
Tabela 4.7 Sintonia de Controladores
50
Tabela 4.8 Adaptação do Ganho do Controlador de Temperatura
53
Tabela 4.9 Condições Iniciais do Sistema Exemplo 3
56
Tabela 4.10 Variação dos “Set-Points” durante a Simulação
57
Tabela 4.11 Condições de Entrada do Sistema Exemplo 3
59
xiii
NOMENCLATURA
LETRAS ARÁBICAS
Símbolo
Descrição
ai, bi
parâmetros do componente i puro na equação de estado
am, bm
parâmetros da mistura
Af
energia livre de Helmhotz na fase f
DLL
)
fi
“dynamic link library”
fugacidade do componente i
Fij
vazão do composto i na corrente j
Hj
entalpia (por unidade de tempo) da corrente j
kij
parâmetro de interação binária
Kc
parâmetro de sintonia do controlador
nc
número de componentes
Nij
número de moles do componente i na fase j
Ni
número total de moles do componente i no interior do
“flash”
np
número de fases
Pf
pressão da fase f
Q
carga térmica no sistema
QF
função objetivo do “flash” UVN
R
constante universal dos gases
S
entropia
T
temperatura
U
energia interna
Up
energia interna na fase p
Vf
volume da f
VTOTAL
volume total do “flash”
Ws
potência de eixo
xij
fração molar do componente i na fase j
xiv
LETRAS GREGAS
Símbolo
Descrição
µij
potencial químico do componente i na fase j
δjp
função delta de Kronecker
τi
parâmetro de sintonia para ação integral do controlador
τD
parâmetro de sintonia para ação derivativa do controlador
α
parâmetro de sintonia para filtro na ação derivativa do
controlador
Ω
função do teste global de estabilidade de fases
SOBRESCRITOS
Símbolo
Descrição
ESPEC
valor especificado
SUBSCRITOS
Símbolo
Descrição
ci
propriedade crítica do componente i
ri
propriedade reduzida do componente i
xv
1
Introdução
A simulação dinâmica de operações da indústria de petróleo e gás natural
tem destaque no cenário de controle de processos, treinamento de operadores,
e diagnóstico e análise de falhas. Neste tipo de aplicação, a etapa de
modelagem termodinâmica é, freqüentemente, feita de forma bastante
simplificada a fim de reduzir o esforço computacional das simulações.
Entretanto, a disponibilidade de computadores cada vez mais velozes está
ampliando a possibilidade de realizar simulações dinâmicas com modelos
termodinâmicos mais rigorosos. Tal esforço se justifica pois, potencialmente,
permite o estabelecimento de condições operacionais e metas de controle
adequadas, fundamentadas no uso de modelos termodinâmicos capazes de
correlacionar e predizer o comportamento de misturas de hidrocarbonetos. Em
tais simulações, várias propriedades termodinâmicas são necessárias, tais
como, entalpia, fugacidade de compostos na mistura, capacidade calorífica,
energia interna, além de diversas de suas derivadas.
Para a obtenção das expressões das diversas propriedades físicas
necessárias a uma simulação dinâmica, a computação algébrica pode ser um
instrumento auxiliar poderoso. Poder-se-ia mesmo pensar em realizar toda a
simulação dinâmica de equipamentos de separação usando programas
comerciais de computação algébrica, tais como Maple® ou Mathematica®.
Entretanto, programas de álgebra computacional são raramente empregados
neste contexto devido ao longo tempo computacional envolvido. Contudo, uma
abordagem híbrida que empregue recursos algébricos para a geração de
códigos em linguagens de cálculos numéricos rápidos, a exemplo de
FORTRAN e C, surge com forte atrativo. Adicionalmente, integrar estes
códigos compilados com um ambiente de computação numérica e visualização
gráfica representa um diferencial significativo como agente facilitador no
desenvolvimento de simuladores de processo.
2
INTRODUÇÃO
O objetivo deste trabalho foi desenvolver a simulação dinâmica de
tambores de “flash”, usando-se este termo de forma ampla para denotar
equipamentos para o armazenamento e a separação de fluidos, típicos da
indústria de refino de petróleo e gás natural. Tal como em diversos modelos
existentes na literatura para aplicações similares, a hipótese central no
desenvolvimento do modelo do equipamento é a de que o equilíbrio
termodinâmico é atingido instantaneamente no interior do tambor de “flash”.
Uma distinção importante, contudo, é que o procedimento de simulação
dinâmica
foi
formulado
usando-se
propriedades
físicas
rigorosamente
calculadas. Especial atenção foi dada à realização do teste de estabilidade de
fases, o que permite que o modelo do “flash” se altere automaticamente, ao
longo da simulação, para incorporar, se necessário, novas fases ao sistema.
De igual modo, testa-se o sistema quanto à eventual necessidade de remoção
de fases ao longo da simulação.
O modelo resultante constitui um sistema de equações algébricodiferenciais. Há diversas estratégias possíveis para a solução de tais sistemas.
Explorá-las sistematicamente e comparar seus desempenhos, entretanto, não
constituiu o escopo central desta tese. Sendo assim, resolveram-se,
separadamente, as equações algébricas do modelo a cada passo de
integração das equações diferenciais.
Do ponto de vista das ferramentas computacionais utilizadas, o
desenvolvimento do procedimento de simulação resultou da integração da
solução proveniente da computação algébrica, com códigos em linguagem
compiláveis e ambiente de resolução numérica e visualização gráfica. Em
particular, usou-se o “software” Mathematica® para a parte de computação
algébrica e a geração automática de códigos em FORTRAN. A resolução das
equações algébricas do modelo foi feita usando-se um código compilado em
FORTRAN, que foi acoplado ao software MATLAB®, responsável pela
integração numérica das equações diferenciais, implementação de estratégias
de controle do processo e pela visualização dos resultados.
Os exemplos desenvolvidos visam a ilustrar os recursos de cálculos
termodinâmicos e a robustez da formulação da resolução de problemas de
INTRODUÇÃO
3
adição e remoção de fases. Entretanto, foram identificadas algumas limitações
na formulação desenvolvida nesta tese, cuja superação deverá ser objeto de
pesquisas adicionais.
Este trabalho está organizado sob a forma de capítulos. O Capítulo 2
contém uma breve Revisão Bibliográfica sobre o processo abordado,
modelagem termodinâmica e resolução numérica.
O Capítulo 3 apresenta a Formulação do problema, destacando o “flash”
com especificação de energia interna, volume e números de mols (“flash”
UVN), a geração de condição inicial e a abordagem numérica adotada.
Exemplos que empregam hidrocarbonetos, e exploram situações de
surgimento e desaparecimento de fases, são apresentados e discutidos em
Resultados no Capítulo 4.
No Capítulo 5, encontram-se as Conclusões e Sugestões.
No Apêndice 1, está disposto o desenvolvimento da função objetivo para
teste de estabilidade de fases. O desenvolvimento computacional do modelo
termodinâmico é apresentado no Apêndice 2.
2
Revisão Bibliográfica
2.1. Simulação Dinâmica de Processos
A simulação de processos é uma ferramenta valiosa no projeto, análise
e operação de processos, e tem se desenvolvido sob o impulso de
especificações rígidas de qualidade de produtos, segurança de processos e
regulamentações ambientais. Um enfoque especial tem sido dedicado à
simulação dinâmica de processos, que depende de uma ferramenta de
resolução de equações que resista a situações como condições iniciais
inconsistentes, restrições de desigualdade e mau condicionamento, para citar
algumas. Adicionalmente, são frequentes as situações em que variáveis de
processos estão sujeitas a restrições, como, por exemplo, de não-negatividade
em temperatura e volume, dificultando a convergência na solução dos modelos.
Outra dificuldade encontrada em simulação de processos é a ocorrência
de comportamentos abruptos, com descontinuidades de derivadas. Marquardt
(1991) aponta que a maioria dos códigos de resolução de equações está apta a
enfrentar modelos de comportamentos suaves. A simulação de processos, no
entanto, apresenta ações discretas e restrições lógicas que podem conduzir a
descontinuidades nas equações do modelo (Gopal e Biegler, 1997). Um
exemplo típico de descontinuidade é encontrado em transição de fases. Esta
situação pode levar a mudanças nas equações do modelo ou mesmo uma
mudança na sua estrutura, exigindo uma representação consistente da
descontinuidade e uma ferramenta de resolução capaz de contornar tal
dificuldade.
Mudanças discretas nas equações de um modelo são marcadas pela
ocorrência de um evento. Os métodos empregados para solução de sistemas
dinâmicos devem ser capazes de localizá-lo.
REVISÃO BIBLIOGRÁFICA
5
Os problemas de simulação dinâmica em um intervalo contínuo podem
ser descritos por um sistema de equações algébrico-diferenciais (EAD): as
equações diferenciais descrevem o comportamento dinâmico enquanto que as
equações algébricas impõem relações físicas do sistema. Uma vez localizado
um evento, um novo conjunto de EAD’s e/ou variáveis passa a descrever o
sistema, e um novo conjunto de condições iniciais deve ser gerado para
continuação da simulação, que tipicamente resultam da condição final do
intervalo anterior e de condições decorrentes do próprio evento (Gopal e
Biegler, 1997). Três aspectos importantes devem ser considerados quando se
trata de problemas dinâmicos não suaves:
1. Localização de eventos nas variáveis de estado, no intervalo de
integração;
2. Identificação do instante de ocorrência do evento (emprego de rotina de
interpolação);
3. Reinicialização da integração a partir do instante da ocorrência do
evento;
Quando a integração é conduzida, o problema fica “preso” no caso
contínuo e nenhuma descontinuidade é permitida, até que uma função
comutadora (switching function) localiza um evento. Estas funções são
projetadas de tal forma que os seus zeros correspondam aos pontos de
descontinuidade. Uma vez localizado o evento, o tempo correspondente é
encontrado por alguma rotina de interpolação. A integração prossegue a partir
deste ponto, com a nova estrutura do modelo e com novas condições iniciais. A
questão crucial é se o conjunto de EAD’s anteriores ao evento pode fornecer a
condição inicial para o intervalo posterior ao evento.
Na simulação de processos químicos, é muito comum a ocorrência de
descontinuidades implícitas decorrentes do cálculo de equilíbrio de fases, como
na simulação de “flashes”. O problema básico é que o número de fases não é
conhecido a priori.
REVISÃO BIBLIOGRÁFICA
6
2.2. Simulação de “Flashes”
Vasos para o armazenamento e separação de fluidos são dos
equipamentos mais estudados no que diz respeito à modelagem do
comportamento dinâmico. A formulação deste tipo de problema aparece em
diversos livros-textos (Luyben, 1989; Raman, 1985). Entretanto, observa-se
também que, na maioria dos estudos de simulação dinâmica destes
equipamentos, a modelagem termodinâmica é muito simplificada, empregandose aproximações tais como a de comportamento ideal das fases. Tais
aproximações reduzem, em larga medida, a complexidade dos sistemas de
equações do modelo e, via de regra, simplificam sua resolução, distanciando,
muitas vezes, os valores preditos dos valores reais. Neste contexto, optou-se
por delimitar a revisão bibliográfica dos textos de simulação dinâmica apenas
aos trabalhos que empregam uma formulação termodinâmica mais rigorosa
para a simulação de tambores de “flash”.
No presente trabalho, conforme está apresentado no Capítulo 3,
referente à formulação do problema, a parte algébrica do sistema de EAD’s
abordado corresponde às equações de um “flash” com valores especificados
para a energia interna (U), volume (V) e número de mols (N) no equipamento
(“flash” UVN). O equilíbrio termodinâmico corresponde à configuração
macroscópica que maximiza a entropia do sistema, sobre os estados que têm
os valores especificados de U, V e N (Callen, 1985). Entretanto, o uso direto de
uma formulação baseada na maximização da entropia é inconveniente pois, em
geral, as propriedades termodinâmicas necessárias à simulação não são
funções explícitas de U, V e N. Isto ocorre porque a maior parte das equações
de estado em uso corrente em Engenharia Química permite obter as
propriedades residuais em termos de temperatura, volume e frações molares
ou, alternativamente, temperatura, pressão e frações molares. Por este motivo,
este tipo de “flash” foi estudado em três trabalhos recentes (Saha e Carroll
(1997); Müller e Marquardt (1997); Michelsen (1999)).
No trabalho de Saha e Carroll (1997) o “flash” isoenergético-isocórico
(UVN) é utilizado para simulação dinâmica de vasos onde é ressaltado que
nenhuma das variáveis intensivas (T e P) são conhecidas. Este fato produz
REVISÃO BIBLIOGRÁFICA
7
maior dificuldade de resolução do que em “flashes” de outros tipos. No modelo,
além das equações de balanço são necessárias as relações de equilíbrio. Os
autores empregaram um método de convergência “super-linear” para resolver o
sistema de equações não-lineares resultantes da formulação, que, criticamente
requer uma boa condição inicial. Como nenhuma das variáveis intensivas são
conhecidas, um esquema de geração de condições iniciais é crucial para a
solução do problema. Os autores ressaltam que uma solução trivial para este
problema de “flash” é bastante usual, particularmente quando uma equação de
estado é utilizada para cálculos de propriedades tanto na fase líquida quanto
na fase vapor. O uso de análise de estabilidade é apresentado para evitar este
problema. Observa-se, ainda, que a formulação apresentada neste artigo está
limitada a problemas bifásicos e que a existe a necessidade de resolver a
equação de estado para cada fase, em cada iteração e passo de integração.
Müller e Marquardt (1997) estudaram o efeito dinâmico de “flashes”
multifásicos, comparando o efeito de usar esquemas simplificados para
detectar o surgimento de fases versus a utilização do teste global de
estabilidade de fases. Dois tipos de especificações foram consideradas: UVN e
TPN, mas escassos detalhes são apresentados a respeito da formulação
adotada e da metodologia de resolução. Todos os exemplos apresentados no
artigo modelam as fases líquidas usando expressões para a energia livre de
Gibbs em excesso, ao contrário do trabalho desenvolvido nesta tese, que utiliza
uma equação de estado para modelar as fases fluidas presentes. Esta é uma
distinção importante porque o uso de equações de estado introduz dificuldades
adicionais devido à grande dependência da pressão com o volume molar de
fases líquidas.
Michelsen (1999) apresenta uma variedade de formulações para a
resolução de “flashes” de importância prática, cujas especificações são,
contudo, pouco usuais ou ainda escassamente estudadas. Conforme é bem
conhecido na literatura de termodinâmica, o autor ressalta que a resolução de
“flashes” TP corresponde a localizar o mínimo global da energia livre de Gibbs
da mistura. Adicionalmente, especificações de PH, PS, TV ou SV também
permitem a seleção de funções termodinâmicas de estado para as quais o
REVISÃO BIBLIOGRÁFICA
8
mínimo global pode ser localizado. São apontadas como vantagens da
abordagem de minimização: (a) a solução desejada é única, (b) análise de
estabilidade pode ser usada para verificar a consistência da solução e
determinar o número de fases em equilíbrio. Exceto para o “flash” TP, as
demais especificações citadas recaem em um problema de minimização (da
função termodinâmica de estado) com restrições não lineares. Duas soluções
são apresentadas. Na primeira, uma malha externa de otimização tem como
variáveis de decisão T e P, e em malha interna é resolvido o problema de
“flash” TP. Na segunda abordagem, uma função objetivo modificada é proposta
onde as restrições são removidas, mas o problema de ponto de mínimo é
substituído por um problema de ponto de sela. As equações resultantes são
resolvidas através de um método de Newton global. Nesta tese, utiliza-se a
segunda abordagem apresentada. Vale ressaltar que Michelsen (1999)
comenta, nas conclusões do seu trabalho, que a formulação do “flash” UVN
iterando diretamente sobre o volume das fases e sobre os números de mols de
cada composto em cada fase, isto é, sem a necessidade de resolver a equação
de estado, merecia estudos adicionais. Esta é a formulação adotada nesta
tese.
Para ressaltar a importância do cálculo de propriedades termodinâmicas
em um ambiente amigável ao usuário, o produto comercial MULTIFLASH
(InfoChem, 2001) de simulação dinâmica de “flashes”, oferece uma ampla
variedade de opções de especificações: PT, PH, TH, PS, TS, PV, TV, PU, TU,
Pf, Tf, UV e HS, onde: P = pressão, T = temperatura, H = entalpia, S= entropia,
V = volume, U = energia interna e f = fração molar. O software utiliza o
ambiente MATLAB®, e calcula diversas propriedades físicas e de transporte:
composição, fase, temperatura, pressão, volume, entalpia, entropia, energia
interna, energia livre de Gibbs, fator de compressibilidade, peso molecular
médio, capacidades caloríficas, viscosidade, condutividade térmica e tensão
superficial. Várias equações de estado estão disponíveis para seleção.
9
REVISÃO BIBLIOGRÁFICA
2.3. Propriedades Termodinâmicas
Embora haja diversas equações de estado propostas na literatura para a
modelagem de sistemas fluidos de hidrocarbonetos, optou-se por realizar todas
as simulações com a equação de estado de Peng-Robinson (1976), que é uma
das mais empregadas na modelagem de tais sistemas. Sendo assim, não faz
parte do escopo deste trabalho uma análise da formulação e do desempenho
de equações de estado. Por este motivo, a revisão bibliográfica referente ao
cálculo
de
propriedades
termodinâmicas
se
limita
ao
aspecto
de
implementação computacional de modelos com o auxílio de técnicas
computacionais modernas, tais como diferenciação automática e computação
algébrica (CA), com ênfase nesta última.
O cálculo de propriedades físicas para projeto, análise e controle de
processos requer, freqüentemente, a execução de expressões longas,
passíveis de erro de implementação quando realizadas manualmente. Esta
situação gera o incentivo para uso de CA. Um emprego de processamento
algébrico de particular interesse é a obtenção de derivadas analíticas de
modelos termodinâmicos (Bischof et al., 1992). Vieira (2001) recomenda o
emprego de diferenciação automática na formulação de problemas de EAD’s.
O uso de CA na implementação de procedimentos de cálculo de
propriedades termodinâmicas foi anteriormente empregado (Taylor e Monagan,
1993; Taylor, 1994; Adams e Taylor, 1994; Silva, 1993; Silva e Castier, 1993b),
utilizando modelos relativamente simples para cálculo de propriedades de
misturas com número de componentes conhecido a priori. Silva e Castier
(1993a) e Taylor (1997) desenvolveram procedimentos de manipulação
algébrica de modelos termodinâmicos que relaxaram esta limitação.
Os programas de CA como Maple® (Char et al., 1991) e Mathematica®
(Wolfram, 1991) realizam não apenas as tarefas de manipulação algébrica mas
também as de operações numéricas e de visualização gráfica, em um ambiente
computacional completo. Contudo, as computações numéricas realizadas em
ambientes de CA são lentas quando comparadas com códigos em FORTRAN
ou C. Adicionalmente, há um forte interesse em acoplar novos modelos
termodinâmicos a programas pré-existentes para simulação de operações de
10
REVISÃO BIBLIOGRÁFICA
processo químico, criando incentivo para uso de pacotes de CA para escrever
procedimentos que implementem modelos termodinâmicos em linguagem
compilável (Castier, 1999).
Um pacote para implementação automática de modelos termodinâmicos,
Thermath, na linguagem de programação Mathematica® foi desenvolvido por
Castier (1999). O Thermath, a partir de um modelo de energia livre de Gibbs de
excesso ( G E ) ou de uma equação de estado, pode ser empregado para gerar
expressões para várias propriedades termodinâmicas e analisar a estrutura das
expressões resultantes, gerando código que as implementa em linguagem de
computação (FORTRAN77, por exemplo).
2.4. Resolução de Equações Algébrico-Diferenciais
Com a introdução de modelos termodinâmicos rigorosos, envolvendo
testes de estabilidade de fases, a complexidade do sistema de equações
algébrico-diferenciais
(EAD’s)
resultante
requer
códigos
numéricos
especialmente desenvolvidos para esta finalidade.
O código DASSL emprega métodos BDF (Backward Differentiation
Formulas) para resolver um sistema de EAD’s ou de equações diferenciais
ordinárias (EDO’s). Os métodos são de passo e ordem variável. O sistema de
EAD’s no DASSL é escrito na forma:
f (t , y , y ' ) = 0
(2.1)
onde y’ denota a derivada temporal de y. Os métodos BDF empregados no
código DASSL requerem a resolução de sistemas não-lineares de grande
dimensão, em cada passo de tempo:
f (t n , y n ,α n y n + β n ) = 0
(2.2)
onde α n e β n são escalares que dependem do método e do passo. O sistema
da Eq. (2.2) é resolvido por uma iteração de Newton modificada, que requer, a
cada iteração, a solução do problema linear:
11
REVISÃO BIBLIOGRÁFICA
Ay n( k +1) = bn( k )
(2.3)
onde a matriz A é dada por:
A = αn
∂f
∂f
+
∂y ' ∂y
(2.4)
que é resolvida por um método específico para matrizes em faixas (banded
method), sensível à largura da faixa. Ou seja, para problemas com faixa larga
(como para problemas de equações diferenciais parciais bi-dimensionais
resolvidas pelo método das linhas), o resolvedor do sistema apresentado na
Eq. (2.3) se torna bastante ineficiente. Brown et al. (1994) apresentou um novo
resolvedor de EAD’s, o DASPK, que utiliza um método GMRES para resolução
do sistema linear, que não requer a matriz A mas sim o produto Av, aproximado
pela
Eq. (2.5):
Av ≅ f (t , y + δv ,α ( y + δv ) + β ) − f (t , y ,αy + β )δ
Neste
trabalho,
optou-se
por
(2.5)
empregar
o
MATLAB®/SIMULINK® pelas ferramentas nele oferecidas para
ambiente
resolução
numérica e visualização gráfica, incluindo código para a resolução de EAD’s. O
código, ode15s, baseia-se em uma variante do algoritmo BDF, o NDF
(numerical differentiation formula), desenvolvido para integrar EDO’s rígidas da
forma (Shampine e Reichelt, 1997; Shampine et al., 1999):
M (t )y ' = f (t , y )
(2.6)
Quando a matriz de massa M(t) é singular, trata-se de um sistema de
EAD’s. De acordo com Shampine et al. (1999), o cálculo de um novo passo
com NDF ou BDF não requer que a matrix M(t) seja não-singular, como ocorre
nos códigos DASSL, LSODI ou SPRINT.
A maioria dos códigos disponíveis para resolução de EAD’s impõe a
necessidade de uma condição inicial consistente. O ode15s oferece, a exemplo
do DASSL e SPRINT, um código para cálculo de condições iniciais. Petzold
(1992) observa que a crítica mais comum dos usuários do DASSL e SPRINT é
12
REVISÃO BIBLIOGRÁFICA
a falta de robustez destes códigos na geração de condições iniciais. Nesta
direção, o trabalho de Vieira (2001) apresenta diversas sugestões para
inicialização em problemas de engenharia química.
2.4.1. O código ode15s
Os códigos BDF têm emprego bastante difundido na resolução de
problemas rígidos (stiffs). Quando são usados passo de tempo constante (h) e
diferenças para trás (backward differences), a fórmula de ordem k, BDFk, para
um passo de (tn,yn) para (tn+1,yn+1) é dada pela Eq. (2.7).
k
1 m
∇ y n +1 − hF (t n +1; y n +1 ) = 0
m =1 m
∑
A equação algébrica para yn+1
(2.7)
é resolvida pelo método iterativo de
Newton, inicializado como:
k
y n(0+)1 = ∑ ∇ m y n
(2.8)
m =1
Por outro lado, as fórmulas NDF (numerical differentiation formula)
(Shampine, 1997) são métodos da forma:
k
1 m
∇ y n +1 − hF (t n +1; y n +1 ) + κγ k ( y n +1 − y n(0+)1 ) = 0
m
m =1
∑
(2.9)
onde
k
1
j =1 j
γk = ∑
(2.10)
O código ode15s é uma implantação com passo de integração quase
constante do método NDF em termos de diferenças para trás. É fornecida a
opção de integração com a fórmula BDF, e ordens inferiores a 5 (o default).
Muitas das táticas adotadas no código se assemelham aos códigos como o
REVISÃO BIBLIOGRÁFICA
13
DIFSUB, DDRIV2, LSODE e VODE. O método de seleção do passo inicial
decorre da observação que é possível formar estimativas de um passo inicial
ótimo a partir das derivadas de F(t,y) em t = t0 (Shampine, 1997).
Os métodos para resolução de problemas rígidos envolvem as derivadas
parciais das funções que definem as equações diferenciais. Os códigos mais
popularizados permitem ao usuário fornecer as derivadas analíticas. O código
ode15s permite, alternativamente, a geração numérica da matriz Jacobiana.
Um aspecto da montagem das derivadas parciais é a opção de vetorizar,
quando a função para cálculo das derivadas é codificada de forma a devolver
uma matriz de valores da função. Uma característica do código é calcular
Jacobianas esparsas e cheias.
O código retém cópia da matriz Jacobiana, como explorado pelo VODE.
A taxa de convergência é monitorada e a iteração é interrompida se for predito
que a convergência não será atingida em 4 iterações. Neste evento, uma nova
matriz Jacobiana é formada, caso contrário, o passo é reduzido. Neste
esquema, o código formará poucas matrizes Jacobianas quando aplicado a um
problema
não-rígido, o que o torna competitivo com os códigos para problemas nãorígidos, em parte, também, devido à eficiente álgebra linear do MATLAB®
(Shampine, 1997).
Adicionalmente, Shampine e Reichelt (1997) ressaltam que o código
ode15s do MATLAB® apresenta-se com o objetivo de tornar a resolução de
uma EAD de índice 1 o mais semelhante possível à resolução de uma EDO.
Desta forma, o código ode15s reconhece uma EAD automaticamente, e aceita
o vetor y0 como uma estimativa inicial, computando um conjunto de condições
y ) consistentes, resolvendo a equação algébrica não-linear:
iniciais ( ~
M (t 0 )(
y~ − y 0
) = f (t 0 , y~)
h
(2.11)
onde h é um passo no tempo suficientemente pequeno. ~
y é um conjunto
consistente de condições iniciais com
14
REVISÃO BIBLIOGRÁFICA
y~ − y 0
y~' = (
)
h
(2.12)
No cômputo da condição inicial, apenas os valores algébricos são
alterados.
Adicionalmente, dado o objetivo de manter semelhança de uso entre as
aplicações EDO/EAD, o MATLAB® estende a opção de “localização de
eventos” para a solução de equações algébrico-diferenciais. Ressalta-se que
este recurso é de grande valia no tratamento de descontinuidades, como
apresentado por Gopal e Biegler (1997).
2.4.2. Resolução de EAD’s no SIMULINK®
O ambiente SIMULINK® é uma ferramenta largamente empregada na
resolução de modelos físicos especificados em uma linguagem direcionada por
diagramas de blocos. Um “loop” algébrico em um diagrama de blocos é uma
série de blocos conectados de tal forma que as saídas afetam as entradas. Em
termos matemáticos, estes “loops” definem uma EAD. As EAD’s formuladas em
SIMULINK devem ter a forma semi-explícita (de índice 1):
u ' = f1(t , u, v )
(2.13)
0 = f2 (t , u, v )
A abordagem indireta (semi-explícita) é chamada por Shampine et al.
(1999) de “Abordagem EDO”. As equações são resolvidas no intervalo [t0,tf ]
com condição inicial u(t0) = u0
e uma estimativa inicial para v(t0) = v0 .
Assumindo que a equação algébrica 0 = f2 (t 0 , u0 , v ) tem solução V próxima a
v0 , o requisito básico é que a matriz Jacobiana
∂f2
seja não-singular na região
∂v
contendo (t0,u0,V), o que torna o problema de índice 1. A resolução implícita
requer a existência de uma função R(t,u) tal que R(t0,u0) = V e
0 = f2 (t 0 , u0 , R(t , u )) . As equações diferenciais são então:
15
REVISÃO BIBLIOGRÁFICA
u ' = f1(t , u, R(t , u )) = f (t , u )
(2.14)
que é uma EDO. A função f(t,u) envolve a resolução de um sistema de
equações algébricas não-lineares, que é um aspecto computacional que dá
complexidade à abordagem. Por outro lado, permite acesso a qualquer rotina
de resolução de EDO’s do MATLAB® e do SIMULINK®, que incluem, entre
outros, dois códigos Runge-Kutta explícitos e um Adams-Bashforth-Moulton,
para problemas não rígidos, e um código NDF, um código Rosenbrock
modificado e um código de regra trapezoidal para sistemas rígidos.
A “Abordagem EDO” é a solução “natural” no ambiente SIMULINK®, que
só está apto a resolver EAD’s semi-explícitas na forma acima (Shampine et al.,
1999), representado na Fig. 2.1, onde ui são saídas de blocos integradores e vi
saídas de blocos não integradores. Como o fluxo de sinais em um diagrama de
blocos
SIMULINK®
é
direcional,
os
blocos
podem
ser
ordenados
topologicamente. Os componentes extremamente conectados representam as
restrições algébricas. A primeira etapa ao resolver 0 = f2 (t , u,v ) para a variável
v
é
particionar
em
sub-conjuntos,
que
podem
ser
resolvidos
independentemente em ordem topológica, por corte nos nós. Cada subconjunto é resolvido com o código HYBRID1 do MINPACK.
Figura 2.1. Diagrama de Blocos de Simulador SIMULINK®
REVISÃO BIBLIOGRÁFICA
16
É importante ressaltar que testes feitos no ambiente MATLAB®
(Shampine et al., 1999) indicam que a “Abordagem EDO” não é competitiva
quando comparada com a abordagem direta. Contudo, não é desprezível a
vantagem apontada anteriormente de que a abordagem indireta abre a
possibilidade de uso de uma grande variedade de códigos de resolução de
EDO’s.
2.5. Da Computação Algébrica à Simulação Dinâmica Rigorosa de
Processos: O Problema de “Flash” Multifásico
Neste capítulo, pretendeu-se abordar os diversos aspectos envolvidos
na simulação dinâmica de processos químicos. O problema abordado neste
estudo aproxima as questões expostas com a proposta de unir a modelagem
termodinâmica rigorosa às ferramentas de resolução numérica disponíveis para
o tratamento numérico em ambiente de programação amigável.
Dentre os aspectos relevantes do trabalho, destacam-se:
1. O equipamento de “flash” é abordado com especificação UVN. A
modelagem envolve uma formulação na qual a condição de equilíbrio
corresponde a um ponto de sela. A formulação adotada nesta tese foi
recentemente (Michelsen, 1999) apontada como merecedora de
estudos;
2. O problema de modelagem termodinâmica é abordado com rigor,
envolvendo teste de estabilidade de fases. A possibilidade de
surgimento ou desaparecimento de fases ao longo da simulação origina
uma descontinuidade implícita que requer eficiência do código de
resolução numérica;
3. Emprego de resultados de computação algébrica para desenvolvimento
de códigos FORTRAN para cálculo de propriedades termodinâmicas.
Nesta tarefa, emprega-se o programa Thermath (Castier, 1999);
4. O código FORTRAN gerado pelo Thermath é integrado ao ambiente
MATLAB® através de procedimento de geração de DLL executável por
REVISÃO BIBLIOGRÁFICA
17
programa em linguagem de programação MATLAB®. Este recurso é
apresentado no Capítulo 3 – Formulação;
5. O resolvedor de EDA’s/EDO’s, ode15s, do MATLAB®. O uso do código
DASSL não foi explorado neste trabalho já que havia opção de
resolução com as próprias ferramentas internas do MATLAB®. A opção
pelo ambiente MATLAB® é que orientou, neste trabalho, a seleção pelo
código ode15s e não o reconhecimento de qualquer superioridade
numérica deste sobre o DASSL;
6. O ambiente SIMULINK® foi explorado como interface amigável para o
problema proposto. Neste, o bloco S-function permite invocar as rotinas
geradas pelo Thermath e o código ode15s, entre outros.
3
Formulação
Neste capítulo, apresenta-se o modelo do equipamento de “flash” com
energia interna (U), volume (V) e número de mols (N) especificados (“flash”
UVN). As propriedades termodinâmicas necessárias para a simulação rigorosa
do equipamento são obtidas por computação algébrica (CA) através do
programa
Thermath
(Castier,
1999)
no
ambiente
Mathematica®. Tal
procedimento gera um código em FORTRAN que pode ser utilizado pelo
simulador. A dedução das expressões algébricas responsáveis pela geração do
código referente à modelagem do “flash” UVN encontra-se no Apêndice 1
(“Montagem das Equações Algébricas de Equilíbrio de Fases”). No Apêndice 2
(“Cálculo das Propriedades Termodinâmicas”), está apresentada a geração das
equações algébricas e das propriedades necessárias para o cálculo do
equilíbrio termodinâmico e teste de estabilidade de fases.
Neste
trabalho,
explora-se
a
possibilidade
de
migrar
códigos
desenvolvidos em linguagem FORTRAN para um ambiente de computação
numérica e visualização gráfica, como o MATLAB®, através de uma interface
desenvolvida.
3.1. “Flash” UVN
Considera-se um tambor de “flash” com s f correntes de entrada e s w
correntes de saída. Para tambores com apenas correntes de entrada (ou
saída), s w (ou s f ) é zero e a mesma formulação é mantida. O volume do
tambor é constante e conhecido, igual a V. Os balanços de massa e energia
são descritos nas Eqs. (3.1) e (3.2):
20
FORMULAÇÃO
s
sw
f
dU
= ∑ H jf • − ∑ H kw• + Q − WS
dt
k =1
j =1
(3.1)
s
sw
f
dN i
= ∑ f ij• − ∑ wik•
dt
k =1
j =1
i = 1,..., nc
(3.2)
Nestas equações, t é o tempo, U é a energia interna, e Ni é o número
de mols do componente i no tanque. A entalpia (por unidade de tempo) da
corrente j de entrada é representada por H jf • e f ij• é a vazão molar do
componente i na corrente de entrada j. Os símbolos H kw• e wik• denotam
valores análogos para a corrente de saída k, e Ws é a potência de eixo, nulo no
caso de um tambor de “flash”. A carga térmica adicionada ao tambor é
representada por Q e nc é o número total de componentes no sistema.
O sistema resultante é composto de
(nc + 1)
equações diferenciais
ordinárias (EDO’s) representadas pelas Eqs. (3.1) e (3.2) que, por integração,
fornecem a trajetória de U e N . Assim, possíveis mudanças no número de
fases em equilíbrio no interior do “flash” não afetam em momento algum o
sistema
de
equações
diferenciais,
ou
seja,
o
número
de
estados
correspondentes às equações diferenciais permanece constante ao longo da
integração.
Sob integração, U , V e N serão conhecidos a qualquer passo de
tempo. Admite-se que, em qualquer intervalo, o material no interior do tambor
esteja em equilíbrio termodinâmico. Com esta hipótese, as condições de
equilíbrio podem ser determinadas pela maximização da entropia para um dado
conjunto de valores UVN .
As condições de máxima entropia fornecem um conjunto de equações
algébricas que devem ser resolvidas a cada passo de integração. No entanto,
na maioria das aplicações práticas, as equações de estado (EDE) são da forma
P = P (T ,V , N ) . Com isso, a formulação de entropia torna-se pouco conveniente
21
FORMULAÇÃO
dado que as propriedades residuais obtidas a partir das EDE são também
funções de T , V e N.
Recentemente, Michelsen (1999) formulou vários problemas de “flash”
menos usuais. No caso de especificações UVN, ele sugeriu o uso de uma
função QF que, convenientemente diferenciada, fornece as condições de
equilíbrio, conforme será apresentado na seqüência. A função QF é dada por:
QF =
A − U ESPEC
RT
(3.3)
onde A é a energia livre de Helmholtz, U ESPEC é o valor especificado de energia
interna e R é a constante universal dos gases. No caso de problemas de
equilíbrio líquido-vapor em “flash”, Michelsen (1999) sugere a diferenciação de
QF em relação ao logaritmo natural da temperatura e do volume de uma das
fases (a fase vapor, por exemplo), e também em relação ao número de mols de
cada componente em uma das fases. Nesta tese, essa formulação foi tentada
mas mostrou-se de difícil convergência quando uma das fases apresentava-se
em pequena quantidade, isto é, muito próxima ao ponto de saturação, tornando
inconveniente a manipulação de funções logarítmicas.
Optou-se, portanto, por trabalhar com as derivadas diretamente em
relação à temperatura, aos volumes de cada fase, e manteve-se a sugestão
quanto ao número de mols. A dedução das equações de modelagem do
equipamento encontra-se no Apêndice 1.
Observou-se que esta segunda abordagem teve mais sucesso. Contudo,
longe do ponto de saturação, as duas opções apresentaram desempenhos
semelhantes e convergiram com o mesmo número de iterações. Assim, para
um sistema contendo np fases, foi obtido o seguinte conjunto de equações
algébricas constituídas das derivadas de QF :
np
⎛ ∂Q ⎞
rT = ⎜ F ⎟ =
⎝ ∂T ⎠V ,N
U
ESPEC
− ∑ Uk
k =1
RT 2
=0
(3.4)
22
FORMULAÇÃO
⎛ ∂Q
rV = ⎜⎜ F
⎝ ∂Vk
⎞
P − Pk
⎟⎟
= K
=0
RT
⎠TN,Vm ≠ k ,K
⎛ ∂Q
rN = ⎜ F
⎜ ∂n
⎝ ij
⎞
µ − µ iJ
⎟
= ij
= 0 , i = 1,..., n c
⎟T ,V
RT
⎠ nkm,k ≠ i ,m ≠ j ,J
k = 1,..., n p
k ≠K
(3.5)
j = 1,..., n p ; j ≠ J
(3.6)
Nestas equações, U k e Pk são a energia interna e a pressão na fase k,
µ ij é o potencial químico do componente i na fase j . O problema foi
formulado admitindo que volume de uma das fases (fase K) é uma variável
dependente, calculada a partir da equação de conservação de volume. Por
conveniência numérica, a fase com maior volume, em cada iteração, foi
admitida como a fase dependente. De forma semelhante, a fase com o maior
número de mols para cada componente foi admitida como a variável
dependente. As Eqs. (3.4 - 3.6) representam um conjunto de (n p + n c (n p − 1))
equações algébricas não-lineares. No caso de existir uma única fase no
sistema, o conjunto reduz-se a uma única equação (Eq. 3.4). As derivadas
parciais necessárias para o cálculo da matriz Jacobiana do conjunto de
equações, que estão deduzidas no Apêndice 1, estão apresentadas a seguir.
Deve-se ressaltar que uma das vantagens desta formulação é que a matriz
Jacobiana do conjunto de equações é a matriz Hessiana de QF ; assim, esta
matriz
é
simétrica
e
é
necessário
calcular
apenas
[(n p + n c (n p − 1)) ∗ (n p + 1 + n c (n p − 1))] / 2 de seus termos, reduzindo o esforço
computacional. O sistema resultante de equações algébricas não lineares foi
resolvido usando-se o método de Newton-Raphson, cujo passo completo, em
cada iteração, é dado pela resolução do seguinte sistema de equações
lineares:
⎡r
⎢ TT
⎢ rVT
⎢
⎢⎣rNT
T
rTV
rVV
rNV
T ⎤
⎡ rT ⎤
rTN
⎡ ∆T ⎤
⎥ ⎢ ⎥
⎢ ⎥
rVN ⎥ ⋅ ⎢∆V ⎥ = − ⎢ rV ⎥
⎥
⎢r N ⎥
rNN ⎥ ⎢⎣∆N ⎥⎦
⎣ ⎦
⎦
(3.7)
23
FORMULAÇÃO
Eventualmente, o passo completo de Newton-Raphson calculado pela
resolução da Eq. (3.7) leva uma ou mais variáveis a valores não-físicos como,
por exemplo, volumes ou números de mols negativos. Na implementação desta
tese, quando isto ocorreu, reduziu-se o tamanho de passo, preservando-se a
direção de busca.
No caso da dinâmica do “flash”, devido à integração das equações
diferenciais, resolve-se em cada instante de tempo, um problema de “flash”
UVN com especificações similares aos “flashes” precedentes. Observou-se que
o uso da solução do “flash” precedente constitui, em geral, uma ótima
estimativa inicial para um novo problema de “flash”. Extrapolações lineares das
duas últimas soluções de problemas precedentes também foram tentadas, com
bons resultados.
Os termos que aparecem na Eq. (3.7) são dados por:
2(U − U ESPEC )
1
⎛ ∂r ⎞
−
rTT = ⎜ T ⎟ =
3
RT
RT 2
⎝ ∂T ⎠V ,N
rVV
⎛ ∂r
=⎜ V
⎜ ∂V
⎝ p
⎞
1
⎟
=−
⎟T ,V ,m ≠ p,P ,
RT
⎠N m
⎛ ∂r
rNN = ⎜⎜ N
⎝ ∂N rs
⎡⎛
⎢⎜ ∂PF
⎢⎜⎝ ∂VF
⎣
⎞
1 ⎡⎛⎜ ∂µ ij
⎟⎟
=
⎢⎜
⎠Tnab,V,,a ≠ r , RT ⎢⎣⎝ ∂N rj
b ≠ s ,S
⎛ ∂U ⎞
⎜
⎟
⎝ ∂T ⎠V ,N
⎞
⎛ ∂P
⎟⎟
− δ f ,p ⎜⎜ f
⎠TN,Vm ,m ≠ p,P ,
⎝ ∂Vf
⎞
⎛
⎟(δ js − δ jS ) + ⎜ ∂µ iJ
⎜ ∂N
⎟
⎝ rJ
⎠
⎡⎛
⎢⎜ ∂U p
⎢⎜ ∂Vp
⎣⎢⎝
rTV
⎛ ∂ rT
=⎜
⎜ ∂V
⎝ p
⎞
1
⎟
=−
⎟T ,V ,m ≠ p,P ,
RT 2
⎠N m
rTN
⎛ ∂r
= ⎜⎜ T
⎝ ∂N rs
⎞
1
⎟⎟
=−
RT 2
⎠TN,abV ,,a ≠ r ,b ≠s,S
rVN
⎛ ∂r
= ⎜⎜ V
⎝ ∂N rs
⎡
⎞
1 ⎢⎛ ∂Pf
⎟⎟
⎜⎜
=
⎠TN,abV ,,a ≠ r ,b ≠s,S RT ⎢⎣⎝ ∂N rf
(3.8)
⎤
⎞
⎟⎟(δ Js − δ JS )⎥
⎥⎦
⎠
⎞
⎛ ∂U
⎟
− ⎜⎜ P
⎟T ,V ,m ≠ p,P ,
⎝ ∂VP
⎠N m
⎡⎛
⎢⎜ ∂U s
⎢⎜⎝ ∂N rs
⎣
⎤
⎞
⎥
⎟⎟
⎥
T
,
V
,
m
≠
p
,
P
,
m
⎠N
⎦
(3.10)
⎤
⎞
⎥
⎟⎟
⎠TN,Vm ,m ≠ p,P , ⎥⎥
⎦
⎛ ∂U S
⎞
⎟⎟
− ⎜⎜
⎠TN,abV ,,a ≠ r ,b ≠ s,S ⎝ ∂N rS
⎤
⎞
⎥
⎟⎟
T
,
V
,
⎠ Nab ,a ≠ r ,b ≠s,S ⎥⎦
⎞
⎛
⎟⎟
(δ fs − δ fS ) − ⎜⎜ ∂PF
⎠TN,abV ,,a ≠ r ,b ≠s,S
⎝ ∂N rF
(3.9)
(3.11)
(3.12)
⎤
⎞
⎟⎟
(δ Fs − δ FS )⎥
⎥
T
,
V
,
⎠ Nab ,a ≠ r ,b ≠s,S
⎦
(3.13)
24
FORMULAÇÃO
1
⎛ ∂r ⎞
(PF − Pf ) + 1
rVT = ⎜ V ⎟ = −
2
RT
RT
⎝ ∂T ⎠V ,N
(µ iJ − µ ij ) 1
⎛ ∂r ⎞
rNT = ⎜ N ⎟ =
+
RT
RT
⎝ ∂T ⎠V ,N
rNV
⎛ ∂r
=⎜ N
⎜ ∂V
⎝ p
⎡
⎞
1 ⎢⎛⎜ ∂µ ij
⎟
=
⎟T ,V ,m ≠ p,P , RT ⎢⎜ ∂V
⎠N m
⎢⎣⎝ j
⎡⎛ ∂µ ij
⎢⎜⎜
⎢⎣⎝ ∂T
⎡⎛ ∂PF ⎞
⎛ ∂P ⎞ ⎤
⎟ −⎜ f ⎟ ⎥
⎢⎜
⎣⎢⎝ ∂T ⎠V ,N ⎝ ∂T ⎠V ,N ⎦⎥
(3.14)
⎞
⎛ ∂µ ⎞ ⎤
⎟⎟ − ⎜ iJ ⎟ ⎥
⎠V ,N ⎝ ∂T ⎠V ,N ⎥⎦
(3.15)
⎞
⎟
(δ − δ jP ) − ⎛⎜⎜ ∂µ iJ
⎟T ,V ,m ≠ p,P , jp
⎝ ∂VJ
⎠N m
⎤
⎞
⎟⎟
(δ Jp − δ JP )⎥⎥
,
,
≠
,
,
T
V
m
p
P
⎠N m
⎥⎦
(3.16)
Uma vez que as Eqs. (3.8) a (3.16) constituem os elementos da matriz
Hessiana de QF deve existir simetria entre diversos termos. Assim, as
submatrizes rNN e rVV devem ser simétricas. Além disso, rTV (Eq. (3.11)) e rVT
(Eq. (3.14)), bem como rTN (Eq. (3.12)) e rNT (Eq. (3.15)) também devem
apresentar simetria. Todas estas simetrias foram testadas durante a
implementação e corretamente obtidas pelo programa computacional.
Como já mencionado, admite-se que o tambor esteja em equilíbrio
termodinâmico a qualquer tempo. Para detectar o possível surgimento de uma
nova fase no tambor, o teste de estabilidade global de fase é empregado
(Michelsen, 1982a). Neste caso, o sistema de EAD’s é modificado para levar
em consideração as equações algébricas relacionadas à existência desta nova
fase. Admite-se o desaparecimento de uma fase sempre que a quantidade
desta, durante a integração do sistema de EAD’s, situa-se abaixo de um
patamar pré-estabelecido. Neste caso, o sistema de EAD’s é modificado pela
remoção das equações algébricas correspondentes. Ressalta-se que apesar
do surgimento ou desaparecimento de fases não modificar o número de
equações
diferenciais,
este
tem
impacto
no
número
de
estados
correspondentes às equações algébricas.
Para o cálculo da condição inicial do tambor, emprega-se o
procedimento desenvolvido por Espósito (1999) e Espósito et al. (2000). Este
FORMULAÇÃO
25
procedimento corresponde à solução de um problema de “flash” multifásico em
estado estacionário, no qual as especificações são a temperatura (T), o volume
(V) e o número de mols (N) de cada componente na mistura (“flash” TVN). Para
o cálculo das propriedades físicas das correntes de entrada, que pode conter
uma ou mais fases, usou-se um procedimento de “flash” (Michelsen, 1982b)
com especificação de temperatura (T), pressão (P) e número de mols (N) de
cada componente.
As propriedades físicas foram calculadas com a EDE de Peng-Robinson
(1976) com as regras de mistura de um fluido de van der Waals (com
parâmetros de interação binária fixados em zero para todos os exemplos). As
expressões para cálculo de propriedades foram obtidas automaticamente com
o programa de computação algébrica Thermath (Castier, 1999) e estão
apresentadas no Apêndice 2. As propriedades de compostos puros
necessárias para caracterizar cada substância foram obtidas em Reid et al.
(1987).
3.2. Implementação Computacional
O fluxograma da Fig. (3.1) apresenta o algoritmo de resolução do “flash”
UVN.
A rotina desenvolvida por Espósito (1999) foi adaptada para interface do
MATLAB® (desenvolvimento de gateway, com pequenas modificações de
código para adequação à finalidade do problema). A DLL gerada é chamada
em código MATLAB® para obter as condições iniciais do problema (“flash” TVN,
código eqfuvn.for). A cada passo de integração, as equações algébricas são
avaliadas em código FORTRAN, que realiza o teste de estabilidade de fases, e
resolve as equações algébricas pertinentes, utilizando as rotinas de cálculo de
propriedades termodinâmicas desenvolvidas pelo procedimento de computação
algébrica do Thermath (Castier, 1999) (“flash” UVN, código uvnslog.for). A DLL
gerada para esta rotina é invocada pelo código de resolução de EAD/EDO´s
ode15s, no ambiente MATLAB®.
A primeira restrição quanto ao tratamento dos sistemas como EAD’s é a
determinação do índice. Por definição, o índice diferencial é o número mínimo
26
FORMULAÇÃO
de vezes que um subgrupo de sistema de EAD’s (ou equações derivadas dele)
precisa ser diferenciado, em relação à variável independente até ser
transformado em um sistema de EDO´s (Costa Jr., 2001).
O código ode15s do MATLAB® permite resolução simultânea de
equações algébricas e diferenciais de índice 1. Utilizando a ferramenta que
possibilita a resolução das EAD’s, todas as rotinas (gateways), necessárias
para a comunicação do MATLAB® com o FORTRAN foram construídas e
testadas. No entanto, o código ode15s diagnosticou o problema como de índice
superior a 1, restando, no escopo deste estudo, a “Abordagem EDO” descrita
no Capítulo 2, e incluída no Fluxograma da Fig. (3.1).
Adicionalmente, como mencionado no item 3.1, o surgimento ou
desaparecimento de fases não modifica o número de estados correspondentes
a
equações
diferenciais
mas
tem
impacto
no
número
de
estados
correspondentes a equações algébricas. Sob este aspecto, a “Abordagem
EDO” para resolução do conjunto de EAD’s apresenta uma vantagem
comparativa, tendo em vista que alterar as dimensões do vetor de estado e
gerar condições para a nova fase pode atribuir rigidez ao sistema em estudo.
Para inicialização do problema, recorreu-se ao “flash” TVN que necessita
das seguintes especificações:
1. número de componentes
2. propriedades críticas (temperatura e pressão) e fator acêntrico
3. densidade global (g/cm3)
4. peso molecular dos componentes presentes
5. coeficientes da equação (Reid et al., 1987) para cálculo da
capacidade calorífica a pressão constante no estado de gás ideal
6. temperatura (K)
7. número de mols iniciais por componente
27
FORMULAÇÃO
Figura 3.1: Algoritmo para Resolução Numérica
FORMULAÇÃO
28
Note-se que o conhecimento dos números de mols, dos pesos
moleculares e da densidade mássica global, permite determinar o volume.
3.2.1. Integração de Códigos FORTRAN ao Ambiente MATLAB®
A utilização do ambiente computacional MATLAB® oferece como
maiores atrativos a disponibilidade de rotinas numéricas para resolução das
EDO’s e um forte elenco de rotinas de visualização gráfica. Contudo, é
indiscutível que a experiência prévia no desenvolvimento de códigos
computacionalmente eficientes para cálculo de propriedades termodinâmicas, e
simulação de operações típicas de processos químicos, apoia-se no uso de
linguagens de programação compiláveis, com destaque, na Engenharia
Química, para o FORTRAN.
Na direção de integrar desenvolvimentos em FORTRAN ao ambiente
MATLAB®, são empregados códigos EXECUTÁVEIS em MATLAB (arquivos
MEX), isto é, que podem ser invocados no ambiente MATLAB®. Os
componentes de um arquivo MEX FORTRAN, ilustrados na Fig. (3.2), incluem:
a) uma rotina computacional que contém o código que se deseja
implementar com o arquivo MEX
b) uma rotina gateway que realiza a interface entre a rotina computacional
e o MATLAB® tendo como ponto de entrada a função mexFunction e
seus parâmetros (nle, nld, ple, pld), onde pld é uma matriz de
argumentos do lado direito do comando, nld é o número de argumentos
de entrada, ple é uma matriz de argumentos de saída do lado esquerdo,
e nple é o número de saídas do lado esquerdo. A gateway chama a
rotina computacional como subrotina.
29
FORMULAÇÃO
(Arquivo MEX)
MATLAB
Rotina Gateway
arq_2() arq_3()
arq_1()
Rotinas FORTRAN
>> [a, b, c] = func(d, e, f, g);
nle = 3
nld = 4
a pld
ple
b
c
d
e
f
g
número de lado esquerdo
Ponteiro para lado esquerdo
número de lado direito
ponteiro para lado direito
Figura 3.2. Estrutura MEX FORTRAN
FORMULAÇÃO
30
3.2.2. Código em ambiente MATLAB®
As rotinas DLL’s desenvolvidas foram utilizadas em código MATLAB®. A
ocorrência de uma nova fase (ou desaparecimento de uma fase) constitui-se
em uma descontinuidade e pode gerar dificuldade de convergência na
resolução das equações algébricas. Para proteger o sistema, desenvolveu-se
uma lógica através do recurso localizador de eventos disponível no código
ode15s. Sempre que é identificada a dificuldade de convergência (normalmente
associada ao evento de mudança de fase), a integração é interrompida
automaticamente e retrocede-se no tempo para, então, retomar a integração
com redução no passo máximo permitido. Após a normalização de
convergência, relaxa-se o passo máximo, de acordo com Tabela 3.1.
3.2.3. Código em ambiente SIMULINK®
Alternativamente, visando à obtenção de uma interface amigável e
recursos disponíveis para simples configuração (como controladores, seletores,
saídas gráficas e mostradores digitais, entre outros) através da construção
modular, desenvolveu-se uma versão para o ambiente SIMULINK®, com o
algoritmo da Fig. (3.1), conforme apresentado na Fig. (3.3).
O
“Tambor de Flash” utiliza o bloco “S-function” através da qual é
possível chamar programas escritos em linguagem MATLAB® ou Funções MEX
em FORTRAN (ou C). Ressalta-se que esta alternativa supera as limitações da
construção apresentada na Fig. (2.1). A Fig. (3.3) apresenta o recurso de
geração de máscaras para blocos, tornando a interface do programa mais
amigável.
31
FORMULAÇÃO
Tabela 3.1: Ajuste de Passo Máximo para Auxílio de Convergência das
Equações Algébricas
tspan=[t0;tfinal];
maxt=10;
opt=odeset('Events',@eventos,'MaxStep',maxt,'InitialStep',0.00001);
while t0 < tfinal
tspan = [t0; tfinal];
[t,y,te,ye,ie] = ode15s(@f,tspan,y0,opt);
y0=y(end,:)';
t0=t(end,:);
if ie= =1
opt=odeset('Events',@eventos,'MaxStep',maxt/100,…
'InitialStep',0.00001);
y0=y(end-1,:)';
t0=t(end-1,:);
M
else
opt=odeset('Events',@eventos,'MaxStep',maxt,…
'InitialStep',0.001);
y0=y(end,:)';
t0=t(end,:);
M
end
end
%*************************************************************
% FUNÇÃO LOCALIZADORA DO EVENTO DE NÃO CONVERGÊNCIA
%*************************************************************
function [value,isterminal,direction]=eventos(t,y)
global flg saidoloop
value = [(flg= =1) saidoloop];
isterminal=[1 1];
% Evento interrompe a integração
direction=[0 0];
% Evento em qualquer direção
Figura 3.3. Modelo em Ambiente SIMULINK®
FORMULAÇÃO
33
Como exemplo do elenco de recursos disponíveis para a construção de
modelos dinâmicos no ambiente SIMULINK®, foi elaborado um modelo,
apresentado na Fig. (3.4), que utiliza:
a)
bloco controlador PID (com janela de diálogo aberta para ilustrar a
interface de configuração);
b)
bloco seletor multi-portas, através do qual foi configurada uma lógica
que modifica o ganho do controlador de acordo com o surgimento de
fases;
c)
bloco Produto para promover a adaptação do ganho do controlador;
d)
bloco gerador de pulsos para pulsar o “setpoint” de temperatura do
tambor em torno de um valor médio, fixado através de um bloco
constante.
3.3. Cálculo de Propriedades Termodinâmicas
As expressões de propriedades termodinâmicas necessárias para a
formulação do problema em estudo, para uma equação de estado tal como a
de Peng-Robinson (1976), são longas e complexas.
O Thermath (Castier, 1999) permite a obtenção de expressões para o
cálculo de propriedades físicas como entalpia, entropia, coeficiente de
fugacidade, capacidade calorífica, etc., a partir de modelos termodinâmicos
definidos pelo usuário, através da geração de procedimentos automáticos em
linguagem FORTRAN. Embora a versão atual do programa esteja limitada à
geração de códigos nesta linguagem, possivelmente, poucas modificações
seriam necessárias para que tais procedimentos sejam gerados em outras
linguagens de programação como C ou em código para MATLAB®.
Figura 3.4. “Flash” com Controle de Temperatura
35
FORMULAÇÃO
A possibilidade de geração de expressões genéricas torna a utilização
do Thermath bastante atrativa. As equações independem, por exemplo, do
número de componentes presentes no sistema. O Thermath analisa todas as
expressões definidas identificando a existência de subexpressões interrelacionadas ou simétricas com a finalidade de estabelecer uma ordem
hierárquica para a resolução destas. Isto é bastante importante quando existe a
preocupação quanto ao esforço computacional, como por exemplo, quando se
deseja calcular a derivada do logaritmo da fugacidade em relação ao número
de mols do componente i em uma mistura de nc componentes. Conforme a Eq.
(3.17), esta derivada é simétrica:
⎛ ∂ ln fˆi
⎜
⎜ ∂n
j
⎝
⎛ ∂ ln fˆj
⎞
⎟
=⎜
⎟T ,P , ⎜ ∂n
i
⎠ nk ≠ j ⎝
⎞
⎟
i=1,...,nc j=1,..,nc
⎟T ,P ,
⎠ nk ≠ i
(3.17)
Sem ser considerada tal igualdade, seria necessário calcular nc2
equações. Ainda que o algoritmo implementado no programa Thermath
ocasionalmente falhe na identificação de tais simetrias em expressões muito
longas, ele permite, nos casos em que funciona perfeitamente, gerar o código
FORTRAN explorando a simetria, o que, neste caso, significa que o programa
gerado avaliaria somente (nc (nc + 1) 2) expressões.
A implementação dos modelos termodinâmicos através do Thermath,
basicamente, pode ser divida em duas partes. A primeira, que compreende
aproximadamente 17% do código total, realiza a dedução das expressões,
utilizando
funções
implícitas
do
Mathematica®
e
procedimentos
complementares. Na segunda parte do código, a análise quanto à eficiência da
implementação, descrita anteriormente, é realizada. Tendo em vista que foi
possível aplicar a versão atual do programa Thermath diretamente nesta tese,
sem necessidade de alterações, não se apresenta uma discussão sobre o
modo de funcionamento do programa. Tal discussão pode ser encontrada em
Castier (1999).
FORMULAÇÃO
36
3.3.1. Propriedades Termodinâmicas
Geralmente, as propriedades termodinâmicas são escritas como a soma
da contribuição de um sistema de referência (geralmente consideram-se gases
ou soluções ideais) e do desvio em relação ao sistema de referência
(chamadas de propriedades residuais ou de excesso, dependendo da
referência adotada). As expressões para contribuições ideais são bem
conhecidas e facilmente determinadas. Porém, a contribuição relativa ao desvio
nem sempre é calculada com facilidade. No caso desta tese, utilizou-se uma
equação de estado para modelar as fases do sistema. Assim sendo, a
contribuição de referência é a de gás ideal que, por sua simplicidade, foi
implementada manualmente em linguagem computacional. O desvio é uma
contribuição usualmente denominada de residual, cuja implementação
computacional foi realizada utilizando-se o programa Thermath.
Os programas de “flash” TPN e TVN necessários para a determinação,
respectivamente, das condições das correntes de entrada e da condição inicial
no tambor de “flash” para a simulação dinâmica encontravam-se praticamente
prontos (Gandur, 1990; Espósito, 1999). Entretanto, foi necessário deduzir e
implementar automaticamente as diversas propriedades termodinâmicas
necessárias para a formulação do “flash” UVN que é resolvido em cada passo
de integração. As propriedades assim implementadas foram:
a) a própria equação de estado, ou seja, a expressão para a pressão;
b) a primeira e a segunda derivadas da pressão em relação ao volume molar,
a temperatura e a frações molares constantes, necessárias para a
resolução da equação de estado pelo algoritmo de Topliss et al. (1988)
durante o teste de estabilidade de fases;
c) a energia interna residual e o logaritmo da fugacidade de cada composto na
mistura e as derivadas destas duas propriedades a volume total constante:
em relação à temperatura (que é a capacidade calorífica residual a volume
constante) e em relação aos números de mols de cada composto, mantidos
FORMULAÇÃO
37
constantes os números de mols dos demais compostos, a temperatura e o
volume total;
d) a entalpia residual.
A entrada de dados para a dedução destas expressões e a rotina em
FORTRAN escrita automaticamente encontram-se no Apêndice 2.
3.3.2. Teste da Estabilidade Global de Fases
O teste de estabilidade global de fases de um sistema contendo nc
componentes, seguindo a formulação proposta por Michelsen (1982), requer a
minimização da função objetivo:
nc
Ω = ∑ xiJ ( µ iJ − µ ij )
(3.18)
i =1
com J ≠ j , onde J é a fase denominada de incipiente e j é uma das fases
presentes no sistema. µ iJ e x iJ representam, respectivamente, o potencial
químico e a fração molar do componente i na fase incipiente, e µ ij corresponde
ao potencial químico do componente i nas demais fases presentes no sistema.
Deve-se notar que o teste de estabilidade é realizado para detectar se uma
solução das equações de equilíbrio termodinâmico, determinada supondo-se
um determinado número de fases, é estável ou não. Portanto, o teste de
estabilidade de fases é realizado a partir de uma configuração do sistema em
que os valores de µ ij são iguais em todas as fases.
Se a função Ω for positiva, o sistema é estável; em caso contrário, o
sistema é instável e, portanto, pelo menos mais uma fase deve ser agregada
para a determinação das condições de equilíbrio termodinâmico. Este teste de
estabilidade de fases tem sido amplamente utilizado na literatura de equilíbrio
de fases desde que foi proposto. A implementação utilizada nesta tese foi a
realizada originalmente por Michelsen (1982), sem modificações substanciais.
Por este motivo, não se faz uma discussão a respeito dos detalhes de
implementação.
38
FORMULAÇÃO
3.3.3. Equação de Estado
A equação de estado utilizada na modelagem do sistema foi a equação
de Peng-Robinson (1976), por ser de uso freqüente na modelagem de
processos químicos e petroquímicos. Esta equação de estado está
representada pela Eq.(3.19):
P=
am
RT
− 2
v − bm v − 2bmv − bm2
(3.19)
onde P representa a pressão do sistema, T a temperatura, R a constante
universal dos gases, v o volume molar e am e bm são os parâmetros de
interação da mistura, calculados segundo as regras de mistura de van der
Walls. O cálculo dos parâmetros a e b para a mistura, bem como os valores
para os compostos puros, estão apresentados a seguir:
nc
nc
i =1
j =1
am = ∑ x i ∑ x j aij
(3.20)
nc
bm = ∑ x i bi
(3.21)
aij = ai a j (1 − k ij )
(3.22)
i =1
[
(
ai =
0,45724R 2Tci2
1 - 1,43754 1 - Tri0,5
Pci
bi =
RTci
8Pci
)]
2
(3.23)
(3.24)
onde Tci e Pci são, respectivamente, a temperatura e pressão críticas, xi é a
fração molar do componente i e kij é o parâmetro de interação binária.
As propriedades termodinâmicas necessárias à simulação foram
deduzidas rigorosamente, usando o programa Thermath, tornando possível
realizar simulações em condições afastadas da idealidade.
FORMULAÇÃO
39
3.4. “Flash” TVN
As condições iniciais necessárias para a resolução das equações do
modelo dinâmico são obtidas através das especificações de um “flash” com
temperatura (T), volume (V) e número de mols (N) conhecidos, ou seja, um
“flash” TVN.
Espósito (1999) desenvolveu um programa, em linguagem FORTRAN,
que consiste na determinação de equilíbrio termodinâmico em sistemas sob
influência de campos gravitacionais. O algoritmo é análogo ao procedimento
desenvolvido por Castier (1988) para sistemas com T, P e N especificados, que
consiste na minimização da energia livre de Gibbs. Porém, este caso possui
especificações de T, V e N e a função objetivo a ser minimizada é a energia
livre de Helmholtz.
O trabalho de Espósito (1999) visava ao estudo do efeito do campo
gravitacional em sistemas em equilíbrio. A proposta do seu trabalho é a
utilização de potenciais termodinâmicos modificados que acrescentam a
influência da gravidade no sistema. No presente trabalho de tese, tal influência
não foi incluída como fonte de interferência significativa nos resultados da
dinâmica dos equipamentos simulados. Por este motivo, algumas modificações
foram realizadas no programa de Espósito (1999), de modo a negligenciar o
efeito do campo gravitacional, tornando-o um “flash TVN” convencional. Deste
modo, os resultados gerados pela versão modificada do programa de Espósito
(1999), sem levar em conta o efeito do campo gravitacional, são consistentes
com as hipóteses adotadas na formulação da dinâmica do tambor de “flash”,
garantindo-se assim a consistência das condições iniciais determinadas.
Exceto pelo negligenciamento do campo gravitacional, o algoritmo e o
programa de Espósito (1999) foram utilizados sem modificações adicionais. Por
este motivo, considera-se que uma discussão detalhada do algoritmo de
Espósito (1999) e de sua implementação computacional fuja ao escopo do
texto desta tese.
FORMULAÇÃO
40
3.5. “Flash” TPN
No balanço de energia, torna-se necessário o cálculo da entalpia das
correntes de entrada e saída. Os valores das correntes de saída são
calculados pela própria rotina que resolve a parte algébrica do modelo
dinâmico, ou seja, as equações do “flash” UVN. Entretanto, os valores
referentes às correntes de entrada são calculados através de um “flash” com
pressão, temperatura e número de mols especificados (“flash” TPN).
Enfocando, neste item, as propriedades das correntes de entrada, devese notar que suas entalpias podem variar com o tempo, caso tais correntes
estejam sujeitas a perturbações. Neste caso, a resolução das equações
diferenciais requer a avaliação de “flashes” TPN para avaliar as propriedades
das correntes de entrada em cada passo de integração.
Adotou-se o algoritmo de “flash” TPN proposto por Michelsen (1982a, b),
que consiste na minimização da energia livre de Gibbs. Do ponto de vista de
implementação computacional, foi utilizada aquela realizada por Gandur
(1990). Em sua tese de mestrado, Gandur (1990) concentrou-se no
desenvolvimento de um algoritmo para “flash” reativo com especificação de
pressão, entalpia e número de mols. Entretanto, o programa computacional
implementado por ele permite também resolver “flashes” TPN em sistemas não
reativos, calculando-se a entalpia como uma das variáveis de saída. Portanto,
o programa de Gandur (1990) pode ser utilizado sem necessidade de nenhuma
modificação importante e, por este motivo, seus detalhes não são discutidos
neste texto de tese.
4
Resultados
A capacidade de simulação do sistema desenvolvido é demonstrada a
partir de exemplos, apresentados a seguir. Em todos os casos, a integração
das equações diferenciais foi realizada utilizando o código ode15s do
MATLAB®.
Em todos os gráficos apresentados neste capítulo, as temperaturas são
em Kelvin; as pressões, em bar; as energias internas, em cm3.bar; as cargas
térmicas, em cm3.bar/min; as vazões, em mol/min; os volumes, em cm3 e o
tempo, em min.
4.1. Sistema Exemplo 1: BENZENO-ETILBENZENO
Uma mistura de benzeno (i = 1) e etilbenzeno (i = 2), inicialmente nas
condições apresentadas na Tabela 4.1 (utilizada para todas as rodadas com
este Sistema Exemplo), é submetida a padrões temporais dispostos na Tabela
4.2.
Tabela 4.1: Condições Iniciais do Sistema Exemplo 1
Temperatura Inicial (K)
373
Volume do Tambor (m3)
1,1
Número de Mols Inicial
i = 1 0,33
i = 2 0,67
A Fig. (4.1) reproduz uma composição de telas de especificação e
acompanhamento da simulação. As séries temporais são reproduzidas na Fig.
(4.2). Observa-se que em t=250 min surge uma nova fase (np=2).
Figura 4.1. Sistema Exemplo 1 Submetido a Padrão Temporal Descrito
na Tabela 4.2: Telas de Especificação e Acompanhamento SIMULINK®
Figura 4.2. Sistema Exemplo 1 sem Controle: Séries Temporais
44
RESULTADOS
Tabela 4.2: Condições de Entrada 1 do Sistema Exemplo 1
Carga Térmica (cm3.bar /min)
1000+100sen(t)
Pressão da Alimentação (bar)
30
Temperatura da Alimentação (K)
271+20 sen(0,005t)
Vazão da Alimentação (mol/min)
0,5
Vazão de Saída (mol/min)
0
Composição da Alimentação
i=1
0,33
i=2
0,67
Na Fig. (4.3), o mesmo sistema descrito na Tabela 4.1 é submetido a
novo padrão temporal de acordo com a Tabela 4.3, e simulado em ambiente
MATLAB®. Os pontos dos gráficos são gerados a cada vez que a rotina de
cálculo das derivadas é chamada pelo integrador ode15s, observando-se a
adaptação do passo a cada variação abrupta da entrada temperatura. Verificase surgimento de fase em tempo de aproximadamente 50 minutos.
Tabela 4.3: Condições de Entrada 2 do Sistema Exemplo 1
Carga Térmica (cm3.bar/min)
10000 + 1000 sen(0,01 t)
Pressão da Alimentação (bar)
5
Temperatura da Alimentação (K)
250, sujeito a pulsos de
amplitude 50 e período
20π minutos
Vazão da Alimentação (mol/min)
2
Vazão de Saída (mol/min)
0
Composição da Alimentação
i=1
0,33
i=2
0,67
45
RESULTADOS
4.2. Sistema Exemplo 2: METANO
A simulação reportada utiliza metano puro, inicialmente na condição
apresentada na Tabela 4.4. O tambor, em um primeiro caso, é alimentado de
acordo com o padrão temporal descrito na Tabela 4.5.
Tabela 4.4: Condições Iniciais do Sistema Exemplo 2
Temperatura Inicial (K)
298,15
Volume do Tambor (m3)
30
Número de Mols Inicial
922,73
As séries temporais de entrada (Q) e resposta (U, P, T e np) estão
apresentadas na Fig. (4.4).
Adicionalmente, para o Sistema Exemplo 2, visando a explorar
surgimento de fase, impôs-se ao sistema novo padrão temporal na
alimentação, conforme mostrado na Tabela 4.6, obtendo-se as séries
temporais apresentadas na Fig. (4.5).
Tabela 4.5: Condições de Entrada 1 do Sistema Exemplo 2
Carga Térmica (cm3.bar/min)
Para t < 500
Q=0
Para t ≥ 500
Q = -100 * (t - 500)
Pressão da Alimentação (bar)
3
Temperatura da Alimentação (K)
298,15
Vazão da Alimentação (mol/min)
1
Vazão de Saída (mol/min)
0
46
RESULTADOS
Tabela 4.6: Condições de Entrada 2 do Sistema Exemplo 2
Carga Térmica (cm3.bar/min)
Para t < 500
Q=0
Para t ≥ 500
Q = -10000 (t - 500)
Pressão da Alimentação(bar)
10
Temperatura da Alimentação (K)
301,62
Vazão da Alimentação (mol/min)
10
Vazão de Saída (mol/min)
0
Observa-se a adaptação do passo em t = 500min, quando é imposta
rampa de remoção de calor no tambor. Em t = 800min, surge uma segunda
fase, líquida, que ocasiona mudança no padrão temporal de P e T. Este fato,
em termos de controle de processo, implicaria na necessidade de adaptação
do ganho do controlador à nova situação.
4.3. Controle para Sistema Exemplo 1
Os recursos do SIMULINK® foram explorados na operação do Sistema
Exemplo 1, com as condições iniciais da Tabela 4.1, exceto pela temperatura
que foi inicializada em 400 K. A Fig. (4.6) mostra a tela do modelo em
SIMULINK®, e algumas telas gráficas. As séries temporais da simulação estão
apresentadas na Fig. (4.7). A pressão e a composição da corrente de
alimentação são aquelas apresentadas na Tabela 4.2. A temperatura da
alimentação foi fixada em 380K. A carga térmica do “flash” foi manipulada por
um controlador de temperatura, que teve seu “set-point” alterado de através de
pulsos conforme Fig. (4.7). A pressão de operação do “flash” foi controlada em
0,45 bar, através de manipulação da vazão de entrada. A vazão de retirada foi
fixada em 0 mol/min.
Figura 4.3. Sistema Exemplo 1 Submetido a Padrão Temporal Descrito na Tabela 4.3
(No diagrama da temperatura, os círculos e cruzes indicam,
respectivamente, os valores no tanque e na corrente de alimentação)
Figura 4.4. Simulação de Injeção de Gás Metano Puro no Tambor de “Flash”
Submetido a Padrão Temporal Descrito na Tabela 4.5
Figura 4.5. Simulação de Injeção de Gás Metano Puro no Tambor de “Flash” Submetido a
Padrão Temporal Descrito na Tabela 4.6
50
RESULTADOS
Na operação reportada, a sintonia do controlador de temperatura seguiu
ajuste mais rígido, relaxando-se a sintonia do controlador de pressão para
minimizar interação entre as malhas, conforme apresentado na Tabela 4.9.
A função de transferência para os controladores está apresentada na
Eq. (4.1):
Gc ( s ) = K c (1 +
τ s
1
+ D )
τ I s αs + 1
(4.1)
Tabela 4.7: Sintonia de Controladores
Neste
Parâmetros de
Controlador de
Controlador de
Sintonia
Temperatura
Pressão
Kc
1000
0,01
τI
0,001
10
τD
0,0005
0,05
α
0,01
0,01
exemplo,
não
ocorreu
surgimento
de
nova
fase
e,
conseqüentemente, não houve dificuldades de controle. Para explorar esta
situação, fixou-se o valor de referência em T = 370K, com pulsos de –2K.
Devido às dificuldades de controle, empregou-se esquema de aceleração do
controlador de temperatura apresentado na Tabela 4.8, A Fig. (4.8) apresenta a
comfiguração de simulação no SIMULINK®.
As séries temporais de resposta de temperatura a mudançcas do “SetPoint” do controlador de temperatura, e o movimento da variável manipulada Q
estão apresentadas na Fig. (4.9).
Figura 4.6. Sistema Exemplo 1 sob Controle de T e P
Figura 4.7. Sistema Exemplo 1 sob Controle de T e P: Séries Temporais
(F representa a vazão da alimentação)
53
RESULTADOS
Tabela 4.8: Adaptação do Ganho do Controlador de Temperatura
Número de Fases
1
2
Sintonia de Controlador (*)
K c = K c ,1 ;
K c ,1
Kc
=
Kc
= 10
K c = 10 K c ,1 ;
τI
τI
τ I ,1
;τ D = τ D ,1
K c ,1
τ I ,1
;τ D = 10τ D ,1
(*) Kc,1 = 1000; τI,1 = 100; τD,1 = 0,005, α = 0,01 (mantido constante)
Apesar da adaptação, verifica-se uma dinâmica mais lenta do processo
(o ciclo dos pulsos precisou ser aumentado em relação à Fig. (4.7) para que os
novos “set-points” fossem alcançados, apesar de ser empregada uma sintonia
muito mais rápida nos controladores para nP = 1). O desempenho insatisfatório
da malha de controle sugere o possível risco de projeto de estratégias de
controladores a partir de modelos dinâmicos que não incorporam o rigor
necessário nos modelos termodinâmicos empregados.
4.4. Sistema Exemplo 3: GLP
Um terceiro Sistema Exemplo foi empregado para teste do programa de
simulação de “flash” multifásico. A Fig. (4.10) representa o esquema do
processo proposto.
O sistema é composto por 6 componentes e as condições iniciais estão
dispostas na Tabela 4.9. O padrão de operação temporal está apresentado na
Tabela 4.10.
Figura 4.8. Adaptação de Ganho do Controlador sob Evento de Surgimento de Fase
Figura 4.9. Desempenho da Estratégia de Controle com Adaptação de Ganho na Ocorrência de Nova Fase
(F representa a vazão de alimentação)
Figura 4.10. Sistema Exemplo 3: Flash
Na Tabela 4.10, está apresentada a variação dos ”set-points” dos
controladores no intervalo de simulação e, na Tabela 4.11 apresentam-se as
condições operacionais da corrente de entrada. O valor inicial da carga térmica
(Q) foi especificado igual a –20000 cm3.bar/min.
Tabela 4.9: Condições Iniciais do Sistema Exemplo 3
Temperatura Inicial (K)
315,0
Volume do Tambor (m3)
Número de Mols Inicial
4,4
Etano(i=1)
10,8
Propeno(i=2) 360,8
Propano(i=3) 146,5
i-Butano(i=4) 233,0
n-Butano(i=5) 233,0
n-Pentano(i=6)
15,9
Na Fig. (4.11) é apresentado o esquema de simulação em ambiente
SIMULINK®. Realiza-se o controle na temperatura do processo através da
carga térmica (Q) adicionada ou removida do processo, da pressão através da
vazão de retirada de vapor presente no sistema e há ainda um controlador de
nível que manipula a vazão da fase líquida que somente entra em operação
quando o volume atingir 200000 cm3.
Tabela 4.10: Variação dos “set-points” durante a Simulação
t < 20
Temperatura (K)
300
Pressão (bar)
5
Nível do tanque (cm3)
10000
20 ≤ t < 30
Temperatura (K)
300
Pressão (bar)
3
Nível do tanque (cm3)
100000
30 ≤ t < 40
Temperatura (K)
305
Pressão (bar)
3
Nível do tanque (cm3)
10000
Temperatura (K)
300
Pressão (bar)
3
Nível do tanque (cm3)
100000
t ≥ 40
Figura 4.11. Tela de Especificação para o Sistema Exemplo 3
em Ambiente SIMULINK®
A Fig. (4.12) apresenta a série temporal referente à simulação. Observase que logo ocorre o surgimento de fase, aproximadamente em 1 min. Para
melhor visualização, a Fig. (4.13) apresenta a evolução da temperatura do
sistema e da carga térmica a ele submetida.
Em t=20 min ocorre uma mudança no “Set-Point” da pressão. Com isto,
observa-se a tentativa de ajuste do controlador, porém na série temporal que
precede tal evento, o “off-set” estava estável dificultando a ação integral do
controlador.
Tabela 4.11: Condições de Entrada do Sistema Exemplo 3
Pressão da Alimentação (bar)
5
Temperatura da Alimentação (K)
270
Vazão da Alimentação (mol/min)
1
Composição da Alimentação
i=1
0,0108
i=2
0,3608
i=3
0,1465
i=4
0,2330
i=5
0,2330
i=6
0,1590
Com a alteração no “set-point” do volume da fase líquido, as equações
algébricas não convergiram para a variável pressão, conforme pode ser
observado na Fig. (4.14). Este fato é reconhecido pela rotins através uma
função evento que pára a integração, retorna ao passo anterior, prosseguindo o
avanço no tempo com passos 100 vezes menor. Isto faz com que o sistema se
recupere e retome a convergência.
Figura 4.12. Sistema Exemplo 1: Série Temporais
(Vap e Liq são, respectivamente, as vazões retirada de vapor
e líquido do sistema e Vol é o nível do tanque)
Figura 4.13. Sistema Exemplo 3: Temperatura do Sistema e Carga Térmica
P
Vap
Figura 4.14. Sistema Exemplo 3: Resposta da Pressão à
modificações no “Set-Point”
Vol
Liq
Figura 4.15 - Resposta do
Controlador de Nível ao Surgimento da Fase
5
Conclusões e Sugestões
A simulação dinâmica de “flashes” bifásicos foi abordada com modelos
rigorosos para cálculo de propriedades termodinâmicas e de equilíbrio de
fases. Utilizou-se o ambiente de computação algébrica Mathematica® para
geração automática de códigos em FORTRAN para cálculo de propriedades
termodinâmicas e o ambiente de resolução numérica e visualização gráfica
MATLAB®/SIMULINK® para a simulação do equipamento.
A implementação envolveu o desenvolvimento das equações para “flash”
com especificação de energia interna, volume e número de mols (“flash” UVN),
resultando em expressões complexas, e o código foi desenvolvido por
computação algébrica. O emprego da ferramenta de desenvolvimento e a
integração com o ambiente MATLAB® permitiu a simulação de processo
envolvendo misturas de hidrocarbonetos, e gás liquefeito de petróleo, de
interesse do Setor de Petróleo e Gás.
Uma conclusão importante, do ponto de vista numérico, diz respeito às
variáveis independentes adotadas para a resolução das equações algébricas
que constituem o modelo do “flash” UVN. Além dos números de mols dos
compostos nas fases, Michelsen (1999) sugeriu que se utilizassem os
logaritmos da temperatura e dos volumes das fases. Observou-se que,
freqüentemente, não se obtém convergência com esta escolha, na vizinhança
dos pontos de saturação, embora, longe destes, o esquema proposto por
Michelsen (1999) funcione bem. Obteve-se uma formulação numericamente
bem mais confiável ao se utilizar a temperatura e os volumes das fases
diretamente, sem uso de logaritmos.
Identificou-se, entretanto, uma limitação na formulação adotada. Apesar
do modelo detectar corretamente o aparecimento de uma terceira fase para
certos conjuntos de especificação, não foi possível obter a convergência das
equações algébricas do modelo nos problemas trifásicos testados. Extensos
CONCLUSÕES E SUGESTÕES
66
testes foram realizados em busca de algum possível erro de programação, que
não foi encontrado. Especula-se que a possível causa para a limitação de
desempenho encontrada esteja na formulação, que caracteriza a solução
através de um ponto de sela. Michelsen (1987), por exemplo, detectou
dificuldades na resolução de “flashes” isentálpicos quando utilizou uma
formulação baseada em um ponto de sela.
Vale também ressaltar que a formulação adotada permite desacoplar,
em certos casos, os subconjuntos das equações algébricas e diferenciais. Em
função do interesse em desenvolver um programa geral de simulação, este não
foi um aspecto explorado nesta tese. Mas, dependendo do tipo de
especificação, pode ser possível obter soluções analíticas das equações
diferenciais, mesmo quando o subconjunto algébrico apresenta grande
complexidade.
Como sugestões para trabalhos futuros, destacam-se:
1. investigar a possibilidade de formular uma função objetivo algébrica (ou
seja, a função cujas derivadas constituem as equações algébricas do
modelo de “flash”) que incorpore um termo de penalidade. Tal
penalidade deve depender do quadrado do desvio entre os valores
especificados e calculado da energia interna, de modo análogo ao que
foi empregado, com êxito, por Michelsen (1987) e por Gandur (1990), no
caso mais geral, com a existência de reações químicas. A expectativa é
que essa formulação permitiria caracterizar a condição de equilíbrio
termodinâmico como um ponto de mínimo da nova função objetivo, o
que apresenta vantagens numéricas em relação à caracterização como
ponto de sela.
2. investigar a resolução das equações algébrico-diferenciais (EAD) pelo
código odes15s do MATLAB® ou, através de uma gateway, pelo código
DASSL. Os exemplos testados, em primeira etapa neste trabalho,
obtiverem do código ode15s a indicação de não se tratar de um sistema
de EAD´s de índice 1. Entretanto, tal indicação foi interpretada como
resultado preliminar, sujeito a verificação adicional. Métodos de redução
CONCLUSÕES E SUGESTÕES
67
de índice deverão ser investigados para que os benefícios numéricos da
resolução simultânea do conjunto de equações algébrico-diferenciais
sejam utilizados para a simulação de operações unitárias mais
complexas como colunas de destilação reativas.
3. existe um interesse crescente no desenvolvimento de ferramentas de
detecção e diagnóstico de falhas baseadas em modelos do processo,
estimação de parâmetros e estimação de estados (Iserman, 1997). A
detecção de falha é o reconhecimento de um comportamento onde pelo
menos uma propriedade característica ou parâmetro do processo se
desvia da situação de normalidade. Nesta linha, a disponibilidade de
modelos rigorosos não só auxiliam na identificação destes eventos como
também no diagnóstico. O impacto de simplificações normalmente
adotadas frente ao rigor termodinâmico apresentado neste trabalho
requer investigações futuras.
4. incorporar, no programa desenvolvido, chaveamentos lógicos que
interrompam a execução no caso de situações não-físicas, como por
exemplo, o escoamento de regiões de pressões mais baixas para
aquelas de pressões mais altas, tal como ocorreu nos exemplos
apresentados nas Seções 4.1 e 4.3.
6
Referências
Adams, S. e Taylor, R., Thermodynamics with Maple V. III - Thermodynamic
Property Relations and the Maxwell Equations, Maple Tech. 1, 68-81,
1994.
Alfradique, M.F., Espósito, R.O. e Castier, M., Automatic Generation of
Procedures for the Simulation of Multistage Separators Using Computer
Algebra,
aceito
para
publicação
em
Chemical
Engineering
Communications, 2001.
Bischof, C.H., Carle, A., Corliss, G.F., Griewank, A. e Hovland, P., ADIFOR:
Generating
Derivative
Codes
from
Fortran
Programs,
Scientific
Programming 1, 1-29, 1992.
Brown, P.N., Hindmarsh, A.C. e Petzold, L.R., Using Krylov Methods in the
Solution of Large-Scale Differential-Algebraic Systems, SIAM J. Sci.
Comput., 15, 1467-1488, 1994.
Callen, H.B., Thermodynamics and an Introduction to Thermostatistics,
2aedição, Wiley, New York, 1985.
Castier, M., Automatic Implementation Of Thermodynamic Models Using
Computer Algebra, Computers and Chemical Engineering, 23, 1229,
1999.
Costa Jr, E.F., Caracterização Estrutural Automática de Sistemas AlgébricoDiferenciais, Exame de Qualificação ao Doutorado, PEQ/COPPE,
Universidade Federal do Rio de Janeiro, 2001.
Espósito, R.O., “Cálculo de Equilíbrio Termodinâmico em Sistemas sob a
Influência de Campos Gravitacionais”, Tese de Mestrado, Escola de
Química, Universidade Federal do Rio de Janeiro, 1999.
Espósito R.O, Castier M. e Tavares, F.W., Calculations of thermodynamic
equilibrium
in
systems
subject
to
gravitational
Engineering Science, 55, 3495-3504, 2000.
fields,
Chemical
70
REFERÊNCIAS
Gandur, M.C., “Equilíbrio Químico e de Fases com Especificação de Pressão e
Entalpia em Sistemas não Ideiais”, PEQ/COPPE, Universidade Federal
do Rio de Janeiro, 1990.
Gopal, V. e Biegler, L.T., Nonsmooth Dynamic Simulation with Linear
Programming Based Methods, Computers and Chemical Engineering,
21(7), 675-689, 1997.
Gopal, V. e Biegler, L.T., Smoothing Methods for the Treatment of
Complementarity Conditions and Nested Discontinuities, AIChE Journal,
45(7), 1535-1547, 1999.
Gopal, V. e. Biegler, L.T, An Optimization Approach to Consistent Initialization
and
Reinitialization
after
Discontinuities
of
Differential
Algebraic
Equations, SIAM J. Cient. Comp., 20(2), 447-467, 1999b.
Isermann, R. e Ballé, P., Trends in the Application of Model Based Fault
Detection and Diagnosis of Technical Process, Control Eng. Pratice, 5,
709-719, 1997.
Marquardt, W., Dynamic Process Simulation Recent Progress and Future
Challenges, in “Chemical Process Control” Eds. Y. Arkun e W.H. Ray,
CACHE – AIChE Publications, 130, 1991.
Michelsen, M.L., The Isothermal Flash Problem. I. Stability Analysis, Fluid
Phase Equilibria 9, 1-19, 1982a.
Michelsen, M.L., The Isothermal Flash Problem. II. Phase Split Calculation,
Fluid Phase Equilibria 9, 21–40, 1982b.
Michelsen, M.L., Multiphase Isenthalpic and Isentropic Flash Algorithms, Fluid
Phase Equilibria, 33, 13-27, 1987.
Michelsen, M.L., State Function Based Flash Specifications, Fluid Phase
Equilibria, 158-160, 617-626, 1999.
Müller, D. e Marquardt, W., Dynamic Multiple-Phase Flash Simulation: Global
Stability Analysis Versus Quick Phase Determination, Computers and
Chemical Engineering, 97, S817-S822, 1997.
Peng, D.Y. e Robinson, D.B., A simple two-constant equation of state, Ind. Eng.
Chem. Fundamentals 15, 59-64, 1976.
REFERÊNCIAS
71
Petzold, L., “Numerical Methods for Differential-Algebraic Equations. Current
Status and Future Directions” , Computational Ordinary Differential
Equations, Eds. J. Cash e I. Gladwell, Clarendon Press, Oxford, 259273, 1992.
Reid, R.C., Prausnitz, J.M. e Poling, B.E., The Properties of Gases and Liquids,
McGraw-Hill Book Company, 4a edição, 1987.
Saha, S. e Carroll, J.J., The Isoenergetic-Isochoric Flash, Fluid Phase
Equilibria, 138, 23-41, 1997.
Shampine, L.F. e Reichelt, M.W., The MATLAB ODE Suite, SIAM Journal on
Scientific Computing, 18, 1-22, 1997.
Shampine, L.F., Reichelt, M.W., e Kierzenka, J.A., Solving Index-1 DAEs in
MATLAB and Simulink, SIAM Review, 41, 538-552, 1999.
Silva, A.S., Diferenciação e Implementação de Modelos Termodinâmicos
Usando Computação Algébrica, Tese de Mestrado, Escola de Química,
Universidade Federal do Rio de Janeiro, 1993.
Silva, A.S. e Castier, M., Automatic Differentiation and Implementation of
Thermodynamic Models Using a Computer Algebra System, Computers
and Chemical Engineering 17, S473-S478, 1993a.
Silva, A.S. e Castier, M., Computação Algébrica: uma Revisão com Aplicações
à Engenharia Química, Revista Brasileira de Engenharia, Caderno de
Engenharia Química 10, 65-94, 1993b.
Taylor, R. e. Monagan, M.B, Thermodynamics with Maple V. I - Equations of
state, Maple Tech., 10, 50-59, 1993.
Taylor, R., Thermodynamics with Maple V. II - Phase equilibria in binary
systems, Maple Tech. 1, 83-92, 1994.
Taylor, R., Automatic Derivation of Thermodynamic Property Functions using
Computer Algebra, Fluid Phase Equilibria 129, 37-47, 1997.
Topliss, R.J., Dimitrelis, D. e Prausnitz, J.M., Computational Aspects of a NonCubic Equation of State for Phase Equilibrium Calculations. Effect of
Density Dependent Mixing Rules. Comput. Chem. Eng., 12, 483-489,
1988.
REFERÊNCIAS
72
Vieira, R.C., Iniciação de Sistemas Algébrico-Diferenciais Empregando
Algoritmos Heurísticos de Otimização, D.Sc. Exame de Qualificação ao
Doutorado, Universidade Federal do Rio de Janeiro, 1999.
Wilson, G., AIChE National Meeting, May 4-7, 1969, Cleveland, OH.
Wolfram, S., Mathematica®: a System for Doing Mathematics by Computer,
Addison-Wesley, 1991, Redwood City, CA.
Apêndice 1
Montagem das Equações Algébricas
de Equilíbrio de Fases
O objetivo deste Apêndice é apresentar a dedução das derivadas da
Eq. (A1.1), que compõem a matriz Jacobiana o cálculo do equilíbrio de fases
do sistema. Essas expressões foram deduzidas com o auxílio do programa
Thermath. Vale ressaltar, contudo, de que a dedução foi menos automatizada
do que a das propriedades termodinâmicas a partir da equação de estado. Isto
ocorrreu porque se considerou que o volume de uma das fases era uma
variável dependente, bem como que, para cada composto, havia uma fase na
qual o número de mols era dependente. De qualquer modo, uma vez deduzidas
as expressões, o código FORTRAN correspondente foi gerado pelo programa
Thermath e está apresentado na seção A2.3 do Apêndice 2.
A função QF proposta por Michelsen (1999) e que constitui o ponto de
partida da dedução,.é dada por:
A − U ESPEC
QF =
RT
(A1.1)
e será derivada em relação à temperatura (T), ao volume das fases (V) e ao
número de mols (N) dos componentes nas fases independentes produzindo as
funções rT, rV e rN respectivamente as Eqs. (A1.5), (A1.9) e (A1.14). As fases
dependentes estão representadas, por conveniência, por letras maiúsculas.
Para a derivada em relação à temperatura, tem-se:
1
1 ⎛ ∂A ⎞
⎛ ∂Q ⎞
rT = ⎜ F ⎟
A − U ESPEC +
=−
⎜
⎟
2
RT ⎝ ∂T ⎠V ,N
RT
⎝ ∂T ⎠V ,N
(
mas:
)
(A1.2)
73
APÊNDICE 1
⎛ ∂A ⎞
= −S
⎜
⎟
⎝ ∂T ⎠V ,N
(A1.3)
e
S=
U−A
T
(A1.4)
Substituindo as Eqs. (A1.3) e (A1.4) na Eq. (A1.2):
rT
(U − U
=−
ESPEC
RT
)
(A1.5)
2
Para a derivada em relação aos volumes de fase independentes, tem-se:
⎛ ∂Q
rv = ⎜⎜ F
⎝ ∂Vf
⎞
1 ⎛ ∂A
⎟⎟
⎜⎜
=
⎠ TN,Vg ,g ≠ f ,F , RT ⎝ ∂Vf
⎡
⎞
1 ⎢⎛ ∂Af
⎟⎟
⎜⎜
=
⎠TN,Vg ,g ≠ f ,F , RT ⎢⎝ ∂Vf
⎣
⎞
⎛ ∂A
⎟⎟
− ⎜⎜ F
⎠ TN,Vg ,g ≠ f ,F , ⎝ ∂VF
⎞
⎛ ∂VF
⎟⎟
⎜⎜
⎠TN,Vg ,g ≠ f ,F , ⎝ ∂Vf
(A1.6)
mas
⎛ ∂A ⎞
⎟ = −P
⎜
⎝ ∂V ⎠
(A1.7)
e a derivada do volume da fase dependente (F) em relação à qualquer outra
fase independente, será sempre igual a –1 pois:
VF = VTOTAL −
np −1
∑V
f =1
f
(A1.8)
com isto,
rV =
1
(PF − Pf )
RT
(A1.9)
Para a derivada em relação aos números de mols independentes, temse:
⎤
⎞⎥
⎟⎟
⎠⎥
⎦
74
APÊNDICE 1
⎛ ∂Q
rN = ⎜ F
⎜ ∂N
ij
⎝
⎞
1 ⎛⎜ ∂A
⎟
=
⎜
⎟ T ,V ,
⎠ Nmn ,m ≠ i ,n ≠ j ,J RT ⎝ ∂N ij
⎞
⎟
⎟ T ,V ,
⎠ N mn ,m ≠ i ,n ≠ j ,J
(A1.10)
mas a derivada da Eq. (A1.10) pode ser dividida:
⎛ ∂A
⎜
⎜ ∂N
⎝ ij
⎛ ∂A j
⎞
⎟
=⎜
⎜
⎟ T ,V ,
⎠ N mn ,m ≠ i ,n ≠ j ,J ⎝ ∂N ij
⎞
⎛ ∂A
⎟
+ ⎜⎜ J
⎟ T ,V ,
⎠ Nmn ,m ≠ i ,n ≠ j ,J ⎝ ∂N iJ
⎛ ∂N iJ
⎞
⎜
⎟⎟
⎜
,
,
T
V
⎠ Nmn ,m ≠ i ,n ≠ j ,J ⎝ ∂N ij
⎞
⎟
⎟
⎠
(A1.11)
seguindo o mesmo raciocínio do apresentado para o volume, é possível
escrever o número de mols do componente i segundo a Eq.(A1.12). A derivada
do número de mols na fase dependente (J) em relação às fases independentes
(j) é portanto, -1.
np −1
N iJ = N i − ∑ N ij
(A1.12)
⎛ ∂A ⎞
mas ⎜
⎟ =µ
⎝ ∂N ⎠ T ,V
(A1.13)
j =1
Com isto, a Eq. (A1.10) pode ser escrita conforme a Eq. (A1.14):
rN =
µ ij − µ iJ
(A1.14)
RT
Derivando as funções rT, rV e rN em relação à temperatura obtem-se as
Eqs.(A1.15), (A1.16) e (A1.17).
(
)
2 U − U ESPEC
1 ⎛ ∂U ⎞
⎛ ∂rT ⎞
−
⎜
⎟ =
⎜
⎟
3
RT
RT 2 ⎝ ∂T ⎠V ,N
⎝ ∂T ⎠V ,N
1
⎛ ∂rv ⎞
=−
(PF − Pf ) + 1
⎜
⎟
2
RT
RT
⎝ ∂T ⎠TN,abV ,,b ≠ s,S
1 ⎡ (µ iJ − µ ij ) ⎛ ∂µ ij
⎛ ∂rN ⎞
+ ⎜⎜
⎢
⎜
⎟ =
T
⎝ ∂T ⎠V ,N RT ⎢⎣
⎝ ∂T
⎤
⎡⎛ ∂P ⎞
⎛ ∂P ⎞
⎥
⎢⎜ F ⎟
−⎜ f ⎟
⎢⎝ ∂T ⎠TN,V ,,b ≠ s,S ⎝ ∂T ⎠TN,V ,,b ≠ s,S ⎥
ab
ab
⎦
⎣
⎞
⎛ ∂µ ⎞ ⎤
⎟⎟ − ⎜ iJ ⎟ ⎥
⎠V ,N ⎝ ∂T ⎠V ,N ⎥⎦
(A1.15)
(A1.16)
(A1.17)
Derivando agora em relação ao volume Vp das fases independentes,
ressaltando que Vp não necessariamente é igual a Vf.
75
APÊNDICE 1
⎛ ∂rT
⎜
⎜ ∂V
⎝ p
⎞
1 ⎛⎜ ∂U
⎟
=−
⎟ T ,V ,m ≠ p,P ,
RT 2 ⎜⎝ ∂V p
⎠N m
⎡
⎞
1
⎢⎛⎜ ∂U P
⎟
=−
2
⎟ T ,V ,m ≠ p,P ,
RT ⎢⎜⎝ ∂VP
⎠N m
⎢⎣
⎞
⎟⎟
⎠ TN,Vm ,m ≠ p,P ,
⎛ ∂V
⋅⎜ P
⎜ ∂V
⎝ p
⎞ ⎛ ∂U p
⎟+⎜
⎟ ⎜ ∂V
⎠ ⎝ p
(A1.18)
mas conforme discutido anteriormente,
⎛ ∂VP
⎜
⎜ ∂V
⎝ p
⎞
⎟
= −1
⎟T ,V ,m ≠ p,P ,
m
⎠N
(A1.19)
portanto,
⎛ ∂rT ⎞
1
⎜
⎟
⎜ ∂V ⎟T ,V ,m ≠ p,P , = − RT 2
⎝ p ⎠N m
⎡⎛
⎛ ∂UP
⎢⎜ ∂U p ⎞⎟
⎜
−
⎢⎜⎝ ∂Vp ⎟⎠T ,Vm ,m ≠ p,P , ⎜⎝ ∂Vp
⎢⎣
N
⎤
⎞
⎥
⎟
⎟T ,V ,m ≠ p,P , ⎥
⎠N m
⎥⎦
(A1.20)
Derivando a função rV em relação ao volume Vp:
⎡
⎛ ∂rv ⎞
⎛ ∂Pf
1 ⎢⎛⎜ ∂PF ⎞⎟
⎜
⎟
⎜
⎜ ∂V ⎟T ,V ,m ≠ p,P , = RT ⎢⎜ ∂V ⎟T ,V ,m ≠ p,P , − ⎜ ∂V
p ⎠
m
⎝ P
⎝ p ⎠N m
⎝
N
⎣⎢
⎤
⎞
⎥
⎟⎟
⎥
,
,
≠
,
,
T
V
m
p
P
m
⎠N
⎦⎥
(A1.21)
mas a Eq. (A1.21) pode ser escrita da forma:
⎡
⎛ ∂rV ⎞
1 ⎢⎛ ∂PF
⎟
⎜
⎜
=
⎜ ∂V ⎟T ,V ,m ≠ p,P , RT ⎢⎜ ∂V
⎝ p ⎠N m
⎣⎝ F
⎛ ∂V ⎞ ⎛ ∂P
⎞
⎟⎟
⋅ ⎜ F ⎟ − ⎜⎜ f
⎟
⎜
⎠TN,Vm ,m ≠ p,P , ⎝ ∂Vp ⎠ ⎝ ∂Vf
⎛ ∂Vf ⎞⎤
⎞
⎟⎥
⎟⎟
⋅⎜
⎜ ∂V ⎟⎥
T
V
m
p
P
,
,
≠
,
,
m
⎠N
⎝ p ⎠⎦
(A1.22)
As derivadas apresentadas na Eq. (A1.22) devem ser analisadas
⎛ ∂V
separadamente. O termo ⎜ F
⎜ ∂V
⎝ p
analisado com maior cuidado:
f =p
⇒
⎛ ∂Vf ⎞
⎜
⎟
=1
⎜ ∂V ⎟T ,V ,m ≠ p,P ,
⎝ p ⎠N m
f ≠p
⇒
⎛ ∂Vf ⎞
⎜
⎟
=0
⎜ ∂V ⎟T ,V ,m ≠ p,P ,
⎝ p ⎠N m
⎞
⎞
⎛
⎟ = −1, porém o termo ⎜ ∂Vf ⎟ deve ser
⎟
⎜
⎟
⎝ ∂Vp ⎠
⎠
⎤
⎞
⎥
⎟
⎟ T ,V ,m ≠ p,P , ⎥
⎠N m
⎦⎥
76
APÊNDICE 1
Com isto, é possível reescrever a Eq. (A1.22):
⎛ ∂rv ⎞
1
⎜
⎟
=−
⎜ ∂V ⎟T ,V ,m ≠ p,P ,
RT
⎝ p ⎠N m
⎡⎛
⎢⎜ ∂PF
⎢⎜⎝ ∂VF
⎣
⎛ ∂P
⎞
⎟⎟
− δ f , p ⋅ ⎜⎜ f
⎝ ∂Vf
⎠TN,Vm ,m ≠ p,P ,
⎤
⎞
⎥
⎟⎟
⎥
T
,
V
,
m
≠
p
,
P
,
⎠N m
⎦
(A1.23)
Derivando rN em relação ao volume da fase p:
⎡
⎛ ∂µ
⎛ ∂rN ⎞
1 ⎢⎛⎜ ∂µ ij ⎞⎟
⎟
⎜
− ⎜ iJ
=
⎜ ∂V ⎟T ,V ,m ≠ p,P , RT ⎢⎜ ∂V ⎟T ,V ,m ≠ p,P , ⎜ ∂V
⎝ p
⎝ p ⎠N m
⎢⎣⎝ p ⎠N m
⎤
⎞
⎥
⎟
⎟T ,V ,m ≠ p,P , ⎥
⎠N m
⎥⎦
(A1.22)
A Eq. (A1.22) pode ser escrita na forma:
⎡
⎛ ∂rN ⎞
1 ⎢⎛⎜ ∂µ ij
⎜
⎟
=
⎜ ∂V ⎟T ,V ,m ≠ p,P , RT ⎢⎜ ∂V
⎝ p ⎠N m
⎢⎣⎝ j
⎞
⎛ ∂Vi ⎞ ⎛ ∂µ iJ
⎟
⎜
⎟−⎜
⎟T ,V ,m ≠ p,P , ⎜ ∂V ⎟ ⎜ ∂V
⎠N m
⎝ p⎠ ⎝ J
⎛ ∂V j
Analisando o termo ⎜
⎜ ∂V
⎝ p
⎛ ∂VJ ⎞⎤⎥
⎞
⎜
⎟
⎟⎟
⎜ ∂V ⎟⎥
,
,
≠
,
,
T
V
m
p
P
m
⎠N
⎝ p ⎠⎥⎦
(A1.23)
⎞
⎟ quando j for a fase dependente, ou seja, j =
⎟
⎠
P:
⎛ ∂V j
⎜
⎜ ∂V
⎝ p
⎞
⎟ = −1
⎟
⎠
(A1.24)
caso contrário, deve-se considerar as possibilidades:
⎛ ∂V j
para j = p: ⎜
⎜ ∂V
⎝ p
⎞
⎟ =1
⎟
⎠
(A1.25)
⎛ ∂V j
para j ≠ p: ⎜
⎜ ∂V
⎝ p
⎞
⎟=0
⎟
⎠
(A1.26)
Com isso, pode-se escrever que:
⎛ ∂V j
⎜
⎜ ∂V
⎝ p
⎞
⎟ = δ jp − δ jP
⎟
⎠
(A1.27)
⎛ ∂V ⎞
Analisando o termo ⎜ J ⎟ seguindo o mesmo raciocínio, obtem-se que:
⎜ ∂V ⎟
⎝ p⎠
⎛ ∂VJ ⎞
⎟ = δ Jp − δ JP
⎜
⎜ ∂V ⎟
⎝ p⎠
(A1.28)
77
APÊNDICE 1
Substituindo as Eqs. (A1.27) e (A1.28) na Eq. (A1.23), chega-se à Eq.
(A1.29):
⎡
⎛ ∂rN ⎞
1 ⎢⎛⎜ ∂µ ij
⎟
⎜
=
⎜ ∂V ⎟T ,V ,m ≠ p,P , RT ⎢⎜ ∂V
j
⎝ p ⎠N m
⎣⎢⎝
⎞
⎛ ∂µ
⎟
(
δ
− δ jP ) − ⎜⎜ iJ
jp
⎟T ,V ,m ≠ p,P ,
⎝ ∂VJ
⎠N m
⎤
⎞
⎟⎟
(δ Jp − δ JP )⎥⎥ (A1.29)
T
V
m
p
P
,
,
≠
,
,
m
⎠N
⎦⎥
Tomando as derivadas agora em relação ao número de mols em relação
ao componente r e a fase s (Nrs):
⎛ ∂rT
⎜⎜
⎝ ∂N rs
⎞
1 ⎛ ∂U
⎜
⎟⎟
=−
2 ⎜
,
,
T
V
RT
⎝ ∂N rs
⎠ Nab ,b ≠ s,S
⎡
⎞
1 ⎢⎛ ∂U S
⎜
⎟⎟
=−
2 ⎜
⎢
,
,
T
V
∂N
RT
⎠ Nab ,b ≠ s,S
⎣⎝ rS
⎞
⎟⎟
⎠ TN,abV ,,b ≠ s,S
⎛ ∂N
⋅ ⎜⎜ rS
⎝ ∂N rs
⎞ ⎛ ∂U s
⎟⎟ + ⎜⎜
⎠ ⎝ ∂N rs
⎤
⎞
⎥
⎟⎟
T
,
V
,
⎠ Nab ,b ≠ s,S ⎥⎦
(A1.30)
A derivada do número de mols do componente r na fase dependente S
em relação à fase independente s será sempre igual a –1. Com isto, obtem-se
a
Eq. (A1.31):
⎛ ∂rT ⎞
1
⎜⎜
⎟⎟
=−
RT 2
⎝ ∂Nrs ⎠TN,abV ,,b ≠ s,S
⎛ ∂U ⎞
1
⎜⎜
⎟⎟
=−
RT 2
⎝ ∂Nrs ⎠TN,abV ,,b ≠ s,S
⎤
⎡⎛
⎛ ∂US ⎞
⎥ (A1.31)
⎢⎜ ∂Us ⎞⎟
⎜
⎟
−
⎜ ∂N ⎟T ,V ,
⎢⎜⎝ ∂Nrs ⎟⎠T ,V ,
⎝ rS ⎠Nab ,b ≠ s,S ⎥⎦
N ab ,b ≠ s ,S
⎣
A derivada da função rV em relação ao número de mols do componente r
na fase s é:
⎛ ∂rv
⎜⎜
⎝ ∂N rs
⎡
⎞
1 ⎢⎛ ∂PF
⎜⎜
⎟⎟
=
⎢
T
V
,
,
RT
∂N
⎠ Nab ,b ≠ s,S
⎣⎝ rs
⎛ ∂P
O termo ⎜⎜ F
⎝ ∂N rs
⎛ ∂P
⎞
⎟⎟
− ⎜⎜ f
⎠ TN,abV ,,b ≠ s,S ⎝ ∂N rs
⎤
⎞
⎥
⎟⎟
T
,
V
,
⎠ Nab ,b ≠ s,S ⎥⎦
⎞
⎟⎟
pode ser escrito conforme a Eq. (A1.33):
T
,
V
,
⎠ Nab ,b ≠ s,S
(A1.32)
78
APÊNDICE 1
⎛ ∂Pf ⎞
⎛ ∂P
⎟⎟
⎜⎜
= ⎜⎜ f
⎝ ∂Nrs ⎠Tnab,V,,b ≠ s,S ⎝ ∂Nrf
⎛ ∂N ⎞
⎞
⎟⎟
⋅ ⎜⎜ rf ⎟⎟
⎠Tnab,V,,b ≠ s,S ⎝ ∂Nrs ⎠Tnab,V,,b ≠ s,S
(A1.33)
A derivada do número de mols do componente r na fase f em relação à
⎛ ∂N ⎞
deve ser analisada cuidadosamente:
fase s ⎜⎜ rf ⎟⎟
⎝ ∂Nrs ⎠Tnab,V,,b ≠ s,S
se f = S :
⎛ ∂N rS
⎜⎜
⎝ ∂N rs
⎞
⎟⎟
= −1
⎠Tnab,V,,b ≠ s,S
(A1.34)
caso contrário, ou seja, se f ≠ S :
se f = s
⎛ ∂N rf
⎜⎜
⎝ ∂N rs
⎞
⎟⎟
=1
T
,
V
,
⎠ nab ,b ≠ s,S
(A1.35)
⎛ ∂N rf
⎜⎜
⎝ ∂N rs
⎞
⎟⎟
=0
⎠ Tnab,V,,b ≠ s,S
(A1.36)
logo,
⎛ ∂N rf
⎜⎜
⎝ ∂N rs
⎞
⎟⎟
= (δ fs − δ fS )
⎠ Tnab,V,,b ≠ s,S
(A1.37)
da mesma forma:
⎛ ∂NrF
⎜⎜
⎝ ∂Nrs
⎞
⎟⎟
= (δ Fs − δ FS )
⎠Tnab,V,,b ≠ s,S
(A1.38)
com isto, a Eq. (A1.32) pode ser escrita conforme a Eq. (A1.39):
⎛ ∂rV
⎜⎜
⎝ ∂N rs
⎡
⎞
1 ⎢⎛ ∂Pf
⎟⎟
⎜⎜
=
⎢
,
,
T
V
∂N
RT
⎠ Nab ,b ≠s,S
⎣⎝ rf
⎞
⎛ ∂P
⎟⎟
⋅ (δ fs − δ fS ) − ⎜⎜ F
⎠TN,abV ,,b ≠s,S
⎝ ∂N rF
⎤
⎞
⎟⎟
⋅ (δ Fs − δ FS )⎥
⎥
⎠TN,abV ,,b ≠s,S
⎦
(A1.39)
79
APÊNDICE 1
Derivando agora o função rN em relação ao número de mols Nrs:
⎛ ∂rN
⎜⎜
⎝ ∂N rs
⎡
⎞
1 ⎢⎛ ∂µ ij
⎜⎜
⎟⎟
=
⎠TN,abV ,,b ≠s,S RT ⎢⎣⎝ ∂N rs
⎞
⎛ ∂µ
⎟⎟
− ⎜⎜ iJ
⎠TN,abV ,,b ≠s,S ⎝ ∂N rs
⎛ ∂µ ij
Analisando o termo ⎜⎜
⎝ ∂N rs
⎤
⎞
⎥
⎟⎟
T
,
V
,
⎠ Nab ,b ≠s,S ⎥⎦
(A1.40)
⎞
⎟⎟
:
⎠TN,abV ,,b ≠s,S
caso j = S, pode-se escrever que:
⎛ ∂µ ij
⎜⎜
⎝ ∂N rs
nc ⎛
⎞
∂µij
⎟⎟
= ∑⎜
⎜
⎠TN,abV ,,b ≠s,S z =1 ⎝ ∂N zj
⎞ ⎛ ∂N zj
⎟⋅⎜
⎟ ⎜ ∂N
⎠ ⎝ rs
⎞
⎟⎟
⎠
(A1.41)
Caso o termo relativo à derivada do número de mols do componente z
na fase j em relação ao componente r na fase s só tem sentido se z=r e será
igual a –1, pois trata-se de fases independentes. Deve-se atentar que neste
caso tem-se um somatório, logo somente este termo não nulo contribuirá para
a definição da Eq. (A1.42).
⎛ ∂µ ij
⎜⎜
⎝ ∂N rs
⎛ ∂µij
⎞
⎟⎟
= −⎜
⎜ ∂N
⎠TN,abV ,,b ≠ s,S
⎝ zj
⎞
⎟
⎟T ,V ,
⎠ Nab ,b ≠s,S
(A1.42)
caso j ≠ S, somente quando se tratar da mesma fase independente, ou seja
quando j = s, e a derivada será igual a 1.
Com isto, pode-se escrever o termo em análise segundo a Eq. (A1.43):
⎛ ∂µ ij
⎜⎜
⎝ ∂N rs
⎛ ∂µ ij
⎞
⎟⎟
=⎜
⎜
⎠TN,abV ,,b ≠ s,S ⎝ ∂N rj
⎞
⎟
⋅ (δ is − δ jS )
⎟T ,V ,
⎠ Nab ,b ≠s,S
(A1.43)
O mesmo raciocínio deve ser aplicado para o outro termo derivativo da
Eq. (A1.40) obtendo-se então a Eq. (A1.44):
⎛ ∂r
rNN = ⎜⎜ N
⎝ ∂N rs
⎞
1 ⎡⎛⎜ ∂µ ij
⎟⎟
=
⎢⎜
⎠Tnab,V,,a ≠ r , RT ⎢⎣⎝ ∂N rj
b ≠ s ,S
⎞
⎛
⎟(δ js − δ jS ) + ⎜ ∂µ iJ
⎜ ∂N
⎟
⎝ rJ
⎠
⎤
⎞
⎟⎟(δ Js − δ JS )⎥
⎥⎦
⎠
(A1.44)
Apêndice 2
Cálculo das Propriedades
Termodinâmicas
A seguir, está apresentada a forma de definição no Mathematica® da
função objetivo e das propriedades termodinâmicas necessárias para os
cálculos da resolução do sistema algébrico apresentado no Apêndice 1.
Este apêndice está dividido em quatro seções. As seções A2.1 e A2.2
apresentam, respectivamente, a entrada e a saída (subrotina em FORTRAN)
do programa Thermath das equações do modelo de “flash” UVN. As seções
A2.3 e A2.4 apresentam a entrada e a saída (subrotina em FORTRAN) do
programa Thermath das propriedades termodinâmicas calculadas a partir da
equação de estado de Peng-Robinson (1976).
A2.1. Entrada dos dados para geração da rotina flsuvn:
deplis2 = {{Aph, {t, V, rn}}, {Uph, {t, V, rn}}, {Sph, {t, V, rn}}, {Cvph, {t, V, rn}}, {pph, {t, V, rn}}, {rmi, {t, V,
rn}}};
{deplis, deplis2} = updatedeplis[deplis, deplis2];
A = Sum[Aph[ie], {ie, 1, np}]
q = (A - Uspec)/(r*t);
Rt0 = D[q, t, NonConstants -> deplis] //. rule1 //. rule2 //. rule3 //. rule4 //. d$Aph$t[phase_] -> -Sph[phase];
Rt = Simplify[Rt0 //. Sph[phase_] -> (Uph[phase] - Aph[phase])/t //. rule1 //. rule2 //. rule3 //. rule4];
dRtdt = Simplify[D[Rt, t, NonConstants -> deplis] //. d$Uph$t[phase_] -> Cvph[phase]];
Rtnew = Rt /. Sum[Uph[ie], {ie, 1, np}] -> Sum[Uph[indepv[indv]], {indv, 1, npm1}] + Uph[kdepv[nzv]];
APÊNDICE 2
81
dRtdV0 = D[Rtnew, V[indepv[k]], NonConstants -> deplis]
dRtdV1 = dRtdV0 //. d$Uph$V[kdepv[nzv], indepv[k]] -> -d$Uph$V[kdepv[nzv]];
dRtdV = dRtdV1 //. Sum[d$Uph$V[indepv[indv], indepv[k]], {indv, 1, npm1}] -> d$Uph$V[indepv[k]];
Rtforn = Rt /. Sum[Uph[ie], {ie, 1, np}] -> Sum[Uph[indepn[i, m]], {m, 1, npm1}] + Uph[kdepn[i]];
dRtdn = D[Rtforn, rn[i, indepn[i, j]], NonConstants -> deplis] //. d$Uph$rn[kdepn[i], i, indepn[i, j]] -> d$Uph$rn[i, kdepn[i]] //. Sum[d$Uph$rn[indepn[i, m], i, indepn[i, j]], {m, 1, npm1}] -> d$Uph$rn[i, indepn[i,
j]];
qforRv = q /. Sum[Aph[ie], {ie, 1, np}] -> Sum[Aph[indepv[indv]], {indv, 1, npm1}] + Aph[kdepv[nzv]];
Rv0 = D[qforRv, V[indepv[k]], NonConstants -> deplis] //. rule1 //. rule2 //. rule3 //. rule4;
Rv1 = Rv0 //. d$Aph$V[kdepv[nzv], indepv[k]] -> pph[kdepv[nzv]];
Rv = Rv1 //. Sum[d$Aph$V[indepv[indv], indepv[k]], {indv, 1, npm1}] -> -pph[indepv[k]];
dRvdt = Simplify[D[Rv, t, NonConstants -> deplis]]//. rule4;
dRvdV0 = D[Rv, V[indepv[m]], NonConstants -> deplis];
dRvdV1 = dRvdV0 //. d$pph$V[kdepv[nzv], indepv[m]] -> -d$pph$V[kdepv[nzv]];
dRvdV = Simplify[dRvdV1 //. d$pph$V[indepv[k], indepv[m]] -> d$pph$V[indepv[k]]*rkr[indepv[k],
indepv[m]]];
dRvdn = D[Rv, rn[i, indepn[i, j]], NonConstants -> deplis] //. d$pph$rn[kdepv[nzv], i, indepn[i, j]] ->
d$pph$rn[i, kdepv[nzv]]*(rkr[kdepv[nzv], indepn[i, j]] - rkr[kdepv[nzv], kdepn[i]])//. d$pph$rn[indepv[k], i,
indepn[i, j]] -> d$pph$rn[i, indepv[k]]*(rkr[indepv[k], indepn[i, j]] - rkr[indepv[k], kdepn[i]]);
qforRnn = q /. Sum[Aph[ie], {ie, 1, np}] -> Sum[Aph[indepn[i, m]], {m, 1, npm1}] + Aph[kdepn[i]];
Rnn0 = D[qforRnn, rn[i, indepn[i, j]], NonConstants -> deplis] //. rule1 //. rule2 //. rule3 //. rule4;
Rnn1 = Rnn0 /. d$Aph$rn[kdepn[i], i, indepn[i, j]] -> -rmi[i, kdepn[i]];
Rnn = Rnn1 /. Sum[d$Aph$rn[indepn[i, m], i, indepn[i, j]], {m, 1, npm1}] -> rmi[i, indepn[i, j]];
dRnndt = Simplify[D[Rnn, t, NonConstants -> deplis]] //. rule4;
82
APÊNDICE 2
dRnndV0 = D[Rnn, V[indepv[m]], NonConstants -> deplis] //. d$rmi$V[i, indepn[i, j], indepv[m]] ->
d$rmi$V[i, indepn[i, j]]*(rkr[indepn[i, j], indepv[m]] - rkr[indepn[i, j], kdepv[nzv]]);
dRnndV = dRnndV0 //. d$rmi$V[i, kdepn[i], indepv[m]]-> d$rmi$V[i, kdepn[i]]*(rkr[kdepn[i], indepv[m]] rkr[kdepn[i], kdepv[nzv]]);
dRnndn0 = D[Rnn, rn[i1, indepn[i1, j1]], NonConstants -> deplis];
dRnndn=dRnndn0 //. (d$rmi$rn[i,indepn[i,j],i1,indepn[i1,j1]]->d$rmi$rn[i,i1,indepn[i,j]]*
(rkr[indepn[i,j],indepn[i1,j1]]-rkr[indepn[i,j],kdepn[i1]]))//.(d$rmi$rn[i,kdepn[i],i1,indepn[i1,j1]]->
d$rmi$rn[i,i1,kdepn[i]]*(rkr[kdepn[i],indepn[i1,j1]]-rkr[kdepn[i],kdepn[i1]]));
equations = {Rt, dRtdt, dRtdV, dRtdn, Rv, dRvdt, dRvdV, dRvdn, Rnn, dRnndt, dRnndV, dRnndn};
prontas = ordeq[equations, {}];
createroutine[flsuvn, prontas];
A2.2. Rotina flsuvn
subroutine flsuvn(
&
green, i$b, i$e,
&
ie$b, ie$e, indepn$b,
&
indepn$e, indepv$b, indepv$e,
&
i1$b, i1$e, j$b,
&
j$e, j1$b, j1$e,
&
k$b, kdepn$b, kdepn$e,
&
kdepv$b, kdepv$e, k$e,
&
m$b, m$e, nzv$b,
&
nzv$e, Cvph, d$rmi$rn,
&
d$rmi$t, d$rmi$V, i,
&
indepn, indepv, i1,
&
kdepn, kdepv, np,
&
pph, r, rkr,
&
rmi, t, Uph,
&
Uspec, o$27, o$28,
&
o$29, o$30, o$31,
&
o$32, o$33, o$34,
&
o$35, o$36, o$37,
&
o$38
&
)
C
C
C
C
C
C
C
C
Parameter statements used to
dimension internal variables
Please check if these values are
adequate for your application; if
they are not, change them manually.
parameter
parameter
parameter
parameter
parameter
parameter
(i_1=1)
(i_2=20)
(i1_1=1)
(i1_2=20)
(j_1=1)
(j_2=20)
83
APÊNDICE 2
parameter
parameter
parameter
parameter
parameter
parameter
parameter
parameter
C
C
C
(j1_1=1)
(j1_2=20)
(k_1=1)
(k_2=20)
(m_1=1)
(m_2=20)
(nzv_1=1)
(nzv_2=20)
Type declarations and dimension statements
implicit Double precision (a-h,o-z)
C
logical
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
&
dimension
dimension
dimension
dimension
dimension
&
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
&
dimension
dimension
dimension
&
dimension
&
C
C
C
green
green (1:12)
Uph(ie$b:ie$e)
kdepv(nzv$b:nzv$e)
kdepn(i1$b:i1$e)
indepv(m$b:m$e)
Cvph(ie$b:ie$e)
indepn(i1$b:i1$e,j1$b:j1$e)
pph(kdepv$b:kdepv$e)
d$Uph$V(kdepv$b:kdepv$e)
d$pph$V(kdepv$b:kdepv$e)
d$pph$t(kdepv$b:kdepv$e)
rmi(i$b:i$e,kdepn$b:kdepn$e)
rkr(kdepn$b:kdepn$e,kdepv$b:kdepv
$e)
d$Uph$rn(i$b:i$e,kdepn$b:kdepn$e)
d$rmi$V(i$b:i$e,kdepn$b:kdepn$e)
d$rmi$t(i$b:i$e,kdepn$b:kdepn$e)
d$pph$rn(i$b:i$e,kdepv$b:kdepv$e)
d$rmi$rn(i$b:i$e,i1$b:i1$e,kdepn$
b:kdepn$e)
w$4(nzv_1:nzv_2)
w$5(k_1:k_2)
w$7(i_1:i_2)
w$8(i_1:i_2,j_1:j_2)
w$10(i_1:i_2,j_1:j_2,nzv_1:nzv_2)
w$13(i_1:i_2,k_1:k_2)
w$14(i1_1:i1_2,i_1:i_2,j1_1:j1_2)
w$15(i_1:i_2,j_1:j_2,k_1:k_2)
w$25(i_1:i_2,m_1:m_2,nzv_1:nzv_2)
w$26(i_1:i_2,i1_1:i1_2,j1_1:j1_2)
o$29(k$b:k$e,nzv$b:nzv$e)
o$30(i$b:i$e,j$b:j$e)
o$31(k$b:k$e,nzv$b:nzv$e)
o$32(k$b:k$e,nzv$b:nzv$e)
o$33(k$b:k$e,m$b:m$e,nzv$b:nzv$e)
o$34(i$b:i$e,j$b:j$e,k$b:k$e,nzv$
b:nzv$e)
o$35(i$b:i$e,j$b:j$e)
o$36(i$b:i$e,j$b:j$e)
o$37(i$b:i$e,j$b:j$e,m$b:m$e,nzv$
b:nzv$e)
o$38(i$b:i$e,i1$b:i1$e,j$b:j$e,j1
$b:j1$e)
Internal function definition
Poneg$(a$dummy,b$dummy)=a$dummy**b$dummy
C
C
C
C
Start executable commands
84
APÊNDICE 2
if (
&
green(1)
&
green(2)
&
green(3)
&
green(4)
&
green(5)
&
green(6)
&
green(7)
&
green(8)
&
green(9)
&
green(10)
&
green(11)
&
green(12)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
C
w$3=1/r
C
end if
C
C
if (
&
green(5)
&
green(7)
&
green(8)
&
green(9)
&
green(11)
&
green(12)
&
) then
.or.
.or.
.or.
.or.
.or.
C
w$12=w$3/t
C
end if
C
C
if (
&
green(1)
&
green(3)
&
green(4)
&
green(6)
&
green(10)
&
) then
.or.
.or.
.or.
.or.
C
w$11=w$3/t**2
C
end if
C
C
if (
&
green(1)
&
green(2)
&
) then
.or.
C
s$43=dfloat(0.)
do 10043 ie=1,np
s$43=s$43 + Uph(ie)
10043 continue
w$6=s$43
C
end if
C
C
if (
&
green(5)
.or.
&
green(6)
&
) then
APÊNDICE 2
C
do 10045 nzv=nzv$b,nzv$e
w$4(nzv)=pph(kdepv(nzv))
10045 continue
do 10047 k=k$b,k$e
w$5(k)=pph(indepv(k))
10047 continue
C
end if
C
C
if (
&
green(9)
.or.
&
green(10)
&
) then
C
do 10049 i=i$b,i$e
w$7(i)=rmi(i,kdepn(i))
10049 continue
do 10051 i=i$b,i$e
do 10052 j=j$b,j$e
w$8(i,j)=rmi(i,indepn(i,j))
10052 continue
10051 continue
C
end if
C
C
if (
&
green(8)
.or.
&
green(11)
&
) then
C
do 10054 i=i$b,i$e
do 10055 j=j$b,j$e
do 10056 nzv=nzv$b,nzv$e
w$10(i,j,nzv)=rkr(indepn(i,j),kdepv(nzv))
10056 continue
10055 continue
10054 continue
do 10058 i=i$b,i$e
do 10059 k=k$b,k$e
w$13(i,k)=rkr(indepv(k),kdepn(i))
10059 continue
10058 continue
do 10061 i=i$b,i$e
do 10062 j=j$b,j$e
do 10063 k=k$b,k$e
w$15(i,j,k)=rkr(indepn(i,j),indepv(k))
10063 continue
10062 continue
10061 continue
C
end if
C
C
if (
&
green(1)
&
) then
C
o$27=w$11*(Uspec - w$6)
C
end if
C
C
85
APÊNDICE 2
if (
&
green(2)
&
) then
C
s$66=dfloat(0.)
do 10066 ie=1,np
s$66=s$66 + Cvph(ie)
10066 continue
o$28=w$3*(-(s$66*t) - 2*Uspec + 2*w$6)/t**
&
3
C
end if
C
C
if (
&
green(3)
&
) then
C
do 10068 k=k$b,k$e
do 10069 nzv=nzv$b,nzv$e
o$29(k,nzv)=w$11*(-d$Uph$V(indepv(k)) + d$
&
Uph$V(kdepv(nzv)))
10069 continue
10068 continue
C
end if
C
C
if (
&
green(4)
&
) then
C
do 10071 i=i$b,i$e
do 10072 j=j$b,j$e
o$30(i,j)=w$11*(-d$Uph$rn(i,indepn(i,j)) +
&
d$Uph$rn(i,kdepn(i)))
10072 continue
10071 continue
C
end if
C
C
if (
&
green(5)
&
) then
C
do 10074 k=k$b,k$e
do 10075 nzv=nzv$b,nzv$e
o$31(k,nzv)=w$12*(w$4(nzv) - w$5(k))
10075 continue
10074 continue
C
end if
C
C
if (
&
green(6)
&
) then
C
do 10077 k=k$b,k$e
do 10078 nzv=nzv$b,nzv$e
o$32(k,nzv)=w$11*(t*(-d$pph$t(indepv(k)) +
&
d$pph$t(kdepv(nzv))) - w$4(nzv) + w$5
&
(k))
10078 continue
86
APÊNDICE 2
10077 continue
C
end if
C
C
if (
&
green(7)
&
) then
C
do 10080 k=k$b,k$e
do 10081 m=m$b,m$e
do 10082 nzv=nzv$b,nzv$e
o$33(k,m,nzv)=-(w$12*(d$pph$V(kdepv(nzv))
&
+ d$pph$V(indepv(k))*rkr(indepv(k),ind
&
epv(m))))
10082 continue
10081 continue
10080 continue
C
end if
C
C
if (
&
green(8)
&
) then
C
do 10084 i=i$b,i$e
do 10085 j=j$b,j$e
do 10086 k=k$b,k$e
do 10087 nzv=nzv$b,nzv$e
o$34(i,j,k,nzv)=w$12*(d$pph$rn(i,kdepv(nzv
&
))*(-rkr(kdepn(i),kdepv(nzv)) + w$10(i
&
,j,nzv)) - d$pph$rn(i,indepv(k))*(-w$1
&
3(i,k) + w$15(i,j,k)))
10087 continue
10086 continue
10085 continue
10084 continue
C
end if
C
C
if (
&
green(9)
&
) then
C
do 10089 i=i$b,i$e
do 10090 j=j$b,j$e
o$35(i,j)=w$12*(-w$7(i) + w$8(i,j))
10090 continue
10089 continue
C
end if
C
C
if (
&
green(10)
&
) then
C
do 10092 i=i$b,i$e
do 10093 j=j$b,j$e
o$36(i,j)=w$11*(t*(d$rmi$t(i,indepn(i,j))
&
- d$rmi$t(i,kdepn(i))) + w$7(i) - w$8(
&
i,j))
10093 continue
87
APÊNDICE 2
10092 continue
C
end if
C
C
if (
&
green(11)
&
) then
C
do 10095 i=i$b,i$e
do 10096 m=m$b,m$e
do 10097 nzv=nzv$b,nzv$e
w$25(i,m,nzv)=-(d$rmi$V(i,kdepn(i))*(-rkr(
&
kdepn(i),kdepv(nzv)) + w$13(i,m)))
10097 continue
10096 continue
10095 continue
do 10099 i=i$b,i$e
do 10100 j=j$b,j$e
do 10101 m=m$b,m$e
do 10102 nzv=nzv$b,nzv$e
o$37(i,j,m,nzv)=w$12*(d$rmi$V(i,indepn(i,j
&
))*(-w$10(i,j,nzv) + w$15(i,j,m)) + w$
&
25(i,m,nzv))
10102 continue
10101 continue
10100 continue
10099 continue
C
end if
C
C
if (
&
green(12)
&
) then
C
do 10104 i1=i1$b,i1$e
do 10105 i=i$b,i$e
do 10106 j1=j1$b,j1$e
w$14(i1,i,j1)=rkr(indepn(i1,j1),kdepn(i))
10106 continue
10105 continue
10104 continue
do 10108 i=i$b,i$e
do 10109 i1=i1$b,i1$e
do 10110 j1=j1$b,j1$e
w$26(i,i1,j1)=-(d$rmi$rn(i,i1,kdepn(i))*(&
rkr(kdepn(i),kdepn(i1)) + w$14(i1,i,j1
&
)))
10110 continue
10109 continue
10108 continue
do 10112 i=i$b,i$e
do 10113 i1=i1$b,i1$e
do 10114 j=j$b,j$e
do 10115 j1=j1$b,j1$e
o$38(i,i1,j,j1)=w$12*(d$rmi$rn(i,i1,indepn
&
(i,j))*(rkr(indepn(i,j),indepn(i1,j1))
&
- w$14(i,i1,j)) + w$26(i,i1,j1))
10115 continue
10114 continue
10113 continue
10112 continue
C
end if
88
APÊNDICE 2
C
C
return
end
C
C
The subroutine was successfully created
A2.3. Entrada dos dados para geração da rotina penrob:
deplis2={{am,{t,x}},{bm,{x}}};
{deplis,deplis2}=updatedeplis[deplis,deplis2];
(* Peng-Robinson EOS *)
p=r*t/(v-bm)-am/(v^2+2*bm*v-bm^2) //.rule1 //.rule2 //.rule4
dpdv=Expand[D[p,v]] //.rule1 //.rule2 //.rule4
dpdrho=Expand[-v^2*(dpdv)] //.rule1 //.rule2 //.rule4
d2pdv2=Expand[D[dpdv,v]] //.rule1 //.rule2 //.rule4
d2pdrho2=Expand[v^3*(2*dpdv+v*d2pdv2)] //.rule1 //.rule2 //.rule4
dpdt=Expand[D[p,t,NonConstants->deplis]] //.rule1 //.rule2 //.rule4
lnfimix=Simplify[resid[g]/(r*t) //.rule1 //.rule2 //.rule4
(* Residual enthalpy *)
hr=Simplify[resid[h] //.rule1 //.rule2 //.rule4 //.rule4
dhrdt=Expand[D[hr,t,NonConstants->deplis] //.rule1 //.rule2 //.rule4
(* Residual internal energy *)
ur=Simplify[resid[u] //.rule1 //.rule2 //.rule4 //.rule4
cvr=D[ur,t,NonConstants->deplis]
lnfii=Simplify[Pmp[lnfimix] //.rule1 //.rule2 //.rule3]//.rule4
rmi = Log[x[k$0]] + lnfii + Log[p])],
(* derivative of rmi with respect to t at constant molar volume and mole fraction *)
drmidt=Simplify[D[rmi,t,NonConstants->deplis]] //.rule1 //.rule2 //.rule3//.rule4
(* derivative of rmi with respect to molar volume at constant temperature and mole fraction *)
89
90
APÊNDICE 2
drmidv=Simplify[D[rmi,v,NonConstants->deplis]] //.rule1 //.rule2 //.rule3 //.rule4
(* derivative of ln fii with respect to mole numbers at constant temperature and pressure *)
dlnfiidnjatp=Dtransf[lnfii,1,p,v] //.rule1 //.rule2 //.rule3 //.rule4
(* derivative of rmi with respect to mole numbers at constant temperature and total volume *)
drmidnjatV=Dtransf[rmi,1,V,v] //.rule1 //.rule2 //.rule3 //.rule4\
(* derivative of p with respect to mole numbers at constant temperature and total volume
*)
dpdnatV=Dtransf[p,1,V,v] //.rule1 //.rule2 //.rule3 //.rule4
durdv=Simplify[D[ur, v, NonConstants->deplis]]//. rule1//.rule2 //. rule3//. rule4
expressions={p, dpdv, d2pdv2, dpdrho, d2pdrho2, dpdt, lnfimix, hr, dhrdt, ur, cvr, lnfii, rmi, drmidt, drmidv,
dlnfiidnjatp, drmidnjatV, durdv, dpdnatV};
analyzedexpressions=ordeq[expressions,{t,x,v},deplis2];
createroutine[penrob,analyzedexpressions];
A2.4. Rotina penrob
subroutine penrob(
&
hierch, green, k$0$b,
&
k$0$e, k$1$b, k$1$e,
&
am, bm, d$1am,
&
d$1bm, d$2am, d$2bm,
&
d$am$t, d$d$1am$t, d$d$am$t$t,
&
r, rkr, t,
&
v, x, o$313,
&
o$314, o$315, o$316,
&
o$317, o$318, o$319,
&
o$320, o$321, o$322,
&
o$323, o$324, o$325,
&
o$326, o$327, o$328,
&
o$329, o$330, o$331
&
)
C
C
C
C
C
C
C
C
Parameter statements used to
dimension internal variables
Please check if these values are
adequate for your application; if
they are not, change them manually.
parameter
parameter
parameter
parameter
C
C
C
(k$0_1=1)
(k$0_2=20)
(k$1_1=1)
(k$1_2=20)
Type declarations and dimension statements
91
APÊNDICE 2
implicit Double precision (a-h,o-z)
C
logical
logical
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
dimension
C
C
C
green
hierch
green (1:19)
hierch(0:3)
x(k$0$b:k$0$e)
d$1bm(k$1$b:k$1$e)
d$1am(k$1$b:k$1$e)
d$d$1am$t(k$0$b:k$0$e)
d$2bm(k$0$b:k$0$e,k$1$b:k$1$e)
d$2am(k$0$b:k$0$e,k$1$b:k$1$e)
rkr(k$0$b:k$0$e,k$1$b:k$1$e)
w$50(k$0_1:k$0_2,k$1_1:k$1_2)
w$80(k$0_1:k$0_2)
w$99(k$1_1:k$1_2)
w$100(k$0_1:k$0_2)
w$125(k$0_1:k$0_2,k$1_1:k$1_2)
w$126(k$0_1:k$0_2,k$1_1:k$1_2)
w$231(k$1_1:k$1_2)
w$234(k$0_1:k$0_2)
w$266(k$1_1:k$1_2)
w$282(k$0_1:k$0_2)
w$291(k$1_1:k$1_2)
w$294(k$0_1:k$0_2)
w$297(k$1_1:k$1_2)
w$298(k$0_1:k$0_2)
w$302(k$0_1:k$0_2)
w$307(k$0_1:k$0_2)
w$310(k$0_1:k$0_2)
w$311(k$0_1:k$0_2,k$1_1:k$1_2)
w$312(k$0_1:k$0_2)
o$324(k$0$b:k$0$e)
o$325(k$0$b:k$0$e)
o$326(k$0$b:k$0$e)
o$327(k$0$b:k$0$e)
o$328(k$0$b:k$0$e,k$1$b:k$1$e)
o$329(k$0$b:k$0$e,k$1$b:k$1$e)
o$331(k$1$b:k$1$e)
Save automatically created variables
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
w$1
w$2
w$3
w$7
w$9
w$10
w$11
w$23
w$24
w$25
w$27
w$28
w$32
w$33
w$34
w$35
w$36
w$38
w$43
w$47
w$48
w$49
w$50
92
APÊNDICE 2
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
save
C
C
C
w$51
w$59
w$60
w$61
w$72
w$73
w$78
w$80
w$83
w$85
w$87
w$89
w$91
w$92
w$93
w$95
w$97
w$98
w$99
w$100
w$112
w$124
w$125
w$126
w$129
w$131
w$152
Internal function definition
Poneg$(a$dummy,b$dummy)=a$dummy**b$dummy
C
C
C
C
C
C
Start executable commands
Open hierarchy level 0
if (hierch(0)) then
C
if (
&
green(7)
&
green(8)
&
green(9)
&
green(10)
&
green(11)
&
green(12)
&
green(13)
&
green(14)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
C
w$3=dble(1.4142135623730951)
w$60=1 + w$3
w$124=1 - w$3
C
end if
C
C
if (
&
green(7)
&
green(12)
&
green(13)
&
green(14)
&
green(15)
.or.
.or.
.or.
.or.
.or.
93
APÊNDICE 2
&
&
&
green(16)
green(17)
) then
.or.
C
w$28=1/r
C
end if
C
C
if (
&
green(7)
&
green(8)
&
green(9)
&
green(10)
&
green(11)
&
) then
.or.
.or.
.or.
.or.
C
w$32=dble(0.7071067811865476)
C
end if
C
end if
C
C
C
C
C
C
Closed hierarchy level 0
Open hierarchy level 1
if (hierch(1)) then
C
if (
&
green(1)
&
green(2)
&
green(3)
&
green(4)
&
green(5)
&
green(7)
&
green(8)
&
green(12)
&
green(13)
&
green(14)
&
green(15)
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
C
w$2=r*t
C
end if
C
C
if (
&
green(7)
&
green(12)
&
green(13)
&
green(15)
&
green(16)
&
green(17)
&
) then
C
w$27=1/t
C
end if
C
.or.
.or.
.or.
.or.
.or.
94
APÊNDICE 2
C
if (
&
green(7)
&
green(15)
&
) then
.or.
C
w$83=w$27*w$28
C
end if
C
C
if (
&
green(16)
&
green(17)
&
) then
.or.
C
w$38=-2*w$2
C
end if
C
C
if (
&
green(14)
&
) then
C
w$34=4*t
C
end if
C
end if
C
C
C
C
C
C
Closed hierarchy level 1
Open hierarchy level 2
if (hierch(2)) then
C
if (
&
green(7)
&
green(8)
&
green(9)
&
green(10)
&
green(11)
&
green(12)
&
green(13)
&
green(14)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
C
w$112=bm*w$60
w$152=bm*w$124
C
end if
C
C
if (
&
green(7)
&
green(8)
&
green(9)
&
green(10)
&
green(11)
&
green(16)
&
green(17)
.or.
.or.
.or.
.or.
.or.
.or.
95
APÊNDICE 2
&
) then
C
w$30=1/bm
C
end if
C
C
if (
&
green(9)
&
green(12)
&
green(13)
&
green(14)
&
green(15)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
C
w$29=bm**2
C
end if
C
C
if (
&
green(5)
&
green(12)
&
green(13)
&
green(14)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
.or.
C
w$1=4*bm
C
end if
C
C
if (
&
green(7)
&
green(8)
&
green(10)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
C
w$59=w$30*dble(-0.5)
C
end if
C
C
if (
&
green(2)
&
green(15)
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
.or.
.or.
C
w$7=2*am
C
end if
C
C
if (
&
green(12)
&
green(13)
.or.
.or.
96
APÊNDICE 2
&
&
&
&
green(14)
green(16)
green(17)
) then
.or.
.or.
C
w$11=2*bm
w$51=am*w$3
w$82=w$28/bm**2
C
end if
C
C
if (
&
green(8)
&
green(10)
&
green(14)
&
green(18)
&
) then
.or.
.or.
.or.
C
w$78=am - d$am$t*t
C
end if
C
C
if (
&
green(4)
&
green(15)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
C
w$25=-2*am
C
end if
C
C
if (
&
green(12)
&
green(13)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
C
w$36=bm*w$2
w$89=w$2*w$29
w$95=-(w$2*w$29)
do 10358 k$0=k$0$b,k$0$e
w$100(k$0)=-am - d$1am(k$0)
10358 continue
C
end if
C
C
if (
&
green(3)
.or.
&
green(5)
.or.
&
green(14)
&
) then
C
w$24=-8*bm
C
end if
C
C
if (
97
APÊNDICE 2
&
&
&
&
green(7)
green(8)
green(10)
) then
.or.
.or.
C
w$87=w$32*w$59
C
end if
C
C
if (
&
green(12)
&
green(13)
&
green(14)
&
) then
.or.
.or.
C
w$86=w$82*dble(0.25)
C
end if
C
C
if (
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
C
w$9=-2*bm
C
end if
C
C
if (
&
green(8)
&
green(10)
&
) then
.or.
C
w$131=w$78*w$87
C
end if
C
C
if (
&
green(9)
&
green(11)
&
) then
.or.
C
w$98=d$d$am$t$t*t*w$30*w$32*dble(0.5)
C
end if
C
C
if (
&
green(9)
&
green(14)
&
) then
.or.
C
w$61=r*w$29
C
end if
C
C
if (
&
green(12)
&
green(13)
.or.
98
APÊNDICE 2
&
) then
C
w$93=w$27*w$86
C
end if
C
C
if (
&
green(16)
&
green(17)
&
) then
.or.
C
w$10=-3*bm
w$23=3*bm
w$46=w$2*w$9
do 10371 k$0=k$0$b,k$0$e
do 10372 k$1=k$1$b,k$1$e
w$50(k$0,k$1)=d$1bm(k$0) + d$2bm(k$0,k$1)
10372 continue
10371 continue
w$72=-am + w$36
w$73=w$25 + w$36
w$91=w$27*w$82
do 10377 k$1=k$1$b,k$1$e
w$99(k$1)=w$46*d$1bm(k$1)
10377 continue
do 10379 k$0=k$0$b,k$0$e
do 10380 k$1=k$1$b,k$1$e
w$126(k$0,k$1)=-d$1am(k$0) - d$1am(k$1) &
d$2am(k$0,k$1)
10380 continue
10379 continue
w$129=-(bm*w$60)
C
end if
C
C
if (
&
green(9)
&
) then
C
w$33=bm*d$am$t
w$43=2*bm*r
C
end if
C
C
if (
&
green(14)
&
) then
C
w$47=w$2*w$24
w$48=d$am$t*t*w$11
w$49=-(bm*d$am$t*t)
w$92=bm*w$51
w$97=w$86/t**2
C
end if
C
C
if (
&
green(15)
&
) then
C
w$35=bm**3
APÊNDICE 2
w$85=-4*w$29
C
end if
C
C
if (
&
green(17)
&
) then
C
do 10392 k$0=k$0$b,k$0$e
w$80(k$0)=1/x(k$0)
10392 continue
do 10394 k$0=k$0$b,k$0$e
do 10395 k$1=k$1$b,k$1$e
w$125(k$0,k$1)=rkr(k$0,k$1) - x(k$0)
10395 continue
10394 continue
C
end if
C
end if
C
C
Closed hierarchy level 2
C
C
C
Open hierarchy level 3
C
if (hierch(3)) then
C
if (
&
green(1)
.or.
&
green(2)
.or.
&
green(3)
.or.
&
green(4)
.or.
&
green(5)
.or.
&
green(6)
.or.
&
green(7)
.or.
&
green(8)
.or.
&
green(9)
.or.
&
green(12)
.or.
&
green(13)
.or.
&
green(14)
.or.
&
green(15)
.or.
&
green(16)
.or.
&
green(17)
.or.
&
green(18)
.or.
&
green(19)
&
) then
C
w$16=v**2
C
end if
C
C
if (
&
green(1)
.or.
&
green(2)
.or.
&
green(3)
.or.
&
green(4)
.or.
&
green(5)
.or.
&
green(6)
.or.
&
green(8)
.or.
&
green(12)
.or.
&
green(13)
.or.
&
green(14)
.or.
99
100
APÊNDICE 2
&
&
&
&
&
&
green(15)
green(16)
green(17)
green(18)
green(19)
) then
.or.
.or.
.or.
.or.
C
w$18=2*v
w$62=-bm + w$18
C
end if
C
C
if (
&
green(1)
&
green(2)
&
green(3)
&
green(4)
&
green(5)
&
green(6)
&
green(8)
&
green(13)
&
green(14)
&
green(15)
&
green(16)
&
green(17)
&
green(18)
&
green(19)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
C
w$141=w$16 + bm*w$62
C
end if
C
C
if (
&
green(1)
&
green(2)
&
green(3)
&
green(4)
&
green(5)
&
green(6)
&
green(7)
&
green(12)
&
green(13)
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
C
w$56=-bm + v
C
end if
C
C
if (
&
green(7)
&
green(8)
&
green(9)
&
green(10)
&
green(11)
&
green(12)
&
green(13)
&
green(14)
&
green(16)
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
.or.
101
APÊNDICE 2
&
&
green(17)
) then
C
w$151=v + w$112
w$183=1/(v + w$152)
w$198=w$151*w$183
w$204=Log(w$198)
C
end if
C
C
if (
&
green(5)
&
green(8)
&
green(9)
&
green(12)
&
green(13)
&
green(15)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
.or.
C
w$15=v**3
C
end if
C
C
if (
&
green(1)
&
green(5)
&
green(6)
&
green(13)
&
green(16)
&
green(17)
&
green(18)
&
green(19)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
.or.
C
w$160=1/w$141
C
end if
C
C
if (
&
green(2)
&
green(3)
&
green(4)
&
green(5)
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
.or.
.or.
.or.
.or.
C
w$166=w$141**(-2)
C
end if
C
C
if (
&
green(8)
&
green(9)
&
green(12)
&
green(13)
&
green(15)
&
green(16)
.or.
.or.
.or.
.or.
.or.
.or.
102
APÊNDICE 2
&
&
green(17)
) then
C
w$177=w$15 + bm*(bm*(bm - 3*v) + w$16)
C
end if
C
C
if (
&
green(2)
&
green(4)
&
green(5)
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
.or.
.or.
.or.
C
w$116=w$56**(-2)
C
end if
C
C
if (
&
green(2)
&
green(4)
&
green(15)
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
.or.
.or.
.or.
C
w$13=bm + v
C
end if
C
C
if (
&
green(8)
&
green(9)
&
green(12)
&
green(13)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
.or.
C
w$196=1/w$177
C
end if
C
C
if (
&
green(7)
&
green(12)
&
green(13)
&
green(14)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
.or.
C
w$19=-2*v
w$148=-w$16 + bm*(bm + w$19)
C
end if
C
C
103
APÊNDICE 2
if (
&
green(1)
&
green(5)
&
green(6)
&
green(13)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
C
w$115=1/w$56
C
end if
C
C
if (
&
green(4)
&
green(5)
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
.or.
.or.
C
w$155=w$116*w$2
C
end if
C
C
if (
&
green(7)
&
green(8)
&
green(15)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
C
w$58=bm - v
C
end if
C
C
if (
&
green(7)
&
green(12)
&
green(13)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
.or.
C
w$174=w$148*w$2
w$197=1/(w$174 + am*w$56)
w$224=Log(w$174*w$197)
C
end if
C
C
if (
&
green(12)
&
green(13)
&
green(14)
&
green(16)
&
green(17)
&
) then
C
w$17=am*v
w$41=-(am*v)
w$54=am + v*w$2
.or.
.or.
.or.
.or.
104
APÊNDICE 2
C
end if
C
C
if (
&
green(7)
&
green(14)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
C
w$170=1/w$148
C
end if
C
C
if (
&
green(2)
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
.or.
C
w$195=-(w$116*w$2) + w$13*w$166*w$7
C
end if
C
C
if (
&
green(12)
&
green(13)
&
green(16)
&
green(17)
&
) then
.or.
.or.
.or.
C
w$67=v + w$11
w$103=w$17 + w$95
w$180=v*(w$41 + bm*w$54) + w$62*w$89
w$226=w$177*w$204
w$240=w$226*w$3
w$254=w$1*w$180 + w$226*w$51
w$269=bm*w$103 + w$177*w$2*w$224 + v*(w$41
&
+ w$36*w$67)
w$279=w$1*w$269
do 10434 k$0=k$0$b,k$0$e
w$282(k$0)=w$279 + w$240*w$100(k$0)
10434 continue
do 10436 k$0=k$0$b,k$0$e
w$294(k$0)=w$254*d$1bm(k$0) + bm*w$282(k$0
&
)
10436 continue
C
end if
C
C
if (
&
green(1)
.or.
&
green(13)
.or.
&
green(17)
&
) then
C
w$194=-(am*w$160) + w$115*w$2
C
end if
C
105
APÊNDICE 2
C
if (
&
green(8)
&
green(16)
&
green(17)
&
) then
.or.
.or.
C
w$119=am*w$58
C
end if
C
C
if (
&
green(16)
&
green(17)
&
green(19)
&
) then
.or.
.or.
C
w$185=am*w$160*(w$18 + w$9)
do 10441 k$1=k$1$b,k$1$e
w$231(k$1)=w$155*d$1bm(k$1) + w$160*(-d$1a
&
m(k$1) + w$185*d$1bm(k$1))
10441 continue
C
end if
C
C
if (
&
green(3)
.or.
&
green(5)
&
) then
C
w$149=-8*w$16 + bm*(-16*v + w$24)
w$184=am*w$166
C
end if
C
C
if (
&
green(8)
.or.
&
green(9)
&
) then
C
w$217=v*w$196
C
end if
C
C
if (
&
green(8)
.or.
&
green(10)
&
) then
C
w$249=w$131*w$204
C
end if
C
C
if (
&
green(4)
.or.
&
green(15)
&
) then
C
w$159=w$13*w$25
C
106
APÊNDICE 2
end if
C
C
if (
&
green(9)
&
green(11)
&
) then
.or.
C
w$250=w$204*w$98
C
end if
C
C
if (
&
green(12)
&
green(13)
&
) then
.or.
C
w$223=w$196*w$93
do 10450 k$0=k$0$b,k$0$e
w$302(k$0)=w$223*w$294(k$0)
10450 continue
C
end if
C
C
if (
&
green(16)
.or.
&
green(17)
&
) then
C
w$55=w$11 + w$19
w$127=(v + w$1)*w$2
w$139=-am + w$2*w$55
w$143=3*w$16 + bm*(w$10 + w$18)
w$147=w$16 + bm*(-6*v + w$23)
w$171=1/w$151
w$175=w$36*(-bm + 2*w$62)
w$191=w$170*w$177
w$209=w$196*dble(-0.25)
w$210=w$143*w$209
w$211=-(w$148*w$197)
w$221=w$196*w$91
w$225=w$147*w$209 + w$59
w$229=w$147*w$204
w$239=am*w$229
w$241=w$171*(-(w$124*w$198) + w$60)
w$242=am*w$241
w$251=w$147*w$224
w$255=w$171*w$177*(1 + (-v + w$129)*w$183)
&
+ w$143*w$204
w$256=w$229 + w$177*w$241
w$265=w$255*w$3
do 10473 k$1=k$1$b,k$1$e
w$266(k$1)=w$256*d$1bm(k$1)
10473 continue
w$277=w$255*w$51 + w$1*(bm*(am + w$2*w$67)
&
+ v*w$73)
do 10476 k$1=k$1$b,k$1$e
w$291(k$1)=w$3*(w$239*d$1bm(k$1) + w$177*(
&
w$204*d$1am(k$1) + w$242*d$1bm(k$1)))
&
+ 4*(w$180*d$1bm(k$1) + bm*(w$175*d$1b
&
m(k$1) + v*(w$58*d$1am(k$1) + w$54*d$1
&
bm(k$1))))
10476 continue
APÊNDICE 2
w$292=w$1*(w$119 + w$2*(w$143*w$224 + w$19
&
1*(-2*w$13 + w$211*(am + w$13*w$38)) +
&
bm*w$67) + v*w$72)
do 10479 k$1=k$1$b,k$1$e
w$297(k$1)=4*(w$269*d$1bm(k$1) + bm*(w$103
&
*d$1bm(k$1) + v*(-(v*d$1am(k$1)) + w$1
&
27*d$1bm(k$1)) + w$2*(w$251*d$1bm(k$1)
&
+ w$191*(w$55*d$1bm(k$1) + w$211*(w$5
&
6*d$1am(k$1) + w$139*d$1bm(k$1)))) + b
&
m*(v*d$1am(k$1) + w$99(k$1))))
10479 continue
do 10481 k$0=k$0$b,k$0$e
w$298(k$0)=w$225*w$294(k$0)
10481 continue
do 10483 k$0=k$0$b,k$0$e
w$307(k$0)=dble(0.25)*(w$277*d$1bm(k$0) +
&
bm*(w$292 + w$265*w$100(k$0))) + w$210
&
*w$294(k$0)
10483 continue
do 10485 k$0=k$0$b,k$0$e
do 10486 k$1=k$1$b,k$1$e
w$311(k$0,k$1)=d$1bm(k$1)*w$298(k$0) + dbl
&
e(0.25)*(d$1bm(k$1)*w$282(k$0) + d$1bm
&
(k$0)*w$291(k$1) + bm*(w$3*(w$226*w$12
&
6(k$0,k$1) + w$100(k$0)*w$266(k$1)) +
&
w$297(k$1)) + w$254*w$50(k$0,k$1))
10486 continue
10485 continue
C
end if
C
C
if (
&
green(1)
&
) then
C
o$313=w$194
C
end if
C
C
if (
&
green(2)
&
) then
C
o$314=w$195
C
end if
C
C
if (
&
green(3)
&
) then
C
o$315=(am*w$149)/w$141**3 + 2*(w$184 + w$2
&
/w$56**3)
C
end if
C
C
if (
&
green(4)
&
) then
C
o$316=w$16*(w$155 + w$159*w$166)
107
APÊNDICE 2
C
end if
C
C
if (
&
green(5)
&
) then
C
o$317=w$15*(w$155*(-2 + w$115*w$18) + (w$1
&
+ v*(6 + w$149*w$160))*w$184)
C
end if
C
C
if (
&
green(6)
&
) then
C
o$318=r*w$115 - d$am$t*w$160
C
end if
C
C
if (
&
green(7)
&
) then
C
o$319=w$83*((w$2*(-bm + w$224*w$58))/w$58
&
+ am*(v*w$170 + w$204*w$87))
C
end if
C
C
if (
&
green(8)
&
) then
C
o$320=-w$2 + (w$119 + w$141*w$2)*w$217 + w
&
$249
C
end if
C
C
if (
&
green(9)
&
) then
C
w$222=w$16*w$196
o$321=-r + r*w$15*w$196 - d$am$t*w$222 + w
&
$250 + w$217*w$33 + w$222*w$43 - v*w$1
&
96*w$61
C
end if
C
C
if (
&
green(10)
&
) then
C
o$322=w$249
C
end if
C
C
if (
108
APÊNDICE 2
&
&
green(11)
) then
C
o$323=w$250
C
end if
C
C
if (
&
green(12)
&
) then
C
do 10500 k$0=k$0$b,k$0$e
o$324(k$0)=w$302(k$0)
10500 continue
C
end if
C
C
if (
&
green(13)
&
) then
C
do 10502 k$0=k$0$b,k$0$e
o$325(k$0)=Log(w$194*x(k$0)) + w$302(k$0)
10502 continue
C
end if
C
C
if (
&
green(14)
&
) then
C
w$84=d$am$t*w$16
w$154=v*(w$41 + w$48) + bm*(-2*w$17 + w$49
&
)
w$189=w$170*w$97
w$237=w$204*w$3
w$268=bm*(v*(w$47 - 4*w$54) + w$34*(d$am$t
&
*v + w$61) + w$204*w$92)
w$274=-((v*w$1 + w$148*w$237)*w$78)
do 10510 k$0=k$0$b,k$0$e
o$326(k$0)=w$189*(w$274*d$1bm(k$0) + bm*(w
&
$268 + w$237*(w$154 + w$148*d$1am(k$0)
&
+ t*(w$84 + w$141*d$d$1am$t(k$0)))))
10510 continue
C
end if
C
C
if (
&
green(15)
&
) then
C
w$165=w$141**2
w$173=(-bm + 4*v)*w$35 + w$16*(w$16 + w$85
&
)
w$200=-(w$165*w$2) + w$58**3*w$7
w$215=w$83/w$177**2
w$216=-(v*(w$165*w$2 + w$159*w$58**2))
do 10517 k$0=k$0$b,k$0$e
o$327(k$0)=w$215*(w$216 + w$173*d$1am(k$0)
&
+ w$200*d$1bm(k$0))
10517 continue
109
APÊNDICE 2
C
end if
C
C
if (
&
green(16)
&
) then
C
do 10519 k$0=k$0$b,k$0$e
w$310(k$0)=-(w$307(k$0)/w$195)
10519 continue
do 10521 k$0=k$0$b,k$0$e
do 10522 k$1=k$1$b,k$1$e
o$328(k$0,k$1)=w$221*(w$231(k$1)*w$310(k$0
&
) + w$311(k$0,k$1))
10522 continue
10521 continue
C
end if
C
C
if (
&
green(17)
&
) then
C
w$207=1/w$194
do 10525 k$0=k$0$b,k$0$e
w$234(k$0)=w$207*w$80(k$0)
10525 continue
w$235=w$195*w$207
do 10528 k$0=k$0$b,k$0$e
w$312(k$0)=-(v*(w$235 + w$221*w$307(k$0)))
10528 continue
do 10530 k$0=k$0$b,k$0$e
do 10531 k$1=k$1$b,k$1$e
o$329(k$0,k$1)=w$221*w$311(k$0,k$1) + w$31
&
2(k$0) + w$234(k$0)*(w$194*w$125(k$0,k
&
$1) + w$231(k$1)*x(k$0))
10531 continue
10530 continue
C
end if
C
C
if (
&
green(18)
&
) then
C
o$330=w$160*w$78
C
end if
C
C
if (
&
green(19)
&
) then
C
w$214=-(v*w$195)
do 10535 k$1=k$1$b,k$1$e
o$331(k$1)=w$214 + w$231(k$1)
10535 continue
C
end if
C
end if
110
APÊNDICE 2
C
C
C
C
Closed hierarchy level 3
return
end
C
C
The subroutin was successfully created
111