Análise e Simulação de processos Químicos (ASPQ-2012)

3_CSTR.sce

// Exercício (Aula 2012-11-05)
//
// Programa para calcular a concentração do reagente para a reação de ordem
// alpha (r = kCa^alpha) A -> B em n reatores CSTR em série.
//
// @date     2012-11-05
// @author   Christian Cândido <[email protected]>
// @license  GPL-v3 <http://www.gnu.org/licenses/gpl-3.0.txt>
// @encoding UTF-8

clear;
getd;

// Entradas
V   = [1.0 2.0 4.0 4.0 3.0]; // Volumes dos reatores
ka  = 0.1; // Constante cinética
F0  = 1.0; // Fluxo inicial
CA0 = 1.0; // Concentração de A na entrada
alpha = 1.5; // Ordem da reação

r_no = 3; // número de reatores

// Chute inicial
C = [CA0; 0.8; 0.7; 0.6; 0.5; 0.4];
C_ant = C + 1;

//Critérios de convergência
iter = 0;
iter_max = 15;
h   = 1.0e-5; // Delta para a derivada
tol = 1.0e-5; // Tolerancia no resultado

// Newton-Raphson: xk+1 = xk - J^-1 F(xk)
while norm( C_ant - C ) > tol & iter < iter_max do
         
    C_ant = C;
    
    for i = 1:r_no
        F(i) = fi(i, C, V, F0, ka);
    end

    for k = 1:r_no // Colunas
        for j = 1:r_no // Linhas
            
            C_pert = C;
            C_pert(k+1) = C_pert(k+1) + h;
            
            J(j, k) = ( fi(j, C_pert, V, F0, ka) - F(j) ) / h;
            
        end
    end
   
    if( norm(F) > tol) then
        C(2:r_no+1) = C(2:r_no+1) - inv(J) * F;
    else
        break;
    end;
    
    iter = iter + 1;
       
end;

printf("CA = ");
disp(C(2:r_no+1));

Xa = (C(1) - C(r_no+1)) / C(1);

printf("\nXA = %f", Xa);

printf("\n\nEm %d iterações", iter);

CSTR_serie.sce

// Exercício (Aula 2012-10-22)
//
// 4 Reatores CSTR em série com corrente de reciclo
// Reação: A -> B, r = kCa
//
// Calcula a concentração de A no interior de cada tanque
//
// @date     2012-10-22
// @author   Christian Cândido <[email protected]>
// @license  GPL-v3 <http://www.gnu.org/licenses/gpl-3.0.txt>
// @encoding UTF-8

// Limpa variaveis
clear()

// Dados de entrada (unidades consistentes)

// Volumes de cada tanque
printf("Digite os valores de V1 a V4") 
V = mscanf(4, "%f")
//V = [1, 2, 0.5, 3]

// Vazão volumétrica de entrada
printf("Digite F_0")
F_0 = scanf("%f")
//F_0 = 1

// Concentração do reagente A na entrada
printf("Digite CA_0")
CA_0 = scanf("%f")
//CA_0 = 1

// Fração de reciclo
printf("Digite r")
r = scanf("%f")
//r = 0.2

// Constante cinética
printf("Digite k")
k = scanf("%f")
// k = 0.1

// Vazao com reciclo
F_l = (1+r)*F_0

// Sistema linear
A = [k*V(1)+F_l, 0, 0, -r*F_0; -F_l, k*V(2)+F_l, 0, 0; 0, -F_l, k*V(3)+F_l, 0; 0, 0, -F_l, k*V(4)+F_l ]
b = [F_0 * CA_0; 0; 0; 0]

// Resolvendo o sistema
//CA = inv(A) * b
CA = linsolve(A, -b)

// Mostra resultado
printf("Ca = ")
disp(CA)

CSTR_trans.sce

// Exercício (Aula 2012-11-12)
//
// Sistema de 3 reatores CSTR em série em regime transiente
//
// @date     2012-11-12
// @author   Christian Cândido <[email protected]>
// @license  GPL-v3 <http://www.gnu.org/licenses/gpl-3.0.txt>
// @encoding UTF-8

clear;
getd;

// Entradas
global Ca0 F0 k V;

Ca0 = 1.0; // mol L^-1
F0  = 1.0; // mol L^-1 s^-1
k   = 0.1;
V   = [1.0, 2.0, 4.0];

// Condições iniciais e de contorno
t0 = 0.0; // s
tf = 20.0; // s
timestep = 0.1; // s
y0 = [0.0; 0.0; 0.0];

// Discretização
t  = [t0:timestep:tf];

// Resolve o sistema de EDOs
[y] = ode(y0, t0, t, sistema_EDO);

// Exibe os gráficos de concentração por tempo em cada reator
clf();
plot(t, y(1,:), 'k');
plot(t, y(2,:), 'kr');
plot(t, y(3,:), 'kb');
xtitle("Ca_i x t");

printf("No final:\n Ca1 = %f\n Ca2 = %f\n Ca3 = %f\n", y(1, $), y(2, $), y(3, $) );

CSTR_trans.sci

// Exercício (Aula 2012-11-12)
//
// Sistema de 3 reatores CSTR em série em regime transiente
//
// Variação da concentração de A
//
// @date     2012-11-12
// @author   Christian Cândido <[email protected]>
// @license  GPL-v3 <http://www.gnu.org/licenses/gpl-3.0.txt>
// @encoding UTF-8

function dCa=CSTR_trans(Ca0, Ca, F0, k, V)
    
    dCa = F0/V * ( Ca0 - Ca ) - k * Ca;

endfunction

f_obj.sci

// Exercício (Aula 2012-10-29)
//
// Estima as temperaturas de saída de um trocador de calor
// pelo método de Newton-Raphson
//
// Função objetivo
//
// @date     2012-10-29
// @author   Christian Cândido <[email protected]>
// @license  GPL-v3 <http://www.gnu.org/licenses/gpl-3.0.txt>
// @encoding UTF-8

function [f]=f_obj(Ts_q)
    
    global Cp_q Cp_f Te_q te_f m_q m_f U A
    
    MLDT = ((Te_q - ts_f(Ts_q)) - (Ts_q - te_f) ) / log((Te_q - ts_f(Ts_q))/(Ts_q - te_f))
    
    f = U * A * MLDT - m_q * Cp_q * (Te_q - Ts_q)
    
endfunction

fi.sci

// Exercício (Aula 2012-11-05)
//
// Programa para calcular a concentração do reagente para a reação de ordem
// alpha (r = kCa^alpha) A -> B em n reatores CSTR em série.
//
// Funções a serem zeradas
//
// @date     2012-11-05
// @author   Christian Cândido <[email protected]>
// @license  GPL-v3 <http://www.gnu.org/licenses/gpl-3.0.txt>
// @encoding UTF-8

function fi = fi(reat, C, V, F0, k, alpha)

    fi = F0 * (C(reat) - C(reat + 1)) - k * C(reat + 1)^alpha*V(reat)

endfunction

sistema_EDO.sci

// Exercício (Aula 2012-11-12)
//
// Sistema de n reatores CSTR em série em regime transiente
//
// Sistema de EDOs
//
// @date     2012-11-12
// @author   Christian Cândido <[email protected]>
// @license  GPL-v3 <http://www.gnu.org/licenses/gpl-3.0.txt>
// @encoding UTF-8

function [dy]=sistema_EDO(t, y)

// Usage: CSTR_tran(Ca0, Ca, F0, k, V)

    global Ca0 F0 k V

    dy1 = CSTR_trans( Ca0,  y(1), F0, k, V(1) )
    dy2 = CSTR_trans( y(1), y(2), F0, k, V(2) )
    dy3 = CSTR_trans( y(2), y(3), F0, k, V(3) )
    
    [dy] = [dy1; dy2; dy3]

endfunction

Trabalho-3_func.sci

// Análise e Simulação de Processos Químicos (10512-0)
// Prof. Antônio José Gonçalves da Cruz
//
// Grupo:
// Christian Carlos Cândido da Silva (389234)
// Fernando Toyoichi Yamane          (389161)
// Jorge Augusto Rodrigues Fanton    (388807)
// 
// Trabalho 3
// Solução de três reatores em série em estado transiente
// desempenhando a reação de ordem alpha:
//       A -> Produtos, r = k Ca^alpha
//
// ** Funções necessárias **

function dCa=CSTR_trans(Ca0, Ca, V)
    
    global F0 k alpha
    
    dCa = F0/V * ( Ca0 - Ca ) - k * Ca^alpha;

endfunction

function [dy]=sistema_EDO(t, y)

// Usage: CSTR_tran(Ca0, Ca, V)

    global Ca0 V

    dy1 = CSTR_trans( Ca0,  y(1), V(1) )
    dy2 = CSTR_trans( y(1), y(2), V(2) )
    dy3 = CSTR_trans( y(2), y(3), V(3) )
    
    [dy] = [dy1; dy2; dy3]

endfunction

Trabalho-3_rotina.sce

// Análise e Simulação de Processos Químicos (10512-0)
// Prof. Antônio José Gonçalves da Cruz
//
// Grupo:
// Christian Carlos Cândido da Silva (389234)
// Fernando Toyoichi Yamane          (389161)
// Jorge Augusto Rodrigues Fanton    (388807)
// 
// Trabalho 3
// Solução de três reatores em série em estado transiente
// desempenhando a reação de ordem alpha:
//       A -> Produtos, r = k Ca^alpha
//
// ** Rotina de resolução **

// Limpeza de variáveis e carregamento de funções
clear;
getd;

// Entradas
global Ca0 F0 k V alpha;

Ca0 = 1.0; // Concentração de entrada do reagente A (mol/L^3)
F0  = 1.0; // Vazão volumétrica (L^3/T)
k   = 1.0; // Constante cinética (mol^(1-alpha)/L^(3(1 - alpha))/T)
V   = [1.0, 2.0, 4.0]; // Volumes dos tanques (L^3)
alpha = 1.0; // 1.4
//alpha = scanf('%f');

// Condições iniciais e de contorno
t0 = 0.0;  // Tempo de início da solução (T)
tf = 8.0; // Tempo final da solução (T)
timestep = 0.1; // Intervalo de exibição (T)
Ca = [0.0; 0.0; 0.0]; // Aproximação inicial da solução (mol/L^3)

// Intervalos de exibição do resultado
t  = [t0:timestep:tf]; // (T)

// Resolve o sistema de EDOs
[Ca] = ode(Ca, t0, t, sistema_EDO); // (mol/L^3)

// Exibe os gráficos de concentração por tempo em cada reator
clf();
title_entity.font_size = 14;
plot(t, Ca(1,:), 'k');
plot(t, Ca(2,:), 'kr');
plot(t, Ca(3,:), 'kb');
xtitle("Ca_i x t");

printf("No final do intervalo considerado:\n Ca1 = %f\n Ca2 = %f\n Ca3 = %f\n", Ca(1, $), Ca(2, $), Ca(3, $) );

//write_csv([t; Ca]', 'relatorios/Trab3_dados_' + string(alpha) +'.csv', ',', '.');

Trabalho-4_func.sci

// Análise e Simulação de Processos Químicos (10512-0)
// Prof. Antônio José Gonçalves da Cruz
//
// Grupo:
// Christian Carlos Cândido da Silva (389234)
// Fernando Toyoichi Yamane          (389161)
// Jorge Augusto Rodrigues Fanton    (388807)
// 
// Trabalho 4
//
// ** Funções necessárias **

function dCa=CSTR_trans(Ca0, Ca, V)
    
    global F0 k alpha
    
    dCa = F0/V * ( Ca0 - Ca ) - k * Ca^alpha;

endfunction

// Sistema de EDOs para o Processo 1
function [dy]=EDOs_P1(t, y)

    global k, alpha

    dy1 = -k(1) * y(1)^alpha(1) - k(2) * y(1)^alpha(2); // dCa/dt
    dy2 = k(1) * y(1)^alpha(1); // dCq/dt
    dy3 = k(2) * y(1)^alpha(2); // dCs/dt
    
    [dy] = [dy1; dy2; dy3];

endfunction

// Sistema de EDOs para o Processo 2
function [dy]=EDOs_P2(t, y)

    global k, alpha

    dy1 = -k(1) * y(1)^alpha(1); // dCa/dt
    dy2 = k(1) * y(1)^alpha(1) - k(2) * y(2)^alpha(2); // dCq/dt
    dy3 = k(2) * y(2)^alpha(2); // dCs/dt
    
    [dy] = [dy1; dy2; dy3];

endfunction

// Sistema de EDOs para o Processo 3
function [dy]=EDOs_P3(t, y)

    global k

    dy1 = -k(1) * y(1); // dCa/dt
    dy2 = k(1) * y(1) - k(2) * y(2) * y(3)^2 - k(3) * y(2); // dCq/dt
    dy3 = 3 * k(2) * y(2) * y(3)^2 + k(3) * y(2) - k(4) * y(3); // dCs/dt
    dy4 = k(4) * y(3); // dCt/dt
    
    [dy] = [dy1; dy2; dy3; dy4];

endfunction

Trabalho-4_rotina.sce

// Análise e Simulação de Processos Químicos (10512-0)
// Prof. Antônio José Gonçalves da Cruz
//
// Grupo:
// Christian Carlos Cândido da Silva (389234)
// Fernando Toyoichi Yamane          (389161)
// Jorge Augusto Rodrigues Fanton    (388807)
// 
// Trabalho 4
// Planta multipropósito executando os processos
// P1: A -> Q, k1
//     A -> S, k2
//
// P2: A -> Q -> S, k1, k2
//
// P3: A -> Q -> S, k0, k1
//            -> S, k2
//     S -> T, k3
//
// ** Rotina de resolução **

// Limpeza de variáveis e carregamento de funções
clear;
getd;

// Entradas
global Ca0 F0 k V alpha;

proc = 1;
caso = 1;

 // Tempo final da solução (T)
timestep = 0.1; // Intervalo de exibição (T)

select caso

case 1 then
    k = [1e-1, 5e-2];
    alpha = [1, 1];
    tf = 70.0;
case 2 then
    k = [1e-2, 2.5e-3];
    alpha = [2, 1];
    tf = 1000.0;
case 2 then
    k = [1e-2, 2.5e-3];
    alpha = [0.5, 1];
    tf = 100.0;
end; 
    
// Condições iniciais e de contorno
t0 = 0;  // Tempo de início da solução (T)
Ca0 = 0.1;
Ca = [Ca0; 0.0; 0.0]; // Aproximação inicial da solução (mol/L^3)

// Intervalos de exibição do resultado
t  = [t0:timestep:tf]; // (T)

select proc

case 1 then
    // Resolve o sistema de EDOs para o processo P1
    [Ca] = ode(Ca, t0, t, EDOs_P1); // (mol/L^3)
    
    // Exibe os gráficos de concentração por tempo de cada espécie

    clf();
    plot(t, Ca(1,:), 'k');
    plot(t, Ca(2,:), 'b');
    plot(t, Ca(3,:), 'r');
    xtitle("Ca_i x t");
    
case 2 then

    // Resolve o sistema de EDOs para o processo P2
    [Ca] = ode(Ca, t0, t, EDOs_P2); // (mol/L^3)

    // Exibe os gráficos de concentração por tempo de cada espécie

    clf();
    plot(t, Ca(1,:), 'k');
    plot(t, Ca(2,:), 'b');
    plot(t, Ca(3,:), 'r');
    xtitle("Ca_i x t");
    
case 3 then

    // Resolve o sistema de EDOs para o processo 3
    k = [1e-3, 2.5e9, 1e-2, 1.0];
    tf = 10000; // Tempo final da solução (T)
    timestep = 10; // Intervalo de exibição (T)
    t  = [t0:timestep:tf]; // (T)
    Ca = [Ca0; 0.0; 0.0; 0.0]; // Aproximação inicial da solução (mol/L^3)
    [Ca] = ode(Ca, t0, t, EDOs_P3); // (mol/L^3)

    // Exibe os gráficos de concentração por tempo de cada espécie

    clf();
    //plot(t, Ca(1,:), 'k');
    plot(t, Ca(2,:), 'k');
    //plot(t, Ca(3,:), 'b');
    //plot(t, Ca(4,:), 'g');
    xtitle("Ca_i x t");

end;
    
// Verificação da conservação da massa no final da solução
printf("A massa residual foi de: %f", Ca0 - sum(Ca(:,$)))

//write_csv([t; Ca]', 'relatorios/Trab4_dados_P' + string(proc) + '_' + string(caso) + '.csv', ',', '.');

trocador.sce

// Exercício (Aula 2012-10-29)
//
// Estima as temperaturas de saída de um trocador de calor
// pelo método de Newton-Raphson
//
// @date     2012-10-29
// @author   Christian Cândido <[email protected]>
// @license  GPL-v3 <http://www.gnu.org/licenses/gpl-3.0.txt>
// @encoding UTF-8

clear;

global Cp_q Cp_f Te_q te_f m_q m_f U A

getd; 

Cp_q = 0.3; // Calor específico do fluido quente
Cp_f = 0.3; // Calor específico do fluido frio
Te_q = 100; // Temperatura de entrada do fluido quente
te_f = 20;  // Temperatura de entrada do fluido frio
m_q = 3e5;  // Vazão mássica do fluido quente
m_f = 12e4; // Vazão mássica do fluid frio

U = 110; // Coeficiente de troca térmica global
A = 250; // Área de troca

epsilon = 1e-4; // Resíduo aceito
iter_limit = 1000; // Limite de iterações
iter = 0;

h = 1e-3; // Delta para cálculo da derivada

Ts_q = 70 // Chute inicial

quoc = 1; // Para iniciar a iteração

while (abs(quoc) > epsilon & iter < iter_limit) do

    df_obj = (f_obj(Ts_q + h) - f_obj(Ts_q - h) ) / 2 / h; // f'(x) = (f(x+h) - f(x-h) )/ 2h
    quoc = f_obj(Ts_q)/df_obj; // f(x)/f'(x)

    Ts_q = Ts_q - quoc // xk+1 = xk - f(x)/f'(x)
    
    iter = iter + 1;

end

printf("----------------------------------\n")
printf("O resultado obtido para Ts_q é: %f\n", Ts_q);
printf("O resultado obtido para ts_f é: %f", ts_f(Ts_q));
printf("\nEm %d iterações", iter);
printf("\n----------------------------------")

ts_f.sci

// Exercício (Aula 2012-10-29)
//
// Estima as temperaturas de saída de um trocador de calor
// pelo método de Newton-Raphson
//
// Temperatura de saída do frio, dada a temperatura de saída do quente
//
// @date     2012-10-29
// @author   Christian Cândido <[email protected]>
// @license  GPL-v3 <http://www.gnu.org/licenses/gpl-3.0.txt>
// @encoding UTF-8

function [t]=ts_f(Ts_q)
    
    global Cp_q Cp_f Te_q te_f m_q m_f U A
    
    t = te_f + (m_q*Cp_q)/(m_f*Cp_f)*(Te_q - Ts_q)
    
endfunction