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?

[Torsten]

To avoid numerical problems that result from self-written code, use a professional MATLAB integrator, e.g. ode15s.

[MeepMeep]

Thank you very much for your quick reply, Torsten!

I tried to redo it as you said, and the solver still has issues solving it. The ode15s solver gives me the following warning:

Warning: Failure at t=0.000000e+00.  Unable to meet
integration tolerances without reducing the step size
below the smallest value allowed (7.905050e-323) at time
t. 
> In ode15s (line 662)
In Separation_gases (line 53

I tried to see if the initial values I provide the solver somehow gives a division by zero, but it does not seem so. I tried other initial conditions, but they keep giving me the same error. I can not change the initial conditions for x_i and F, as these are coming from another equipment. Do you know what could have gone wrong? This is my code:

%% 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 = 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;

par.x0 = [par.F_feed ; par.x_H2_F ; par.x_N2_F ;par.x_NH3_F ;par.x_O2_F ; par.y_H2_init; par.y_N2_init ; par.y_NH3_init ; par.y_O2_init]
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-4,'AbsTol',[1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4]);
[par.Z, par.X] = ode15s(@(z,x) dae(z,x,par),par.tspan,par.x0, par.options);

disp('Z:');
    disp(par.Z);
    disp('X:');
    disp(par.X);
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))))
        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) ))
        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) ))
        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) ))
        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

[Torsten]

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]

Hi Torsten

Once again a big thank you for responding so quickly!

I changed the last equations and implemented fsolve (see my code below). I was wondering if there was also a way in which I could specify how to make sure these initial values are positive? Right now the program works, but along the membrane length the y_i values become negative, which should not physically be possible.

Thank you once again!

% 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

% 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.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',1.e-6*ones(9,1));
[par.Z, par.X] = ode15s(@(z,x) dae(z,x,par),par.tspan,par.x0, par.options);

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)
    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 * ((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

Answer 1

[Torsten]

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);

                                             Norm of      First-order   Trust-region
 Iteration  Func-count     ||f(x)||^2           step       optimality         radius
     0          5             340.232                             131              1
     1         10             136.345              1             59.4              1
     2         15             8.35211        2.02452              7.4            2.5
     3         20             0.49738        1.05013            0.918           5.06
     4         25            0.028113       0.556721            0.113           5.06
     5         30          0.00165936       0.300089           0.0166           5.06
     6         35         0.000106594       0.155651          0.00514           5.06
     7         40         5.30496e-06      0.0637521          0.00054           5.06
     8         45         2.50833e-07      0.0226726         2.91e-05           5.06
     9         50         5.92523e-09     0.00864321         6.21e-06           5.06
    10         55         1.41591e-11      0.0019109         3.04e-07           5.06
    11         60         1.22041e-16    0.000103533         8.92e-10           5.06
    12         65         1.01353e-26    3.05755e-07         8.32e-15           5.06

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

par.initial_values=par.initial_values.^2;
par.initial_values

ans = 1x4

0.3089 0.5876 0.1031 0.0003

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

sum(par.initial_values)

ans = 1.0000

% 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);

Warning: Failure at t=7.224374e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t.

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

[MeepMeep]

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

[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?

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.

[MeepMeep]

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()

                                             Norm of      First-order   Trust-region
 Iteration  Func-count     ||f(x)||^2           step       optimality         radius
     0          5             40.1484                            6.65              1
     1         10             25.8796              1             5.37              1
     2         15             3.64956            2.5             1.91            2.5
     3         20           0.0108228        1.67017           0.0875           6.25
     4         25          0.00204699       0.906251           0.0397           6.25
     5         30         0.000448611       0.395039           0.0206           6.25
     6         35         1.91812e-05       0.143139          0.00431           6.25
     7         40         2.62624e-08      0.0258282          0.00016           6.25
     8         45         3.47022e-14    0.000873717         1.84e-07           6.25

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

ans = 1x4

0.4329 0.5601 0.0069 0.0001

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = struct with fields:
        thickness_membrane: 1.0000e-05
                   Perm_H2: 8.3750e-14
                   Perm_N2: 3.0150e-15
                  Perm_NH3: 2.7805e-13
                   Perm_O2: 1.2395e-14
                         R: 8.3140
                         T: 298.1500
                    Per_H2: 5.1000e-08
                    Per_N2: 1.0000e-09
                   Per_NH3: 2.5000e-09
                    Per_O2: 1.0000e-09
                    F_feed: 3.1499
                    x_H2_F: 0.0268
                    x_N2_F: 0.9682
                   x_NH3_F: 0.0049
                    x_O2_F: 1.0000e-04
                       P_f: 3500000
                       P_P: 100000
                   A_range: 15
                  A_module: 47.2479
                         l: 4
                         r: 3.7500e-04
                         n: 5.0157e+03
                   visc_H2: 8.8000e-05
                   visc_N2: 1.7820e-05
                  visc_NH3: 9.9000e-06
                   visc_O2: 2.0400e-05
                      BETA: 0.0286
             initial_guess: [1 2 3 4]
            initial_values: [0.4329 0.5601 0.0069 5.7850e-05]
                        x0: [9x1 double]
                     zspan: [0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 0.1300 0.1400 0.1500 0.1600 0.1700 0.1800 ... ] (1x401 double)
                         M: [9x9 double]
                   options: [1x1 struct]
                         Z: [401x1 double]
                         X: [401x9 double]
    mass_balance_retentate: [401x1 double]
     mass_balance_permeate: [401x1 double]
                    F_H2_R: 4.4689e-04
                    F_N2_R: 2.4115
                   F_NH3_R: 0.0088
                    F_O2_R: 2.4907e-04
                    F_H2_P: 0.0028
                    F_N2_P: 0.7197
                   F_NH3_P: 0.0063
                    F_O2_P: 7.4336e-05
                 ret_flows: [4x1 double]
                perm_flows: [4x1 double]

%% 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]

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]

Hi!

I used this paper as reference. Hope this will help you!

[Sachin kumar Ramayampet]

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.

1. 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.