Multicomponent gas separation in hollow fiber membranes
Mostra commenti meno recenti
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?
4 Commenti
Torsten
il 8 Giu 2024
To avoid numerical problems that result from self-written code, use a professional MATLAB integrator, e.g. ode15s.
MeepMeep
il 8 Giu 2024
Torsten
il 8 Giu 2024
I think the last four equations should start as
x(6) - (par.Per_H2*x(2)... )
x(7) - (par.Per_N2*x(3)... )
x(8) - (par.Per_NH3*x(4) ...)
x(9) - (par.Per_O2*x(5) ...)
because these are algebraic equations of the form
0*dyi/dt = y(i) - rhs(i) (i=6,...,9)
Try to give initial guesses for the yi that satisfy the algebraic equations - otherwise you might get the problem that ode15s cannot compute consistent initial values. You might want to use "fsolve" to compute these consistent initial values for the yi before calling ode15s.
MeepMeep
il 8 Giu 2024
Risposta accettata
par.thickness_membrane = 1e-6; % thickness membrane [m]
par.Perm_H2 = 250*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
par.Perm_N2 = 9*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
par.Perm_NH3 = 830*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
par.Perm_O2 = 37*3.35e-16;
par.R=8.314;
par.T = 273.15+25; % correlation temPeratures [K]
par.Per_H2 = par.Perm_H2/par.thickness_membrane; % Permeance of H2 [mol/m2 s Pa]
par.Per_N2 = par.Perm_N2/par.thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
par.Per_NH3 = par.Perm_NH3/par.thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
par.Per_O2 = par.Perm_O2/par.thickness_membrane;
%T_g = 273.15+155; % glass transition temPerature of polymer membrane [K]
%% Input parameters
par.F_feed = 1133.95/3.6; % feed [mol/s]
par.x_H2_F = 0.0268322; % [-]
par.x_N2_F =0.96828; % [-]
par.x_NH3_F = 0.00475768; % [-]
par.x_O2_F = 0.000129826;
%par.y_H2_init = par.Per_H2*par.x_H2_F/ (par.Per_H2*par.x_H2_F + par.Per_N2*par.x_N2_F + par.Per_NH3*par.x_NH3_F + par.Per_O2*par.x_O2_F);
%par.y_N2_init = par.Per_N2*par.x_N2_F/ (par.Per_H2*par.x_H2_F + par.Per_N2*par.x_N2_F + par.Per_NH3*par.x_NH3_F + par.Per_O2*par.x_O2_F);
%par.y_NH3_init = par.Per_NH3*par.x_NH3_F/ (par.Per_H2*par.x_H2_F + par.Per_N2*par.x_N2_F + par.Per_NH3*par.x_NH3_F + par.Per_O2*par.x_O2_F);
%par.y_O2_init = par.Per_O2*par.x_O2_F/ (par.Per_H2*par.x_H2_F + par.Per_N2*par.x_N2_F + par.Per_NH3*par.x_NH3_F + par.Per_O2*par.x_O2_F);
par.P_f = 35e5; % pressure feed [Pa]
par.P_P = 1e5; % pressure Permeate [Pa]
%% Assumptions membrane
par.L_fiber = 15;
par.r_fiber = (0.75e-3)/2;
par.n_fiber = 4000;
par.visc_H2 = 88e-6;
par.visc_N2 = 17.82e-6; % [Pa s]
par.visc_NH3 = 9.9e-6;
par.visc_O2 = 20.4e-6;
par.mu = par.visc_N2; % most present
par.BETA = par.P_P/par.P_f;
% Initial guess for x(6) to x(9)
%par.initial_guess = [0.3 ;0.1 ;0.5 ;0.1]; % Example guess, you may need to adjust
%par.initial_guess = [0.1;0.015;0.1;0.015]; % Example guess, you may need to adjust
par.initial_guess = [1 2 3 4];
% Use fsolve to find consistent initial values
options = optimoptions('fsolve', 'Display', 'iter','TolX',1e-10,'TolFun',1e-10);
[par.initial_values, fval, exitflag] = fsolve(@(x) algebraic_constraints(x, par), par.initial_guess, options);
par.initial_values=par.initial_values.^2;
par.initial_values
sum(par.initial_values)
% Check if fsolve was successful
if exitflag <= 0
error('fsolve did not converge to a solution');
end
par.x0 = [par.F_feed ; par.x_H2_F ; par.x_N2_F ;par.x_NH3_F ;par.x_O2_F ;par.initial_values(1);par.initial_values(2); par.initial_values(3); par.initial_values(4) ];
par.tspan = 0:0.05:par.L_fiber;
par.M = [1 0 0 0 0 0 0 0 0 ; 0 1 0 0 0 0 0 0 0 ; 0 0 1 0 0 0 0 0 0 ; 0 0 0 1 0 0 0 0 0 ; 0 0 0 0 1 0 0 0 0 ; 0 0 0 0 0 0 0 0 0 ; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 ;0 0 0 0 0 0 0 0 0];
par.options = odeset('Mass',par.M,'RelTol',1e-6,'AbsTol',1e-6*ones(9,1));
[par.Z, par.X] = ode15s(@(z,x) dae(z,x,par),par.tspan,par.x0, par.options);
plot(par.Z,par.X(:,2:end))
function out = dae(z,x,par)
out = [- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P)))
1/x(1) *2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P) + x(2)*(- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P))))
1/x(1) *2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P) + x(3)*(- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P))))
1/x(1) *2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P) + x(4)*(- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P))))
1/x(1) *2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P) + x(5)*(- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P))))
x(6) - (par.Per_H2*x(2) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_H2*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))
x(7) - (par.Per_N2*x(3) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_N2*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))
x(8) - (par.Per_NH3*x(4) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_NH3*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))
x(9) - (par.Per_O2*x(5) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_O2*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))];
%x(1) = F
%x(2) = xH2
%x(3) = xN2
%x(4) = xHNH3
%x(5) = xO2
%x(6) = yH2
%x(7) = yN2
%x(8) = yNH3
%x(9) = yO2
end
function F = algebraic_constraints(x, par)
x=x.^2;
sum_terms = (x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2);
F = [
x(1) - par.Per_H2*par.x_H2_F * sum_terms / (1-par.BETA + par.BETA*par.Per_H2*sum_terms)
x(2) - par.Per_N2*par.x_N2_F * sum_terms / (1-par.BETA + par.BETA*par.Per_N2*sum_terms)
x(3) - par.Per_NH3*par.x_NH3_F * sum_terms / (1-par.BETA + par.BETA*par.Per_NH3*sum_terms)
x(4) - par.Per_O2*par.x_O2_F * sum_terms / (1-par.BETA + par.BETA*par.Per_O2*sum_terms)
];
end
6 Commenti
MeepMeep
il 9 Giu 2024
Hi Torsten!
I was wondering if there was a way to prevent the y_i values from becoming negative, or if this is simply a consequence of the program? I found a way to keep them positive by setting the membrane length (L_fiber) to maximum 0.3.
%% Gas separation in a cocurrent hollow fiber membrane
function par = Separation_gases()
par.thickness_membrane = 1e-6; % thickness membrane [m]
par.Perm_H2 = 250*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
par.Perm_N2 = 9*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
par.Perm_NH3 = 830*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
par.Perm_O2 = 37*3.35e-16;
par.R=8.314;
par.T = 273.15+25; % correlation temPeratures [K]
par.Per_H2 = par.Perm_H2/par.thickness_membrane; % Permeance of H2 [mol/m2 s Pa]
par.Per_N2 = par.Perm_N2/par.thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
par.Per_NH3 = par.Perm_NH3/par.thickness_membrane; % Permeance of N2 [mol/m2 s Pa]
par.Per_O2 = par.Perm_O2/par.thickness_membrane;
%T_g = 273.15+155; % glass transition temPerature of polymer membrane [K]
%% Input parameters
par.F_feed = 1133.95/3.6; % feed [mol/s]
par.x_H2_F = 0.0268322; % [-]
par.x_N2_F =0.96828; % [-]
par.x_NH3_F = 0.00475768; % [-]
par.x_O2_F = 0.000129826;
par.y_H2_init = par.Per_H2*par.x_H2_F/ (par.Per_H2*par.x_H2_F + par.Per_N2*par.x_N2_F + par.Per_NH3*par.x_NH3_F + par.Per_O2*par.x_O2_F);
par.y_N2_init = par.Per_N2*par.x_N2_F/ (par.Per_H2*par.x_H2_F + par.Per_N2*par.x_N2_F + par.Per_NH3*par.x_NH3_F + par.Per_O2*par.x_O2_F);
par.y_NH3_init = par.Per_NH3*par.x_NH3_F/ (par.Per_H2*par.x_H2_F + par.Per_N2*par.x_N2_F + par.Per_NH3*par.x_NH3_F + par.Per_O2*par.x_O2_F);
par.y_O2_init = par.Per_O2*par.x_O2_F/ (par.Per_H2*par.x_H2_F + par.Per_N2*par.x_N2_F + par.Per_NH3*par.x_NH3_F + par.Per_O2*par.x_O2_F);
par.P_f = 35e5; % pressure feed [Pa]
par.P_P = 1e5; % pressure Permeate [Pa]
%% Assumptions membrane
par.L_fiber = 0.3; %[m]
par.r_fiber = (0.75e-3)/2;
par.n_fiber = 4000;
par.visc_H2 = 88e-6;
par.visc_N2 = 17.82e-6; % [Pa s]
par.visc_NH3 = 9.9e-6;
par.visc_O2 = 20.4e-6;
par.mu = par.visc_N2; % most present
par.BETA = par.P_P/par.P_f;
% Initial guess for x(6) to x(9)
par.initial_guess = [par.y_H2_init; par.y_N2_init;par.y_NH3_init; par.y_O2_init]; % Example guess, you may need to adjust
%par.initial_guess = [0.1;0.015;0.1;0.015]; % Example guess, you may need to adjust
% Use fsolve to find consistent initial values
options = optimoptions('fsolve', 'Display', 'iter');
[par.initial_values, fval, exitflag] = fsolve(@(x) algebraic_constraints(x, par), par.initial_guess, options);
par.initial_values
% Check if fsolve was successful
if exitflag <= 0
error('fsolve did not converge to a solution');
end
par.x0 = [par.F_feed ; par.x_H2_F ; par.x_N2_F ;par.x_NH3_F ;par.x_O2_F ;par.initial_values(1);par.initial_values(2); par.initial_values(3); par.initial_values(4) ];
par.tspan = 0:0.01:par.L_fiber;
par.M = [1 0 0 0 0 0 0 0 0 ; 0 1 0 0 0 0 0 0 0 ; 0 0 1 0 0 0 0 0 0 ; 0 0 0 1 0 0 0 0 0 ; 0 0 0 0 1 0 0 0 0 ; 0 0 0 0 0 0 0 0 0 ; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 ;0 0 0 0 0 0 0 0 0 ];
par.options = odeset('Mass',par.M,'RelTol',1e-3,'AbsTol',1.e-3*ones(9,1));
[par.Z, par.X] = ode15s(@(z,x) dae(z,x,par),par.tspan,par.x0, par.options);
plot(par.Z,par.X(:,2:end))
end
function out = dae(z,x,par)
out = [- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P)))
1/x(1) *2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P) + x(2)*(- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P))))
1/x(1) *2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P) + x(3)*(- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P))))
1/x(1) *2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P) + x(4)*(- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P))))
1/x(1) *2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P) + x(5)*(- 2*3.14*par.r_fiber*par.L_fiber*par.n_fiber*((par.Per_H2*(x(2)*par.P_f + x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f + x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f + x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f + x(9)*par.P_P))))
x(6) - (par.Per_H2*x(2) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_H2*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))
x(7) - (par.Per_N2*x(3) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_N2*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))
x(8) - (par.Per_NH3*x(4) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_NH3*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))
x(9) - (par.Per_O2*x(5) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_O2*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))];
1- x(6) - x(7) - x(8) - x(9);
1- x(2) - x(3) - x(4) - x(5);
x(1) = max(x(1), 0);
x(2) = max(x(2), 0);
x(3) = max(x(3), 0);
x(4) = max(x(4), 0);
x(5) = max(x(5), 0);
x(6) = max(x(6), 0);
x(7) = max(x(7), 0);
x(8) = max(x(8), 0);
%x(1) = F
%x(2) = xH2
%x(3) = xN2
%x(4) = xHNH3
%x(5) = xO2
%x(6) = yH2
%x(7) = yN2
%x(8) = yNH3
%x(9) = yO2
end
function F = algebraic_constraints(x, par)
F = [
x(1) - (par.Per_H2*par.x_H2_F * ((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_H2*((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) )))
x(2) - (par.Per_N2*par.x_N2_F * ((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_N2*((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) )))
x(3) - (par.Per_NH3*par.x_NH3_F * ((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_NH3*((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) )))
x(4) - (par.Per_O2*par.x_O2_F * ((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_O2*((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) )))
];
end
I was wondering if there was a way to prevent the y_i values from becoming negative, or if this is simply a consequence of the program?
Given the xi, the yi are solutions of the algebraic system of equations you solved with "fsolve" to get consistent initial values. ode15s tries to follow the solution path for the yi given the xi. So I don't see a way to keep the yi positive because they just follow the xi. It could be that your algebraic system for the yi has several solutions so that the solution should abruptly make a jump once they get negative - that seems quite improbable, but I don't know.
You could try to turn the algebraic equations to differential equations by differentiating them with respect to z - but this will become quite complicated.
You could also try to compute the yi inside your function "dae" using "fsolve" and remove them as algebraic solution variables for ode15s, but I don't know whether this will change anything.
But before making substantial changes, I'd first check the program and the balance equations for possible errors.
Hi Torsten! I indeed found some issues in the code, and now it is working just fine! Thank you very much for all of your help
This is the working code:
Separation_gases()
%% Gas separation in a cocurrent hollow fiber membrane
function par = Separation_gases()
par.thickness_membrane = 10e-6; % thickness membrane [m]
par.Perm_H2 = 250*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
par.Perm_N2 = 9*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
par.Perm_NH3 = 830*3.35e-16; % Permeability H2 [mol m/m2 s Pa]
par.Perm_O2 = 37*3.35e-16;
par.R=8.314;
par.T = 273.15+25; % correlation temPeratures [K]
par.Per_H2 = 510e-10; % Permeance of H2 [mol/m2 s Pa]
par.Per_N2 = 10e-10; % Permeance of N2 [mol/m2 s Pa]
par.Per_NH3 = 25e-10; % Permeance of N2 [mol/m2 s Pa]
par.Per_O2 = 10e-10;
%T_g = 273.15+155; % glass transition temPerature of polymer membrane [K]
%% Input parameters
par.F_feed = (1133.95/3.6)/100; % feed [mol/s]
par.x_H2_F = 0.0268; % [-]
par.x_N2_F =0.9682; % [-]
par.x_NH3_F = 0.0049; % [-]
par.x_O2_F = 0.0001;
par.P_f = 35e5; % pressure feed [Pa]
par.P_P = 1e5; % pressure Permeate [Pa]
%% Assumptions membrane
par.A_range = 15; % [m2/(mol/s)]
par.A_module =par.A_range *par.F_feed; % effective area
par.l=4; %[m]
par.r=0.75e-3/2;
par.n = par.A_module/(2*3.14*par.r*par.l);
par.visc_H2 = 88e-6;
par.visc_N2 = 17.82e-6; % [Pa s]
par.visc_NH3 = 9.9e-6;
par.visc_O2 = 20.4e-6;
par.BETA = par.P_P/par.P_f;
% Initial guess for x(6) to x(9)
par.initial_guess = [1 2 3 4]; % Example guess, you may need to adjust
% Use fsolve to find consistent initial values
options = optimoptions('fsolve', 'Display', 'iter');
[par.initial_values, fval, exitflag] = fsolve(@(x) algebraic_constraints(x, par), par.initial_guess, options);
par.initial_values
% Check if fsolve was successful
if exitflag <= 0
error('fsolve did not converge to a solution');
end
par.x0 = [par.F_feed ; par.x_H2_F ; par.x_N2_F ;par.x_NH3_F ;par.x_O2_F ; par.initial_values(1) ;par.initial_values(2) ; par.initial_values(3);par.initial_values(4)];
par.zspan = 0:0.01:4;
par.M = [1 0 0 0 0 0 0 0 0 ; 0 1 0 0 0 0 0 0 0 ; 0 0 1 0 0 0 0 0 0 ; 0 0 0 1 0 0 0 0 0 ; 0 0 0 0 1 0 0 0 0 ; 0 0 0 0 0 0 0 0 0 ; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 ;0 0 0 0 0 0 0 0 0 ];
par.options = odeset('Mass',par.M,'RelTol', 1e-10, 'AbsTol', 1e-10);
[par.Z, par.X] = ode15s(@(z,x) dae(z,x,par),par.zspan,par.x0, par.options);
% Check mass balances
par.mass_balance_retentate =par.X(:,2) + par.X(:,3) + par.X(:,4)+ par.X(:,5) ;
par.mass_balance_permeate = par.X(:,6) + par.X(:,7) + par.X(:,8)+ par.X(:,9);
%% Retentate flows
par.F_H2_R = par.X(end,2)*par.X(end,1); % [-]
par.F_N2_R =par.X(end,3)*par.X(end,1); % [-]
par.F_NH3_R =par.X(end,4)*par.X(end,1); % [-]
par.F_O2_R = par.X(end,5)*par.X(end,1);
%% Permeate flows
par.F_H2_P = par.X(end,6)*(par.F_feed -par.X(end,1)); % [-]
par.F_N2_P =par.X(end,7)*(par.F_feed -par.X(end,1)); % [-]
par.F_NH3_P =par.X(end,8)*(par.F_feed -par.X(end,1)); % [-]
par.F_O2_P = par.X(end,9)*(par.F_feed -par.X(end,1));
%% Vector retentate
par.ret_flows = [par.F_H2_R; par.F_N2_R; par.F_NH3_R;par.F_O2_R];
%% Vector permeate
par.perm_flows = [par.F_H2_P; par.F_N2_P; par.F_NH3_P;par.F_O2_P];
figure(1)
plot(par.Z,par.X(:,2:5))
colororder([1 0 0; 0 1 0; 0 0 1; 0 1 1])
xlabel('Distance along membrane [m]')
ylabel('Molar fraction in retentate [-]')
legend('xH2','xN2', 'xNH3', 'xO2')
figure(2)
plot(par.Z,par.X(:,5:end))
colororder([1 0 0; 0 1 0; 0 0 1; 0 1 1])
xlabel('Distance along membrane [m]')
ylabel('Molar fraction in permeate [-]')
legend('yH2','yN2', 'yNH3', 'yO2')
%
figure(3)
plot(par.Z, par.X(:,1))
xlabel('Distance along membrane [m]')
ylabel('Flow rate feed [mol/s] [-]')
legend('Feed flow rate [mol/s')
end
function out = dae(z,x,par)
out = [- par.A_module*((par.Per_H2*(x(2)*par.P_f - x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f - x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f - x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f - x(9)*par.P_P)))
-1/x(1) *(par.A_module*par.Per_H2*(x(2)*par.P_f - x(6)*par.P_P) + x(2)*(- par.A_module*((par.Per_H2*(x(2)*par.P_f - x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f - x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f - x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f - x(9)*par.P_P)))))
-1/x(1) *(par.A_module*par.Per_N2*(x(3)*par.P_f - x(7)*par.P_P) + x(3)*(- par.A_module*((par.Per_H2*(x(2)*par.P_f - x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f - x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f - x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f - x(9)*par.P_P)))))
-1/x(1) *(par.A_module*par.Per_NH3*(x(4)*par.P_f - x(8)*par.P_P) + x(4)*(- par.A_module*((par.Per_H2*(x(2)*par.P_f - x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f - x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f - x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f - x(9)*par.P_P)))))
-1/x(1) *(par.A_module*par.Per_O2*(x(5)*par.P_f - x(9)*par.P_P) + x(5)*(- par.A_module*((par.Per_H2*(x(2)*par.P_f - x(6)*par.P_P)) + (par.Per_NH3*(x(4)*par.P_f - x(8)*par.P_P)) + (par.Per_N2*(x(3)*par.P_f - x(7)*par.P_P)) + (par.Per_O2*(x(5)*par.P_f - x(9)*par.P_P)))))
x(6) - (par.Per_H2*x(2) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_H2*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))
x(7) - (par.Per_N2*x(3) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_N2*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))
x(8) - (par.Per_NH3*x(4) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_NH3*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))
x(9) - (par.Per_O2*x(5) * ((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_O2*((x(6)/par.Per_H2) + (x(7)/par.Per_N2) + (x(8)/par.Per_NH3) + (x(9)/par.Per_O2) )))];
%x(1) = F
%x(2) = xH2
%x(3) = xN2
%x(4) = xHNH3
%x(5) = xO2
%x(6) = yH2
%x(7) = yN2
%x(8) = yNH3
%x(9) = yO2
end
function F = algebraic_constraints(x, par)
F = [
x(1) - (par.Per_H2*par.x_H2_F * ((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_H2*((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) )))
x(2) - (par.Per_N2*par.x_N2_F * ((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_N2*((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) )))
x(3) - (par.Per_NH3*par.x_NH3_F * ((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_NH3*((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) )))
x(4) - (par.Per_O2*par.x_O2_F * ((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) ) / (1-par.BETA + par.BETA*par.Per_O2*((x(1)/par.Per_H2) + (x(2)/par.Per_N2) + (x(3)/par.Per_NH3) + (x(4)/par.Per_O2) )))
];
end
Sachin kumar Ramayampet
il 20 Giu 2024
Hi!
Could you please provide the reference paper used for this mathematical modelling. I am also trying to model the binary gas separation using assymetric hollow fiber membrane. As the system of equations are hoighly Non linear and stiff ODEs, I have also faced convergence issues and got the y values as negatives similar to your case.
Thanking you in advance.
MeepMeep
il 20 Giu 2024
Sachin kumar Ramayampet
il 21 Giu 2024
Hi!
Thank you for sharing the reference. I have some other questions regarding the code. If you could answer then it would be a great help.
- While initializing the membrane parameters (see below ) you have used par.A_range and then par.A_module instead of taking 2 * pi * R * L * N (which is given in the reference paper). So would you please explain is their any particular reason for taking par.A_module in this way.
%% Assumptions membrane
par.A_range = 15; % [m2/(mol/s)]
par.A_module =par.A_range *par.F_feed; % effective area
par.l=4; %[m]
par.r=0.75e-3/2;
par.n = par.A_module/(2*3.14*par.r*par.l);
2. In the plot of "molar Fraction in permeate", the mole fraction of the NH3 is increasing over the length of fiber module over the other components eventhough the initial mole fraction of NH3 in feed as well as the permeance of NH3 are considerably low when compared with other components. Could you explain this behaviour? (I am asking this because I am not familiar with this Topic.)
3. If I am correct this code is for the cocurrent configuration right?
I hope you would clarify the above questions and Thank you in advance.
Più risposte (0)
Categorie
Scopri di più su Develop Apps Using App Designer in Centro assistenza e File Exchange
il 8 Giu 2024
il 21 Giu 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
