ATIVIDADE COMPLETA DE SIMULAÇÃO EQUAÇÕES DIFERENCIAIS

Prévia do material em texto

Análise e Simulação de Processos:
Relatório da Atividade II
Florianópolis
Sumário
1 Problema 1 5
1.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Problema 2 9
2.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3 Problema 3 13
3.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2.1 Item a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2.2 Item b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3.1 Item a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3.2 Item b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4 Problema 4 17
4.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
5 Problema 5 20
5.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
5.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
5.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
6 Problema 6 23
6.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
6.2 Resolução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
6.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
7 Problema 7 25
7.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
7.2 Resolução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
7.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
8 Problema 8 28
8.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
8.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
8.2.1 Modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
8.2.2 Resolução - Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
8.2.3 Resolução - Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . 30
8.2.4 Resolução - Solver . . . . . . . . . . . . . . . . . . . . . . . . . . 31
1
8.3 Item c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
8.3.1 Resolução - Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
8.3.2 Resolução - Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . 33
8.3.3 Resolução - Solver . . . . . . . . . . . . . . . . . . . . . . . . . . 34
8.4 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
8.4.1 Item b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
8.4.2 Item c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
9 Problema 9 41
9.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
9.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
9.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
10 Problema 10 43
10.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
10.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
10.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
11 Problema 11 46
11.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
11.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
11.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
12 Problema 12 49
12.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
12.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
12.2.1 Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
12.2.2 Método de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . 50
12.2.3 Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
12.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
13 Problema 13 53
13.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
13.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
13.2.1 Efeito de cada termo na equação . . . . . . . . . . . . . . . . . . 55
13.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
13.3.1 Problema principal . . . . . . . . . . . . . . . . . . . . . . . . . . 56
13.3.2 Fast Fourier Transform (FFT) . . . . . . . . . . . . . . . . . . . . 58
14 Problema 14 60
14.1 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
14.2 Item a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
14.2.1 Solução estacionária . . . . . . . . . . . . . . . . . . . . . . . . . 61
14.2.2 Solução transiente . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
14.2.3 Item d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
14.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
14.3.1 Solução transiente - Sistema I . . . . . . . . . . . . . . . . . . . . 66
14.3.2 Solução transiente - Sistema I . . . . . . . . . . . . . . . . . . . . 67
14.3.3 Solução estacionária - Sistema I . . . . . . . . . . . . . . . . . . . 68
2
14.3.4 Solução estacionária - Sistema II . . . . . . . . . . . . . . . . . . 70
15 Problema 15 72
15.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
15.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
15.2.1 Quantidade inicial: (3 ; 1) . . . . . . . . . . . . . . . . . . . . . . 72
15.2.2 Quantidade inicial: (7 ; 4) . . . . . . . . . . . . . . . . . . . . . . 73
15.2.3 Quantidade inicial: (5 ; 2) . . . . . . . . . . . . . . . . . . . . . . 73
15.2.4 Quantidade inicial: (5,01 ; 2) . . . . . . . . . . . . . . . . . . . . 74
15.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
16 Problema 16 76
16.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
16.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
16.2.1 Item a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
16.2.2 Items b e c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
16.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
17 Problema 17 79
17.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
17.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
17.3 Código associado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
18 Problema 18 86
18.1 Enunciado do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
18.2 Solução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3
Introdução
Na tentativa de descrever os fenômenos f́ısicos presentes no mundo, é muito comum
deparar-se com equações diferenciais. Na verdade, praticamente tudo que observamos
no nosso dia-a-dia ou em situações extraordinárias, em todos os lugarese em todas as
escalas de grandeza, é definido através deste tipo de equação.
Na engenharia de processos não poderia ser diferente. Desta forma, é de fundamental
importância que o engenheiro tenha conhecimento e domı́nio de diversas metodologias de
resolução para estes problemas.
Neste trabalho, são resolvidos dezoito problemas selecionados de equações diferenciais
aplicadas à engenharia de processos. A ferramenta computacional utilizada foi Python
3.8.5 (default, Jul 28 2020, 12:59:40). Os códigos foram executados no sistema ope-
racional Ubuntu 20.04.1 LTS. As bibliotecas Matplotlib 3.3.1 (gráficos), Numpy 1.19.1
(manipulação de arrays) e Scipy 1.5.3 (computação cient́ıfica) foram utilizadas. Como
biblioteca auxiliar, foi utilizada a TikzPlotLib, que converte gráficos do Matplotlib para
LATEX.
4
1 Problema 1
1.1 Enunciado do problema
Um balanço de massa para um composto qúımico em um reator de mistura pode ser
escrito como
V
dc
dt
= F −Qc− kV c2
onde V é o volume (12 m3), c é a concentração (g/m3), F é a taxa de alimentação
(175 g/min), Q é a vazão (1 m3/min) e k é uma taxa de reação de segunda ordem
(0, 15 m3/g/min). Se c(0) = 0, resolva a EDO até que a concentração atinja um ńıvel
estável. Utilize o método de Euler expĺıcito e o método de Runge-Kutta de quarta ordem
para resolver este problema.
Desafio: se ignorarmos o fato de que as concentrações devem ser positivas, encontre
um intervalo de condições iniciais tal que você obteria uma trajetória muito diferente da
que foi obtida com c(0) = 0. Relacione seus resultados com as soluções estacionárias.
1.2 Solução
Para o problema proposto, a resolução da equação diferencial ordinária é apresentada
na Figura 1.1, com uma discretização que considera um espaçamento de 0, 1 min na
dimensão do tempo.
0 0.5 1 1.5 2 2.5
0
2
4
6
8
10
Tempo (min)
C
on
ce
n
tr
aç
ão
(g
/m
3
)
Euler Expĺıcito
Runge Kutta 4ª Ordem
Figura 1.1: Solução do problema 1 através do método de Euler expĺıcito e do método de
Runge-Kutta de quarta ordem.
As curvas estão suficientemente próximas, especialmente nas regiões lineares do gráfico.
Considera-se o resutado obtido pelo método de Runge-Kutta de quarta ordem mais
confiável devido ao número de parâmetros calculados para obtenção de um dado ponto.
Resposta ao desafio:
5
A solução estacionária para o problema é
F −Qc− kV c2 = 0.
Este polinômio em c tem duas soluções. Naturalmente, uma delas é a encontrada na
resolução gráfica. Para encontrar a outra, basta utilizar a fórmula de Baskhara para a
equação, que fornece:
c2 =
−Q+
√
Q2 + 4kFV
−2kV
= −10, 14 g/m3
Utilizando o valor desta solução estacionária como chute inicial e resolvendo a equação
diferencial ordinária com o método de Runge-Kutta de quarta ordem, obtém-se a solução
da Figura 1.2.
0 0.5 1 1.5 2 2.5
−10.6
−10.4
−10.2
−10
−9.8
−9.6
Tempo (min)
C
on
ce
n
tr
aç
ão
(g
/m
3
)
Figura 1.2: Solução do desafio do problema 1.
1.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-06T22:22:00-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T01:48:01-02:00
5
6
7 import numpy as np
8 import matplotlib.pyplot as plt
9 import euler_explicito
10 import RK_4
11
12 # Parâmetros do problema
6
13
14 V = 12 # m^3
15 F = 175 # g/min
16 Q = 1 # m^3/min
17 k = 0.15 # m^3/g/min
18 c0 = 0 # concentraç~ao inicial
19 dt = 0.1 # stepsize no tempo
20
21 # EDO a ser resolvida
22
23 def dcdt(t, c):
24 dcdt = 1 / V * (F - Q * c - k * V * c ** 2)
25 return dcdt
26
27 # Resoluç~ao
28
29 t_RK, c_RK = RK_4.RK_4(0, 2.5, c0, dt, dcdt)
30 t_euler, c_euler = euler_explicito.Euler(0, 2.5, c0, dt, dcdt)
31
32 # Plot
33
34 plt.plot(t_euler, c_euler, color="#E00650", label="Euler Explı́cito")
35 plt.plot(t_RK, c_RK, color="#2D2E7B", label="Runge Kutta 4ª Ordem")
36 plt.xlabel("Tempo (min)")
37 plt.ylabel("Concentraç~ao ($g/m^3$)")
38 plt.xlim(0, 2.5)
39 plt.ylim(0, 10)
40 plt.legend()
41 plt.grid()
42 plt.show()
1 # @Author: eduardo
2 # @Date: 2020-11-06T22:22:00-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T01:10:04-02:00
5
6 import numpy as np
7 import matplotlib.pyplot as plt
8 import euler_explicito
9 import RK_4
10
11 # Parâmetros do problema
12
13 V = 12 # m^3
14 F = 175 # g/min
15 Q = 1 # m^3/min
16 k = 0.15 # m^3/g/min
7
17 c0 = -10.141822724758168 # concentraç~ao inicial
18 dt = 0.5 # stepsize no tempo
19
20
21 def dcdt(t, c):
22 dcdt = 1 / V * (F - Q * c - k * V * c**2)
23 return dcdt
24
25
26 t_RK, c_RK = RK_4.RK_4(0, 2.5, c0, dt, dcdt)
27 t_euler, c_euler = euler_explicito.Euler(0, 2.5, c0, dt, dcdt)
28
29 plt.plot(t_euler, c_euler, color="#E00650", label="Euler Explı́cito")
30 plt.plot(t_RK, c_RK, color="#2D2E7B", label="Runge Kutta 4ª Ordem")
31 plt.xlabel("Tempo (min)")
32 plt.ylabel("Concentraç~ao ($g/m^3$)")
33 plt.xlim(0, 2.5)
34 plt.grid()
35 plt.show()
8
2 Problema 2
2.1 Enunciado do problema
Se cin = cb(1− e−0,12t), calcule a concentração no escoamento de sáıda de uma substância
conservativa (sem reações) para um único tanque de mistura como uma função do tempo.
Utilize o método de Euler expĺıcito, o método de Runge-Kutta de quarta ordem e um
solver adequado para fazer os cálculos. Use os valores de cb = 40 mg/m
3, Q = 6 m3/min,
V = 100 m3 e c0 = 20 mg/m
3. Faça os cálculos de t = 0 a 100min. Trace seus resultados
com a concentração no escoamento de entrada em função do tempo.
2.2 Solução
Para uma discretização com ∆t = 2 s, apresenta-se a solução para o problema na Fi-
gura 2.1. Note como a resposta obtida através do solver e do método de Runge-Kutta
de quarta ordem são bastante próximas entre si. O método de Euler expĺıcito apresenta
boa precisão em relação aos outros métodos em algumas regiões do gráfico, todavia, em
intervalos onde há grandes mudanças na derivada do gráfico, as soluções apresentam
diferenças consideráveis entre si.
0 20 40 60 80 100
20
25
30
35
40
Tempo (min)
C
on
ce
n
tr
aç
ão
(m
g
/m
3
)
Euler Expĺıcito
Runge Kutta 4ª Ordem
Solver
Figura 2.1: Solução do problema 2.
Ao diminuir-se o espaçamento entre os pontos da discretização, percebe-se como as
respostas utilizando diferentes métodos passam a convergir entre si. Para um ∆t = 0, 1
as três curvas se tornam praticamente indistingúıveis no gráfico, de modo que optou-se
por apresentar os gráficos para os três métodos utilizados de forma individual devido à
grande proximidade dos resultados obtidos.
A Figura 2.2 apresenta a solução pelo método de Euler expĺıcito.
9
0 20 40 60 80 100
20
25
30
35
40
Tempo (min)
C
on
ce
n
tr
aç
ão
(m
g
/m
3
)
Euler Expĺıcito
Figura 2.2: Solução do problema 2 através do método de Euler expĺıcito.
A Figura 2.3 apresenta a solução pelo método de Runge-Kutta de quarta ordem.
0 20 40 60 80 100
20
25
30
35
40
Tempo (min)
C
on
ce
n
tr
aç
ão
(m
g
/m
3
)
Runge Kutta 4ª Ordem
Figura 2.3: Solução do problema 2 através do método de Runge-Kutta de quarta ordem.
A Figura 2.4 apresenta a solução pelo solver solve ivp do pacote Scipy de python com
o método de Radau.
10
0 20 40 60 80 100
20
25
30
35
40
Tempo (min)
C
on
ce
n
tr
aç
ão
(m
g
/m
3
)
Solver
Figura 2.4: Solução do problema 2 através de solver.
2.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-07T16:53:42-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T01:48:36-02:00
5
6 import numpy as np
7 import matplotlib.pyplot as plt
8 import euler_explicito
9 import RK_4
10 import scipy.integrate
11
12 # Parâmetros do problema
13
14 V = 100 # m^3
15 Q = 6 # m^3/min
16 cb = 40 # mg/m^3
17 c0 = 20 # concentraç~ao inicial
18 dt = 0.1 # stepsize no tempo
19
20 # EDO a ser resolvida
21
22 def dcdt(t, c):
23 dcdt = Q / V * (cb * (1 - np.exp(-0.12* t)) - c)
24 return dcdt
25
26 # Resoluç~ao
11
27
28 t_euler, c_euler = euler_explicito.Euler(0, 100, c0, dt, dcdt)
29 t_RK, c_RK = RK_4.RK_4(0, 100, c0, dt, dcdt)
30 sol = scipy.integrate.solve_ivp(dcdt, (0, 100), [c0], method="Radau", max_step = dt)
31
32 # Plot
33
34 plt.plot(t_euler, c_euler, color="#E00650", label="Euler Explı́cito")
35 plt.plot(t_RK, c_RK, color="#2D2E7B", label="Runge Kutta 4ª Ordem")
36 plt.plot(sol.t, sol.y[0], color="#F1A91B", label="Solver")
37 plt.xlabel("Tempo (min)")
38 plt.ylabel("Concentraç~ao ($mg/m^3$)")
39 plt.xlim(0, 100)
40 plt.legend()
41 plt.grid()
42 plt.show()
12
3 Problema 3
3.1 Enunciado do problema
A água do mar com uma concentração de 8000 g/m3 é bombeada para um tanque bem
misturado a uma vazão de 0, 6 m3/h. Por causa de um defeito no trabalho do projeto, a
água está evaporando do tanque a uma vazão de 0, 025 m3/h. A solução salina deixa o
tanque a uma vazão de 0, 6 m3/h.
(a) Se o tanque continha originalmente 1 m3 da solução de entrada, quanto tempo
depois de ligar a bomba de sáıda o tanque secará?
(b) Use métodos numéricos para determinar a concentração de sal no tanque como uma
função do tempo.
3.2 Solução
3.2.1 Item a
A variação de volume no tanque é dada simplesmente pela quantidade de água que
evapora, já que as vazões de entrada e de sáıda são iguais.
Assim, a Equação 3.1 representa um bom modelo para o problema:
dV
dt
= −0, 25 (3.1)
Resolvendo com os métodos de Euler e Runge-Kutta de quarta ordem:
0 5 10 15 20 25 30 35 40
0
0.2
0.4
0.6
0.8
1
Tempo (h)
V
ol
u
m
e
d
o
ta
n
q
u
e
(m
3
)
Euler Expĺıcito
Runge Kutta 4ª Ordem
Figura 3.1: Solução do problema 3 através dos métodos de Euler e Runge Kutta de quarta
ordem.
13
3.2.2 Item b
Para este problema, a modelagem obtida baseia-se em um balanço de massa para o sal
no tanque, representado por
dcV
dt
= Q(cin − c).
Aplicando a regra da cadeia obtém-se
V
dc
dt
+ C
dV
dt
= Q(cin − c).
Sabe-se do item a que
dV
dt
= −Qvap. Aplicando e simplificando, é posśıvel chegar a
Equação 3.2, equação diferencial a ser resolvida.
dc
dt
=
Q(cin − c) +Qvapc
V (t)
(3.2)
0 2 4 6 8 10
8,000
8,050
8,100
8,150
8,200
8,250
8,300
8,350
Tempo (h)
C
on
ce
n
tr
aç
ão
(m
g
/m
3
)
Euler Expĺıcito
Runge-Kutta 4ª Ordem
Figura 3.2: Solução do problema 3 (item b) através dos métodos de Euler e Runge Kutta
de quarta ordem.
3.3 Código associado
3.3.1 Item a
1 # @Author: eduardo
2 # @Date: 2020-11-07T18:35:38-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-19T18:43:40-02:00
5
6 import numpy as np
14
7 import matplotlib.pyplot as plt
8 import euler_explicito
9 import RK_4
10
11 Q = 0.6 # m^3/h
12 dt = 0.01 # stepsize no tempo
13 Qv = 0.025 # m^3/h
14 V0 = 1 # m^3
15
16
17 def dVdt(t, V):
18 dVdt = -Qv
19 return dVdt
20
21
22 t_euler, V_euler = euler_explicito.Euler(0, 40, V0, dt, dVdt)
23 t_RK, V_RK = RK_4.RK_4(0, 40, V0, dt, dVdt)
24
25 plt.plot(t_euler, V_euler, color="#E00650", label="Euler Explı́cito")
26 plt.plot(t_RK, V_RK, linestyle=":", color="#2D2E7B", label="Runge Kutta 4ª Ordem")
27 plt.xlabel("Tempo (h)")
28 plt.ylabel("Volume do tanque ($m^3$)")
29 plt.xlim(0, 40)
30 plt.ylim(0, 1)
31 plt.legend()
32 plt.grid()
33
34 plt.show()
3.3.2 Item b
1 # @Author: eduardo
2 # @Date: 2020-11-07T17:27:26-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-19T19:12:42-02:00
5 import numpy as np
6 import matplotlib.pyplot as plt
7 from scipy.integrate import odeint
8 import RK_4
9 import euler_explicito
10
11
12 Q = 0.6 # m^3/min
13 Qv = 0.025 # m^3/min
14 ci = 8000 # concentraç~ao inicial
15 dt = 0.01 # stepsize no tempo
16
15
17
18 def dcdt(t, c):
19 dcdt = (Q * (ci - c) + Qv * c) / (1 - 0.025 * t)
20 return dcdt
21
22
23 t, c = RK_4.RK_4(0, 10, 8000, dt, dcdt)
24 t_euler, V_euler = euler_explicito.Euler(0, 10, 8000, dt, dcdt)
25
26 # c = odeint(dcdt, 8000, t)
27
28 plt.plot(t_euler, V_euler, color="#E00650", label="Euler Explı́cito")
29 plt.plot(t, c, linestyle ='--',color="#2D2E7B", label='Runge-Kutta 4ª Ordem')
30 plt.xlabel("Tempo (h)")
31 plt.ylabel("Concentraç~ao ($mg/m^3$)")
32 plt.ylim(8000, 8350)
33 plt.xlim(0, 10)
34 plt.grid()
35 plt.legend()
36 plt.show()
16
4 Problema 4
4.1 Enunciado do problema
Uma esfera de gelo que tem 6 cm de diâmetro é removida de um congelador a 0 C e é
colocada em uma tela de rede à temperatura ambiente Ta = 20 C. Qual será o diâmetro
da esfera de gelo como uma função do tempo fora do congelador (supondo que toda a
água que derreteu caia imediatamente pela tela)? O coeficiente de transferência de calor
h para uma esfera em um quarto sem vento é cerca de 3 W/(m2K). O fluxo de calor da
esfera de gelo para o ar é dado por
Fluxo =
q
A
= h(Ta − T )
onde q é o calor e A é a área da superf́ıcie de uma esfera. Use um método numérico
para fazer seus cálculos. Observe que o calor latente de fusão é 333 kJ/kg e que a
densidade do gelo é aproximadamente 917 kg/m3.
4.2 Solução
Como a superf́ıcie da esfera está em constante fusão, sabe-se que sua temperatura está
fixa em 0 °C. Além disso, considera-se que todo calor fornecido para a esfera na forma
convectiva é utilizado no processo do derretimento.
Desta forma, a variação do volume da esfera pode ser dado pela Equação 4.1:
dV
dt
=
−A(t)h(Ta − T )
Qfρ
(4.1)
Fazendo uma mudança de variáveis para escrever o volume e a área superficial da
esfera em função de seu raio,
V =
4
3
πr3
dV = 4πr2dr
A = 4πr2
Substituindo,
4πr2dr
dt
=
4πr2h(Ta − T )
Qfρ
Simplificando, chega-se a equação diferencial ordinária que deve ser resolvida, expressa
pela Equação 4.2:
dr
dt
=
h(Ta − T )
Qfρ
(4.2)
A solução numérica para este problema é apresentada na Figura 4.1.
17
0 0.2 0.4 0.6 0.8 1 1.2 1.4
·105
0
1
2
3
4
5
6
·10−2
Tempo (min)
D
iâ
m
et
ro
(m
)
Runge Kutta 4ª Ordem
Figura 4.1: Solução do problema 4 através do método de Runge Kutta de quarta ordem.
4.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-07T18:48:50-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-19T19:18:15-02:00
5
6 import numpy as np
7 import matplotlib.pyplot as plt
8 import euler_explicito
9 import RK_4
10 from scipy.integrate import odeint
11
12 # Parâmetros do problema
13
14 Tinf = 20 # m^3
15 T = 0 # ºC
16 h = 3 # mg/m^3
17 Qfusao = 333000 # J/Kg
18 rho = 917 # Kg/m^3
19 dt = 1000 # stepsize no tempo
20 r0 = 0.03 # Raio no tempo 0
21
22
23 def dcdt(t, r):
24 drdt = -h * (Tinf - T) / (Qfusao * rho)
25 return drdt
18
26
27
28 t_RK, r_RK = RK_4.RK_4(0, 160000, r0, dt, dcdt)
29
30 plt.plot(t_RK, 2 * r_RK, color="#E00650", label="Runge Kutta 4ª Ordem")
31 plt.xlabel("Tempo (min)")
32 plt.ylabel("Diâmetro (m)")
33 plt.xlim(0, 153000)
34 plt.ylim(0,r0*2)
35 plt.legend()
36 plt.grid()
37 plt.show()
19
5 Problema 5
5.1 Enunciado do problema
As seguintes equações definem a concentração de trés reagentes:
dca
dt
= −10cacc + cb
dcb
dt
= 10cacc − cb
dcc
dt
= −10cacc + cb − 2cc
Se as condições iniciais forem ca = 50, cb = 0 e cc = 40, encontre as concentrações
para instantes entre 0 c 3 s.
5.2 Solução
A solução numérica para o problema é apresentada na Figura 5.1.
0 0.5 1 1.5 2 2.5 3
0
10
20
30
40
50
Tempo (s)
C
on
ce
n
tr
aç
ão
(-
)
Ca
Cb
Cc
Figura 5.1: Solução do problema 5 utilizando o solver com método de impĺıcito de Radau.
Para métodos expĺıcitos, este problema exigiu passos de magnitude muito pequena,
de forma que seu custo computacional aumentou de forma gigantesca.
Observe que neste problema tudo acontece em tempos muito próximos de zero. Para
um intervalo em x menor, é posśıvel visuzizar melhor como as concentrações variam antes
de atingir o estado permanente.
20
0 2 · 10−2 4 · 10−2 6 · 10−2 8 · 10−2 0.1
0
10
20
30
40
50
Tempo(s)
C
on
ce
n
tr
aç
ão
(-
)
Ca
Cb
Cc
Figura 5.2: Solução do problema 5 utilizando o solver com método de impĺıcito de Radau.
5.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-07T19:49:54-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-19T19:31:50-02:00
5
6 import numpy as np
7 import matplotlib.pyplot as plt
8 import scipy.integrate
9
10
11 ca0 = 50
12 cb0 = 0
13 cc0 = 40
14 c0 = [ca0, cb0, cc0]
15 dt = 0.5
16
17
18 def solve(t, c):
19 c1, c2, c3 = c
20 dcadt = -10 * c1 * c3 + c2
21 dcbdt = 10 * c1 * c3 - c2
22 dccdt = -10 * c1 * c3 + c2 - 2 * c3
23 return dcadt, dcbdt, dccdt
24
25
21
26 sol = scipy.integrate.solve_ivp(solve, (0, 3), c0, method="Radau", max_step=dt)
27
28 plt.plot(sol.t, sol.y[0], color="#E00650", label="Ca")
29 plt.plot(sol.t, sol.y[1], color="#2D2E7B", label="Cb")
30 plt.plot(sol.t, sol.y[2], color="#E8B038", label="Cc")
31 plt.xlim(0, 3)
32 plt.ylim(0, 50)
33 plt.grid()
34 plt.legend()
35 plt.xlabel('Tempo (s)')
36 plt.ylabel('Concentraç~ao (-)')
37
38 plt.show()
22
6 Problema 6
6.1 Enunciado do problema
A reação A → B ocorre em dois reatores em série. Os reatores estão bem misturados,
mas não cstão no cstado estacionário. O balanço de massa para cada reator de tanque
agitado é mostrado a seguir:
dCA1
dt
=
1
τ
(CA0 − CA1)− kCA1
dCB1
dt
= −1
τ
CB1 + kCA1
dCA2
dt
=
1
τ
(CA1 − CA2)− kCA2
dCB2
dt
=
1
τ
(CB1 − CB2) + kCA2
onde CA0 é a concentração de A na entrada do primeiro reator, CA1 é a concentração
de A na sáıda do primeiro reator (e na entrada do segundo), CA2 é a concentração de A
na sáıda do segundo reator, CB1 é a concentração de B na sáıda do primeiro reator (e
na entrada do segundo), CB2 é a concentração de B no segundo reator, τ é o tempo de
residência para cada reator e k é a taxa constante para a reação de A para produzir B.
Se CA0 for igual a 20, encontre as concentrações de A e B em ambos os reatores durante
os primeiros 10 minutos de operação. Use k = 0, 12/min e τ = 5 min e suponha que as
condições iniciais de todas as variáveis dependentes sejam zero.
6.2 Resolução
A solução numérica utilizando o método de Runge-Kutta de quarta ordem para o pro-
blema é apresentada na Figura 6.1.
0 2 4 6 8 10
0
2
4
6
8
10
12
Tempo (min)
C
on
ce
n
tr
aç
ão
(-
)
Ca1
Cb1
Ca2
Cb2
Figura 6.1: Solução do problema 6 utilizando o método de Runge-Kutta de quarta ordem.
23
6.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-07T20:49:22-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-19T20:53:08-02:00
5
6 import numpy as np
7 import matplotlib.pyplot as plt
8 import RK_4
9 import scipy.integrate
10
11 k = 0.12
12 tau = 5 # minutos
13 c0 = np.zeros(4)
14 ca0 = 20
15 dt = 0.01
16
17
18 def solve(t, c):
19 ca1, cb1, ca2, cb2 = c
20 dca1dt = 1 / tau * (ca0 - ca1) - k * ca1
21 dcb1dt = -1 / tau * cb1 + k * ca1
22 dca2dt = 1 / tau * (ca1 - ca2) - k * ca2
23 dcb2dt = 1 / tau * (cb1 - cb2) + k * ca2
24 return dca1dt, dcb1dt, dca2dt, dcb2dt
25
26
27 t, c = RK_4.RK_4_SE(0, 10, c0, dt, solve)
28
29 sol2 = scipy.integrate.solve_ivp(solve, (0, 10), c0, method="Radau", max_step=dt)
30
31
32 plt.plot(t, c[0], color="#E00650", label="Ca1")
33 plt.plot(t, c[1], color="#2D2E7B", label="Cb1")
34 plt.plot(t, c[2], color="#E8B038", label="Ca2")
35 plt.plot(t, c[3], label="Cb2")
36 # plt.plot(sol2.t, sol2.y[3], label="Solver")
37 plt.xlim(0, 10)
38 plt.ylim(0, 12.5)
39 plt.grid()
40 plt.legend()
41 plt.xlabel('Tempo (min)')
42 plt.ylabel('Concentraç~ao (-)')
43
44 plt.show()
24
7 Problema 7
7.1 Enunciado do problema
Um reator de batelada pode ser descrito pelas seguintes equações:
dC
dt
= −e(−10/(T+273))C
dT
dt
= 1000e(−10/(T+273))C − 10(T − 20)
onde C é a concentração do reagente e T a temperatura do reator. Inicialmente, o
reator está a 15 C e tem uma concentração do reagente C de 1, 0 mol/L. Encontre a
concentração e a temperatura do reator como uma função do tempo.
7.2 Resolução
A solução numérica utilizando o método de Radau para a concentração no problema é
apresentada na Figura 7.1.
0 0.5 1 1.5 2 2.5 3
0
0.2
0.4
0.6
0.8
1
t (-)
C
(m
ol
/L
)
C
Figura 7.1: Solução do problema 7 utilizando o método de Radau.
A solução numérica utilizando o método de Radau para a concentração no problema
é apresentada na Figura 7.1.
25
0 0.5 1 1.5 2 2.5 3
20
40
60
80
t (-)
T
(°
C
)
T
Figura 7.2: Solução do problema 7 utilizando o método de Radau.
7.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-07T23:45:47-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:30:55-02:00
5
6 import numpy as np
7 from numpy import exp
8 import matplotlib.pyplot as plt
9 import RK_4
10 import scipy.integrate
11
12 c0 = 1 # mol/L
13 T0 = 15 # ºC
14 x0 = [c0, T0]
15 dt = 0.5
16
17
18 def solve(t, y):
19 C, T = y
20 dCdt = -exp(-10 / (T + 273)) * C
21 dTdt = 1000 * exp(-10 / (T + 273)) * C - 10 * (T - 20)
22 return dCdt, dTdt
23
24
25 t, c = RK_4.RK_4_SE(0, 3, x0, dt, solve)
26
26
27 sol = scipy.integrate.solve_ivp(solve, (0, 10),
28 x0,
29 method="Radau",
30 max_step=dt)
31
32 fig, axs = plt.subplots(1, 2)
33 axs[0].plot(sol.t, sol.y[0], color="#E00650", label="C")
34 axs[1].plot(sol.t, sol.y[1], color="#2D2E7B", label="T")
35 axs[1].set_ylabel('Temperatura')
36 axs[0].set_xlabel('Tempo')
37 axs[1].set_xlabel('Tempo')
38 plt.show()
27
8 Problema 8
8.1 Enunciado do problema
Considere o reator de mistura perfeita encamisado, conduzindo reaçio quimica exotérmica,
representado ao lado. Pede-se:
(a) Obter o modelo matemático que representa este processo.
(b) Analisar a dinâmica do processo para o cenário em que a temperatura de entrada
é mantida fixa em 350 K e, em t = 0 a concentração no reator é de 0, 1 mol/L a
temperatura do CSTR é 600 K e a temperatura da camisa é 500 K.
(c) Em t = 50 s a temperatura na alimentação do CSTR é alterada para 298 K. Analise
a dinâmica do processo frente a esta modificaç̃ıo.
Realize sua análise utilizando os métodos de Euler e Runge-Kutta de quarta ordem
(RK4), bem como um solver adequado. Comente cada linha do algoritmo e apresente
uma discuss̃ıo sobre o processo em estudo.
Informações adicionais:
Alimentação do reator Camisa Outros Parâmetros
F = 25 L/s Fj = 5L/s CA,0 = 4 mol/L
V = 250 L Vj = 40 L k0 = 800 s
−1
ρ = 1000 kg/m3 ρ = 800 kg/m3 E/R = 4500 K
cp = 2500 J/Kg ·K cp = 5000 J/Kg ·K (−∆H) = 250 kJ/mol
T0 = 350 K Tj,0 = 25 C UA = 20 kW/K
8.2 Solução
8.2.1 Modelo
Considera-se para o problema a reação irreverśıvel de primeira ordem
A −−→ B.
O consumo da espécie A é dado pela Equação 8.1:
dCa
dt
=
Ca0 − Ca
τ
− kCa (8.1)
Sendo τ o tempo espacial do reator, dado pela razão entre V e F. A geração do produto
B, por outro lado, pode ser expressa pela Equação 8.2:
dCb
dt
=
Cb
τ
+ kCa (8.2)
Um balanço energético para o interior do reator permite a obtenção da Equação 8.3:
dT
dt
=
T0 − T
τ
+
∆HkCa
ρCp
− UA(T − Tj)
ρCpV
(8.3)
O balanço para a camisa térmica do reator fornece a Equação 8.4:
28
dTj
dt
=
T0j − Tj
τj
+
UA(T − Tj)
ρjCpjVj
(8.4)
Finalmente, a constante K para a reação é dada pela Equação 8.5:
K = K0 · exp
(
− E
RT
)
(8.5)
8.2.2 Resolução - Euler
A resolução do problema consiste na obtenção na solução para estas cinco equações aco-
pladas. As concentrações em função do tempo são apresentadas na Figura 8.1. As
temperaturas em função do tempo são apresentadas na Figura 8.2.
0 10 20 30 40 50
0
1,000
2,000
3,000
4,000
Tempo (s)
C
on
ce
n
tr
aç
ão
(m
ol
/m
3
)
Ca
Cb
Figura 8.1: Solução do problema 8 utilizando o método de Euler.
29
0 10 20 30 40 50
450
500
550
600
650
Tempo (s)
T
em
p
er
at
u
ra
(K
)
T
Tj
Figura 8.2: Solução do problema 8 utilizando o método de Euler.
8.2.3Resolução - Runge-Kutta
As concentrações em função do tempo são apresentadas na Figura 8.3. As temperaturas
em função do tempo são apresentadas na Figura 8.4.
0 10 20 30 40 50
0
1,000
2,000
3,000
4,000
Tempo (s)
C
on
ce
n
tr
aç
ão
(m
ol
/m
3
)
Ca
Cb
Figura 8.3: Solução do problema 8 utilizando o método de Runge-Kutta de quarta ordem.
30
0 10 20 30 40 50
450
500
550
600
650
Tempo (s)
T
em
p
er
at
u
ra
(K
)
T
Tj
Figura 8.4: Solução do problema 8 utilizando o método de Runge-Kutta de quarta ordem.
8.2.4 Resolução - Solver
As concentrações em função do tempo são apresentadas na Figura 8.5. As temperaturas
em função do tempo são apresentadas na Figura 8.6.
0 10 20 30 40 50
0
1,000
2,000
3,000
4,000
Tempo (s)
C
on
ce
n
tr
aç
ão
(m
ol
/m
3
)
Ca
Cb
Figura 8.5: Solução do problema 8 utilizando o solver com método de Radau.
31
0 10 20 30 40 50
450
500
550
600
650
Tempo (s)
T
em
p
er
at
u
ra
(K
)
T
Tj
Figura 8.6: Solução do problema 8 utilizando o solver com método de Radau.
8.3 Item c
Neste item, avaliamos uma perturbação na temperatura de alimentação do reator a partir
de 50 segundos.
8.3.1 Resolução - Euler
As concentrações em função do tempo são apresentadas na Figura 8.7. As temperaturas
em função do tempo são apresentadas na Figura 8.8.
0 20 40 60 80 100 120
0
1,000
2,000
3,000
4,000
Tempo (s)
C
on
ce
n
tr
aç
ão
(m
ol
/m
3
)
Ca
Cb
Figura 8.7: Solução do problema 8 (item c) utilizando o método de Euler.
32
0 20 40 60 80 100 120
300
400
500
600
Tempo (s)
C
on
ce
n
tr
aç
ão
(m
ol
/m
3
)
T
Tj
Figura 8.8: Solução do problema 8 (item c) utilizando o método de Euler.
8.3.2 Resolução - Runge-Kutta
As concentrações em função do tempo são apresentadas na Figura 8.9. As temperaturas
em função do tempo são apresentadas na Figura 8.10.
0 20 40 60 80 100 120
0
1,000
2,000
3,000
4,000
Tempo (s)
C
on
ce
n
tr
aç
ão
(m
ol
/m
3
)
Ca
Cb
Figura 8.9: Solução do problema 8 (item c) utilizando o método de Runge-Kutta de
quarta ordem.
33
0 20 40 60 80 100 120
300
400
500
600
Tempo (s)
C
on
ce
n
tr
aç
ão
(m
ol
/m
3
)
T
Tj
Figura 8.10: Solução do problema 8 (item c) utilizando o método de Runge-Kutta de
quarta ordem.
8.3.3 Resolução - Solver
As concentrações em função do tempo são apresentadas na Figura 8.11. As temperaturas
em função do tempo são apresentadas na Figura 8.12.
0 20 40 60 80 100 120
0
1,000
2,000
3,000
4,000
Tempo (s)
C
on
ce
n
tr
aç
ão
(m
ol
/m
3
)
Ca
Cb
Figura 8.11: Solução do problema 8 (item c) utilizando o solver com método de Radau.
34
0 20 40 60 80 100 120
300
400
500
600
Tempo (s)
T
em
p
er
at
u
ra
(K
)
T
Tj
Figura 8.12: Solução do problema 8 (item c) utilizando o solver com método de Radau.
Este caso pode ser considerado como um ótimo exemplo de como é importante simu-
lar processos e de uma análise de sensibilidade detalhada. A resolução deste problema
permitiu a visualização de como a variação de um único parâmetro - neste caso a tem-
peratura de entrada de uma corrente no reator - altera completamente sua dinâmica.
Em um caso prático, caso isto acontecesse, a concentração do produto de interesse seria
zerada, gerando sem dúvidas muitos problemas técnicos e econômicos na planta que este
reator está operando.
8.4 Código associado
8.4.1 Item b
1 # @Author: eduardo
2 # @Date: 2020-11-19T13:45:19-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:35:24-02:00
5
6 import matplotlib.pyplot as plt
7 import numpy as np
8 from numpy import exp
9 import scipy.integrate
10 import RK_4
11 import euler_explicito
12
13
14 # Parâmetros do modelo
15 # Alimentaç~ao do reator
16
17 F = 0.025 # m^3/s
35
18 V = 0.250 # m^3
19 rho = 1000 # kg/m^3
20 cp = 2500 # J/kg.k
21 T0 = 350 # K
22 tau = V / F
23
24 # Camisa
25
26 F_j = 0.005 # m^3/s
27 V_j = 0.04 # m^3
28 rho_j = 800 # kg/m^3
29 cp_j = 5000 # J/kg.K
30 T0_j = 25 + 273.15 # K
31 tau_j = V_j / F_j
32
33 # Demais parâmetros
34
35 Ca0 = 4000 # mol/m^3
36 k0 = 800 # s^-1
37 ER = 4500 # K
38 dH = 250000 # J/mol
39 UA = 20000 # W/K
40
41 dt = 0.1
42 x0 = [100, 0, 600, 500]
43
44
45 def solve(t, x):
46 Ca, Cb, T, Tj = x
47 K = k0 * np.exp(-ER / T)
48 dCadt = (Ca0 - Ca) / tau - K * Ca
49 dCbdt = -Cb / tau + K * Ca
50 dTdt = (T0 - T) / tau + dH * K * Ca / (rho * cp) - UA * (T - Tj) / (rho *
51 cp * V)
52 dTjdt = (T0_j - Tj) / tau_j + UA * (T - Tj) / (rho_j * cp_j * V_j)
53 return dCadt, dCbdt, dTdt, dTjdt
54
55
56 t, c = euler_explicito.Euler_SE(0, 50, x0, dt, solve)
57
58 fig, axs = plt.subplots(1, 2)
59 axs[0].plot(t, c[0], color="#E00650", label="Ca")
60 axs[0].plot(t, c[1], color="#E8B038", label="Cb")
61
62 axs[0].set_xlabel('Tempo (s)')
63 axs[0].set_ylabel('Concentraç~ao $(mol/m^3)$')
64 axs[0].legend()
65 axs[0].grid(True)
36
66
67 axs[1].plot(t, c[2], color="#2D2E7B", label="T")
68 axs[1].plot(t, c[3], color="#E8B038", label="Tj")
69
70
71 axs[1].set_ylabel('Temperatura (K)')
72 axs[1].set_xlabel('Tempo (s)')
73 axs[1].legend()
74 axs[1].grid(True)
75 plt.title('Euler')
76 plt.show()
77
78 t, c = RK_4.RK_4_SE(0, 50, x0, dt, solve)
79
80 fig, axs = plt.subplots(1, 2)
81 axs[0].plot(t, c[0], color="#E00650", label="Ca")
82 axs[0].plot(t, c[1], color="#E8B038", label="Cb")
83
84 axs[0].set_xlabel('Tempo (s)')
85 axs[0].set_ylabel('Concentraç~ao $(mol/m^3)$')
86 axs[0].legend()
87 axs[0].grid(True)
88
89 axs[1].plot(t, c[2], color="#2D2E7B", label="T")
90 axs[1].plot(t, c[3], color="#E8B038", label="Tj")
91
92
93 axs[1].set_ylabel('Temperatura (K)')
94 axs[1].set_xlabel('Tempo (s)')
95 axs[1].legend()
96 axs[1].grid(True)
97 plt.title('Runge-Kutta')
98 plt.show()
99
100 sol = scipy.integrate.solve_ivp(solve, (0, 50), x0, method="Radau", max_step=dt)
101
102 fig, axs = plt.subplots(1, 2)
103 axs[0].plot(sol.t, sol.y[0], color="#E00650", label="Ca")
104 axs[0].plot(sol.t, sol.y[1], color="#E8B038", label="Cb")
105
106 axs[0].set_xlabel('Tempo (s)')
107 axs[0].set_ylabel('Concentraç~ao $(mol/m^3)$')
108 axs[0].legend()
109 axs[0].grid(True)
110
111 axs[1].plot(sol.t, sol.y[2], color="#2D2E7B", label="T")
112 axs[1].plot(sol.t, sol.y[3], color="#E8B038", label="Tj")
113
37
114
115 axs[1].set_ylabel('Temperatura (K)')
116 axs[1].set_xlabel('Tempo (s)')
117 axs[1].legend()
118 axs[1].grid(True)
119 plt.title('Solver')
120 plt.show()
8.4.2 Item c
1 # @Author: eduardo
2 # @Date: 2020-11-19T20:29:26-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:36:51-02:00
5
6 import matplotlib.pyplot as plt
7 import numpy as np
8 from numpy import exp
9 import scipy.integrate
10 import RK_4
11 import euler_explicito
12
13 # Parâmetros do modelo
14 # Alimentaç~ao do reator
15
16 F = 0.025 # m^3/s
17 V = 0.250 # m^3
18 rho = 1000 # kg/m^3
19 cp = 2500 # J/kg.k
20 tau = V / F
21
22 # Camisa
23
24 F_j = 0.005 # m^3/s
25 V_j = 0.04 # m^3
26 rho_j = 800 # kg/m^3
27 cp_j = 5000 # J/kg.K
28 T0_j = 25 + 273.15 # K
29 tau_j = V_j / F_j
30
31 # Demais parâmetros
32
33 Ca0 = 4000 # mol/m^3
34 k0 = 800 # s^-1
35 ER = 4500 # K
36 dH = 250000 # J/mol
37 UA = 20000 # W/K
38
38
39 dt = 0.1
40 x0 = [100, 0, 600, 500]
41
42
43 def solve(t, x):
44 if t < 50:
45 T0 = 350
46 else:
47 T0 = 298
48 Ca, Cb, T, Tj = x
49 K = k0 * np.exp(-ER / T)
50 dCadt = (Ca0 - Ca) / tau - K * Ca
51 dCbdt = -Cb / tau + K * Ca
52 dTdt = (T0 - T) / tau + dH * K * Ca / (rho * cp) - UA * (T - Tj) / (rho *
53 cp * V)
54 dTjdt = (T0_j - Tj) / tau_j + UA * (T - Tj) / (rho_j * cp_j * V_j)
55 return dCadt, dCbdt, dTdt, dTjdt
56
57 t, c = euler_explicito.Euler_SE(0, 125, x0, dt, solve)
58
59 fig, axs = plt.subplots(1, 2)
60 axs[0].plot(t, c[0], color="#E00650", label="Ca")
61 axs[0].plot(t, c[1], color="#E8B038", label="Cb")62
63 axs[0].set_xlabel('Tempo (s)')
64 axs[0].set_ylabel('Concentraç~ao $(mol/m^3)$')
65 axs[0].legend()
66 axs[0].grid(True)
67
68 axs[1].plot(t, c[2], color="#2D2E7B", label="T")
69 axs[1].plot(t, c[3], color="#E8B038", label="Tj")
70
71
72 axs[1].set_ylabel('Temperatura (K)')
73 axs[1].set_xlabel('Tempo (s)')
74 axs[1].legend()
75 axs[1].grid(True)
76 plt.title('Euler')
77 plt.show()
78
79 t, c = RK_4.RK_4_SE(0, 125, x0, dt, solve)
80
81 fig, axs = plt.subplots(1, 2)
82 axs[0].plot(t, c[0], color="#E00650", label="Ca")
83 axs[0].plot(t, c[1], color="#E8B038", label="Cb")
84
85 axs[0].set_xlabel('Tempo (s)')
39
86 axs[0].set_ylabel('Concentraç~ao $(mol/m^3)$')
87 axs[0].legend()
88 axs[0].grid(True)
89
90 axs[1].plot(t, c[2], color="#2D2E7B", label="T")
91 axs[1].plot(t, c[3], color="#E8B038", label="Tj")
92
93
94 axs[1].set_ylabel('Temperatura (K)')
95 axs[1].set_xlabel('Tempo (s)')
96 axs[1].legend()
97 axs[1].grid(True)
98 plt.title('Runge-Kutta')
99 plt.show()
100
101 sol = scipy.integrate.solve_ivp(solve, (0, 125), x0, method="Radau", max_step=dt)
102
103 fig, axs = plt.subplots(1, 2)
104 axs[0].plot(sol.t, sol.y[0], color="#E00650", label="Ca")
105 axs[0].plot(sol.t, sol.y[1], color="#E8B038", label="Cb")
106
107 axs[0].set_xlabel('Tempo (s)')
108 axs[0].set_ylabel('Concentraç~ao $(mol/m^3)$')
109 axs[0].legend()
110 axs[0].grid(True)
111
112 axs[1].plot(sol.t, sol.y[2], color="#2D2E7B", label="T")
113 axs[1].plot(sol.t, sol.y[3], color="#E8B038", label="Tj")
114
115
116 axs[1].set_ylabel('Temperatura (K)')
117 axs[1].set_xlabel('Tempo (s)')
118 axs[1].legend()
119 axs[1].grid(True)
120 plt.title('Solver')
121 plt.show()
40
9 Problema 9
9.1 Enunciado do problema
O seguinte sistema é um exemplo clássico de EDOs rigidas que podem ocorrer na cinética
de reações quimicas:
dc1
dt
= −0, 013c1 − 1000c1c3
dc2
dt
= −2500c2c3
dc3
dt
= −0, 013c1 − 1000c1c3 − 2500c2c3
Resolva essas EDOs de t = 0 a 50 com condições iniciais c1(0) = c2(0) = 1 e c3(0) = 0
9.2 Solução
A solução numérica utilizando o método impĺıcito Radau para o problema é apresentada
na Figura 6.1.
0 10 20 30 40 50
0
0.2
0.4
0.6
0.8
1
1.2
1.4
Tempo (-)
C
on
ce
n
tr
aç
ão
c1
c2
c3
Figura 9.1: Solução do problema 9 utilizando o método de Radau com um solver.
9.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-08T00:02:53-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:38:42-02:00
5 import numpy as np
6 from numpy import exp
41
7 import matplotlib.pyplot as plt
8 import scipy.integrate
9
10 c10 = 1
11 c20 = 1
12 c30 = 0
13 c0 = [c10, c20, c30]
14 dt = 0.1
15
16
17 def solve(t, c):
18 c1, c2, c3 = c
19 dc1dt = -0.013 * c1 - 1000 * c1 * c3
20 dc2dt = -2500 * c2 * c3
21 dc3dt = -0.013 * c1 - 1000 * c1 * c3 - 2500 * c2 * c3
22 return dc1dt, dc2dt, dc3dt
23
24 sol = scipy.integrate.solve_ivp(solve, (0, 50),
25 c0,
26 method="Radau",
27 max_step=dt)
28
29 plt.plot(sol.t, sol.y[0], color="#E00650", label="c1")
30 plt.plot(sol.t, sol.y[1], color="#2D2E7B", label="c2")
31 plt.plot(sol.t, sol.y[2], color="#E8B038", label="c3")
32 plt.xlabel('Tempo (-)')
33 plt.ylabel('Concentraç~ao')
34 plt.xlim(0, 50)
35 plt.grid()
36 plt.legend()
37 plt.show()
42
10 Problema 10
10.1 Enunciado do problema
Um biofilme com espessura Lf (cm) cresce na superficie de um sólido. Depois de atra-
vessar uma camada de difusão de espessura L (cm) , um composto quimico A se difunde
no biofilme onde é sujeito a uma reação de primeira ordem irreversivel que o converte em
um produto B. O balanço de massa estacionário pode ser usado para deduzir as seguintes
cquações diferenciais ordinárias para o composto A :
D
d2ca
dx2
= 0, 0 ≤ x < L
Df
d2ca
dx2
− kca = 0, L ≤ x < L+ Lf
em que D é o coeficiente de difusão na camada de difusão (0, 8 cm2/dia), Df é o
coeficiente de difusão no biofilme (0, 64 cm2dia) e k é a taxa de primeira ordem para a
conversão de A em B (0, 1 /dia). As seguintes condições de contorno são válidas:
ca = ca,0, em x = 0
dca
dx
= 0, em x = L+ Lf
em que ca,0 é a concentração de água no volume do liquido (100 mol/L. Faça uma
representação esquemática do processo. Use o método das diferenças finitas para calcular
a distribuição estacionária de A de x = 0 a L+Lf , em que L = 0, 008 cm e Lf = 0, 004 cm.
Utilize diferenças finitas centradas com ∆x = 0, 001 cm.
10.2 Solução
Para 0 ≤ x < L, ou seja, na região do sólido,
D
dca
dx
= 0.
Como, D 6= 0,
dca
dx
= 0.
Uma discretização utilizando o método das diferenças finitas fornece a Equação 10.1:
ci+1a − 2cia + ci−1a
∆x2
= 0 (10.1)
Para o domı́nio do biofilme, L ≤ x < L+ Lf ,
Df
d2ca
d2x
− kca = 0.
A discretização para esta região é dada pela Equação 10.2:
ci+1a +
(
k∆x2
D
− 2
)
cia + c
i−1
a = 0 (10.2)
43
Note como para este problema, para calcular cia são necessários os valores de c
i
a, c
i−1
a e
ci+1a . Destes, para cada i, somente o valor passado é conhecido. Isto sugere a necessidade
do desenvolvimento de um sistema de equações para a solução do problema.
Além disso, é necessário incluir as condições de contorno no problema. Desta forma,
ao discretizar-se o domı́nio de x com n pontos, sabemos da condição de fluxo nulo em
x = L+ Lf que
cna − cn−1a = 0.
Da condição em x = 0,
c1a = c
0
a.
Estas duas condições definem, respectivamente, a primeira e a última linha da matriz
associada ao sistema linear. A matriz do sistema é grande e esparsa (tridiagonal), mesmo
para uma discretização com intervalos grandes, como sugerido no exerćıcio. Solucionando-
se o sistema linear, obtém-se um gráfico como o da Figura 10.1.
0 0.2 0.4 0.6 0.8 1 1.2
·10−2
99.9996
99.9997
99.9998
99.9999
100
x (cm)
C
on
ce
n
tr
aç
ão
(m
ol
/L
)
Figura 10.1: Solução para o problema 10 obtido através do método de diferenças finitas.
10.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-08T10:08:00-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:39:17-02:00
5
6
7 import numpy as np
8 import matplotlib.pyplot as plt
9
44
10 L = 0.008 # cm
11 Lf = 0.004 # cm
12 D = 0.8 # cm^2 / dia
13 Df = 0.64 # cm^2 / dia
14 k = 0.1 # dia^-1
15 dx = 0.001 # cm
16 n = int(L / dx)
17 ca0 = 100 # mol / L
18 m = int((L + Lf) / dx)
19
20 A = np.zeros((m + 1, m + 1))
21 b = np.zeros(m + 1)
22
23 A[0, 0] = 1
24 b[0] = ca0
25
26
27 for i in range(1, n + 1):
28 A[i, i] = -2
29 A[i, i - 1] = 1
30 A[i, i + 1] = 1
31
32 for i in range(n + 1, m):
33 A[i, i] = -k * dx ** 2 / Df - 2
34 A[i, i - 1] = 1
35 A[i, i + 1] = 1
36
37 A[m, m] = -1
38 A[m, m - 1] = 1
39 b[n] = 0
40
41
42 u = np.linalg.solve(A, b)
43
44 x = np.arange(0, (L + Lf) + dx, dx)
45
46 plt.plot(x, u, color="#E00650")
47 plt.xlabel("x (cm)")
48 plt.ylabel("Concentraç~ao (mol/L)")
49 plt.xlim(0, L + Lf)
50 plt.grid()
51 plt.show()
45
11 Problema 11
11.1 Enunciado do problema
A scguinte cquaçio diferencial descreve a concentração estacionária de uma substância
que reage com cinética de primeira ordem em um reator de fluxo em pistão axialmente
disperso:
D
d2c
dx2
− U dc
dx
− kc = 0
em que D é o coeficiente de dispersão (m2/h), c é a concentração (mol/L), x é a distância
(m), U é a velocidade (m/h) e k é a taxa de reação (m/h). As condições de contorno
podem ser formuladas como:
Ucin = Uc(x = 0)−D
dc
dx
(x = 0)
dc
dx
(x = L) = 0
em que cin é a concentração do fluxo de entrada (mol/L) e L é o comprimento do reator
(m). Essas são chamadas de condições de contorno de Danckwerts. Use a abordagem de
diferenças finitas para determinar a concentração como uma função da distância, dados
os scguintes parâmetros: D = 5000 m2/h, U = 100m/h, k = 2 /h , L = 100 m, e
cin = 100 mol/L. Utilize aproximaçõ por diferenças centradas com ∆x = 10 m para
obter sua solução.
11.2 SoluçãoA resolução deste problema exige uma abordagem similar a do problema 10. A discre-
tização na forma da Equação 11.1:(
D
∆x2
+
U
2∆x
)
ci+1 −
(
2D
∆x2
+ k
)
ci +
(
D
∆x2
− U
2∆x
)
ci−1 = 0 (11.1)
Em x = L, a condição de contorno discretizadas fica
cn − cn−1 = 0
Em x = 0, a condição de contorno de Danckwerts fornece(
U +
D
∆x
)
c1 − D
∆x
c2 = Ucin
A solução do problema é apresentada na Figura 11.1.
46
0 20 40 60 80 100
30
40
50
60
x (cm)
C
on
ce
n
tr
aç
ão
(m
ol
/L
)
Figura 11.1: Solução para o problema 11 obtido através do método de diferenças finitas.
11.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-08T12:19:15-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T01:53:03-02:00
5
6 import numpy as np
7 import matplotlib.pyplot as plt
8
9 L = 100 # m
10 D = 5000 # m^2 / h
11 U = 100 # mol/h
12 k = 2 # h^-1
13 dx = 10 # cm
14 n = int(L / dx)
15 cin = 100 # mol/L
16
17 A = np.zeros((n + 1, n + 1))
18 b = np.zeros(n + 1)
19
20 A[0, 0] = U + D / dx
21 A[0, 1] = -D / dx
22 b[0] = U * cin
23
24 A[n, n] = 1
25 A[n, n - 1] = -1
26 b[n] = 0
47
27
28 for i in range(1, n):
29 A[i, i] = -(2 * D / dx ** 2 + k)
30 A[i, i - 1] = D / dx ** 2 + U / (2 * dx)
31 A[i, i + 1] = D / dx ** 2 - U / (2 * dx)
32
33
34 u = np.linalg.solve(A, b)
35
36 x = np.arange(0, L + dx, dx)
37
38
39 plt.plot(x, u, color="#E00650")
40 plt.xlabel("x (cm)")
41 plt.ylabel("Concentraç~ao (mol/L)")
42 plt.xlim(0, L)
43 plt.grid()
44 plt.show()
48
12 Problema 12
12.1 Enunciado do problema
O crescimento de bactéria em um reator descont́ınuo utiliza um substrato solúvel. A
captação do substrato é representada por um modelo loǵıstico com limitação de Michaelis-
Menten. A morte de bactérias produz detritos que são subsequentemente convertidos a
substratos por hidrólise. Além disso, as bactérias também excretam alguns substratos
diretamente. Morte, hidrólise e excreção são todos simulados por reações de primeira
ordem.
Os balanços de massa podem ser escritos como:
dX
dt
= µmax
(
1− X
K
)(
S
Ks + S
)
X − kdX − keX
dC
dt
= kdX − khC
dS
dt
= keX + khC − µmax
(
1− X
K
)(
S
Ks + S
)
X
em que X, C e S são as concentrações (mg/L de bactéria, detritos e substrato,
respectivamente; µmax é a taxa de crescimento máxima (1/dia); K é a capacidade de
suporte loǵıstico (mg/L) ; Ks é a constante de meia saturação de Michaclis-Menten
(mg/L) ; kd é a taxa de mortalidade (1/dia); ke é a taxa de excreção (1/dia); e kh é a
taxa de hidrólise (1/d). Simule as concentrações de t = 0 a 100 dias, dadas as condições
iniciais X(0) = 1 mg/L, S(0) = 100 mg/L, e C(0) = 0 mg/L. Utilize os seguintes
parâmetros nos cálculos: µmax = 10 /dia, K = 10 mg/L, Ks = 10 mg/L kd = 0, 1/ dia,
ke = 0, 1 /dia e kh = 0, 1 /dia. Utilize os métodos numéricos de Euler explicito e Runge-
Kutta de quarta ordem, além de um solver de sua preferència.
49
12.2 Solução
12.2.1 Método de Euler
0 20 40 60 80 100
0
20
40
60
80
100
X
C
S
Figura 12.1: Resolução do problema 12 com o método de Euler.
12.2.2 Método de Runge-Kutta
0 20 40 60 80 100
0
20
40
60
80
100
X
C
S
Figura 12.2: Resolução do problema 12 com o método de Runge-Kutta de quarta ordem.
50
12.2.3 Solver
0 20 40 60 80 100
0
20
40
60
80
100
X
C
S
Figura 12.3: Resolução do problema 12 com método de Radau implementado por solver.
12.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-08T16:07:40-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:41:50-02:00
5
6 import numpy as np
7 from numpy import exp
8 import matplotlib.pyplot as plt
9 import RK_4
10 import euler_explicito
11 import scipy.integrate
12
13 cX0 = 1 # mg/L
14 cC0 = 0 # mg/L
15 cS0 = 100 # mg/L
16 c0 = [cX0, cC0, cS0]
17 dt = 0.1
18
19 umax = 10 # dia^-1
20 K = 10 # mg/L
21 ks = 10 # mg/L
22 kd = 0.1 # dia^-1
23 ke = 0.1 # dia^-1
24 kh = 0.1 # dia ^-1
25
51
26
27 def solve(t, c):
28 X, C, S = c
29 dXdt = umax * (1 - X / K) * (S / (ks + S)) * X - kd * X - ke * X
30 dCdt = kd * X - kh * C
31 dSdt = ke * X + kh * C - umax * (1 - X / K) * (S / (ks + S)) * X
32 return dXdt, dCdt, dSdt
33
34
35 t, y = euler_explicito.Euler_SE(0, 100, c0, dt, solve)
36 plt.plot(t, y[0], color="#E00650", label="X")
37 plt.plot(t, y[1], color="#2D2E7B", label="C")
38 plt.plot(t, y[2], color="#E8B038", label="S")
39 plt.legend()
40 plt.title('Euler')
41 plt.show()
42
43 t, y = RK_4.RK_4_SE(0, 100, c0, dt, solve)
44
45 plt.plot(t, y[0], color="#E00650", label="X")
46 plt.plot(t, y[1], color="#2D2E7B", label="C")
47 plt.plot(t, y[2], color="#E8B038", label="S")
48 plt.legend()
49 plt.title('Runge-Kutta')
50 plt.show()
51
52 sol = scipy.integrate.solve_ivp(solve, (0, 100),
53 c0,
54 method="Radau",
55 max_step=dt)
56
57 plt.plot(sol.t, sol.y[0], color="#E00650", label="X")
58 plt.plot(sol.t, sol.y[1], color="#2D2E7B", label="C")
59 plt.plot(sol.t, sol.y[2], color="#E8B038", label="S")
60 plt.legend()
61 plt.title('Solver')
62 plt.show()
52
13 Problema 13
13.1 Enunciado do problema
A equação de advecção-difusão é usada para calcular a distribuição de concentração ao
longo do comprimento de um biorreator retangular, conforme mostra a figura abaixo:
∂C
∂t
= D
∂2C
∂x2
− U ∂C
∂x
− kc
onde C é a concentração (mg/m3), t é o tempo (min), D é o coeficiente de difusão
m2/min, x é a distância ao longo do eixo longitudinal do tanque (m), em que x = 0
na entrada do tanque, U é a velocidade na direção x (m/min) e k é a taxa de reação
(min−1) pela qual o componente A alimentado ao biorreator é consumido.
Desenvolva um esquema para resolver essa cquação numericamente. Teste o seu es-
quema para k = 0, 15, D = 100 e U = 1 para um reator com comprimento de 10 m
operando em regime cont́ınuo. Use ∆x = 1 e um passo de ∆t = 0, 005 Suponha que a
concentração de entrada seja 100 e que a concentração inicial no bioreator seja zero. Faça
a simulação para o intervalo de t = 0 até 100 min, faça os gráficos necessários (ao longo
de x ) e avalie os resultados obtidos.
13.2 Solução
É importante destacar que ao implementar-se o problema com a discretização sugerida,
os erros ficaram muito maiores que os toleráveis. Desta forma, como é posśıvel reparar
no código associado, foi utilizada uma discretização com intervalos menores de modo a
refinar a solução. Infelizmente, não pode melhorar-se a discretização no espaço para que
não houvesse a violação das condições de estabilidade de Von-Neumann. A discretização
ficou na forma da Equação 13.1, sendo i o ı́ndice que diz respeito ao espaço e j o ı́ndice
que diz respeito ao tempo.
cj+1i − c
j
i
∆t
= D
(
cji+1 − 2c
j
i + c
j
i−1
∆x2
)
−
cj+1i+1 − c
j
i−1
∆x
− kcji (13.1)
A Figura 13.1 mostra a solução do problema após 100 minutos. Para tempos menores,
a solução é mostrada na Figura 13.2. Através desta figura, pode-se notar que a evolução
da concentração com o tempo ocorre de forma abrupta e após os minutos iniciais a solução
do problema pouco varia com o tempo.
53
0 2 4 6 8 10
94
95
96
97
98
99
100
x (m)
c
(k
g
/m
3
)
Após 100 min
Figura 13.1: Solução após 100 minutos para o problema 13.
x (cm)
0
2
4
6
8
10
Te
mp
o (
mi
n)
0
20
40
60
80
100
Co
nc
en
tra
çã
o 
(m
g/
m
3 )
0
20
40
60
80
100
Figura 13.2: Solução transiente para o problema 13.
54
13.2.1 Efeito de cada termo na equação
O pacote Numpy de Python calcula de forma rápida a transformada rápida de Fourier.
Isto é interessante pois permite de maneira fácil analisar o peso de cada parâmetro da
equação diferencial parcial. A Figura 13.3 mostra o efeito do termo difusivo, a Figura 13.4
o efeito do termo avectivo e a Figura 13.5 o efeito do termo reativo da equação. Os efeitos
compostos são mostrados na Figura 13.6. Estas figuras não representamo problema em
questão por terem condições inciais e contornos diferentes. Sua riqueza está, no entanto,
nos insights que elas fornecem acerca de cada parâmetro na equação.
0 200 400 600 800
0
20
40
60
80
Figura 13.3: D = 100. U = 0. k = 0.
0 200 400 600 800
0
20
40
60
80
Figura 13.4: D = 0. U = 1. k = 0.
55
0 200 400 600 800
0
20
40
60
80
Figura 13.5: D = 0. U = 0. k = 0,15.
0 200 400 600 800
0
20
40
60
80
Figura 13.6: D = 100. U = 1, k = 0,15.
13.3 Código associado
13.3.1 Problema principal
1 # @Author: eduardo
2 # @Date: 2020-11-18T00:38:26-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:42:37-02:00
5
56
6 import numpy as np
7 import matplotlib.pyplot as plt
8
9 # Entrar valores dos parâmetros e constantes
10 D = 100
11 U = 1
12 k = 0.15
13 dt = 0.0005 # Discretizaç~ao no tempo
14 dx = 0.5 # Discretizaç~ao no espaço
15 c0 = 100
16 tf = 100 # min
17 xf = 10 # m
18
19 x = np.arange(0, xf + dx, dx)
20
21 c = np.zeros_like(x)
22 c[0] = c0
23 c_novo = c.copy()
24
25 C = []
26 for j in range(int(tf / dt) + 1):
27 c_novo[0] = c0
28 # Iterar sobre os elementos no interior da malha
29 c_novo[1:-1] = dt * ((D / (dx**2) * (c[2:] - 2 * c[1:-1] + c[:-2]) - U /
30 (2 * dx) * (c[2:] - c[:-2]) - k * c[1:-1])) + c[1:-1]
31 # Impor condiç~oes de contorno dc/dt|c=100 = 0
32 c_novo[-1] = c_novo[-2]
33 c = c_novo.copy()
34 C.append(c)
35
36 C = np.array(C)
37
38 plt.plot(x, c, label="Após 100 min", color="#E00650")
39 plt.xlabel(r"$x\ (m)$")
40 plt.ylabel(r"$c\ (kg/m^3)$")
41 plt.legend()
42 plt.grid()
43 plt.xlim(0, 10)
44 plt.show()
45
46 dplot = 10000
47 u_plot = C[0:-1:dplot, :]
48
49 ysplot = []
50 for j in range(u_plot.shape[0]):
51 ys = dplot * np.ones(u_plot.shape[1]) * dt * j
52 ysplot.append(ys)
53
57
54 fig = plt.figure()
55 ax = fig.gca(projection='3d')
56 surf = ax.plot_wireframe(x, ysplot, u_plot, color='black')
57 ax.set_zlabel('Concentraç~ao ($mg/m^3$)')
58 ax.set_zlim(0, 100)
59 ax.set_xlabel('x (cm)')
60 ax.set_xlim(0, 10)
61 ax.set_ylabel('Tempo (min)')
62 ax.set_ylim(0, 100)
63 plt.show()
64
65
66 ax = plt.axes(projection='3d')
67 ax.plot_surface(x, ysplot, u_plot, rstride=1, cstride=1,
68 cmap='viridis', edgecolor='none')
69 plt.show()
13.3.2 Fast Fourier Transform (FFT)
1 # @Author: eduardo
2 # @Date: 2020-11-12T22:23:18-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T03:34:54-02:00
5
6 import numpy as np
7 import matplotlib.pyplot as plt
8 from scipy.integrate import odeint
9
10 L = 100 # Length of domain
11 dx = 0.1
12 N = int(L / dx) # Number of discrete points
13 x = np.arange(0, L, dx)
14
15 # Define discrete wavenumbers
16 kappa = 2 * np.pi * np.fft.fftfreq(N, d=dx)
17
18 # Initial condition
19 u0 = np.zeros_like(x)
20 u0[int((L / 2 - L / 5) / dx):int((L / 2 + L / 5) / dx)] = 1
21 print(u0)
22 dt = 0.05
23 t = np.arange(0, 5, dt)
24
25 D, U, k = 100, 1, 0.15
26
27
28 def rhsConvecAdvec(u, t, kappa):
58
29 uhat = np.fft.fft(u)
30 d_uhat = (1j) * kappa * uhat
31 dd_uhat = -np.power(kappa, 2) * uhat
32 d_u = np.fft.ifft(d_uhat)
33 dd_u = np.fft.ifft(dd_uhat)
34 du_dt = D * dd_u - U * d_u - k * u
35 return du_dt.real
36
37
38 u = odeint(rhsConvecAdvec, u0, t, args=(kappa, ))
39
40 plt.figure()
41 plt.imshow(np.flipud(u), aspect=8)
42 plt.set_cmap("magma")
43 plt.show()
59
14 Problema 14
Luss e Amundson (Stability of batch catalytic fluidized beds, AIChE J., 14 (2), 211,1968)
estudaram um modelo simplificado para a dinâmica de um reator cataĺıtico de leito flui-
dizado no qual uma reação irreverśıvel em fase gasosa A → B ocorre. As equações de
conversão obtidas a partir dos balanços de massa e energia, bem como a constante cinética
para esta reação, são apresentadas abaixo (chamaremos este sistema de equações como
conjunto 1)
dP
dτ
= Pe − P +Hg (Pp − P )
dT
dτ
= Te − T +HT (Tp − T ) +HW (TW − T )
dPp
dτ
=
Hg
A
[P − Pp(1 +K)]
dTp
dτ
=
HT
C
[(T − Tp) + FKPp]
K = 0, 0006 exp
(
20, 7− 15000
Tp
)
onde T é a temperatura na fase fluido (◦R), P é a pressão parcial do reagente na
fase fluido (atm), Tp é a temperatura na superf́ıcie do catalisador (
◦R), Pp é a pressão
parcial do reagente na superficie do catalisador (atm), K é a constante de velocidade
de reação (adimensional), τ é o tempo adimensional, e o subscrito e indica condições
na alimentação. As constantes adimensionais são dadas como Hg = 320, Te = 600,
HT = 266, 67, HW = 1, 6, TW = 720, F = 8000, A = 0, 17142, C = 205, 74 e Pe = 0, 1.
Aiken e Lapidus (An effective integration method for typical stiff systems, AIChE J.,
20 (2), 368,1974) susequentemente usaram o conjunto 1 como teste para um programa que
desenvolveram, mas reescreveram as equações introduzindo valores numéricos no sistema
e arredondando alguns coeficientes. Os autores obtiveram, então, o sistema de equações
abaixo, denominado aqui conjunto 2
dP
dτ
= 0, 1 + 320Pp − 321P
dT
dτ
= 1752− 269T + 267Tp
dPp
dτ
= 1, 88× 103 [P − Pp(1 +K)]
dTp
dτ
= 1, 3 (T − Tp) + 1, 04× 104KPp
K = 0, 0006 exp
(
20, 7− 15000
Tp
)
(a) Introduza os valores apresentados por Luss e Amundson no conjunto 1 e observe as
diferenças entre o conjunto 1 e o conjunto 2 .
(b) Encontre todas as soluções em regime permanente tanto para o conjunto 1 quanto
para o conjunto 2 no intervalo (500 ≤ T ≤ 1300)◦ R.
60
(c) Resolva tanto o conjunto 1 quanto o conjunto 2 utilizando como condições iniciais
P = 0, 1, T = 600, Pp = 0 e Tp = 761. Integre no intervalo 0 ≤ τ ≤ 1500.
(d) Explique as diferenças entre as soluções dinâmicas e em regime estacionário obtidas
para os conjuntos 1 e 2 .
14.1 Solução
14.2 Item a
As equações com os parâmetros de Luss e Amundson são:
dP
dτ
= 0, 1− P + 320 (Pp − P )
dT
dτ
= 600− T + 266, 67 (Tp − T ) + 1, 6 (720− T )
dPp
dτ
=
320
0, 17142
[P − Pp(1 +K)]
dTp
dτ
=
266, 67
205, 74
[(T − Tp) + 8000KPp]
K = 0, 0006 exp
(
20, 7− 15000
Tp
)
Nota-se que as diferenças são pontuais e, percentualmente, moderadas.
14.2.1 Solução estacionária
Para o sistema de equações I, encontrar a solução estacionária para o problema consiste
em encontrar todas as soluções estacionárias para
0 = Pe − P +Hg (Pp − P )
0 = Te − T +HT (Tp − T ) +HW (TW − T )
0 = P − Pp(1 +K)
0 = T − Tp + FKPp
K = 0, 0006 exp
(
20, 7− 15000
Tp
)
dentro do intervalo pedido. Observe que, ao utilizar-se a segunda equação, pode-se
obter o valor de Tp somente em função de T . Ao utilizar-se a quinta equação, pode-se
obter o valor de K somente em função de Tp, já conhecida. Com a quarta equação, é
posśıvel obter o valor de Pp somente em função de valores conhecidos e finalmente, com
a primeira equação, é posśıvel computar P . Para cada T , estes valores são solução para
o sistema não linear se a terceira equação resultar em zero.
Desta forma, um posśıvel algoŕıtmo para calcular os pontos estacionários do problema
seria:
61
1. Com o primeiro T do intervalo, calcular Tp com
Tp = −
Te − T −HtT +Hw ∗ (Tw − T )
Ht
;
2. Com T e TP , calcular K com
K = 0, 0006 exp
(
20, 7− 15000
Tp
)
;
3. Com T , TP e K, calcular PP com
PP =
T − Tp
FK
;
4. Com os valores já conhecidos, calcular P com
P = PP (1 +K);
5. Calcular o valor da função auxiliar:
f = Pe − P +Hg(PP − P );
6. Repetir para todos os Ts do intervalo.
Sabe-se que quando a função auxiliar for zero, passa-se por um ponto que satisfaz o
sistema de equações para o estado estacionário.
A Figura 14.1 mostra o valor da função auxiliar versus a temperatura.
700 750 800 850 900 950 1,000
−0.1
−5 · 10−2
0
5 · 10−2
0.1
Temperatura (ºR)
F
u
n
çã
o
A
u
x
il
ia
r
Figura 14.1: Valor da função auxiliar para o cálculo do estado estacionário na questão
14.
62
Valores menores que 675 e maiores de 1000 estão fora do gráfico pois a função não
muda de sinal neste intervalo. Assim, percebe-seque existem três estados estacionários.
Utilizando os pontos que fazem a função mudar de sinal para refinar a solução, obtém-se:
Estado estacionário I - Conjunto I
T (R) 690,45
TP (R) 690,61
k 0,0002164
P (atm) 0,09353
PP (atm) 0,09351
Tabela 14.1: Valores obtidos para o primeiro estado estacionário.
Estado estacionário II - Conjunto I
T (R) 758,34
TP (R) 759,17
P (atm) 0,06705
PP (atm) 0,06694
k 0,001538
Tabela 14.2: Valores obtidos para o segundo estado estacionário.
Estado estacionário III - Conjunto I
T (R) 912,77
TP (R) 915,10
P (atm) 0,006822
PP (atm) 0,0065304
k 0,004459
Tabela 14.3: Valores obtidos para o terceiro estado estacionário.
Para o conjunto II de equações, o procedimento é exatamente o mesmo. Todavia, para
este sistema de equações, encontrou-se somente uma solução estacionária, apresentada na
Tabela 14.4:
Estado estacionário - Conjunto II
T (R) 1208,28
TP (R) 1210,78
PP (atm) 0,0001268
P (atm) 0,0004386
k 2,442
Tabela 14.4: Valores obtidos para o estado estacionário do conjunto de equações.
14.2.2 Solução transiente
A Figura 14.2 apresenta as soluções transientes para a pressão do conjunto I. A Figura 14.3
apresenta a solução para as temperaturas.
63
0 200 400 600 800 1,000 1,200 1,400
0
2 · 10−2
4 · 10−2
6 · 10−2
8 · 10−2
0.1
τ (s)
P
re
ss
ão
(a
tm
)
P
Pp
Figura 14.2: Solução transiente para as pressões do conjunto I.
0 200 400 600 800 1,000 1,200 1,400
600
650
700
750
800
850
900
τ (s)
T
em
p
er
at
u
ra
(º
R
)
T
Tp
Figura 14.3: Solução transiente para as temperaturas do conjunto I.
Para o conjunto II, as soluções das pressões e temperaturas são apresentadas na Fi-
gura 14.4 e na Figura 14.5, respectivamente.
64
0 200 400 600 800 1,000 1,200 1,400
0
2 · 10−2
4 · 10−2
6 · 10−2
8 · 10−2
0.1
τ (s)
T
em
p
er
at
u
ra
(º
R
)
P
Pp
Figura 14.4: transiente para as pressões do conjunto II
0 200 400 600 800 1,000 1,200 1,400
600
700
800
900
1,000
1,100
1,200
τ (s)
T
em
p
er
at
u
ra
(º
R
)
T
Tp
Figura 14.5: Solução para as temperaturas do conjunto II.
14.2.3 Item d
A diferença entre os estados estacionários se deve às mudanças nos parâmetros das
equações. Isto evidencia que, mesmo com um modelo exatamente igual em termos de
estrutura das equações e partindo do mesmo ponto (condições iniciais iguais), pode-se
obter variações significativas nos resultados de um modelo matemático no fenômeno f́ısico,
caso este apresente sensibilidade aos parâmetros. Esta discussão será mais aprofundada
no problema 18.
65
14.3 Código associado
14.3.1 Solução transiente - Sistema I
1 # @Author: eduardo
2 # @Date: 2020-11-11T23:25:48-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:44:10-02:00
5
6 import numpy as np
7 from numpy import exp
8 import matplotlib.pyplot as plt
9 import RK_4
10 import euler_explicito
11 import scipy.integrate
12
13 Hg = 320
14 Te = 600
15 Ht = 266.67
16 Hw = 1.6
17 Tw = 720
18 F = 8000
19 A = 0.17142
20 C = 205.74
21 Pe = 0.1
22
23
24 def solve(tau, c):
25 P, T, Pp, Tp = c
26 K = 0.0006 * exp(20.7 - 15000 / Tp)
27 dPdtau = Pe - P + Hg * (Pp - P)
28 dTdtau = Te - T + Ht * (Tp - T) + Hw * (Tw - T)
29 dPpdtau = Hg / A * (P - Pp * (1 + K))
30 dTpdtau = Ht / C * ((T - Tp) + F * K * Pp)
31 return dPdtau, dTdtau, dPpdtau, dTpdtau
32
33
34 c0 = [0.1, 600, 0, 761]
35
36 # t, c = euler_explicito.Euler_SE(0, 100, c0, dt, solve)
37
38 dt = 1
39 sol = scipy.integrate.solve_ivp(solve, (0, 1500),
40 c0,
41 method="Radau",
42 max_step=dt)
43
44 plt.plot(sol.t, sol.y[1], color="#2D2E7B", label="T")
66
45 plt.plot(sol.t, sol.y[3], label="$T_p$")
46 plt.xlabel(r'$ \tau $')
47 plt.ylabel("Temperatura (ºR)")
48 plt.grid()
49 plt.legend()
50 plt.show()
51
52 plt.plot(sol.t, sol.y[0], color="#E00650", label="P")
53 plt.plot(sol.t, sol.y[2], color="#E8B038", label="$P_p$")
54 plt.xlabel(r'$ \tau $')
55 plt.ylabel('Press~ao (atm)')
56 plt.grid()
57 plt.legend()
58 plt.show()
14.3.2 Solução transiente - Sistema I
1 # @Author: eduardo
2 # @Date: 2020-11-12T11:45:43-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:44:24-02:00
5 import numpy as np
6 from numpy import exp
7 import matplotlib.pyplot as plt
8 import RK_4
9 import euler_explicito
10 import scipy.integrate
11
12
13 def solve(tau, c):
14 P, T, Pp, Tp = c
15 K = 0.0006 * exp(20.7 - 15000 / Tp)
16 dPdtau = 0.1 + 320 * Pp - 321 * P
17 dTdtau = 1752 - 269 * T + 267 * Tp
18 dPpdtau = 188e3 * (P - Pp * (1 + K))
19 dTpdtau = 1.3 * (T - Tp) + 1.04e4 * K * Pp
20 return dPdtau, dTdtau, dPpdtau, dTpdtau
21
22
23 c0 = [0.1, 600, 0, 761]
24
25 sol = scipy.integrate.solve_ivp(solve, (0, 1500), c0, method="Radau")
26
27
28 plt.plot(sol.t, sol.y[1], color="#2D2E7B", label="T")
29 plt.plot(sol.t, sol.y[3], label="$T_p$")
30 plt.xlabel(r'$ \tau $')
67
31 plt.ylabel("Temperatura (ºR)")
32 plt.grid()
33 plt.legend()
34 plt.show()
35
36 plt.plot(sol.t, sol.y[0], color="#E00650", label="P")
37 plt.plot(sol.t, sol.y[2], color="#E8B038", label="$P_p$")
38 plt.xlabel(r'$ \tau $')
39 plt.ylabel('Press~ao (atm)')
40 plt.grid()
41 plt.legend()
42 plt.show()
14.3.3 Solução estacionária - Sistema I
1 # @Author: eduardo
2 # @Date: 2020-11-20T16:41:42-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-20T20:28:26-02:00
5
6 import numpy as np
7 from numpy import exp
8 from scipy.optimize import fsolve
9 import matplotlib.pyplot as plt
10
11 Hg = 320
12 Te = 600
13 Ht = 266.67
14 Hw = 1.6
15 Tw = 720
16 F = 8000
17 A = 0.17142
18 C = 205.74
19 Pe = 0.1
20
21
22 def f1(T):
23 # Tp
24 return -(Te - T - Ht * T + Hw * (Tw - T)) / Ht
25
26
27 def f2(Tp):
28 # K
29 return 0.0006 * exp(20.7 - 15000 / Tp)
30
31
32 def f3(K, Tp, T):
68
33 # Pp
34 return -(T - Tp) / (F * K)
35
36
37 def f4(Pp, K):
38 # P
39 return Pp * (1 + K)
40
41
42 def f5(Pp, P):
43 return Pe - P + Hg * (Pp - P)
44
45
46 T_guess = np.arange(500, 1301, 0.01)
47
48 ee = []
49 sol = []
50 for guess in T_guess:
51 add = False
52 Tp = f1(guess)
53 K = f2(Tp)
54 Pp = f3(K, Tp, guess)
55 P = f4(Pp, K)
56 sol.append(f5(Pp, P))
57 if len(sol) > 2:
58 if sol[-1] < 0 and sol[-2] > 0:
59 ee.append([guess, Tp, Pp, P, K])
60 if sol[-1] > 0 and sol[-2] < 0:
61 ee.append([guess, Tp, Pp, P, K])
62
63 plt.plot(T_guess, sol)
64 plt.xlabel('Temperatura (ºR)')
65 plt.ylabel('Funç~ao Auxiliar')
66 plt.ylim(-0.1, .1)
67 plt.xlim(675, 1000)
68 plt.grid()
69 plt.show()
70
71
72 def solve(x):
73 T, Tp, Pp, P, K = x
74 f1 = K - 0.0006 * exp(20.7 - 15000 / Tp)
75 f2 = Pe - P + Hg * (Pp - P)
76 f3 = Te - T + Ht * (Tp - T) + Hw * (Tw - T)
77 f4 = Hg / A * (P - Pp * (1 + K))
78 f5 = Ht / C * ((T - Tp) + F * K * Pp)
79 return f1, f2, f3, f4, f5
80
69
81
82 print(fsolve(solve, ee[0], factor=0.1))
83 print(fsolve(solve, ee[1]))
84 print(fsolve(solve, ee[2]))
14.3.4 Solução estacionária - Sistema II
1 # @Author: eduardo
2 # @Date: 2020-11-20T20:13:03-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:15:41-02:00
5
6 import numpy as np
7 from numpy import exp
8 from scipy.optimize import fsolve
9 import matplotlib.pyplot as plt
10
11 Hg = 320
12 Te = 600
13 Ht = 266.67
14 Hw = 1.6
15 Tw = 720
16 F = 8000
17 A = 0.17142
18 C = 205.74
19 Pe = 0.1
20
21
22 def f1(T):
23 # Tp
24 return (1752 - 269 * T) / -267
25
26
27 def f2(Tp):
28 # K
29 return 0.0006 * exp(20.7 - 15000 / Tp)
30
31
32 def f3(K, Tp, T):
33 # Pp
34 return -1.3 * (T - Tp) / (1.04e4 * K)
35
36
37 def f4(Pp, K):
38 # P
39 return Pp * (1 + K)
40
70
41
42 def f5(Pp, P):
43 return 0.1 + 320 * Pp - 321 * P
44
45
46 T_guess = np.arange(500, 1301, 0.01)
47
48 ee = []
49 sol = []
50 for guess in T_guess:
51 add = False
52 Tp = f1(guess)
53 K = f2(Tp)
54 Pp = f3(K, Tp, guess)
55 P = f4(Pp, K)
56 sol.append(f5(Pp,P))
57 if len(sol) > 2:
58 if sol[-1] < 0 and sol[-2] > 0:
59 ee.append([guess, Tp, Pp, P, K])
60 if sol[-1] > 0 and sol[-2] < 0:
61 ee.append([guess, Tp, Pp, P, K])
62 print(ee)
63 plt.plot(T_guess, sol)
64 plt.xlabel('Temperatura (ºR)')
65 plt.ylabel('Funç~ao Auxiliar')
66 plt.ylim(-0.1, .1)
67 # plt.xlim(675, 1000)
68 plt.grid()
69 plt.show()
70
71
72 def solve(x):
73 T, Tp, Pp, P, K = x
74 f1 = K - 0.0006 * exp(20.7 - 15000 / Tp)
75 f2 = 0.1 + 320 * Pp - 321 * P
76 f3 = 1752 - 269 * T + 267 * Tp
77 f4 = P - Pp * (1 + K)
78 f5 = 1.3 * ((T - Tp) + 1.04e4 * K * Pp)
79 return f1, f2, f3, f4, f5
80
81
82 print(fsolve(solve, ee[0], factor=10))
71
15 Problema 15
15.1 Enunciado do problema
Beltrami (2002) e Gray (2006) postularam o seguinte problema para representar a dinâmica
de duas populações de bactérias competitivas. As variáveis são as densidades de po-
pulação de cada uma das espécies. O termo entre parênteses representa o crescimento
devido às limitações do ambiente e o último termo representa os efeitos negativos da
competição entre as espécies. O primeiro termo representa a densidade de população se
apenas uma das espécies estivesse presente. O segundo termo é proporcional ao número
de interações entre as espécies e representa o decĺınio da densidade de população devido
à competição:
dx1
dt
= 9x1
(
1− x1
9
)
− 2x1x2
dx2
dt
= 6x2
(
1− x2
12
)
− x1x2
Utilize o método de Euler expĺıcito para encontrar uma solução para esse sistema de
EDOs quando as seguintes densidades de população estão presentes inicialmente em um
certo meio: (x1, 0;x2, 0) = (3; 1), (7; 4), (5; 2), (5, 01; 2). Apresente os resultados na forma
de gráficos e faça comentários sobre o estado estacionário em cada caso.
15.2 Solução
Este problema apresenta uma grande quantidade de estados permanentes, que variam de
acordo com o input inicial.
15.2.1 Quantidade inicial: (3 ; 1)
0 1 2 3 4 5 6
0
2
4
6
8
x1
x2
Figura 15.1: Solução para quantidade inicial de (3 ; 1) para o problema 15.
Neste caso, o modelo convergiu rapidamente para o ponto estacionário. x2 apresentou
um leve aumento antes de descer.
72
15.2.2 Quantidade inicial: (7 ; 4)
0 1 2 3 4 5 6
0
2
4
6
8
10
12 x1
x2
Figura 15.2: Solução para quantidade inicial de (7 ; 4) para o problema 15.
O modelo convergiu para o ponto estacionário mais rápido do que no item anterior.
Todavia, antes da convergência, a concentração de bactérias apresentou uma dinâmica
bastante interessantes, com fortes variações no intervalo (0,1).
15.2.3 Quantidade inicial: (5 ; 2)
0 1 2 3 4 5 6
2
2.5
3
3.5
4
4.5
5
x1
x2
Figura 15.3: Solução para quantidade inicial de (5 ; 2) para o problema 15.
A estimativa inicial foi um estado estacionário.
73
15.2.4 Quantidade inicial: (5,01 ; 2)
0 1 2 3 4 5 6
0
2
4
6
8
x1
x2
Figura 15.4: Solução para quantidade inicial de (5,01 ; 2) para o problema 15.
A convergência ocorreu de forma sutil e monótona para o estado permanente.
15.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-08T21:39:02-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T01:54:23-02:00
5
6 import euler_explicito
7 import matplotlib.pyplot as plt
8
9
10 def solve(t, x):
11 x1, x2 = x
12 dx1dt = 9 * x1 * (1 - x1 / 9) - 2 * x1 * x2
13 dx2dt = 6 * x2 * (1 - x2 / 12) - x1 * x2
14 return dx1dt, dx2dt
15
16
17 ci = [[3, 1], [7, 4], [5, 2], [5.01, 2]]
18
19 dx = 0.1
20 for i, c0 in enumerate(ci):
21 ti, xi = euler_explicito.Euler_SE(0, 6, c0, dx, solve)
22 plt.plot(ti, xi[0], color="#E00650", label='x1')
23 plt.plot(ti, xi[1], color="#E8B038", label='x2')
74
24 plt.grid()
25 plt.legend()
26 plt.show()
75
16 Problema 16
16.1 Enunciado do problema
Dependendo dos valores das constantes de velocidade, os sistemas reativos podem ser
stiff. Um problema stiff bem conhecido é o exemplo de Robertson, que consiste em EDOs
não lineares caracterizadas por uma grande diferença entre os valores das constantes de
velocidade. As três reações que ocorrem nesse sistema são:
A
k1→ B
B +B
k2→ B + C
B + C
k3→ A+ C
onde as constantes de velocidade são kI = 0, 04, k2 = 3 · 107, k3 = 1 · 104e as
concentrações iniciais das três espécies são CA0 = 1, CB0 = 0 e CCO = 0
(a) Assuma que as reações são conduzidas em um reator batelada operado em condições
isotérmicas e derive o sistema de EDOs.
(b) Utilize um solver adequado para problemas stiff e plote a solução para t ∈ [0, 100].
(c) Avalie o desempenho de solvers não-stiff e stiff (MATLAB ode45 e ode15s, por
exemplo) analisando o número de passos necessários.
16.2 Solução
16.2.1 Item a
Considerando que todas as reações são elementares, o sistema de equações diferenciais
ordinárias fica:
dca
dt
= −k1ca + k3cbcc (16.1)
dcb
dt
= k1ca − k2c2b +
1
2
k2c
2
b − k3cbcc (16.2)
dcc
dt
=
1
2
k2c
2
b (16.3)
16.2.2 Items b e c
A Figura 16.1 apresenta a solução para o problema.
76
0 20 40 60 80 100
0
0.2
0.4
0.6
0.8
1
Tempo
C
on
ce
n
tr
aç
ão
Ca (stiff)
Cb (stiff)
Cc (stiff)
Ca (não-stiff)
Cb (não-stiff)
Cc (não-stiff)
Figura 16.1: Solução stiff e não-stiff para o problema 16.
Note como as soluções praticamente se sobrepõe. Todavia, para chegar a este resultado
foi necessário uma discretização muito menor para o método não-stiff. Para o método
de Radau, que consiste em um método de Runge-Kutta impĺıcito, foram necessários 21
pontos na discretização. Para o solver que implementa o método exṕıcito de Runge-
Kutta de 4ª ordem, foram necessários 100002 pontos na discretização para a solução não
divergir.
16.3 Código associado
1 # @Author: eduardo
2 # @Date: 2020-11-12T18:50:09-02:00
3 # @Last modified by: eduardo
4 # @Last modified time: 2020-11-21T02:46:26-02:00
5
6 import numpy as np
7 from numpy import exp
8 import matplotlib.pyplot as plt
9 import scipy.integrate
10
11 k1 = 0.04
12 k2 = 3e7
13 k3 = 1e4
14
15 Ca0 = 1
16 Cb0 = 0
17 Cc0 = 0
18 C0 = [Ca0, Cb0, Cc0]
19
77
20
21 def solve(t, c):
22 Ca, Cb, Cc = c
23 dCadt = -k1 * Ca + k3 * Cb * Cc
24 dCbdt = k1 * Ca - 0.5 * k2 * Cb**2 - k3 * Cb * Cc
25 dCcdt = 0.5 * k2 * Cb**2
26 return dCadt, dCbdt, dCcdt
27
28
29 solstiff = scipy.integrate.solve_ivp(solve, (0, 100), C0, method="Radau")
30 solnstiff = scipy.integrate.solve_ivp(solve, (0, 100),
31 C0,
32 method="RK45",
33 max_step=0.001)
34
35 plt.plot(solstiff.t, solstiff.y[0], color="#E00650", label="$C_a$ (stiff)")
36 plt.plot(solstiff.t, solstiff.y[1], color="#2D2E7B", label="$C_b$ (stiff)")
37 plt.plot(solstiff.t, solstiff.y[2], color="#E8B038", label="$C_c$ (stiff)")
38 plt.xlabel("Tempo")
39 plt.ylabel("Concentraç~ao")
40 plt.grid()
41 plt.legend()
42 plt.xlim(0, 100)
43 plt.ylim(0, 1)
44 plt.show()
45
46
47
48 plt.plot(solnstiff.t,
49 solnstiff.y[0],
50 linestyle='--',
51 label="$C_a$ (n~ao-stiff)")
52 plt.plot(solnstiff.t,
53 solnstiff.y[1],
54 linestyle='--',
55 label="$C_b$ (n~ao-stiff)")
56 plt.plot(solnstiff.t,
57 solnstiff.y[2],
58 linestyle='--',
59 label="$C_c$ (n~ao-stiff)")
60 plt.xlabel("Tempo")
61 plt.ylabel("Concentraç~ao")
62 plt.grid()
63 plt.legend()
64 plt.xlim(0, 100)
65 plt.ylim(0, 1)
66 plt.show()
78
17 Problema 17
17.1 Enunciado do problema
O seguinte sistema de equações representa um problema de reação com transferência de
calor em um pellet de catalisador poroso:
∂c
∂t
=
1
r2
∂
∂r
(
r2
∂c
∂r
)
− φ2 ·R(c, T ), ∂c
∂r
(0, t) = 0, c(1, t) = 1
Le
∂T
∂t
=
1
r2
∂
∂r
(
r2
∂T
∂r
)
+ β · φ2 ·R(c, T ), ∂T
∂r
(0, t) = 0, T (1, t) = 1
R(c, T ) = c · exp
[
γ
(
1− 1
T
)]
Resolva este problema utilizando o método das diferenças finitas considerando o se-
guinte conjunto de parâmetros: φ = 1, 1, γ = 30, β = 0, 15 e Le = 1050. Apresente o
resultado graficamente.
17.2 Solução
Para discretizar as equações é necessário, inicialmente, aplicar a regra da cadeia.