P2 (1)

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

clear all
clc
% Entrada R1
F1=40; % ft3/h
T1=530; %ºR
Ca1=0.5; % lb mol/ft3
Cb1=0.5; % lb mol/ft3
Cc1=0; % lb mol/ft3
Cd1=0; % lb mol/ft3
Ce1=0;
mVa1= 95.365 ; % ft3/h vazão do vapor F
mVa1Erro=mVa1;
% R1 -> R2
T2=449.68; %390.86
F2=40; % ft3/h
F2ant=F2;
F2erro=F2;
Ca2= 0.2765;
Cb2= 0.2765;
Cc2= 0.2235;
Cd2= 0.2235;
Ce2=0;
Cf2=0;
F3=F2;
Ce3=Cc2;
Cf3=0;
T3=550;
% Saída R2
F4=80; % ft3/h
F4ant=F4;
F4erro=F4;
Ca4=0.096600; % lb mol/ft3
Cb4=0.096600; % lb mol/ft3
Cc4= 1.8513e-03;
Cd4=0.1634;
Ce4=1.8513e-03;
Cf4=0.109915;
T4=620; % ?R TEMPERATURA DE SAIDA DO REATOR V2
% Dimensoes
V1=78.539816; % ft3 VOLUME REATOR 1
Vj2=48.080897; % ft3 VOLUME JAQUETA
A2=145.97; % ft2 ÃREA SUPERFICIAL ENTRE A JAQUETA E O REATOR
Xt=10; % ft COMPRIMENTO DA TUBULAÇÃO QUE LIGA O REATOR 1 AO 2
Rt=0.5; %ft RAIO DO TUBO
V2=78.539816; % ft3 VOLUME REATOR 2
% Temperatura ambiente
Tamb=500;
% Coeficiente global de troca termica da jaqueta, tubulação que liga 1 ao 2 e do vapor
U1=250; % BTU/(h ft2 ?R) TANQUE 1
k= 0.029; %condutividade termica do material
% Propriedades do Fluido reativo
cp=0.75; % BTU/(lb ?R)
ro=50; % lb/ft3
% Propriedades do Fluido refrigerante (da jaqueta)
cpj=1; % BTU/(lb ?R)
roj=50; % lb/ft3 62.3
% Propriedades do Fluido de aquecimento
cpa=0.7; % BTU/(lb ?R)
roa=60 ; % lb/ft3
hv=1800; % BTU/lb mol entalpia vapor saturado
hl=250; % BTU/lb mol entalpia liquido condensado
% Propriedades do trocador de calor
Ut1=300;
Rt1=2.15;
Lt1=10;
Tq=650;
Fq=20;
% Constantes de Arrenius Reator 1
Aa1=7.08e10; % fator pre-exponencial 1/h
Ea1=22000; % energia de ativação reator 1 BTU/lb mol
Ra1=1.99; % constante dos gases BTU/(lb mol ?R)
% Constantes de Arrenius Reator 2
Aa2=7.08e10; % fator pre-exponencial 1/h
Ea2=18000; % energia de ativacao reator 1 BTU/lb mol
Ra2=1.99; % constante dos gases BTU/(lb mol ?R)
% Jaqueta de refrigerante
Fj = -20.870; % ft3/h
Fjerro = Fj;
Tj0 = 530;
Tj = 622.65; % ºR TEMPERATURA DO REFRIGERANTE FINAL 625
% Operacao do R1
Tset1=T2;
Vset1=V1; % ft3
% Condicoes de op do reator 2
Tset2=T4; % T4
Vset2=V2; % ft3
% Controlador do Volume do reator 1
intPIDvl1=0;
erro1Vvl1=0;
erro2Vvl1=0;
Rv1_max = 0.1;
Kfvl1=-7; % PROPORCIONAL
Tdfvl1=0.08; % DERIVATIVO
Tifvl1=0.8; % INTEGRAL
% Controlador de vazao de vapor do reator 1
intPIDva1=0;
erro1Vva1=0;
erro2Vva1=0;
Kva1=20; % PROPORCIONAL
Tdva1=0.01; % DERIVATIVO
Tiva1=80; % INTEGRAL
% Controlador do Volume do reator 2
intPIDvl2=0;
erro1Vvl2=0;
erro2Vvl2=0;
Rv2_max = 0.1;
Kfvl2=-5; % PROPORCIONAL
Tdfvl2=0.15; % DERIVATIVO
Tifvl2=2; % INTEGRAL
% Controlador da vazao da jaqueta do reator 2
intPIDj2=0;
erro1Vj2=0;
erro2Vj2=0;
Kfj2=-10; % PROPORCIONAL
Tdfj2=0; % DERIVATIVO
Tifj2=0.6; %I NTEGRAL
% Calor de reacao reator 1 e reator 2
deltaH1=30000; %BTU/lb mol de A
deltaH2=-30000; %BTU/lb mol de A
% Calculado o produto do reator 1 (Ca2, Cb2, Cc2, Cd2, T2)
V1Ca2=V1*Ca2;
V1Cb2=V1*Cb2;
V1Cc2=V1*Cc2;
V1Cd2=V1*Cd2;
V1Ce2=V1*Ce2;
V1T2=V1*T2;
% Calculado o produto do reator 2 (Ca4, Cb4, Cc4, Cd4, Ce4, Cf4, T4)
V2Ca4=V2*Ca4;
V2Cb4=V2*Cb4;
V2Cc4=V2*Cc4;
V2Cd4=V2*Cd4;
V2Ce4=V2*Ce4;
V2Cf4=V2*Cf4;
V2T4=V2*T4;
% variaveis de tempo
t = 0;
tf = 15;
dt = 0.005;
i = 1;
%lembra que o disturbio coloca antes do while aqui
%Disturbio
%Tset1=Tset1+2;
%Tset2=Tset2+2;
Vset1=Vset1+2;
Vset2=Vset2+2;
while t < tf
 % Controlador PID de vazao de vapor do reator 1
 %----------------------------------------------------------
 % A integral do PID no volume
 intPIDva1 = intPIDva1 + (Tset1-T2)*dt;
 % A derivada do PID no volume
 derPIDva1 = (erro2Vva1-erro1Vva1)/dt;
 erro1Vva1 = erro2Vva1;
 % A parte proporcional do PID no volume
 proPIDva1 = (Tset1-T2);
 % A equacao do PID no volume
 mVa1 = mVa1Erro + Kva1*proPIDva1 + Kva1*Tdva1*derPIDva1 + (Kva1/Tiva1)*intPIDva1;
 %----------------------------------------------------------
 % Controlador PID de volume do reator 1
 %----------------------------------------------------------
 % A integral do PID no volume
 intPIDvl1 = intPIDvl1 + (Vset1-V1)*dt;
 % A derivada do PID no volume
 derPIDvl1 = (erro2Vvl1-erro1Vvl1)/dt;
 erro1Vvl1 = erro2Vvl1;
 % A parte proporcional do PID no volume
 proPIDvl1 = (Vset1-V1);
 % A equacao do PID no volume
 F2 = F2erro + Kfvl1*proPIDvl1 + Kfvl1*Tdfvl1*derPIDvl1 + (Kfvl1/Tifvl1)*intPIDvl1;
 %----------------------------------------------------------
 if (abs(F2 - F2ant) > Rv1_max)
 if (F2 >= F2ant)
 F2 = F2ant + Rv1_max;
 elseif (F2 < F2ant)
 F2 = F2ant - Rv1_max;
 endif
 endif
 F2ant=F2;
 % -------- Reator 1 --------
 % Calculo de calor trocado
 Q1=(mVa1*(hv-hl));
 % k da reacao
 k1=(Aa1)*exp(-Ea1/(Ra1*T2));
 % Calculando as derivadas no reator 1
 dV1dt=F1-F2;
 dV1Cadt=F1*Ca1-(F2)*Ca2-V1*k1*Ca2*Cb2;
 dV1Cbdt=F1*Cb1-(F2)*Cb2-V1*k1*Ca2*Cb2;
 dV1Ccdt=F1*Cc1-(F2)*Cc2+V1*k1*Ca2*Cb2;
 dV1Cddt=F1*Cd1-(F2)*Cd2+V1*k1*Ca2*Cb2;
 dV1T2dt=F1*T1-(F2)*T2-deltaH1*V1*k1*Ca2*Cb2/(cp*ro)+Q1/(cp*ro);
 % Calculando as variaveis do reator 1
 V1=V1+dV1dt*dt;
 V1Ca2=V1Ca2+dV1Cadt*dt;
 V1Cb2=V1Cb2+dV1Cbdt*dt;
 V1Cc2=V1Cc2+dV1Ccdt*dt;
 V1Cd2=V1Cd2+dV1Cddt*dt;
 V1T2=V1T2+dV1T2dt*dt;
 % Concentracoes finais
 Ca2=V1Ca2/V1;
 Cb2=V1Cb2/V1;
 Cc2=V1Cc2/V1;
 Cd2=V1Cd2/V1;
 T2=V1T2/V1;
 % Controlador PID de volume do reator 2
 %----------------------------------------------------------
 % A integral do PID no volume
 intPIDvl2 = intPIDvl2 + (Vset2-V2)*dt;
 % A derivada do PID no volume
 derPIDvl2 = (erro2Vvl2-erro1Vvl2)/dt;
 erro1Vvl2 = erro2Vvl2;
 % A parte proporcional do PID no volume
 proPIDvl2 = (Vset2-V2);
 % A equa??o do PID no volume
 F4 = F4erro + Kfvl2*proPIDvl2 + Kfvl2*Tdfvl2*derPIDvl2 + (Kfvl2/Tifvl2)*intPIDvl2;
 %----------------------------------------------------------
 if (abs(F4 - F4ant) > Rv2_max)
 if (F4 >= F4ant)
 F4 = F4ant + Rv2_max;
 elseif (F4 < F4ant)
 F4 = F4ant - Rv2_max;
 endif
 endif
 F4ant=F4;
 % Controlador PID na Vazao de Refrigerante do reator 2
 %----------------------------------------------------------
 % A integral do PID no volume
 intPIDj2= intPIDj2 + (Tset2-T4)*dt;
 % A derivada do PID no volume
 derPIDj = (erro2Vj2-erro1Vj2)/dt;
 erro1Vj2 = erro2Vj2;
 % A parte proporcional do PID no volume
 proPIDj = (Tset2-T4);
 % A equa??o do PID no volume
 Fj = Fjerro + Kfj2*proPIDj + Kfj2*Tdfj2*derPIDj + (Kfj2/Tifj2)*intPIDj2;
 %----------------------------------------------------------
 % -------- Reator 2 --------
 % Concentracao de E
 Ce3=Cc2;
 % Controle de vazao do E
 F3 = (F2 * Cc2) / Ce3;
 % k2 da reação no reator 2
 k2=(Aa2)*exp(-Ea2/(Ra2*T4));
 % O calor trocado com a camisa
 Q2=U1*A2*(T4-Tj);
 % Calculando as derivadas no reator 2
 dV2dt=F2+F3-F4;
 dV2Ccdt=F2*Cc2-(F4)*Cc4-V2*k2*Cc4*Ce4;
 dV2Cedt=F3*Ce3-(F4)*Ce4-V2*k2*Cc4*Ce4;
 dV2Cfdt=-(F4)*Cf4+V2*k2*Cc4*Ce4;
 dV2T4dt= (F3*T3) + (F2*T2)-F4*T4-deltaH2*V2*k2*Cc4*Ce4/(cp*ro)-Q2/(cp*ro);
 dTjdt=Fj*(Tj0-Tj)/Vj2 + Q2/(cpj*roj*Vj2);
 % Calculando as variaveis
 V2=V2+dV2dt*dt;
 V2Cc4=V2Cc4+dV2Ccdt*dt;
 V2Ce4=V2Ce4+dV2Cedt*dt;
 V2Cf4=V2Cf4+dV2Cfdt*dt;
 V2T4=V2T4+dV2T4dt*dt;
 Tj=Tj+dTjdt*dt;
 %Concentracoes finais
 Cc4=V2Cc4/V2;
 Ce4=V2Ce4/V2;
 Cf4=V2Cf4/V2;
 T4=V2T4/V2;
 % Ajuste dos erros dos volumes
 erro2Vvl1 = Vset1-V1;
 erro2Vva1 = Tset1-T2;
 erro2Vvl2 = Vset2-V2;
 erro2Vj2 = Tset2-T4;
 % Vetores para os grafico
 tempoV(i)=t;
 % Vetores para os grafico do reator 1
 T2V(i)=T2;
 V1V(i)=V1;
 F2V(i)=F2;
 Ca2V(i)=Ca2;
 Cb2V(i)=Cb2;
 Cc2V(i)=Cc2;
 Cd2V(i)=Cd2;
 mVa1V(i)=mVa1;
% Trocador de calor 1
%balanco de energia moodle
alfa1 = F4*ro*cp;
beta1 = 2*(3.1415)*Ut1*Rt1*Lt1;
gama1 = Fq*roa*cpa;
T4_q = (gama1*Tq*(-exp(beta1/alfa1)) + alfa1*T4*exp(beta1/gama1) + gama1*Tq*exp(beta1/gama1) - gama1*T4*exp(beta1/gama1))/(alfa1*exp(beta1/gama1) - gama1*exp(beta1/alfa1));
 % Vetores para os graficos do reator 2
 T4V(i)=T4_q;
 V2V(i)=V2;
 F4V(i)=F4;
 Cc4V(i)=Cc4;
 Ce4V(i)=Ce4;
 Cf4V(i)=Cf4;
 TjV(i)=Tj;
 FjV(i)=Fj;
 t = t + 0.005;
 i = i + 1;
end
% Desenho do grafico do reator 1
figure(1)
%figure('Name','Reator 1');
subplot(3,3,3); plot(tempoV,T2V,'.');xlabel('tempo');ylabel('temperatura'); hold on;
subplot(3,3,2); plot(tempoV,V1V,'.');xlabel('tempo');ylabel('volume'); hold on;
subplot(3,3,1); plot(tempoV,F2V,'.');xlabel('tempo');ylabel('vazao'); hold on;
subplot(3,3,4); plot(tempoV,Ca2V,'.');xlabel('tempo');ylabel('Ca'); hold on;
subplot(3,3,5); plot(tempoV,Cb2V,'.');xlabel('tempo');ylabel('Cb'); hold on;
subplot(3,3,6); plot(tempoV,Cc2V,'.');xlabel('tempo');ylabel('Cc'); hold on;
subplot(3,3,7); plot(tempoV,Cd2V,'.');xlabel('tempo');ylabel('Cd'); hold on;
subplot(3,3,8); plot(tempoV,mVa1V,'.');xlabel('tempo');ylabel('Massa de vapor'); hold on;
% Desenho do grafico do reator 2
figure(2)
% figure('Name','Reator 2');
subplot(3,3,3); plot(tempoV,T4V,'.');xlabel('tempo');ylabel('temperatura'); hold on;
subplot(3,3,2); plot(tempoV,V2V,'.');xlabel('tempo');ylabel('volume'); hold on;
subplot(3,3,1); plot(tempoV,F4V,'.');xlabel('tempo');ylabel('vazão'); hold on;
subplot(3,3,4); plot(tempoV,Cc4V,'.');xlabel('tempo');ylabel('Cc'); hold on;
subplot(3,3,5); plot(tempoV,Ce4V,'.');xlabel('tempo');ylabel('Ce'); hold on;
subplot(3,3,6); plot(tempoV,Cf4V,'.');xlabel('tempo');ylabel('Cf'); hold on;
subplot(3,3,7); plot(tempoV,TjV,'.');xlabel('tempo');ylabel('Temperatura de refrigerante'); hold on;
subplot(3,3,8); plot(tempoV,-FjV,'.');xlabel('tempo');ylabel('vazao do refrigerante'); hold on;