Universidade Federal de São Carlos
Campus São Carlos
Departamento de Engenharia Química
Kenji Urazaki Junior
Treinamento em simulação de processos químicos utilizando o
programa científico Scilab
São Carlos
2013
1
Kenji Urazaki Junior
Treinamento em simulação de processos químicos utilizando o programa
científico Scilab
Apresentação do projeto desenvolvido no
Programa Jovens Talentos para Ciência, ação
fomentada pela Capes e CNPq, como
treinamento para iniciação científica.
Orientador: Prof. Dr. Wu Hong Kwong.
São Carlos
2013
2
Sumário
1 INTRODUÇÃO ....................................................................................................................... 4
1.1 Do programa JTPC ........................................................................................................... 4
1.2 Do projeto ......................................................................................................................... 4
1.3 Do Ten Problems .............................................................................................................. 4
1.4 Do Scilab ........................................................................................................................... 5
2 OBJETIVOS ............................................................................................................................ 5
3 REVISÃO BIBLIOGRÁFICA E TEÓRICA .......................................................................... 5
3.1 Problema 1, ....................................................................................................................... 6
3.2 Problema 2 ........................................................................................................................ 7
3.3 Problema 4 ...................................................................................................................... 10
3.4 Problema 5 ...................................................................................................................... 10
3.5 Problema 6 ...................................................................................................................... 12
3.6 Problema 7 ...................................................................................................................... 13
3.7 Problema 8 ...................................................................................................................... 14
4 MÉTODOS NUMÉRICOS ................................................................................................... 16
4.1 Sobre os Métodos Numéricos ......................................................................................... 16
4.2 Método da bissecção ....................................................................................................... 16
4.3 Método Newton-Raphson ............................................................................................... 18
4.4 Divisão pela esquerda ..................................................................................................... 18
4.5 Método do Gauss-Jacobi ................................................................................................. 18
4.6 Método de Euler .............................................................................................................. 20
4.7 Solucionadores ................................................................................................................ 21
5 RESULTADOS ..................................................................................................................... 21
5.1 Problema 1 ...................................................................................................................... 21
5.2 Problema 2 ...................................................................................................................... 23
5.3 Problema 4 ...................................................................................................................... 23
5.4 Problema 5 ...................................................................................................................... 23
5.5 Problema 6 ...................................................................................................................... 24
5.6 Problema 7 ...................................................................................................................... 25
5.7 Problema 8 ...................................................................................................................... 26
6 CONCLUSÃO ....................................................................................................................... 26
3
7 TRABALHOS FUTUROS .................................................................................................... 27
8 APÊNDICE ........................................................................................................................... 27
8.1 Linhas de código de programação dos problemas .......................................................... 27
Problema 1......................................................................................................................... 27
Problema 2......................................................................................................................... 28
Problema 4......................................................................................................................... 29
Problema 5......................................................................................................................... 30
Problema 6......................................................................................................................... 31
Problema 7......................................................................................................................... 32
Problema 8......................................................................................................................... 33
9 REFERÊNCIA BLIBIOGRÁFICA ....................................................................................... 33
4
1 INTRODUÇÃO
1.1 Do programa JTPC
No ano de 2012 por uma iniciativa da Capes e do CNPq realizou-se um programa de
incentivo à iniciação científica que consistia em inserir os estudantes, que ingressaram no
mesmo ano em Universidades e institutos federais, no meio científico. O programa Jovens
Talentos para Ciência selecionou por uma prova de conhecimentos gerais estudantes de
graduação de todas as áreas do conhecimento, eram ofertadas seis mil bolsas de estudo,
dividas entre todo o Brasil, de acordo do número de inscritos em cada Universidade. Após a
vigência de 12 meses da bolsa, os candidatos estariam preparados para serem bolsistas do
Pibid, do Programa Ciências sem Fronteiras e outros programas de iniciação científica.
1.2 Do projeto
Este projeto teve caráter de treinamento em solução de problemas de Engenharia
Química com auxílio de programação computacional. Ao final saber-se-ia programar e
estruturar soluções para problemas matemáticos com aplicações na Engenharia. O software
científico utilizado foi o Scilab e os problemas foram retirados do artigo Ten problems. A
maioria dos problemas foi resolvida com métodos numéricos, que foram estudados no projeto
(os seus algoritmos deveriam ser inseridos e executados no software). Ao final do projeto, os
resultados dos problemas foram discutidos.
1.3 Do Ten Problems
O Artigo é uma coleção de 10 problemas de Engenharia Química que podem ser
solucionados por um software matemático utilizando cálculo numérico. Os problemas traziam
as informações dos métodos numéricos que deviam ser usados para solucioná-los, os
conceitos de engenharia utilizados e o problema em si.
5
1.4 Do Scilab
O Scilab (Scientific Laboratory) é um programa para cálculo científico gratuito desde
1994 e desenvolvido desde 1990 pelo INRIA (“Institut Nationale de Recherche em
Informatique et en Automatique”) e o ENPC (“Ecole Nationale des Ponts et Chaussée”) na
França. O Scilab é um programa open source, possui uma vasta biblioteca com funções
matemáticas, polinômios, listas, sistemas lineares, solucionadores numéricos; e sua
programação é feita em linha admitindo a incorporação de programas em linguagens com C e
Fortran. O Scilab é muito utilizado internacionalmente em ambientes acadêmicos e industriais
e está em constante atualização. Suas funções são similares ao MATLAB, mas sua utilização
neste projeto foi motivada por ele ser um open source.
2 OBJETIVOS
O projeto teve como objetivo o aprendizado de como a computação numérica auxilia
na resolução de problemas de Engenharia Química, tendo como ponto central o aprendizado
dos métodos numéricos e a familiarização com o software científico.
3 REVISÃO BIBLIOGRÁFICA E TEÓRICA
Os problemas a serem resolvidos eram dez, neste projeto sete foram resolvidos: dois
de introdução à Engenharia Química, um de termodinâmica, um de dinâmica de fluidos, um
de transferência de calor, um de transferência de massa, um de processo de separação. Todos
os problemas envolviam algum conceito de engenharia química, assuntos abordados ao
decorrer do curso de graduação.
A seguir os problemas presentes no Ten problems e resolvidos no projeto serão
enunciados, com numeração semelhante ao artigo, apresentando o título, área na engenharia
química e problema matemático a ser modelado. Depois de enunciado, cada problema
recebera uma breve introdução dos conceitos utilizados, todos retirados dos livros de
engenharia química.
6
3.1 Problema 1, Volume molar e fator de compressibilidade a partir da equação de van
der Waals, introdução à Engenharia Química, solução de uma equação não linear.
No problema o volume molar e fator de compressibilidade de um gás deveriam ser
calculados com o uso da equação de van der Waals.
A lei de gases ideais que relaciona Pressão Volume e Temperatura se aplica bem em
pequenas pressões (aproximadamente a da atmosfera). Para pressões mais altas, uma equação
mais complexa deve ser utilizada, o que normalmente necessita de uma solução numérica.
A equação de van der Waals utilizada é dada por:
(P+a/V²)(V-b) = RT
Onde
a= (27/64)(R²Tc²/Pc)
b= RTc/8Pc
E as variáveis são:
P = Pressão (atm);
V = Volume molar (L/gmol);
T = Temperatura em Kelvin;
R = constante dos gases (0.08206 atmL/gmolK);
Tc = Temperatura crítica (405.5 K para a amônia);
Pc = Pressão crítica (111.3 atm para a amônia).
O termo a/V² é a correção para a atração entre as moléculas do gás (a/V²= 0 quando
trabalhamos com um gás ideal) e o termo b é a correção para o volume ocupado pelas
moléculas de gás. As constantes a e b são determinadas experimentalmente para cada gás.
Pressão reduzida é definida por:
Pr = P/Pc
O fator de compressibilidade de um gás é dado por:
7
Z = PV/RT
No problema teve-se de calcular o volume molar e o fator de compressibilidade para o
gás amônia a pressão de P = 56 atm e temperatura de T = 450 K usando a equação de van der
Waals. Os cálculos deveriam ser repetidos para as pressões reduzidas de Pr = 1, 2, 4,10 e 20, e
depois o fator de compressibilidade deveria ser analisado em função de Pr.
Dever-se-ia resolver equações não lineares de terceiro grau para a variável de volume molar.
A equação era:
PV³ - (RT + Pb) V²+aV-ab=0
O método numérico utilizado foi o de bissecção, explicado na seção materiais e métodos.
3.2 Problema 2, Balanço de materiais em um trem de separação, introdução à
Engenharia Química(Balanço de Massa), solução simultânea de equações lineares.
No problema um sistema com Xileno, Estireno, Tolueno e Benzeno é separado por um
sistema de colunas destiladoras. O esquema abaixo mostra as colunas, as correntes F, D, B,
D1, B1, D2 e B2 com suas composições molares e a vazão molar.
8
Figura 3.2-1 Fonte: A collection of ten numerical problems in chemical engineering solved by various
mathematical software packages.
O conceito de balanço de massa parte do princípio de que a massa é conservada.
Portanto, em qualquer sistema tem-se que para uma determinada substância:
Sai = Entra + Reage – acumula
O termo reage pode ser positivo ou negativo dependendo se dentro do sistema a
substância é produzida ou consumida.
Fazendo o balanço de massa para os quatro componentes para todo o sistema, como
não há reação (há a conservação da quantidade de matéria) e nem acúmulo para eles, tem-se
que Entra = Sai. Então:
Xileno: 0.07D1 + 0.18B1 + 0.15D2 + 0.24B2 = 0.15*70
Estireno: 0.04D1 + 0.24B1 + 0.10D2 + 0.65B2 = 0.25*70
9
Tolueno: 0.54D1 + 0.42B1 + 0.54D2 + 0.10B2 = 0.40*70
Benzeno: 0.35D1 + 0.16B1 + 0.21D2 + 001 = 0.20*70
Fazendo o balanço de massa para a coluna de destilação dois, onde é determinada a
vazão e a fração molar da corrente D, tem-se:
Vazão molar total: D = D1 + B1
Xileno: XDX*D = 0.07D1 + 0.18B1
Estireno: XDS*D = 0.04D1 + 0.24B1
Tolueno: XDT*D = 0.54D1 + 0.42B1
Benzeno: XDB*D = 0.35D1 + 0.16B1
Onde XDX = Fração molar de Xileno em D, XDS = Fração molar de Estireno em D,
XDT = Fração molar de Tolueno em D, XDB = Fração molar de Benzeno em D.
Analogamente, fazendo o balanço de massa na coluna três, onde é determinada a
vazão e a fração molar da corrente B, tem-se:
Vazão molar total: B = D2 + B2
Xileno: XBX*B = 0.15D2 + 0.24B2
Estireno: XBS*B = 0.10D2 + 0.65B2
Tolueno: XBT*B = 0.54D2 + 0.10B2
Benzeno: XBB*B = 0.21D2 + 0.01B2
A partir das quatro primeiras inequações foi possível determinar D1, D2, B1, B2. O
método utilizado para a solução foi a divisão a esquerda, o que é explicada na seção materiais
e métodos, trabalhando com a matriz do sistema linear. A vazão e a composição molar das
correntes B e D foram determinadas com as outras oito equações, apenas substituindo os
valores nelas, no próprio Scilab.
10
3.3 Problema 4, Equilíbrio de múltiplas reações entre gases, Termodinâmica, solução
simultânea de equações não lineares.
Neste problema tem-se 7 gases em um reator de volume constante em que ocorrem 3
reações de equilíbrio:
A + B <-> C + D
B + C <-> X + Y
A + X <-> Z
Um sistema algébrico descreve a situação de equilíbrio das reações. Utilizando as
relações termodinâmicas de equilíbrio e as relações estequiométricas para cada reação,
teremos um conjunto de sete equações:
KC1 = CC*CD/CA*CB KC2 = CX*CY/CB*CC KC3 = CZ/CA*CX
CA = CA0 – CD – CZ CB = CB0 – CD – CY CC = CD – CY
CY = CX + CZ
Nas equações CA, CB, CC, CD, CX, CY e CZ são as concentrações das diversas
espécies envolvidas nas reações.
No problema teve-se que solucionar o sistema de equações quando CA0 = CB0 = 1.5,
KC1 = 1.06, KC2 = 2.63 e KC3 = 5 para as estimativas iniciais de que CD = CX = CZ = 0.
Neste problema utilizou-se o solucionador integrado do Scilab, fsolve.
3.4 Problema 5, velocidade terminal de queda, Dinâmica de fluidos, solução de uma
equação não linear.
O conceito envolvido neste problema é o cálculo da velocidade terminar de partículas
sólidas em queda em fluidos sob a força da gravidade.
Uma partícula acelerará sob a ação da gravidade num fluido até que a força de arrasto
iguale-se à força da gravidade, momento em que ela começa a cair com uma velocidade
constante, terminal, dada por:
11
ut = sqrt(2gmp(Ρp – P)/P*PpApCd)
Onde
g = aceleração da gravidade (9.80665 m/s²);
mp = massa da partícula;
Pp = densidade da partícula(kg/m³);
P = densidade do fluido (kg/m³);
Ap = Área da partícula projetada na direção do movimento (m²);
Cd = coeficiente de arrasto.
Para partículas esféricas homogêneas, que é o caso do problema, a equação fica:
vt = sqrt(4gDp(Ρp – P)/P*Cd)
Se existe um movimento relativo entre uma partícula e o meio que ela está imersa, o
fluido do meio exercerá uma força de arrasto sobre a partícula. O valor do coeficiente de
arrasto é uma função do número de Reynolds. O número de Reynolds é dado por:
Re = DpvtP/u
Onde
u é a viscosidade em Pa*s ou Kg/m*s
Para número de Reynolds pequenos (Re < 0.1), pela lei de Stokes tem-se:
Cd = 24/Re
Para número de Reynolds intermediário (0.1 <= Re <= 1000), tem-se:
Cd = (24/Re)(1 + 0.14Re^0.7)
Para o número de Reynolds (1000 <= Re <= 350000) no intervalo em que se aplica a lei de
Newton, ação e reação, tem-se:
Cd = 0.44
12
Para números de Reynolds muito altos (Re > 350000), onde o coeficiente de arrasto cai
drasticamente, tem-se:
Cd = 0.19 – 8*10^4/Re
No problema temos de calcular a velocidade terminal de partículas de carvão (Pp =
1800 kg/m³ e Dp = 0.208*10^-3 m) afundando na água à temperatura de T = 298.15 K( p =
994.6 kg/m³ e u = 8.931*10^-4kg/ms) sob a ação da gravidade e depois numa centrífuga onde
a aceleração é 30 vezes a da gravidade.
O problema foi resolvido utilizando o solucionador fsolve.
3.5 Problema 6, Troca de calor em uma série de tanques, Transferência de Calor,
solução simultânea de um sistema de equações diferenciais ordinárias.
No problema uma solução de vários componentes era aquecida numa série de 3
tanques antes de ser alimentada numa coluna de destilação. Os tanques eram bem misturados,
admitindo então uma uniformidade de temperatura do óleo dentro de cada tanque e estavam
isolados do meio exterior. O óleo entrava a uma temperatura de T0= 20ºC no primeiro tanque
a uma vazão mássica de W1 = 100 kg/min, vapor saturado a Tsteam = 250ºC condensava nas
serpentinas aquecedoras dentro do tanque. Saindo do primeiro tanque a uma temperatura T1
com a mesma vazão mássica que entrava no tanque 1 entrava num tanque 2, o qual era
também aquecido por uma serpentina onde ocorria a condensação do vapor a 250ºC.
Analogamente, o óleo passava então por um 3 tanque.
O aquecimento do óleo era causado pelo resfriamento do vapor. A taxa de
transferência de calor é dada pela expressão:
Q = UA(Tsteam – T)
Onde
UA = Coeficiente de transferência de calor * Área da serpentina (UA = 10 KJ/min.ºC);
T = Temperatura do óleo no tanque (ºC);
Q = Taxa de transferência de calor (kJ/min).
13
Como a energia era conservada no sistema e energia não é criada, nem destruída,
pode-se realizar o balanço de energia para cada tanque individualmente utilizando o seguinte
balanço:
Energia Acumulada = Energia que entra – Energia que sai
A densidade do óleo é assumida constante, bem como a vazão mássica em todos os
tanques, ou seja, W = W1 = W2 = W3 e M = M1 = M2 = M3.
Os balanços de energia dos tanques ficam:
Tanque 1
dÛ1/dt = Û01 + Q – Ûf1
M*Cp*(dT1/dt) = W*Cp*T0 + UA(Tsteam – T1) – W*Cp*T1
Onde
Cp = Capacidade calorífica do óleo (Cp = 2.0 kJ/kg).
Analogamente, para os outros tanques:
Tanque 2:
M*Cp*(dT2/dt) = W*Cp*T1 + UA(Tsteam – T2) – W*Cp*T2
Tanque 3:
M*Cp*(dT3/dt) = W*Cp*T2 + UA(Tsteam – T3) – W*Cp*T3
Rearranjando as EDOs temos 3 diferenciais que podem ser resolvidas por uma solução
numérica. O problema nos pede as temperaturas de quase equilíbrio nos 3 tanques, fazendo o
tempo decorrido ser bem grande, e o tempo necessário para que a temperatura T3 alcance o
valor de 99% do equilíbrio. O sistema de EDOs foi resolvido utilizando o solucionador ode
do Scilab.
3.6 Problema 7, Difusão com uma reação química em um disco, Fenômeno de
transporte, solução de uma equação diferencial ordinária de segunda ordem com 2
condições de contorno.
14
No problema tem-se uma reação química entre A e B e a difusão dos componentes. A
é o reagente, B é o produto. O reagente A é absorvido pelo conteúdo e é consumido dentro
dele. A equação que modela esse sistema é uma EDO de segunda ordem:
d²Ca/dz² = (k/DAB)*Ca
Onde
Ca = concentração do reagente A (kgmol/m³);
z = variável de distância (m);
k = constante da reação homogênea de consumo de A (1/s);
DAB = coeficiente de difusão em área (m²/s).
A interpretação geométrica do sistema é uma superfície exposta que não permite
difusão pela sua superfície externa. Então se tem as condições de contorno:
Ca = Ca0 z = 0
dCa/dz = 0 z = L
Onde Ca0 é a concentração de Ca na superfície z = 0 e como não há difusão na
superfície externa, a derivada de Ca em relação a z é 0.
O sistema possui uma solução analítica da forma:
Ca =Ca0*(cosh(L*sqrt(k/DAB)*(1-z/L)))/cosh(L*sqrt(k/DAB))
Para a solução numérica, utilizou-se o solucionador ode implementado ao método de
Newton-Raphson, para que a condição de contorno da derivada zero estivesse satisfeita.
3.7 Problema 8, destilação de uma mistura binária, Processo de separação, solução de
um sistema de uma equação diferencial ordinária e uma equação não linear.
Destilação é um processo de separação de componentes de uma mistura líquido/vapor
que coexistem a uma mesma temperatura e pressão. Componentes mais voláteis (temperatura
de ebulição mais baixa) tendem a permanecer no estado de vapor; e os menos voláteis
(temperatura de ebulição mais alta) tendem a permanecer no estado líquido. O resultado é que
com o aquecimento da mistura, a fase de vapor fica mais ricamente concentrada com os
15
componentes mais voláteis; e, a fase líquida, com os componentes menos voláteis. Na
destilação deste problema, a mistura é alimentada e o processo ocorre em batelada, ou seja,
não há uma alimentação contínua. O líquido remanescente L vai diminuindo e a temperatura,
por estar o líquido mais rico com o componente de maior temperatura de ebulição, vai
aumentando. A equação que descreve o valor de L em função da fração molar do componente
menos volátil, x2, é uma diferencial da forma:
dL/dx2 = L/(x2*(k2 – 1))
Onde
k2 = Razão do equilíbrio vapor/líquido.
O valor de ki é uma relação que quantifica a tendência do componente i em se vaporizar. O
seu valor é dado por:
ki = Pi/P
Onde
Pi = Pressão de vapor do componente;
P = Pressão total do sistema.
A pressão de vapor Pi pode ser modelada pela equação de Antoine que utiliza 3 parâmetros
experimentais, A, B, C para um componente i tem-se:
Pi = 10^(Ai – Bi/(T+Ci))
Onde
T = Temperatura do sistema em ºC.
Os valores de A, B e C para o benzeno e para o tolueno são:
Benzeno: A = 6.90565 B = 1211.033 C = 220.79
Tolueno: A = 6.95464 B = 1344.8 C = 219.482
A temperatura do sistema de destilação segue o ponto de bolha da mistura. Quando se
tem uma mistura binária, o ponto de bolha é a temperatura em que a primeira bolha de vapor
16
aparece. O ponto de bolha é definido implicitamente usando as razões do equilíbrio
vapor/líquido das 2 substâncias. A equação algébrica é:
k1*x1 + k2*x2 = 1
No problema uma mistura de benzeno (componente 1) e tolueno (componente 2) eram
destiladas e após a fração de tolueno ter um valor determinado, pedia–se a quantidade de
líquido remanescente. A solução deste problema exigiu a solução de uma sistema com 1
equação algébrica não linear, solucionada por fsolve, e 1 EDO, solucionada por ode.
4 MÉTODOS NUMÉRICOS
4.1 Sobre os Métodos Numéricos
Métodos numéricos são formas de solucionar problemas utilizando os conhecimentos
de Cálculo e funções. Eles diferem da solução analítica por dar uma solução aproximada do
problema, muitas vezes implementá-los é muito mais prático, simples e rápido para solucionar
diversos problemas que envolvem funções. Métodos numéricos são baseados em
aproximações, realizando diversas iterações (repetição do algoritmo do método), por serem
estruturados em algoritmos que devem ser repetidos diversas vezes a sua praticidade se
encontra quando entra em conjunto com alguma linguagem de programação que admite laço,
como por exemplo, o software utilizado neste projeto. Métodos numéricos também podem ser
utilizados para otimização, controle de sistemas, construção de gráficos. Alguns dos
principais métodos de resolução serão apresentados a seguir, todos retirados de livros de
Cálculo Numérico, ver nas referências bibliográficas.
4.2 Método da bissecção
Este método é utilizado no problema 1 para encontrar a raiz de uma função(f(x)=0) e
utiliza o Teorema do Valor Intermediário do cálculo.
17
O Teorema do Valor Intermediário do cálculo diz que se temos uma função f(x)
continua num intervalo [a,b] e f(a) >0 e f(b)<0, temos uma f(c) = 0, em que c pertence ao
intervalo [a,b]. Então a função tem uma raíz real c.
No método de bissecção diversas iterações são feitas para que se encontre a raiz da
função. A figura a seguir ilustra o método e seu algoritmo, representando 3 iterações que
aproximam-se cada vez mais do valor ξ, a raiz da função.
Figura 3.1 Fonte: Cálculo Numérico, Aspectos Teóricos e Computacionais - 2 Edição
Algoritmo do método:
Tendo a0 e b0
f(a0)<0 e f(b0)>0 temos um x0 = (a0 + b0)/2
Se f(x0)<0 então fazemos, a1= x0 e b1=b0. Se f(x0)>0 então fazemos, b1=x0 e a1=a0.
A iteração é repetida: x1 = (a0+b0)/2... Até que xi convirja para a raiz com um erro e.
18
4.3 Método Newton-Raphson
Este método é utilizado também para encontrar a raiz de uma função, ele utiliza a
noção da reta tangente a um ponto de f(x). A função iterada é :
Xk+1 = Xk – f(Xk)/f’(Xk)
No processo de iteração, o valor Xk+1 é onde a reta tangente a f(Xk) intercepta a abcissa.
Conforme iterada, a função convergiria ao valor do zero da função original. A imagem
a seguir ilustra o método descrito.
Figura 4.2 Fonte: http://pt.wikipedia.org/wiki/Ficheiro:Ganzhi001.jpg
4.4 Divisão pela esquerda
A divisão pela esquerda consiste em multiplicar os dois lados da representação
matricial pela inversa da matriz dos coeficientes numéricos. No SciLab a divisão pela
esquerda é representada por \. Quando multiplicamos uma matriz pela sua inversa temos sua
identidade, portanto, na solução pela divisão a esquerda, as variáveis são retornadas isoladas.
4.5 Método do Gauss-Jacobi
19
Este método iterativo é utilizado para solução de sistema de n equações lineares.
Normalmente é utilizado quando o número de equações é grande, em que este método
encontra utilidade.
Dada uma matriz Ax = b do sistema linear de n equações e x = Cx + g. Tendo um sistema:
Figura 4.5-1 Fonte: Cálculo Numérico, Aspectos Teóricos e Computacionais - 2 Edição
Desta forma tem-se x = Cx + g onde:
Figura 4.5-2 Fonte: Cálculo Numérico, Aspectos Teóricos e Computacionais - 2 Edição
E
20
Figura 4.5-3 Fonte: Cálculo Numérico, Aspectos Teóricos e Computacionais - 2 Edição
O método consiste em a partir do valor de x(0) obter uma sequência x(1), ..., x(k) que
convirja para a solução depois de k iterações da relação:
x(a+1) = Cx(a) + g
4.6 Método de Euler
O método de Euler pode ser utilizado para resolver um problema de valor inicial: y’ =
f(x,y), y(x0) = y0. Como se sabe o valor de x0 e y0 tem-se o valor de y’(x0,y0) = f(x0, y0).
Temos a reta tangente ao ponto (x0, y0) da função diferencial conhecida e sua equação é:
R: r0(x) = y(x0) + (x – x0)*y’(x0)
Escolhendo um h = x(k+1) – x(k), y(x1) é aproximadamente y1 = r0(x1), ou seja,
y1 = y0 + h*f(x0,y0)
O raciocínio pode ser repetido para y(x2) e, assim, a relação iterativa para o método fica:
y(k+1) = y(k) + h*f(x(k),y(k))
Interpretando geometricamente o problema:
21
Figura 4.6-1 Fonte: Cálculo Numérico, Aspectos Teóricos e Computacionais - 2 Edição
Quanto menor o intervalo entre x0 e x1 tem-se maior precisão no valor de y, este raciocínio é
usado para métodos mais elaborados que utilizam o mesmo princípio do método de Euler.
4.7 Solucionadores
O Scilab possui solucionadores incorporados na sua programação. Neste projeto, 2
deles foram amplamente utilizados. O solucionador fsolve é utilizado para se encontrar a raíz
de uma função polinomial. O solucionador ode é utilizado para modelagem de uma equação
diferencial ordinária de primeira ordem.
5 RESULTADOS
5.1 Problema 1
22
A-) O volume molar da amônia é de 0.574893 litros/mol e sua compressibilidade é de
0.871828 à pressão de 56 e temperatura 450.
Z= 0.871828 Pr = 0.503145
Se resolvido como um gás ideal (PV = nRT) o valor do volume molar da amônia é 0.6594107
L/mol, um valor maior de volume, com significada diferença. O valor de volume molar
calculado pela relação do gás ideal se aproxima cada vez mais do valor calculado pela
equação de van der Waals quanto mais a pressão se aproxima de 1 atmosfera.
B-)
PR V(litros/mol) Z
1 0.233506 0.7038
2 0.0772709 0.465797
4 0.0606556 0.731277
10 0.0508755 1.53342
20 0.0461735 2.78339
C-)
O fator de compressibilidade apresenta este comportamento de decrescer da pressão
reduzida 1 para 2 por, em pressões que ultrapassam a pressão crítica, a fase líquido/ vapor não
ser mais distinguível. Esta ultrapassagem do ponto crítico do gás é muito bem percebida pelo
volume molar, que diminui bruscamente da pressão reduzida 1 para a pressão reduzida 2.
0
0,5
1
1,5
2
2,5
3
0 5 10 15 20 25
Z x PR
Z x PR
23
5.2 Problema 2
A-) Os valores das taxas molares em cada saída são:
D1=26.25 mol/s
B1=17.5 mol/s
D2=8.75 mol/s
B2=17.5 mol/s
B-) Os valores das taxas molares em cada saída são:
D=43.75 mol/s, em que consistia de 4.9875 mols de xileno, 5.25 mol de estireno, 21.525 mol
de tolueno, 11.9875 mol de benzeno.
B=26.25 mol/s, em que consistia de 5.5125 mols de xileno, 12.25 mol de estireno, 6.475 mol
de tolueno, 2.0125 mol de benzeno.
5.3 Problema 4
A-) A concentração final de cada espécie é :
A = 0.420689
B = 0.242897
C = 0.153565
D = 0.705334
X = 0.177792
Y = 0.551769
Z = 0.373977.
5.4 Problema 5
24
A-) A velocidade terminal para a partícula acelerada pela gravidade é 0.0157816 m/s
B-) A velocidade terminal para a partícula acelerada por 30 vezes a gravidade é 0.206021 m/s
5.5 Problema 6
O gráfico foi desenhado no Scilab e representa a temperatura do óleo no tanque
1(curva azul), tanque 2 (curva verde), tanque 3(curva vermelha) em função do tempo que
varia de 0 a 160 s. A temperatura de equilíbrio nos tanques 1, 2 e 3, respectivamente foram:
30.9524 °C, 41.3832°C, 51.3173°C.
25
5.6 Problema 7
No gráfico do valor da concentração de A em função do valor de z tem-se a solução
analítica (curva verde) e a solução pelo método numérico (curva azul). Nota-se que os dois
métodos possuem soluções muito semelhantes.
26
5.7 Problema 8
A curva acima descreve a quantidade de líquido remanescente na destilação(L) em
função da fração molar do tolueno na mistura (0.8 ao final). A quantidade de líquido
remanescente foi de 13.755 mols.
6 CONCLUSÃO
Modelar e solucionar problemas de Engenharia através de recursos computacionais é
uma forma mais simples e prática considerando cada tipo de problema. Cada solução possui
erros inerentes ao processo de iteração, o qual nunca é exato. Na modelagem dos problemas
erros são incorporados, portanto, no momento de programar um método numérico deve-se
analisar o tempo que se dispões para sua resolução, o erro máximo que é admitido, o poder de
processamento da máquina onde se escreve o código e o tipo de problema. No projeto
aplicaram-se diversos métodos para cada problema, não se preocupou com o tempo para
completar as iterações, pois os problemas eram resolvidos isoladamente e possuíam estrutura
muito simples. O treinamento foi satisfatório, a familiaridade com o processo de programação
foi alcançada e os métodos foram absorvidos em parte. Novas fronteiras para o aprendizado
de métodos numéricos foram abertas. No quinto semestre do curso o aprendizado dos
métodos será aprofundado no curso de Cálculo numérico.
27
7 TRABALHOS FUTUROS
Três problemas da coletânea ficaram a serem resolvidos, por questão de tempo de
projeto e eventualidades ao decorrer do ano de pesquisa. Os outros sete problemas podem ser
analisados quanto a outras formas de solucioná-los, há inúmeras com a utilização de métodos
numéricos. O treinamento deu uma base de conhecimento para o aprendizado de otimização
de processos químicos. A otimização possui diversas técnicas para busca de melhores
soluções para os problemas (encontrar máximos ou mínimos).
8 APÊNDICE
8.1 Linhas de código de programação dos problemas:
Problema 1
1. // Calculo do volume molar do gás amônia e o seu fator de compressibilidade com o
uso da equação de Van der Waals
2. // (P+a/V²)*(V-b)=R*T
3. Tc = 405.5; // Temperatura critica da amônia em Kelvin
4. Pc = 111.3; // Pressão critica da amônia em atm
5. R = 0.08206; // Constante dos gases, atm.L/g.mol.K
6. a= (27/64)*((R*R*Tc*Tc)/Pc);
7. b= (R*Tc)/(8*Pc);
8. printf('O Problema 1 será resolvido para determinar o volume molar e a
compressibilidade de uma amostra do gás de Amônia');
9. T = 450;
10. // K
11. Pp= 56;
12. //Atm
13. e= input('Qual a precisão desejada (e<1 e e>0)? ');
14. // Pressão a ser usada nos casos
15. Pr = [(Pp/Pc) 1 2 4 10 20];
16. c = 1;
17. while c<=6
18. P = Pc * Pr(c);
19. W= P;
20. A = -(R*T + P*b);
21. B = (a);
22. C = -(a*b);
23. // Equação do Volume molar: W*V³+A*V²+B*V+C=0 (1)
24. //Método da bisseção para encontrar a raiz
25. // x=0 satisfaz que a equação 1 seja negativa, -(a*b) é negativo já que a e b são
positivos
26. Ri=0;
27. //Achar x que a equação 1 fique >0
28. //Chute inicial:
29. X=((R*T)/P)*3;
28
30. //teste:
31. if W*X*X*X+A*X*X+B*X+C >0 then
32. Si= X;
33. else printf('Defina o segundo valor por outro método')
34. end
35. //precisão desejada
36. n= ((log(Si-Ri)-log(e))/log(2))-1;
37. //método bisseção
38. i=1;
39. while i<= n
40. xi=(Ri+Si)/2;
41. if W*xi*xi*xi+A*xi*xi+B*xi+C<0 then
42. Ri=xi;
43. Si = Si;
44. else
45. Ri=Ri;
46. Si = xi;
47. end
48. i= i + 1;
49. end
50. Z = (P*xi)/(R*T);
51. printf('O volume molar da amônia é de %g litros/mol e sua compressibilidade é de %g
à pressão de %g e temperatura %g,com precisão de %g\n', xi, Z , P, T , e);
52. //Estudo de Z e Pr
53. printf('Z= %g Pr = %g\n', Z , Pr(c));
54. c = c + 1;
55. end
Problema 2
1. //Problema 2 - Steady state material balances on separation train
2. // Solução de um sistema de equações lineares.
3. //Balanço do material nos componentes individuais no trem de separação, na forma
matricial;
4. A = [0.07 0.18 0.15 0.24; 0.04 0.24 0.10 0.65; 0.54 0.42 0.54 0.10; 0.35 0.16 0.21
0.01];
5. b = [0.15*70; 0.25*70; 0.40*70; 0.20*70];
6. //Parte a do problema
7. //temos que : [xylene; styrene; toluene; benzene] = A*[D1 ; B1; D2; B2]= b
8. //divisão à esquerda(Multiplicação de b pela inversa de A)
9. a = A\b;
10. //atribuição dos valores
11. D1= a(1,1);
12. B1= a(2,1);
13. D2= a(3,1);
14. B2= a(4,1);
15. printf('Os valores das taxas molares em cada saída são:\n D1=%g mol/s\n B1=%g
mol/s\n D2=%g mol/s\n B2=%g mol/s\n', D1, B1, D2, B2);
29
16. //Parte b
17. //Saída D; na sequencia de equações temos: Fração molar de xylene, styrene, toluene e
benzene
18. XdxD = [0.07 0.18]*[D1; B1];
19. XdsD = [0.04 0.24]*[D1; B1];
20. XdtD = [0.54 0.42]*[D1; B1];
21. XdbD = [0.35 0.16]*[D1; B1];
22. D = XdxD+ XdsD+ XdtD + XdbD;
23. //Saída B; mesma sequencia da anterior
24. XbxB = [0.15 0.24]*[D2; B2];
25. XbsB = [0.10 0.65]*[D2; B2];
26. XbtB = [0.54 0.10]*[D2; B2];
27. XbbB = [0.21 0.01]*[D2; B2];
28. B = XbxB+ XbsB+ XbtB + XbbB;
29. printf('Os valores das taxas molares em cada saída são:\n D=%g mol/s, em que
consistia de %g mol de xylene, %g mol de styrene, %g mol de toluene, %g mol de
benzene\n B=%g mol/s, em que consistia de %g mol de xylene, %g mol de styrene,
%g mol de toluene, %g mol de benzene', D, XdxD, XdsD, XdtD, XdbD, B, XbxB,
XbsB, XbtB, XbbB)
Problema 4
1. clear
2. k1=1.06;
3. k2=2.63;
4. k3=5;
5. function f=eq(x)
6. f(1)=x(3)*x(4)-k1*x(1)*x(2);
7. f(2)=x(5)*x(6)-k2*x(2)*x(3);
8. f(3)=x(7)-k3*x(1)*x(5);
9. f(4)=1.5-x(4)-x(7)-x(1);
10. f(5)=1.5-x(4)-x(6)-x(2);
11. f(6)=x(4)-x(6)-x(3);
12. f(7)=x(5)+x(7)-x(6);
13. endfunction
14. xi=[1.5 1.5 0 0 0 0 0];
15. x=fsolve(xi,eq,%eps);
16. printf('Para o estado a as concentração finais são:\n A:%g\n B::%g\n C::%g\n D::%g\n
X::%g\n Y::%g\n Z::%g\n',x(1),x(2),x(3),x(4),x(5),x(6),x(7))
17. function f=eq(x)
18. f(1)=x(3)*(x(4)+1)-k1*x(1)*x(2);
19. f(2)=(x(5)+1)*x(6)-k2*x(2)*x(3);
20. f(3)=(x(7)+1)-k3*x(1)*(x(5)+1);
21. f(4)=1.5-(x(4)+1)-(x(7)+1)-x(1);
22. f(5)=1.5-(x(4)+1)-x(6)-x(2);
23. f(6)=(x(4)+1)-x(6)-x(3);
24. f(7)=(x(5)+1)+(x(7)+1)-x(6);
25. endfunction
30
26. xi=[1.5 1.5 0 1 1 0 1];
27. x=fsolve(xi,eq,%eps);
28. printf('Para o estado B, as concentração finais são:\n A:%g\n B::%g\n C::%g\n
D::%g\n X::%g\n Y::%g\n Z::%g\n',x(1),x(2),x(3),x(4),x(5),x(6),x(7))
29. xi=[3 3 5 15 15 5 15];
30. x=fsolve(xi,eq,%eps);
31. printf('Para o estado c as concentração finais são:\n A:%g\n B::%g\n C::%g\n D::%g\n
X::%g\n Y::%g\n Z::%g\n',x(1),x(2),x(3),x(4),x(5),x(6),x(7))
Problema 5
1. //Resolução do problema 5 a partir do fsolve do scilab
2. Pp=1800;
3. Dp= 0.000208;
4. T= 298.15;
5. p=994.6;
6. u=8.931*(10^(-4));
7. g=9.80665;
8. function F=f(vt)
9. Re=(Dp*vt*p)/u;
10. if Re<0.1 then
11. Cd=(24/Re);
12. elseif Re<=1000 then
13. Cd= (24/Re)*(1+0.14*(Re^(0.7)))
14. elseif Re<=350000 then
15. Cd=0.44
16. else
17. Cd= 0.19-(8*(10^4))/Re;
18. end
19. F=sqrt((4*g*(Pp-p)*Dp)/(3*Cd*p))- vt
20. endfunction
21. [vt]=fsolve(0.1,f,%eps)
22. printf('A velocidade terminal para a partícula acelerada pela gravidade é %g m/s \n',
vt)
23. g=30*g;
24. function F=f(vt)
25. Re=(Dp*vt*p)/u;
26. if Re<0.1 then
27. Cd=(24/Re);
28. elseif Re<=1000 then
29. Cd= (24/Re)*(1+0.14*(Re^(0.7)))
30. elseif Re<=350000 then
31. Cd=0.44
32. else
33. Cd= 0.19 - (8*(10^4))/Re;
34. end
35. F=sqrt((4*g*(Pp-p)*Dp)/(3*Cd*p))-vt
36. endfunction
31
37. [vt]=fsolve(0.1,f,%eps)
38. printf('A velocidade terminal para a partícula acelerada por 30 vezes a gravidade é %g
m/s', vt)
Problema 6
1. //// Problema 6, differential equations...
2. //Constantes do processo
3. //Heat capacity of the oil KJ/kg
4. Cp=2;
5. //KJ/min*ºC
6. UA= 10;
7. //Temperatura inicial ºC
8. To=20;
9. //Vazão mássica kg/min
10. W=100;
11. //Massa de óleo kg
12. M= 1000;
13. //Temperatura do vapor ºC
14. Ts=250;
15. function tdot=f(t, T)
16. tdot(1)=[W*Cp*(To-T(1))+UA*(Ts-T(1))]/(M*Cp)
17. tdot(2)=[W*Cp*(T(1)-T(2))+UA*(Ts-T(2))]/(M*Cp)
18. tdot(3)=[W*Cp*(T(2)-T(3))+UA*(Ts-T(3))]/(M*Cp)
19. endfunction
20. c=0:1:5000;
21. for i = 1:101
22. t(1)= 0
23. t(i)= (i-1)
24. T0=[20 20 20];
25. t0=0;
26. [X]=ode(T0,t0,t,f);
27. T1(1)=20
28. T2(1)=20
29. T3(1)=20
30. a = 1 + 3*t(i);
31. b = 2 + 3*t(i);
32. c = 3 + 3*t(i);
33. T1(i)= X(a);
34. T2(i)= X(b);
35. T3(i)= X(c);
36. end
37. printf('As temperaturas de equilibrio nos tanques 1, 2 e 3, respectivamente foram:%g,
%g, %g', T1(161), T2(161), T3(161))
38. plot(t, [T1 T2 T3])
32
Problema 7
////Solução do problema 7 da solução analítica.
////Concentração inicial:(kgmol/m³)
CA0=0.2;
////Constante s^-1
k=10^(-3);
////coeficiente de difusão binário? m²/s
DAB= 1.2*(10^(-9));
//// m, distância da superfície
L=10^(-3);
for c= 2:1:100
z(1)=0
z(c) = c*0.00001
CAb=CA0*(cosh(L*sqrt(k/DAB)*(1-z/L)))/cosh(L*sqrt(k/DAB));
end
////Solução do problema pela EDO de segunda grau.
////x(1)= Ca, x(2) = dca/dt
p=1
dCadt(100) = 1;
while abs(dCadt(100)) > 10^(-6)
for i = 2:1:100
z(1)=0
z(i) = i*0.00001
function ydot=f(z, x)
ydot(2)= k*x(1)/DAB;
ydot(1)= x(2);
ydot(3) = x(4);
ydot(4) = k*x(3)/DAB;
endfunction
A(1) = 1;
A0 = A(p)
x01=[CA0 A0 0 1];
z0=0;
[X] = ode(x01,z0,z,f)
CA(1) = 0.2;
dCadt(1) = A(p);
a = 1 + 4*(i-2);
b = 2 + 4*(i-2);
d = 4 + 4*(i-2);
dCadt(i)= X(b);
CA(i) = X(a);
DX(b) = X(d);
end
p= p + 1
A(p)= A(p-1) - ((X(b))/(X(d)))
end
plot(z, [CA CAb])
33
Problema 8
////Problema 8 - Binary Batch Destillation
////A mistura a ser destilada é uma mistura binária de benzeno e tolueno.
////Dados:
//constantes do benzeno
A1=6.90565;
B1=1211.033;
C1=220.79;
//constantes do tolueno
A2=6.95464;
B2=1344.8;
C2=219.482;
//Informações do sistema
P= 1.2*760;
//Pressão mm de HG
L0= 100;
//100 mols da mistura inicial
x10=0.60;
//fração molar de benzeno
x20=0.40;
//fração molar de tolueno
for i = 2:1:101
X(1)=0.400
X(i) = 0.400 + 0.004*(i-1)
function F=temp(x)
F= (((10^(A1-(B1/(x+C1))))/P)*(1-X(i))) + ((10^((A2-(B2/(x+C2))))/P)*(X(i))) - 1
endfunction
////Função não linear(Resolvendo-a teremos a temperatura x do sistema)
x(i)=fsolve(80,temp)
k2(i) = (10^((A2 - (B2/(x(i)+C2)))))/P
//Solução pelo método de Euler
h(i) = X(i)- X(i-1)
L(1) = 100
L(i) = L(i-1) + h(i)*L(i-1)/(X(i)*(k2(i)-1))
end
plot(X, L)
9 REFERÊNCIA BLIBIOGRÁFICA
CUTLIP, M. HWALEK, J. et al. A collection of ten numerical problems in chemical
engineering solved by various mathematical software packages. New York: John Wiley &
Sons, 1998.
34
BADINO Jr, Alberto Colli; CRUZ, Antonio José Gonçalves. Fundamentos de Balanços de
Massa e Energia: um texto básico para análise de processos químicos. São Carlos:
EdUFSCar, 2010.
ARENALES, Selma Helena Vasconcelos; DAREZZO, Arthur. Cálculo numérico:
aprendizagem com apoio de software. São Paulo: Thomson Learning, 2008.
RUGGIERO, Márcia A. Gomes; ROCHA, Vera Lúcia da Lopes. Cálculo Numérico,
Aspectos Teóricos e Computacionais. São Paulo: Pearson Education, 1996.
Perry, Robert, et al. Perry's Chemical Engineers' Handbook. New York: McGraw-Hill
Professional, 2007.
Dean, Jean A. . Lange's Handbook of Chemistry. New York: McGraw-Hill
Professional,1999.
Van Ness, H. C., et al. Introduction to chemical engineering thermodynamics. New York:
McGraw-Hill Professional, 2005.
Lacerda, E. G. M. Programando com Scilab: versão: 0.24. UFRN: Departamento de
Engenharia de Computação e Automação (DCA), 2011.