matlab

function Modestino_HW5 % Created by Miguel Modestino, 03/15/2017 % part (ii) Nx=100; Ny=200; L=0.1; %m d=0.003; %m CA0=5e3; %mol/m3 Vmax= 1e-3; %m/s [c,dx,dy]=reactor(Nx,Ny,L,d,CA0,Vmax); c(Nx*Ny) x=ones(Nx,1); y=ones(Ny,1); C=zeros(Ny,Nx); for i=1:Nx for j=1:Ny n=(i-1)*Ny+j; C(j,i)=c(n); x(i)=i*dx; y(j)=j*dy; end end imagesc(x,y,C) xlabel('x [m]') ylabel('y [m]') title('Concentration of A, part (ii)') colorbar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Part (iii) Vmax=logspace(-3,-1,10); CAf=zeros(10,1); for i=1:10 [c,dx,dy]=reactor(Nx,Ny,L,d,CA0,Vmax(i)); j=[1:Ny]; nAf(i)=sum(c((Nx-1)*Ny+j)*dy*dx); end figure [AX,H1,H2]=plotyy(Vmax',(CA0*dx*d-nAf)/(CA0*dx*d),Vmax',CA0*dx*d-nAf); ylabel(AX(1),'Conversion of A') ylabel(AX(2),'Production of B [mol s^-^1 m^-^1]') xlabel('V_m_a_x [m s^-^1]') title('Part (iii)') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Part (iv) x0=[1e-4 0.1 1e-4]'; options = optimoptions('fmincon', 'OptimalityTolerance', 5e-6); [Xopt,FVAL,EXITFLAG,OUTPUT]=fmincon(@(x)CostFunction(x),x0,[],[],[],[],[5e-6;5e-3;0],[inf;1;inf],@(x)ConcCons(x),options) 1 d=Xopt(1) L=Xopt(2) Vmax=Xopt(3) [c,dx,dy]=reactor(Nx,Ny,L,d,CA0,Vmax); x=ones(Nx,1); y=ones(Ny,1); C=zeros(Ny,Nx); for i=1:Nx for j=1:Ny n=(i-1)*Ny+j; C(j,i)=c(n); x(i)=i*dx; y(j)=j*dy; end end figure imagesc(x,y,C) xlabel('x [m]') ylabel('y [m]') title('Concentration of A, part (iv)') colorbar end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Reactor Function function [c,dx,dy]=reactor(Nx,Ny,L,d,c0,Vmax) dx=L/(Nx-1); dy=d/(Ny-1); j=10; F=96500; Di=1e-9; %m2/s ri=-j/F; %mol/s/m2 a3=Di/dx^2; a4=Di/dy^2; a5=a4; A=spalloc(Nx*Ny,Nx*Ny,5*(Nx-2)*(Ny-2)+3*Ny+2*(Nx-2)); b=zeros(Nx*Ny,1); for i=2:Nx-1 for j=2:Ny-1 y=j*dy; Vx=4*Vmax/d*(y-y^2/d); a1=Vx/dx+Di/dx^2; a2=-Vx/dx-2*Di/dx^2-2*Di/dy^2; n1=((i-1)-1)*Ny+j; n2=(i-1)*Ny+j; n3=((i+1)-1)*Ny+j; n4=(i-1)*Ny+(j+1); n5=(i-1)*Ny+(j-1); A(n2,n1)=a1; A(n2,n2)=a2; A(n2,n3)=a3; A(n2,n4)=a4; A(n2,n5)=a5; end end for j=1:Ny i=1; n=(i-1)*Ny+j; A(n,n)=1; b(n)=c0; end for j=1:Ny i=Nx; n=(i-1)*Ny+j; n_prev=(i-2)*Ny+j; A(n,n)=1; A(n,n_prev)=-1; end for i=2:Nx-1 j=1; n0=(i-1)*Ny+j; n1=(i-1)*Ny+j+1; A(n0,n0)=-1; A(n0,n1)=1; end for i=2:Nx-1 j=Ny; n0=(i-1)*Ny+j; n1=(i-1)*Ny+j-1; A(n0,n0)=-1; A(n0,n1)=1; b(n0)=-ri/Di*dx; end c=A\b; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Cost Function function E_B=CostFunction(x) d=x(1); L=x(2); Vmax=x(3); CA0=5000; Nx=100; Ny=200; sigma= 0.1;% 1/(ohm-m) mu = 1e-3;% Kg/(m-s) j = 10; %A/m2 Q=2/3*Vmax*d; % m2/s Pohm = j^2*d*L/sigma; Phyd = 12*mu*Q^2*L/d^3; Ptot=Pohm+Phyd; [c,dx,dy]=reactor(Nx,Ny,L,d,CA0,Vmax); j=[1:Ny]; nAf = sum(c((Nx-1)*Ny+j)*dy*dx); nBf = CA0*dx*d-nAf; E_B=Ptot/nBf; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [C,Ceq]=ConcCons(x) d=x(1); L=x(2); Vmax=x(3); CA0=5000; Nx=100; Ny=200; [c,dx,dy]=reactor(Nx,Ny,L,d,CA0,Vmax); Ceq=[]; C = - c(Nx*Ny); end