Code

The following is a Matlab function for evaluating the Example Problem for Experimental Design.
~~M file available here: http://files.charlesmartinreid.com/ExperimentalDesign/ToyProblem_cmr.m~~
function [y_x1 y_x2 y_x3 y_out] = ToyProblem_cmr(m, k, Lmix, z1, z2, z3)
% This is designed for the VUQIUMap.pdf sent to me by Alex Abboud

% inputs:
%   k = reaction rate constant
%   Lmix = mixing length
%   m1 = mass flow rate with M1
%   m2 = mass flow rate with M2
% outputs:
%   y_x1 = vector mass frac of 3 components at location x1
%   y_x2 = vector mass frac of 3 components at location x2
%   y_x3 = vector mass frac of 3 components at location x3
%   y_out = vector mass frac of 3 components at end of pipe
%   Y = 3D vector space time & components

% Initial and Boundary Conditions on y:
y0 = [1; 0; 0] ;  % Initial mass fractions
y1 = [1; 0; 0] ;  % Mass fractions at inlet 1
y2 = [0; 1; 0] ;  % Mass fractions along inlet 2

% Specify Physical Parameters:
rho = 1 ;  % Density
%now an input k   = 50 ;  % Reaction Rate constant
M1  = m ;  % Mass flow rate of inlet 1  %M1=M2 for all
M2  = m ;  % Mass flow rate of inlet 2
A1  = 1 ;  % Area at inlet 1
L   = 3 ;  % Length of reactor
%now input Lmix= 1 ;  % Mixing length
U = M1/(rho*A1) ;  % Velocity (constant) in reactor
T = 1 * L/U ;      % Simulation time (as multiple of residence time)

% Specify Numerical Parameters:
Nz = 501 ;             % Number of points across length (including inlet)
%Nz = 4001;
dz = L/(Nz-1) ;        % Spatial step length
dt = 0.5 * dz/U ;      % Time step length (may be reduced slightly)
%dt = 0.25 * dz/U;
Nt = ceil(T/dt) + 1 ;  % Number of points in time (including initial)
Nzmix = round(Nz*Lmix/L) ;  % Z index where mixing length ends

z = linspace(0,L,Nz) ;
t = linspace(0,T,Nt) ;
dt = t(2) ;

% Area at each location along the length
A = zeros(Nz) ;
A(1:Nzmix) = A1 + M2/(rho*U*Lmix)*z(1:Nzmix) ;  % Area in mixing section
A(Nzmix:Nz)= A1 + M2/(rho*U) ;  % Area after mixing section

Y = zeros(3,Nz,Nt) ;
% Initial Condition,  t=0:
Y(:,:,1) = repmat(y0,[1,Nz]) ;
% Boundary Condition, z=0:
Y(:,1,1) = y1 ;

% Time Stepping
for j = 1:Nt-1
   % Euler Time Step:
   %Y(:,:,j+1) = Y(:,:,j) + dt*ToyProblem_RHS(Y(:,:,j),U,rho,A,L,M2,y2,k,dz,Nzmix) ;

   % Heun's RK2 Time Step (TVD):
   %k1 = dt*ToyProblem_RHS(Y(:,:,j),   U,rho,A,L,Lmix,M2,y2,k,dz,Nzmix) ;
   %k2 = dt*ToyProblem_RHS(Y(:,:,j)+k1,U,rho,A,L,Lmix,M2,y2,k,dz,Nzmix) ;
   %Y(:,:,j+1) = Y(:,:,j) + (k1 + k2)/2 ;

   % Shu and Osher's RK3 Time Step (TVD):
   k1 = dt*ToyProblem_RHS(Y(:,:,j),          U,rho,A,L,Lmix,M2,y2,k,dz) ;
   k2 = dt*ToyProblem_RHS(Y(:,:,j)+k1,       U,rho,A,L,Lmix,M2,y2,k,dz) ;
   k3 = dt*ToyProblem_RHS(Y(:,:,j)+(k1+k2)/4,U,rho,A,L,Lmix,M2,y2,k,dz) ;
   Y(:,:,j+1) = Y(:,:,j) + (k1 + k2 + 4*k3)/6 ;

   % RK4 Time Step (?TVD?):
   %k1 = dt*ToyProblem_RHS(Y(:,:,j),     U,rho,A,L,Lmix,M2,y2,k,dz,Nzmix) ;
   %k2 = dt*ToyProblem_RHS(Y(:,:,j)+k1/2,U,rho,A,L,Lmix,M2,y2,k,dz,Nzmix) ;
   %k3 = dt*ToyProblem_RHS(Y(:,:,j)+k2/2,U,rho,A,L,Lmix,M2,y2,k,dz,Nzmix) ;
   %k4 = dt*ToyProblem_RHS(Y(:,:,j)+k3,  U,rho,A,L,Lmix,M2,y2,k,dz,Nzmix) ;
   %Y(:,:,j+1) = Y(:,:,j) + (k1 + 2*k2 + 2*k3 + k4)/6 ;
end

dx = 1/(L*Nz);

x1_loc = round(z1/L*Nz);
x2_loc = round(z2/L*Nz);
x3_loc = round(z3/L*Nz);

%function outputs for comparison
y_x1 = Y(:,x1_loc,Nt);
y_x2 = Y(:,x2_loc,Nt);
y_x3 = Y(:,x3_loc,Nt);

y_out = Y(:,Nz,Nt);



% % Plot Results:
% figure('Units','Normalized', 'Position',[1/8,1/8,6/8,6/8])
% plot(z,Y(1,:,1), z,Y(1,:,ceil(Nt/4)), z,Y(1,:,ceil(Nt/2)), z,Y(1,:,ceil(3*Nt/4)), z,Y(1,:,Nt), 'LineSmoothing','On')
% xlabel('Position, z','FontName','Times','FontSize',22)
% ylabel('Mass Fraction of Species A, y_A','FontName','Times','FontSize',22)
% legend('Initial Time', '1/4 Final Time', '1/2 Final Time', '3/4 Final Time', 'Final Time')
% set(gca, 'FontName','Times', 'FontSize',18)
% 
% figure('Units','Normalized', 'Position',[1/8,1/8,6/8,6/8])
% plot(z,Y(2,:,1), z,Y(2,:,ceil(Nt/4)), z,Y(2,:,ceil(Nt/2)), z,Y(2,:,ceil(3*Nt/4)), z,Y(2,:,Nt), 'LineSmoothing','On')
% xlabel('Position, z','FontName','Times','FontSize',22)
% ylabel('Mass Fraction of Species B, y_B','FontName','Times','FontSize',22)
% set(gca, 'FontName','Times', 'FontSize',18)
% 
% figure('Units','Normalized', 'Position',[1/8,1/8,6/8,6/8])
% plot(z,Y(3,:,1), z,Y(3,:,ceil(Nt/4)), z,Y(3,:,ceil(Nt/2)), z,Y(3,:,ceil(3*Nt/4)), z,Y(3,:,Nt), 'LineSmoothing','On')
% xlabel('Position, z','FontName','Times','FontSize',22)
% ylabel('Mass Fraction of Species P, y_P','FontName','Times','FontSize',22)
% set(gca, 'FontName','Times', 'FontSize',18)
ToyProblem cmr.m

From charlesreid1

Code

See Also
