650 likes | 837 Views
MATLAB EM SIMULAÇÃO E CONTROLE DE PROCESSOS. PROF a . OFÉLIA DE Q.F. ARAÚJO ofelia@eq.ufrj.br www.eq.ufrj.br/links/h2cin/eqe487. Desafio. Desafio. Desafio. Caso I. Introdução à Engenharia Química: rápido aprendizado na solução de problemas introdutórios. Caso I. Caso I. Caso I.
E N D
MATLAB EM SIMULAÇÃO E CONTROLE DE PROCESSOS PROFa. OFÉLIA DE Q.F. ARAÚJO ofelia@eq.ufrj.br www.eq.ufrj.br/links/h2cin/eqe487
Caso I Introdução à Engenharia Química: rápido aprendizado na solução de problemas introdutórios
Caso I % Definição das constantes do modelo R = 1; % h/m2 A = 2; % m2 Fe = 10; % m3/h % Tempo de simulação t = 0.0:0.01:10.0; % h % Simulação da altura de líquido h = R*Fe*(1 - exp(-t/(R*A))); % m % Visualização da simulação plot(t,h) title('Simulação do tanque de nível') xlabel('Tempo (h)') ylabel('Altura (m)')
Caso II Modelagem e Dinâmica de Processos: como rapidamente um aluno pode simular um processo, avaliando comportamentos dinâmicos complexos.
CasoIIa 2 Tanques em série, volumes V1 e V2, vazão volumétrica F e concentração C
CasoIIa demo
CasoIIa xsi=[0.2 0.6 0.9 1. 2] K=2;tau=10; for i=1:length(xsi) G(i)=tf([K],[tau^2 2*xsi(i)*tau 1]); end step(G)
CasoIIb Parâmetros function [sys,x0] = reator(t,x,u,flag,U,A,DeltaH,ro,Cp,E,R,k0) % % Simula um reator CSTR (mistura perfeita) no qual se conduz uma % reação exotérmica (A->B), resfriado por serpentina % % U = 150 BTU/(h.ft2.R), coeficiente de troca térmica % A = 250 ft2, área de troca térmica % DeltaH = -30000 BTU/lbm, calor de reação % ro = 50 lb/ft3, densidade % Cp = 0.75 BTU/(lbm.R), calor específico % E = 30000 BTU/lbm, energia de ativação % R = 1.99 BTU/(lbm.R), constante dos gases % k0 = 7.08e10 1/h, termo pré-exponencial da constante de reação % switch flag case 0 % Dimensiona o sistema e inicializa os estados % sys=[estados,0,saídas,entradas,0,0] sys = [3,0,3,5,0,0]; % Condições iniciais cr = 0.1315; %lbm/ft3, concentração inicial no reator T = 500; %R, temperatura do reator V = 200; %ft3, volume do reator x0 = [cr T V]; Blá-blá-blá Condiçõesiniciais
CasoIIb case 1 % Calcula as derivadas % Atualiza entradas Cre = u(1); %lbm/ft3, concentração da alimentação Fe = u(2); %ft3/hr, vazão de alimentação F = u(3); %vazão de retirada Tc = u(4); %R, temperatura do fluido de refrigeração Te = u(5); %R, temperatura da alimentação % Cálculo das derivadas Cr = x(1); T = x(2); V = x(3); k = k0*exp(-E/(R*T)); dCr = (Fe*Cre-F*Cr)/V) - k*Cr; dV = Fe-F; dT = (Fe*Cp*ro*(Te-T) + DeltaH*k*Cr*V - U*A*(T-Tc)) /(V*ro*Cp); sys = [dCr; dT; dV]; case 3 % Calcula as saídas sys = [x(1) x(2) x(3)]; otherwise sys = []; end Entradas O MODELO!!
CasoIIc demo
CasoIIc G1=tf(2,[5 1]); G2=tf(-1,[1 1]); step(G1+G2)
CasoIId “Experimento” para geracão de dados experimentais: Temperatura do Tanque 1 Temperatura do Tanquer 2 %Roda Laboratorio de EQ sim('ldeq'); dados=[T1 T2(:,2)]; plot(T1(:,1),T1(:,2),'*b'), xlabel('Tempo'), ylabel('T1') save dados dados demo
CasoIId demo
CasoIId %Obtenção de Funções de Transferência a partir de respostas %temporais a perturbação degrau. %Carrega séries temporais load dados tempo=dados(:,1); y=dados(:,2); degrau=2; %Valor inicial para parâmetros Tau=10; Ganho=10; Tmorto=10; X0=[Tau Ganho Tmorto]; opcao=optimset('Display','iter'); Xotimo=fminsearch('fobj',X0,opcao,degrau,y,tempo)
CasoIId function[f] = fobj(X,degrau,y,tempo) Tau = abs(X(1)); Ganho = abs(X(2)); TM = abs(X(3)); Tempo = max(0,(tempo-TM)); ycalc = degrau*Ganho*(1-exp(-tempo/Tau)); %Somatório do quadrado dos erros f=(y-ycalc)'*(y-ycalc); %DESENHA PARA ACOMPANHAR A OTIMIZAÇÃO plot(tempo,y,'*',tempo,ycalc,'m') %ESCREVE NO GRÁFICO OS VALORES DOS PARÂMETROS title(num2str(X)) %A rotina a seguir força o MATLAb a desenhar o gráfico imediatamente drawnow
Caso III Dinâmica e Controle de Processos: simular estratégias de controle e sintonizar controladores com MATLAB, SIMULINK e Toolboxes de Controle e Otimização.
CasoIIIa Expansão em Frações Parciais ....
CasoIIIa Expansão em Frações Parciais .... % Coeficientes de P(s) e Q(s) em ordem decrescente de potências de s Den = [1 9 26 24 0]; Num = [1 1 1]; [R,P,K] = residue(Num,Den)
CasoIIIa Expansão em Frações Parciais .... E inversão para o domínio do Tempo
CasoIIIa % Simulação Den = [1 9 26 24 0]; Num = [1 1 1]; [R,P,K] = residue(Num,Den) Tempo=0:0.1:10; X=zeros(size(Tempo)); for i=1:length(R) X=X+R(i)*exp(P(i)*Tempo); end plot(Tempo,X)
CasoIIIa % Estabilidade Den = [1 9 26 24 0]; Num = [1 1 1]; roots(Den) ans = 0 -4.0000 -3.0000 -2.0000 ESTÁVEL
CasoIIIa U(S) X(S) demo
CasoIIIb Estabilidade de Malha de Controle demo
kc=[3 4 5 6 7 8]; simb=['c', 'm', 'g', 'b', 'r', 'k']; Texto=[]; kp=8; den=conv(conv([1 2],[1 2]),[1,2]); num=kp; Gp=tf(num,den); Gv=tf(1,1); Gm=tf(1,1); t=linspace(1,10,100); hold on for i=1:length(kc) Gc=tf(kc(i),1); G=Gc*Gv*Gp/(1+Gc*Gv*Gp*Gm); [y,t]=step(G,t); Texto=[Texto; ['K_C = ' num2str(kc(i))]]; plot(t,y,simb(i),'LineWidth',2) end legend(Texto) xlabel('Tempo') ylabel('Variável Controlada')