Dear Sir, Madam,
I’m having difficulties employing this simple model to perform analysis with UQlab. By any chance Could you identify any issues with the composition of this function?
Best regards,
Simon
function [H2tot, H2Flux, Recovery,xfeed_H2_all] = Script_function(X)
% Initialize constants that do not depend on Pfeed and T
% X(:, 1) corresponds to the feed pressure 'Pfeed' [Pa]
% X(:, 2) corresponds to the temperature 'T' [K]
Qin = 0.5; % Total volumetric flow rate entering the reactor[l/min]
Vm = 24.464; % Molar Volume[l/mol]
Pperm = 1E5; % Total pressure in the permeate[Pa]
xin_H2 = 1; % Molar fraction of H2 entering the reactor[-]
n = 0.749; % Membrane characteristic pressure exponent[-]
Ea = 7805.3; % Membrane activation energy[J/mol]
P0 = 4.373e-10; % Membrane permeability coefficient (per unit membrane thickness) [mol/s/Pa/m]
t = 2e-6; % Membrane total thickness (protective+selective layers)[m]
d = 0.014; % Membrane external diameter[m]
l = 0.08; % Membrane length[m]
R = 8.314; % Universal gas constant[J/mol/K]
lm = 0:0.00001:l; % Discretized membrane length [m]
% Preallocate arrays for results
H2tot = zeros(size(X, 1), 1);
H2Flux = zeros(size(X, 1), 1);
Recovery = zeros(size(X, 1), 1);
xfeed_H2_all = zeros(size(X, 1), length(lm));
for idx = 1:size(X, 1)
% Extract Pfeed and T from the input matrix
Pfeed = X(idx, 1);
T = X(idx, 2);
% Calculate dependent variables
Fin_tot = Qin / (Vm * 60); % Total molar flow rate entering the reactor [mol/s]
xin_i = 1 - xin_H2; % Molar fraction of i-species entering the reactor[-]
Fin_H2 = Fin_tot * xin_H2; % Molar flow rate of H2 entering the reactor[mol/s]
Fin_i = Fin_tot * xin_i; % Molar flow rate of i-species entering the reactor[mol/s]
% Calculate the membrane permeance
t_PD = t; % Membrane thickness of the selective layer [m]
Permeance_H2 = (P0 / t_PD) * exp(-Ea / (R * T)); % H2 permeance [mol/s/m/Pa^n]
% Discretization of the membrane length
da = (3.14 * d * l) / length(lm); % Discretized membrane area [m2]
% Initialize variables for the discretized membrane calculations
F_i = zeros(1, length(lm));
F_H2 = zeros(1, length(lm));
Fperm_H2 = zeros(1, length(lm));
Pperm_H2 = zeros(1, length(lm));
Pfeed_i = zeros(1, length(lm));
Pfeed_H2 = zeros(1, length(lm));
xfeed_i = zeros(1, length(lm));
xfeed_H2 = zeros(1, length(lm));
% Initial conditions at the first discretization point
F_i(1) = Fin_i;
F_H2(1) = Fin_H2;
xfeed_i(1) = xin_i;
xfeed_H2(1) = xin_H2;
% Perform calculations along the length of the membrane
for i = 1:length(lm)
xfeed_i(i) = F_i(i) / Fin_tot;
xfeed_H2(i) = (F_H2(i) - Fperm_H2(i)) / (Fin_tot - Fperm_H2(i));
Pfeed_i(i) = Pfeed * xfeed_i(i);
Pfeed_H2(i) = Pfeed * xfeed_H2(i);
if (Pfeed_H2(i)^n - Pperm_H2(i)^n) <= 0
Fperm_H2(i+1) = Fperm_H2(i);
else
Fperm_H2(i+1) = Permeance_H2 * da * (Pfeed_H2(i)^n - Pperm_H2(i)^n);
end
Pperm_H2(i+1) = Pperm;
F_i(i+1) = F_i(i); % Assuming no permeation of i-species
F_H2(i+1) = F_H2(i) - Fperm_H2(i+1);
end
% Calculate the total mole flow permeated for H2, flux, and recovery
H2tot(idx) = sum(Fperm_H2); % Total mole flow permeated [mol/s]
H2Flux(idx) = H2tot(idx) / (3.14 * d * l); % Total flux permeated [mol/s/m^2]
Recovery(idx) = sum(Fperm_H2) / (xin_H2 * Fin_tot); % Recovery rate [-]
xfeed_H2_all(idx, :) = xfeed_H2;
end
end