ERI

 
Cinética                            Reactores Químico

           

Cinética de Reacções Homogéneas

Imagine que tem um vidro de petri e um bom microscópio.
Dentro do petri estão uns insectos pequeninos que têm comida.
Vamos vê-los a crescer e morrer?

% Simulação de Monte-Carlo de um reactor descontínuo
% com mecanismo irreversível série paralelo
%  A + B → 2 A,  r1=k1.NA.NB, Comendo multiplicam-se
%  A + B→ A + C, r2=k2.NA.NB2, ou produzem lixo
%  A + C→ 2 C,   r3=k3.NA.NC, se comerem lixo morrem
%  A → C,        r4=k4.NA, e morrem simplesmente
% NA-Número de insectos vivos
% NB-Número de unidades de comida
% NC-Número de unidade de lixo
clf;hold on;
% repetir 11 vezes a experiência
for i=0:1:10;
clear;
% Vector das velocidades das constantes de velocidade que podem ser vistas como proporcionais às probabilidades
% r1/(r1 + r2 + r3 + r4) de se reproduzirem
% r2/(r1 + r2 + r3 + r4) de comerem e produzirem lixo
% r3/(r1 + r2 + r3 + r4) de comerem lixo e morrerem
% r4/(r1 + r2 + r3 + r4) de morrerem
% ks=[k1 k2 k3 k4]
ks=[0.1 0.001 0.01 0.001];
% Vector do nº inicial de
% Nos=[insectos comida lixo]
Nos=[100 500 0];
% espaçamento dos resultados
pt=0.001;
% tempo total da experiência
ft=2;
% mandar fazer a simulação i
A=montecarlo(ks,Nos,pt,ft);
% A=[t,NA(t),NB(t),NC(t)
plot(A(:,1),A(:,2),A(:,1),A(:,3),A(:,1),A(:,4));
end
legend('vivos','comida','lixo');

Published with MATLAB® 7.5

function Ns=montecarlo(ks,Nos,pt,ft)
% Program to do a kinetic Monte Carlo simulation of the situation
% http://www.eng.buffalo.edu/Courses/ce561/LecNotes.html
%described in CE 561, Homework 2, Problem 2.
%
% Inputs:
% ks = vector of rate constants, [k1 k2 k3 k4]
% Nos = vector of initial numbers of each species, [NAo NBo NCo]
% pt = time interval at which to save/print solution
% ft = final time, the time to end the simulation
%
% Output:
% Ns = matrix containing the solution.
%        the 1st column is time, 2nd through 4th column are NA, NB, NC
%
t=0;
tp=0;
Ns=[t Nos];
NA=Nos(1);
NB=Nos(2);
NC=Nos(3);
while (t<ft)
% compute the reaction rates
  r1=ks(1)*NA*NB;
  r2=ks(2)*NA*NB*NB;
  r3=ks(3)*NA*NC;
  r4=ks(4)*NA;
  rt=r1+r2+r3+r4;
% generate two random numbers
  rns=rand(1,2);
% pick which reaction happened
% based on one of the random numbers
  if (rns(2)<=(r1/rt))
%    reaction 1
     NA=NA+1;
     NB=NB-1;
  elseif (rns(2)<=(r1+r2)/rt)
%    reaction 2
     NB=NB-1;
     NC=NC+1;
  elseif (rns(2)<=(r1+r2+r3)/rt)
%    reaction 3
     NA=NA-1;
     NC=NC+1;
  else
%    reaction 4
     NA=NA-1;
     NC=NC+1;
  end
% advance the time, based on the other
% random number
  dt=-1/rt*log(rns(1));
  t=t+dt;
  if ((t-tp)>=pt)
     Ns=[Ns; t NA NB NC];
     tp=tp+pt;
  end
% If all the bugs are gone, all the rates go to zero, and
% we need to stop.
  if (NA==0)
     disp(['All the bugs are gone at t = ' num2str(t)])
     t=ft;
  end
end
%disp('[ ]')
%disp(['t      NA     NB       NC'])
%sz=size(Ns);
%for i=1:sz(1)
%  disp(Ns(i,:))
end

Published with MATLAB® 7.5