Multicomponent gas separation in hollow fiber membranes

Hi!
I’m trying to model gas separation of a gas with 4 components (H2, N2, N3H and O2) using a hollow fiber membrane. In papers I found the following equations:

(the change in the retentate flow rate)
( the change in the mole fraction of component i in the retentate)
( the change in the mole fraction of component i in the permeate)
I am solving these equations by first using the backwards finite difference method which turns the first two equations into:

Next, I am trying to use Newton’s method to solve all 9 equations together. This is my code:
%% Gas separation in a cocurrent hollow fiber membrane
thickness_membrane = 1e-6; % thickness membrane [m]
Perm_H2 = 250*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_N2 = 9*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_NH3 = 830*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
Perm_O2 = 37*3.35e-16;
R=8.314;
T = 273.15+25; % correlation temperature [K]
Per_H2 = Perm_H2/thickness_membrane; % Permeance of H2 [mol/m2 s Pa]
Per_N2 = Perm_N2/thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
Per_NH3 = Perm_NH3/thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
Per_O2 = Perm_O2/thickness_membrane;

%% Input parameters
F_feed = 1133.95/3.6; % feed [mol/s]
x_H2_F = 0.0268; % [-]
x_N2_F =0.9682; % [-]
x_NH3_F = 0.0049; % [-]
x_O2_F = 0.0001;

P_F = 35e5; % pressure feed [Pa]
P_P = 1e5; % pressure Permeate [Pa]

%% Assumptions membrane
L_fiber = 10;
r_fiber = (0.75e-3)/2;
n_fiber = 4000;

% Define the mesh
step= 80;
mesh = 0:(L_fiber/step):L_fiber; % linear space vector
nmesh= numel(mesh);
h = mesh(2)-mesh(1);
visc_H2 = 88e-6;
visc_N2 = 17.82e-6; % [Pa s]
visc_NH3 = 9.9e-6;
visc_O2 = 20.4e-6;

mu = visc_N2; % most present

syms xH2 yH2 Ff xN2 yN2 xNH3 yNH3 xO2 yO2
variables = [ Ff; xH2; xN2 ; xNH3; xO2; yH2; yN2; yNH3; yO2];
Final_results = zeros(10,10,numel(mesh));
xH2prev = x_H2_F;
xN2prev = x_N2_F;
xNH3prev = x_NH3_F;
xO2prev = x_O2_F;
%P_Pprev= P_P;
Ffprev = F_feed;

BETA = P_P/P_F;
for z = 1:nmesh
%% Equations to solve making use of the backward difference method
eq1 = (Ff-Ffprev)/h + 2*3.14*r_fiber*L_fiber*n_fiber*((Per_H2*(xH2*P_F + yH2*P_P)) + (Per_NH3*(xNH3*P_F + yNH3*P_P)) + (Per_N2*(xN2*P_F + yN2*P_P)) + (Per_O2*(xO2*P_F + yO2*P_P)) );

eq2 = (xH2- xH2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_H2*(xH2*P_F + yH2*P_P) + xH2*(Ff-Ffprev)/h);
eq3 = (xN2- xN2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_N2*(xN2*P_F + yN2*P_P) + xN2*(Ff-Ffprev)/h);
eq4 = (xNH3- xNH3prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_NH3*(xNH3*P_F + yNH3*P_P) + xNH3*(Ff-Ffprev)/h);
eq5 = (xO2- xO2prev)/h – 1/Ff * (2*3.14*r_fiber*L_fiber*n_fiber*Per_O2*(xO2*P_F + yO2*P_P) + xO2*(Ff-Ffprev)/h);

eq6 = yH2 – (Per_H2*xH2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_H2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq7 = yN2 – (Per_N2*xN2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_N2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq8 = yNH3 – (Per_NH3*xNH3 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_NH3*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));
eq9 = yO2 – (Per_O2*xO2 * ((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) ) / (1-BETA + BETA*Per_O2*((yH2/Per_H2) + (yN2/Per_N2) + (yNH3/Per_NH3) + (yO2/Per_O2) )));

%eq10 = ((Pp – Ppprev)/h) -( (8*R*T*mu*(F_feed-Ff))/ (3.14*r_fiber^4*n_fiber*Pp));

F = [eq1;eq2;eq3;eq4;eq5;eq6;eq7;eq8;eq9];
J=jacobian([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9],variables);
x_0 = [10 ;0.5; 0.001; 0.1; 0.0001; 0.3; 0.01; 1; 0.001];

Final_results(1,1:9,1) = x_0′;
iterations = 100;
for iter=1:iterations
%% Newton
% Substitute the initial values into the Jacobian matrix
J_subs = subs(J, variables, x_0);
F_subs = subs(F,variables,x_0);
aug_matrix = [J_subs, -F_subs];
n= 9;
for i = 1:n-1
for j = i+1:n
factor = aug_matrix(j, i) / aug_matrix(i, i);
aug_matrix(j, 🙂 = aug_matrix(j, 🙂 – factor * aug_matrix(i, :);
end
end

y_0 = zeros(n, 1);
y_0(n) = aug_matrix(n, n+1) / aug_matrix(n,n);
for i = n-1:-1:1
y_0(i) = (aug_matrix(i, n+1) – aug_matrix(i, i+1:n) * y_0(i+1:n)) / aug_matrix(i, i);
end

x_result = y_0 + x_0;
Final_results(iter+1,1:9,z) = x_result’;
err = norm(x_result – x_0);
Final_results(iter+1,10,z) = err;
x_0 = x_result;
if err <= 1e-8

%P_Pprev = x_result(10);

disp([‘Converged after ‘, num2str(iter), ‘ iterations’]);
break;
end
end
xH2prev = x_0(2);
xN2prev = x_0(3);
xNH3prev = x_0(4);
xO2prev = x_0(5);
Ffprev = x_0(1);

end

solution_Ff_Per_step = [];
solution_xH2_Per_step = [];
solution_xN2_Per_step = [];
solution_xNH3_Per_step = [];
solution_xO2_Per_step = [];
solution_yH2_Per_step = [];
solution_yN2_Per_step = [];
solution_yNH3_Per_step = [];
solution_yO2_Per_step = [];
solution_P_P_Per_step = [];

for m=1:numel(mesh)
solution_Ff_Per_step = [Final_results(end,1,m) solution_Ff_Per_step];
solution_xH2_Per_step = [Final_results(end,2,m) solution_xH2_Per_step];
solution_xN2_Per_step = [Final_results(end,3,m) solution_xN2_Per_step];
solution_xNH3_Per_step = [Final_results(end,4,m) solution_xNH3_Per_step];
solution_xO2_Per_step = [Final_results(end,5,m) solution_xO2_Per_step];
solution_yH2_Per_step = [Final_results(end,6,m) solution_yH2_Per_step];
solution_yN2_Per_step = [Final_results(end,7,m) solution_yN2_Per_step];
solution_yNH3_Per_step = [Final_results(end,8,m) solution_yNH3_Per_step];
solution_yO2_Per_step = [Final_results(end,9,m) solution_yO2_Per_step];
%solution_P_P_Per_step = [Final_results(end,10,m) solution_P_P_Per_step];

end
figure(1)
plot(mesh,solution_Ff_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Flow rate retentate [m/s]’)
figure(2)
plot(mesh,solution_xH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in retentate’)
figure(3)
plot(mesh,solution_xN2_Per_step )
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in retentate’)
figure(4)
plot(mesh,solution_xNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in retentate’)
figure(5)
plot(mesh,solution_xO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in retentate’)
figure(6)
plot(mesh,solution_yH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in Permeate’)
figure(7)
plot(mesh,solution_yN2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in Permeate’)
figure(8)
plot(mesh,solution_yNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in Permeate’)
figure(9)
plot(mesh,solution_yO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in Permeate’)
figure(10)
% plot(mesh,solution_P_P_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Permeate pressure [bar]’)
For some reason, when the method converges it will always return the initial values of xH2prev, xNH3prev, xO2prev. xN2prev and Ffprev, meaning that there is no change in mole fractions along the membrane. When I try to change the initial conditions, the graphs show spikes in random positions whilst being 0 in other positions along the membrane. This should not be the case. Could you please help me figure out what went wrong in my code?Hi!
I’m trying to model gas separation of a gas with 4 components (H2, N2, N3H and O2) using a hollow fiber membrane. In papers I found the following equations:

%% Input parameters
F_feed = 1133.95/3.6; % feed [mol/s]
x_H2_F = 0.0268; % [-]
x_N2_F =0.9682; % [-]
x_NH3_F = 0.0049; % [-]
x_O2_F = 0.0001;

P_F = 35e5; % pressure feed [Pa]
P_P = 1e5; % pressure Permeate [Pa]

%% Assumptions membrane
L_fiber = 10;
r_fiber = (0.75e-3)/2;
n_fiber = 4000;

mu = visc_N2; % most present

%eq10 = ((Pp – Ppprev)/h) -( (8*R*T*mu*(F_feed-Ff))/ (3.14*r_fiber^4*n_fiber*Pp));

F = [eq1;eq2;eq3;eq4;eq5;eq6;eq7;eq8;eq9];
J=jacobian([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9],variables);
x_0 = [10 ;0.5; 0.001; 0.1; 0.0001; 0.3; 0.01; 1; 0.001];

y_0 = zeros(n, 1);
y_0(n) = aug_matrix(n, n+1) / aug_matrix(n,n);
for i = n-1:-1:1
y_0(i) = (aug_matrix(i, n+1) – aug_matrix(i, i+1:n) * y_0(i+1:n)) / aug_matrix(i, i);
end

x_result = y_0 + x_0;
Final_results(iter+1,1:9,z) = x_result’;
err = norm(x_result – x_0);
Final_results(iter+1,10,z) = err;
x_0 = x_result;
if err <= 1e-8

%P_Pprev = x_result(10);

disp([‘Converged after ‘, num2str(iter), ‘ iterations’]);
break;
end
end
xH2prev = x_0(2);
xN2prev = x_0(3);
xNH3prev = x_0(4);
xO2prev = x_0(5);
Ffprev = x_0(1);

end

end
figure(1)
plot(mesh,solution_Ff_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Flow rate retentate [m/s]’)
figure(2)
plot(mesh,solution_xH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in retentate’)
figure(3)
plot(mesh,solution_xN2_Per_step )
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in retentate’)
figure(4)
plot(mesh,solution_xNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in retentate’)
figure(5)
plot(mesh,solution_xO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in retentate’)
figure(6)
plot(mesh,solution_yH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in Permeate’)
figure(7)
plot(mesh,solution_yN2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in Permeate’)
figure(8)
plot(mesh,solution_yNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in Permeate’)
figure(9)
plot(mesh,solution_yO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in Permeate’)
figure(10)
% plot(mesh,solution_P_P_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Permeate pressure [bar]’)
For some reason, when the method converges it will always return the initial values of xH2prev, xNH3prev, xO2prev. xN2prev and Ffprev, meaning that there is no change in mole fractions along the membrane. When I try to change the initial conditions, the graphs show spikes in random positions whilst being 0 in other positions along the membrane. This should not be the case. Could you please help me figure out what went wrong in my code? Hi!
I’m trying to model gas separation of a gas with 4 components (H2, N2, N3H and O2) using a hollow fiber membrane. In papers I found the following equations:

%% Input parameters
F_feed = 1133.95/3.6; % feed [mol/s]
x_H2_F = 0.0268; % [-]
x_N2_F =0.9682; % [-]
x_NH3_F = 0.0049; % [-]
x_O2_F = 0.0001;

P_F = 35e5; % pressure feed [Pa]
P_P = 1e5; % pressure Permeate [Pa]

%% Assumptions membrane
L_fiber = 10;
r_fiber = (0.75e-3)/2;
n_fiber = 4000;

mu = visc_N2; % most present

%eq10 = ((Pp – Ppprev)/h) -( (8*R*T*mu*(F_feed-Ff))/ (3.14*r_fiber^4*n_fiber*Pp));

F = [eq1;eq2;eq3;eq4;eq5;eq6;eq7;eq8;eq9];
J=jacobian([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9],variables);
x_0 = [10 ;0.5; 0.001; 0.1; 0.0001; 0.3; 0.01; 1; 0.001];

y_0 = zeros(n, 1);
y_0(n) = aug_matrix(n, n+1) / aug_matrix(n,n);
for i = n-1:-1:1
y_0(i) = (aug_matrix(i, n+1) – aug_matrix(i, i+1:n) * y_0(i+1:n)) / aug_matrix(i, i);
end

x_result = y_0 + x_0;
Final_results(iter+1,1:9,z) = x_result’;
err = norm(x_result – x_0);
Final_results(iter+1,10,z) = err;
x_0 = x_result;
if err <= 1e-8

%P_Pprev = x_result(10);

disp([‘Converged after ‘, num2str(iter), ‘ iterations’]);
break;
end
end
xH2prev = x_0(2);
xN2prev = x_0(3);
xNH3prev = x_0(4);
xO2prev = x_0(5);
Ffprev = x_0(1);

end

end
figure(1)
plot(mesh,solution_Ff_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Flow rate retentate [m/s]’)
figure(2)
plot(mesh,solution_xH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in retentate’)
figure(3)
plot(mesh,solution_xN2_Per_step )
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in retentate’)
figure(4)
plot(mesh,solution_xNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in retentate’)
figure(5)
plot(mesh,solution_xO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in retentate’)
figure(6)
plot(mesh,solution_yH2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction H_{2} in Permeate’)
figure(7)
plot(mesh,solution_yN2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction N_{2} in Permeate’)
figure(8)
plot(mesh,solution_yNH3_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction NH_{3} in Permeate’)
figure(9)
plot(mesh,solution_yO2_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Mole fraction O_{2} in Permeate’)
figure(10)
% plot(mesh,solution_P_P_Per_step)
xlabel(‘Distance along membrane module [m]’)
ylabel(‘Permeate pressure [bar]’)
For some reason, when the method converges it will always return the initial values of xH2prev, xNH3prev, xO2prev. xN2prev and Ffprev, meaning that there is no change in mole fractions along the membrane. When I try to change the initial conditions, the graphs show spikes in random positions whilst being 0 in other positions along the membrane. This should not be the case. Could you please help me figure out what went wrong in my code? iterative, newton, backwards difference, finite MATLAB Answers — New Questions

Cart

Cart