Lusanda Hlela Assignment 4 HLLLUS001

28/09/2018

Question 1
Code Output:

Observations:
Case 1: In the absence of diffusion limitations, whereby the reaction rate is
dominated by the rate of transport of reactants through the reaction medium, the
provided design equations are manipulated such that their rate of reaction is only
dependent of the species concentration.
Modelling the above scenario using the Eigen matrix system required uncoupling the
diffusion limited design equation using a diagonal matrix that transforms the coupled
species concentrations to uncoupled species concentrations. Translation between
these concentrations resulted in the figure above, showing the species
concentrations across the reactor volume.
It is apparent from this figure that species concentrations do not change after 2 m3,
implying that the reactions rates are minimized or close to zero at this reactor
volume.
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

Observations:
Case 2: This case explores the scenario whereby diffusion is controlled in the
reaction system. Same as in case, the Eigen matrix system is used to manipulate the
dependence of the design equations.
In this case two design equations were manipulated, one describing the mole
balance using averaged reaction rates across the catalyst volume, and the other
describing the mole balance due to reacting species diffusion in spherical particles
and the reaction mechanisms in the catalyst pores.
It is apparent from the figure above that the reaction rates are minimized, or no
reaction occurs at approximately 17 m3. This may be attributed to the reaction
system being dominated by the diffusion of the reacting species across the catalyst
pores rather than the transport of these species.
The main observation is that for this specific type of reaction, higher reaction rates
translating to low reactor volumes are achieved by diffusion unlimited systems.
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

Observations:
Case 2 was also explored in this scenario. The figure above only shows product
species concentrations because the objective is to obtain the changes in specie 2
concentration (A2) when its diffusivity is decreased by a factor of 10.
Firstly, this reaction system is assumed to be favored by no diffusion limitations, thus
by logic changing species diffusivities should not affect or have minimal effects on
the reaction conversion and A2 formed.
Diffusion controlled: question 1b Diffusion controlled: question 1 c
A1 converted: 0.8316441 A1 converted: 0.823114
A2 formed: 0.6661682 A2 formed: 0.6491079

The above code output observations confirm that this type of reaction system is not
affected by the diffusion mechanism, but the rate of transport of the reacting species
across the catalyst pores because there is minimum apparent change in the
maximum conversion and the formed species 2 concentration.
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

Code Input:

//Reaction rates
k12 = 4
k21 = 1
k23 = 1
k32 = 4

//Initial concentrations
a0 = [1 0 0]'

printf('\n \n (a) Diffusion un-controlled\n')

//Reaction rate matrix

for n = 1:3
b(:, n) = b0(n)*exp(Y(n,n)*V);
end

a = (X*b')'

scf(1);clf;
plot(V, abs(a))
xlabel('Reactor volume [m^3]','fontsize',4)
ylabel('Species concentrations','fontsize',4)
title('1a) No diffusion limitations','fontsize',4)
legend('A1','A2','A3')

printf('\n \n (b) Diffusion controlled\n')

D = diag([0.01,0.01,0.01])*1E-4;
R = 2*1E-2;
[Y,P] = spec(-inv(D)*K);
PH = real(sqrt(P)*R);
i = find(diag(PH)<%eps);
PH(i,i) = %eps;
Ks = K*(3*Y/(PH^2))*(PH*cothm(PH)-eye(3,3))*inv(Y)
[X,Y] = spec(Ks);
b0 = inv(real(X))*a0;

for n = 1:3
b(:,n) = b0(n)*exp(real(Y(n,n))*V);
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

end

a = (X*b')';

printf('\n \n A1 converted\n')
disp(1-a(101,1))
printf('\n \n A2 formed\n')
disp(a(101,2))

scf(2);clf;
plot(V, abs(a));
xlabel('Reactor volume [m^3]','fontsize',4)
ylabel('Species concentrations','fontsize',4)
title('1b) Diffusion limited','fontsize',4)
legend('A1','A2','A3')

printf('\n \n (c) Diffusion controlled\n')

D(2,2) = 0.001*1E-4;
[Y,P] = spec(-inv(D)*K);
PH = real(sqrt(P)*R);
i = find(diag(PH)<%eps);
PH(i,i) = %eps;
Ks = K*(3*Y/(PH^2))*(PH*cothm(PH)-eye(3,3))*inv(Y)
[X,Y] = spec(Ks);
b0 = inv(real(X))*a0;

for n = 1:3
b(:,n) = b0(n)*exp(real(Y(n,n))*V);
end

a = (X*b')';
x = 1 - a(:,1);

printf('\n \n A1 converted\n')
disp(1-a(101,1))
printf('\n \n A2 formed\n')
disp(a(101,2))
scf(3);clf;
plot(abs(x), abs(a(:,2:3)));
xlabel('Conversion x','fontsize',4)
ylabel('Species concentrations','fontsize',4)
title('1c) Diffusion limited','fontsize',4)
legend('A2','A3',2)
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

Question 2
Code Output:
Library routine
A1:
11.401909 -8.6836275
-8.6836275 11.401909

A2:
-0.7357585 0.551819
-1.4715172 1.1036381

Series approximation
A1:
11.401909 -8.6836275
-8.6836275 11.401909

A2:
-0.7349346 0.551407
-1.4704187 1.1030888

Eigen Decomposition
A1:
11.401909 -8.6836275
-8.6836275 11.401909

A2:
-0.7357588 0.5518191
-1.4715176 1.1036382
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

Observations:
Three methods were employed to evaluate a matrix exponential function that is
defined by an infinite series.
Library routine: Scilab function expm (square matrix exponential) was used. The
approximated values of matrix A1 and A2 are obtained as the series terms converge.
Series approximation: This method involved using any series number of terms where
by the square matrix approximations converge to the same value as the library
routine. It was observed that Matrix A1 converged at lower number of series terms
relative to Matrix A2. This might be due to the propagation of magnitude of the series
terms producing a slower converging effect for Matrix A2 than A1.
Eigen decomposition: The eigen matrix systems is explored here, and it is specified
in the question it is an alternative way of series approximation. Matrix approximation
in this method is produced by the product of matrix eigen vector(A), diagonal eigen
values matrix(B) and the inverse eigen vector matrix(A-1). The eigen vectors and
values were evaluated using the spec function, hence it returns these vector and
value matrices for a given square matrix.
The main observation is that all the explored methods produced the same results,
even though the terms of approximation might be different.
Code Input:

printf('\n Library routine')

printf('\n A1:')
A1 = expm(A1)
disp(A1)

printf('\n A2:')
A2 = expm(A2)
disp(A2)

printf('\n Series approximation')

printf('\n A1:')
A = diag([1,1]);
B = diag([1,1]);
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

for i = 1:50
A = A + (1/factorial(i))*(A1^i)
B = B + (1/factorial(i))*(A2^i)
end

disp(A)

printf('\n A2:')
disp(B)

printf('\n Eigen Decomposition')

printf('\n A1:')
[A,B] = spec(A1);
A1 = A*B*inv(A);
a1 = exp(diag(B));
A1 = A*diag(a1)*inv(A);
disp(A1)

printf('\n A2:')
[A,B] = spec(A2);
A2 = A*B*inv(A);
a2 = exp(diag(B));
A2 = A*diag(a2)*inv(A);
disp(A2)
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

Question 3
Code output:

Observations:
Analyzing the concentration change in each bulb as time progresses due to diffusion
along the capillary between bulbs for each gas was performed under the principles of
eigen matrix systems and pseudo steady state for a binary case.

Hydrogen:
Appears to be in bulb 1 before diffusion occurs. Seem to have the highest diffusion
rate and reaches equilibrium faster relative to the other gases.

Carbon dioxide:
Appears to be in bulb 2 before diffusion occurs. Seem to have the lowest diffusion
rate, reaching equilibrium slower relative to the other gases.
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

Nitrogen:
A strange behavior is observed relative to the behavior observed from the
other gases. This may be attributed to Nitrogen being in both bulbs in equal amount
before diffusion. The immediate increase of Nitrogen concentration in bulb 2 might
have resulted in the high diffusion rate of Hydrogen leaving bulb 2 into bulb 1. The
immediate decrease of Nitrogen in bulb 1 must have also been caused by the high
rate of diffusion of Hydrogen into bulb 1. This conclusion is based on maintaining the
pseudo steady state condition.

Code input:

V0 = 77.99*1E-6;
VL = 78.63*1E-6;
dtube = 2.08*1E-3;
L = 85.9*10E-3;
T = 35.2+273.15;
P = 101.325;
AL = 25810; //effective value from binary data

x10 = 0
x20 = 0.50086
x30 = 0.49914
x1L = 0.50121
x2L = 0.49879
x3L = 0

X0 = [x10 x20 x30];

D12=83.3; D13=68; D23=16.8

n=3

for i=1:n-1
B(i,1:n-1)=-xf(i)*((1)./Dij(i,1:n-1)-1/Dij(i,n))
B(i,i)=B(i,i)+sum(xf./Dij(i,:))
end

D = inv(B)

[P,d] = spec(D)
Lusanda Hlela Assignment 4 HLLLUS001
28/09/2018

X0 = inv(P)*X0(1,1:2)'
XL = inv(P)*xL(1,1:2)'
Xf = inv(P)*xf(1,1:2)'

t=linspace(1,290000,51)';

for i=1:n-1
x0(:,i)=exp(-bb*d(i,i)*t)*(X0(i)-Xf(i))+Xf(i);
xl(:,i)=exp(-bb*d(i,i)*t)*(XL(i)-Xf(i))+Xf(i);
end

X1 = (P*x0')';
X2 = (P*xl')';

X1(:,3)=ones(1,51)'-X1(:,1)-X1(:,2);
X2(:,3)=ones(1,51)'-X2(:,1)-X2(:,2);

scf(0);clf
plot(t,X1,'-')
plot(t,X2,'--')
legend('H2 - B1','N2 - B1','CO2 - B1','H2 - B2','N2 - B2','CO2- B2',4)
xlabel('time [s]','fontsize',4)
ylabel('Gas concentration','fontsize',4)
title('Pseudo steady state diffusion between Bulbs 1 & 2','fontsize',4)