Passing out additions parameters after ODE solver.

I am trying to pass pass out the velocity array at the end of adsorption model. I would like to create a matrix of time and velocity to plot a graph. with the other graphs where their parameters are from the ODE solver. Here is the code with some of the irrelavent parts taken out:

clear,clc
%% MAIN
% Physical Parameters
L   = 1;            % Column length m 
r   = 0.25;         % Bed radius m 

t0  = 0;            % Initial Time
tf  = 2;           % Final time
dt  = 0.1;          % Time step
t   = t0:dt:tf;     % Time vector
dz  = 0.05;          % Mesh size
z   = 0:dz:L;       % Mesh generation
n   = numel(z);

%%%%%%% blended gas parameters %%%%%%%%%%%
%%% calc u input

TPF = [1 300 8e5 2e5 1e5];   % Feed: Velocity (m/s), tempurature (K) and pressure (Bar) [Vfeed Tfeed PH PI PL] 

yF  = 0.7;              % Mole fraction for methane 

% y is an array size n*5 of y1 = 1:n, q1 = n+1:2*n,  
% q2 = 2*n+1:3*n, T = 3*n+1:4*n, P = 4*n+1:5*n
sol = adsorptionSolver(t,z,yF,TPF);

Unrecognized function or variable 'uhalf'.

Error in solution>adsorptionModel (line 115)
Velocity=cat(1,uhalf,u);

Error in solution>adsorptionModelODE (line 106)
[dydt, ~] = adsorptionModel(t, y, h, n, yiFeed, TPFeed);

Error in solution>@(t,y)adsorptionModelODE(t,y,h,n,yFeed,TPFeed) (line 102)
out = ode15s(@(t,y) adsorptionModelODE(t,y,h,n,yFeed,TPFeed),t,y0);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 148)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in solution>adsorptionSolver (line 102)
out = ode15s(@(t,y) adsorptionModelODE(t,y,h,n,yFeed,TPFeed),t,y0);

% sol.x is time steps

purityh = 100*Moley2(n,end) / (sol.y(n,end) + Moley2(n,end));
fprintf('Purity of hydrogen is %f%% after %4.2fs\n', purityh, tf)

    figure(1)
    
    subplot(1,2,1)
    mesh(sol.x,z,sol.y(1:n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('mole fraction y1')
    title('mole fraction of methane')
    
    subplot(1,2,2)
    mesh(sol.x,z,Moley2)
    xlabel('t')
    ylabel('bed length')
    zlabel('mole fraction y2')
    title('mole fraction of hydrogen')
    
     figure(2)
    
    subplot(1,2,1)
    mesh(sol.x,z,sol.y(n+1:2*n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('loading q1 mol/kg')
    title('loading of methane')
    
    subplot(1,2,2)
    mesh(sol.x,z,sol.y(2*n+1:3*n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('loading q2 mol/kg')
    title('loading of hydrogen')
    
    figure(3)
    
    subplot(1,2,1)
    mesh(sol.x,z,sol.y(3*n+1:4*n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('Temp')
    title('temp of system')
    subplot(1,2,2)
    mesh(sol.x,z,sol.y(4*n+1:5*n,1:end))
    xlabel('time')
    ylabel('bed length')
    zlabel('Pressure')
    title('Pressure of system')
    
    %% Solver function
function out = adsorptionSolver(t,z,yFeed,TPFeed)
n = numel(z);      % Size of mesh grid
h = diff(z);
h = h(1);
% start up condidtions
Tw = 300;                   % Ambient Tempurature K
Pw = 8e5;                % Ambient Pressure Pa
y1w = 0.7;

% Start up conditions Conditions
y1_0       = zeros(n,1) + y1w;     

q1_0        = zeros(n,1);
q2_0        = zeros(n,1);

T_0         = zeros(n,1) + Tw;

P_0         = zeros(n,1) + Pw;

% y is an array size n*4 of c1 = 1:n, q1 = n+1:2*n, c2 = 2*n+1:3*n, q2 = 3*n+1:4*n
y0      = [y1_0; q1_0; q2_0; T_0; P_0];

% % Solving
out = ode15s(@(t,y) adsorptionModelODE(t,y,h,n,yFeed,TPFeed),t,y0);

end

Ive tried to use another function to separate the two adsorption model fucntion outputs but i havent got it to work and it seems inefficient

function dydt = adsorptionModelODE(t, y, h, n, yiFeed, TPFeed)
    [dydt, ~] = adsorptionModel(t, y, h, n, yiFeed, TPFeed);
end

%% Adsorption model
function [dydt, Velocity] = adsorptionModel(t,y,h,n,yiFeed,TPFeed)

% t is used to calculate boundary pressure during a few cycle steps 
% rest of code and odes

Velocity=cat(1,uhalf,u);
Velocity([1:2:end,2:2:end])=Velocity;

%%%%%%%%%%%%%%%% Concatenate vectors of time derivatives
dydt = [dy1dt;dq1dt;dq2dt;dTdt;dPdt];
end

Ive seen various forum pages but being a beginner they have confused me. I think persistent is the best way to implement it but ive tried and failed.

Thanks in advance,

Tom

[Dyuman Joshi]

There is an undefined parameter in your code, see the edit above.

For the function adsorptionModel the only the input u is used to define the output, so why include other inputs in the function definition?

On a quick glance at your code, I don't see any need to use peristent variables.

[Thomas Stone-Wigg]

Sorry i didnt want to overwhelm this page with the whole code:

clear,clc
%% MAIN
% Physical Parameters
L   = 1;            % Column length m 
r   = 0.25;         % Bed radius m 

t0  = 0;            % Initial Time
tf  = 2;           % Final time
dt  = 0.1;          % Time step
t   = t0:dt:tf;     % Time vector
dz  = 0.05;          % Mesh size
z   = 0:dz:L;       % Mesh generation
n   = numel(z);

%%%%%%% blended gas parameters %%%%%%%%%%%
%%% calc u input

TPF = [1 300 8e5 2e5 1e5];   % Feed: Velocity (m/s), tempurature (K) and pressure (Bar) [Vfeed Tfeed PH PI PL] 

yF  = 0.7;              % Mole fraction for methane 

%%%%%%% Preliminary calculations %%%%%%%%%%%

A = pi*r^2;         % Bed cross-sectional area m2
V = A*L;            % Bed volume m^3
mBed = A*L*426.7;   % Mass of adsorbent (in bed) 426.7 = bed density
f=1;
%%%%%%%%% loop around here to work out p and n and mole fraction for each
%%%%%%%%% loop

% y is an array size n*5 of y1 = 1:n, q1 = n+1:2*n,  
% q2 = 2*n+1:3*n, T = 3*n+1:4*n, P = 4*n+1:5*n
sol = adsorptionSolver(t,z,yF,TPF)
% sol.x is time steps
% sol.y is 
Moley2 = 1-sol.y(1:n,1:end);

format long

purityh = 100*Moley2(n,end) / (sol.y(n,end) + Moley2(n,end));
fprintf('Purity of hydrogen is %f%% after %4.2fs\n', purityh, tf)

% purityh = sol.y(3*n,:) ./(sol.y(n,:) + sol.y(3*n,:));
%purityh = (sol.y(n,:) + sol.y(3*n,:)); % mole fraction
% figure
% plot(sol.x, purityh)
% xlabel('time')
% ylabel('purity of hydrogen')
 
 
 
    figure(1)
    
    subplot(1,2,1)
    mesh(sol.x,z,sol.y(1:n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('mole fraction y1')
    title('mole fraction of methane')
    
    subplot(1,2,2)
    mesh(sol.x,z,Moley2)
    xlabel('t')
    ylabel('bed length')
    zlabel('mole fraction y2')
    title('mole fraction of hydrogen')
    
     figure(2)
    
    subplot(1,2,1)
    mesh(sol.x,z,sol.y(n+1:2*n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('loading q1 mol/kg')
    title('loading of methane')
    
    subplot(1,2,2)
    mesh(sol.x,z,sol.y(2*n+1:3*n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('loading q2 mol/kg')
    title('loading of hydrogen')
    
    figure(3)
    
    subplot(1,2,1)
    mesh(sol.x,z,sol.y(3*n+1:4*n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('Temp')
    title('temp of system')
    subplot(1,2,2)
    mesh(sol.x,z,sol.y(4*n+1:5*n,1:end))
    xlabel('time')
    ylabel('bed length')
    zlabel('Pressure')
    title('Pressure of system')
    

%% Solver function
function out = adsorptionSolver(t,z,yFeed,TPFeed)
n = numel(z);      % Size of mesh grid
h = diff(z);
h = h(1);
% start up condidtions
Tw = 300;                   % Ambient Tempurature K
Pw = 8e5;                % Ambient Pressure Pa
y1w = 0.7;

% Start up conditions Conditions
y1_0       = zeros(n,1) + y1w;     

q1_0        = zeros(n,1);
q2_0        = zeros(n,1);

T_0         = zeros(n,1) + Tw;

P_0         = zeros(n,1) + Pw;

% y is an array size n*4 of c1 = 1:n, q1 = n+1:2*n, c2 = 2*n+1:3*n, q2 = 3*n+1:4*n
y0      = [y1_0; q1_0; q2_0; T_0; P_0];

% % % Creating mass matrix
% M               = eye(5*n);                         % Mass matrix for temporal domain integrator
% M(1,1)          = 0;                                % Algebraic equation for yi1 at z = 0;
% M(n,n)          = 0;                                % Algebraic equation for yi1 at z = L
% M(2*n+1,2*n+1)  = 0;                                % Algebraic equation for yi2 at z = 0;
% M(3*n,3*n)      = 0;                                % Algebraic equation for yi2 at z = L
% M(4*n+1,4*n+1)  = 0;
% opts = odeset('Mass',M,'MassSingular','yes');       % creates mass matrix for the function ode15s

% % Solving


out = ode15s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0);

%[Velocity, ~] = adsorptionModel(t,y,h,n,yFeed,TPFeed)
% Solving
%out = ode23s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0);
end
%% Adsorption model
function [dydt, Velocity] = adsorptionModel(t,y,h,n,yiFeed,TPFeed)

% Variables allocation
% y is an array size n*5 of y1 = 1:n, q1 = n+1:2*n,  
% q2 = 2*n+1:3*n, T = 3*n+1:4*n, P = 4*n+1:5*n
y1 = y(1:n);
y2 = 1 - y1;

q1 = y(n+1:2*n);
q2 = y(2*n+1:3*n);

T  = y(3*n+1:4*n);
P  = y(4*n+1:5*n);

% Z half points
yhalf = zeros(n+1,1);
Thalf = zeros(n+1,1);
Phalf = zeros(n+1,1);
uhalf = zeros(n+1,1);
u = zeros(n,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Constant physical properties and basic simulation parameters
R = 8.3145;                 % ideal gas constant J/mol/K 
D   = 1.3e-5;               % Axial Dispersion coefficient m2/s
epl = 0.4043;               % Void fraction of bed
eplp = 0.546;               % Void fraction of particle
eplt = epl + (1-epl)*eplp;  % Total void fraction 
k   = [0.136;0.259];        % (lumped) Mass Transfer Coefficient 1/s
rhop = 716.3;               % Particle density kg/m3
% rhob = 426.7;               % Bed density kg/m3 = (1-epl)*rhop
Cps = 1046.7;               % heat capacity of solid J/kg/K 
Kl = 1.2e-6;                % Thermal diffusity J/m/s/K
Kz = 0.09;                   % Effective gas thermal conductivity 
deltaH = [24124; 8420];     % heat of adsorption J/mol
% Vfeed = TPFeed(1);
Tfeed = TPFeed(2);
PH = TPFeed(3); 
% PI = TPFeed(4);
PL = TPFeed(5);
lambda = 0.5;                % rate of pressure change (1/s)

% Ergun equation parameters
visc = 3.73e-8;     % gas viscosity kg/m/s
Rp = 5.41e-3;        % particle radius m

% Constants for CH4
c1a = -0.703029; c1b = 108.4773; c1c = -42.52157; c1d = 5.862788; c1e = 0.678565;
% Constants for H2
c2a = 33.066178; c2b = -11.363417; c2c = 11.432816; c2d = -2.772874; c2e = -0.158558;


% Cpg = heat capacity (J/mol*K)
Cpg1 = c1a + c1b*(T/1e3) + c1c*(T/1e3).^2 + c1d*(T/1e3).^3 + c1e./((T/1e3).^2);
Cpg2 = c2a + c2b*(T/1e3) + c2c*(T/1e3).^2 + c2d*(T/1e3).^3 + c2e./((T/1e3).^2);
Cpg = (y1.*Cpg1 + y2.*Cpg2);
% cool prop for viscosity eqs

%%%%%%%% Boundary conditions: j = 0.5 and n + 0.5 %%%%%%%%%%%%%%%%%%%%%%%%    
    %% Pressurization
    
% Phalf(1) = PH - (PH - PL)*exp(-lambda*t);
% Phalf(n+1) = P(n);
% 
% uhalf(1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((P(1)-Phalf(1)));
% %vhalf(n+1) = -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*(Phalf(n+1)-P(n));
% uhalf(n+1) = 0;
% 
% %yhalf(1) = 0.7;
% yhalf(1) = (y1(1)+yiFeed*uhalf(1)*h/D/2) / (1+uhalf(1)*h/D/2);
% %yhalf(1) = uhalf(1)*h/(2*D)*(yiFeed-y1(1)) + y1(1)
% 
% yhalf(n+1) = y1(n);
% 
%     % Inlet gas condtions
%     rhogInlet = ((Phalf(1)) .* (yhalf(1)*16.04e-3 + (1-yhalf(1))*2.02e-3)) ./ (R*Tfeed);
%     Cpg1in = c1a + c1b*(Tfeed/1e3) + c1c*(Tfeed/1e3).^2 + c1d*(Tfeed/1e3).^3 + c1e./((Tfeed/1e3).^2);
%     Cpg2in = c2a + c2b*(Tfeed/1e3) + c2c*(Tfeed/1e3).^2 + c2d*(Tfeed/1e3).^3 + c2e./((Tfeed/1e3).^2);
%     CpgInlet = (yhalf(1).*Cpg1in + (1-yhalf(1)).*Cpg2in);
%     %
% 
% Thalf(1) = (T(1) + Tfeed*uhalf(1)*rhogInlet*CpgInlet*h/Kl/2)/(1 + uhalf(1)*rhogInlet*CpgInlet*h/Kl/2);
% Thalf(n+1) = T(n);
    
    %% Adsorption

    uhalf(1) = TPFeed(1);
    %Phalf(1) = PH;
    Phalf(1)  =  P(1) + (uhalf(1)*h*visc/2) / (4/150*(epl/(1-epl))^2*Rp^2);
    Phalf(n+1) = PH;
    uhalf(n+1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((Phalf(n+1)-P(n)));
    yhalf(1) = (y1(1)+yiFeed*uhalf(1)*h/D/2) / (1+uhalf(1)*h/D/2);
    yhalf(n+1) = y1(n);
    
     % Inlet gas condtions
        rhogInlet = ((Phalf(1)) .* (yhalf(1)*16.04e-3 + (1-yhalf(1))*2.02e-3)) ./ (R*Tfeed);
        Cpg1in = c1a + c1b*(Tfeed/1e3) + c1c*(Tfeed/1e3).^2 + c1d*(Tfeed/1e3).^3 + c1e./((Tfeed/1e3).^2);
        Cpg2in = c2a + c2b*(Tfeed/1e3) + c2c*(Tfeed/1e3).^2 + c2d*(Tfeed/1e3).^3 + c2e./((Tfeed/1e3).^2);
        CpgInlet = (yhalf(1).*Cpg1in + (1-yhalf(1)).*Cpg2in);
        %

    Thalf(1) = (T(1) + Tfeed*uhalf(1)*rhogInlet*CpgInlet*h/Kl/2)/(1 + uhalf(1)*rhogInlet*CpgInlet*h/Kl/2);
    Thalf(n+1) = T(n);

    %% Depres

%     Phalf(1) = P(1);
%     Phalf(n+1) = PI + (PH - PI)*exp(-lambda*t);
%     
%     uhalf(1) = 0;
%     %uhalf(1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((P(1)-Phalf(1)));
%     uhalf(n+1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((Phalf(n+1)-P(n)));
%     
%     yhalf(1) = y1(1);
%     yhalf(n+1) = y1(n);
%     
%     Thalf(1) = T(1);
%     Thalf(n+1) = T(n);
    %% Purge
%     Phalf(1) = PL + (PI - PL)*exp(-lambda*t);
%     Phalf(n+1) = P(n);
%     
%     uhalf(1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((P(1)-Phalf(1)));
%     uhalf(n+1) = 0; 
%     %uhalf(n) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((Phalf(n+1)-P(n)));
%     
%     yhalf(1) = y1(1);
%     yhalf(n+1) = y1(n);
%     
%     Thalf(1) = T(1);
%     Thalf(n+1) = T(n);
  
%%%%%%%% Velocity wall values: j 1.5 to n - 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
uhalf(2:n) =   -(1/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((P(2:n) - P(1:n-1)));

%%%%%%%% Wall values: j = 1.5 %%%%%%%%%%%%%%%%%%%%%%%%
alpha0P = (2/3) / (P(2)-P(1) + 1e-10)^4;
alpha1P = (1/3) / (2*(P(1)-Phalf(1)) + 1e-10)^4;
Phalf(2) = alpha0P/(alpha0P+alpha1P)*(0.5*(P(1)+P(2))) + alpha1P/(alpha0P+alpha1P)*(2*P(1)-Phalf(1));

alpha0y = (2/3) / (y1(2)-y1(1) + 1e-10)^4;
alpha1y = (1/3) / (2*(y1(1)-yhalf(1)) + 1e-10)^4;
yhalf(2) = alpha0y/(alpha0y+alpha1y)*(0.5*(y1(1)+y1(2))) + alpha1y/(alpha0y+alpha1y)*(2*y1(1)-yhalf(1));

alpha0T = (2/3) / (T(2)-T(1) + 1e-10)^4;
alpha1T = (1/3) / (2*(T(1)-Thalf(1)) + 1e-10)^4;
Thalf(2) = alpha0T/(alpha0T+alpha1T)*(0.5*(T(1)+T(2))) + alpha1T/(alpha0T+alpha1T)*(2*T(1)-Thalf(1));

%%%%%%%% Wall values: j = 2.5 to n - 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
alpha0P = (2/3) ./ (P(3:n)-P(2:n-1) + 1e-10).^4;
alpha1P = (1/3) ./ (P(2:n-1)-P(1:n-2) + 1e-10).^4;
Phalf(3:n) = alpha0P./(alpha0P+alpha1P).*(0.5.*(P(2:n-1)+P(3:n))) + alpha1P./(alpha0P+alpha1P).*(1.5.*P(2:n-1)- 0.5* P(1:n-2));

alpha0y = (2/3) ./ (y1(3:n)-y1(2:n-1) + 1e-10).^4;
alpha1y = (1/3) ./ (y1(2:n-1)-y1(1:n-2) + 1e-10).^4;
yhalf(3:n) = alpha0y./(alpha0y+alpha1y).*(0.5*(y1(2:n-1)+y1(3:n))) + alpha1y./(alpha0y+alpha1y).*(1.5.*y1(2:n-1)- 0.5* y1(1:n-2));

alpha0T = (2/3) ./ (T(3:n)-T(2:n-1) + 1e-10).^4;
alpha1T = (1/3) ./ (T(2:n-1)-T(1:n-2)+ 1e-10).^4;
Thalf(3:n) = alpha0T./(alpha0T+alpha1T).*(0.5*(T(2:n-1)+T(3:n))) + alpha1T./(alpha0T+alpha1T).*(1.5.*T(2:n-1)- 0.5* T(1:n-2));



% rhog = (P .* (y1*16.04e-3 + y2*2.02e-3)) ./ (R*T);  % ρ = PM/RT, M: H2 = 2.02 g/mol CH4 = 16.04 g/mol Air= 28.97 g/mol 50/50 H2 and CH4 = 9.03
% % Central velocity eqautions
u(1:n) =   -(1/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((Phalf(2:n+1) - Phalf(1:n)));
rhog = (P .* (y1*16.04e-3 + y2*2.02e-3)) ./ (R*T);  % ρ = PM/RT, M: H2 = 2.02 g/mol CH4 = 16.04 g/mol Air= 28.97 g/mol 50/50 H2 and CH4 = 9.03
%%%%%%%%%%%%%%% Langmuir Equation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % CH4
a11 = -9.5592;	a21 = 4638.5;	b11 = 3.725e-4/1e5;	b21 = 1.443e4;
% % H2
a12 = -23.131;	a22 = 8069.1;	b12 = 2.248/1e5;	b22 = -1.435e4;
qs1 = a11 + (a21./T);
qs2 = a12 + (a22./T);
B1 = b11.*exp(b21./(8.3145*T));
B2 = b12.*exp(b22./(8.3145*T));
qstar1 = (qs1.*B1.*P.*y1) ./ (1+((P.*y1.*B1)+(P.*(1-y1).*B2))); %mols/kg
qstar2 = (qs2.*B2.*P.*(1-y1)) ./ (1+((P.*y1.*B1)+(P.*(1-y1).*B2)));

LT1 = qstar1-q1;
LT2 = qstar2-q2;

%%%%%%%%%%%%%%% Eqautions for loading of component i %%%%%%%%%%%%%%
dq1dt = k(1)*LT1;
dq2dt = k(2)*LT2;

%%%%%%%%%%%%%% Balance equations

dPdt = zeros(n,1);
dy1dt = zeros(n,1);
dTdt = zeros(n,1);

% component mass balance
dy1dt(1) = D/h.*(((y1(2)-y1(1))/h) - 2/h*(y1(1)-yhalf(1))) - u(1)/h.*(yhalf(2)-yhalf(1)) - R*T(1)./P(1).*((1-epl)/epl).*rhop.*(dq1dt(1) - y1(1).*(dq1dt(1)+dq2dt(1)));
dy1dt(2:n-1) = D/h.*((y1(3:n)-y1(2:n-1))/h - (y1(2:n-1)-y1(1:n-2))/h) - u(2:n-1)/h.*(yhalf(3:n)-yhalf(2:n-1)) - R*T(2:n-1)./P(2:n-1).*((1-epl)/epl).*rhop.*(dq1dt(2:n-1) - y1(2:n-1).*(dq1dt(2:n-1)+dq2dt(2:n-1)));
dy1dt(n) = D/h.*(2/h*(yhalf(n+1)-y1(n))-(y1(n)-y1(n-1))/h) - u(n)/h.*(yhalf(n+1)-yhalf(n)) - R*T(n)./P(n).*((1-epl)/epl).*rhop.*(dq1dt(n) - y1(n).*(dq1dt(n)+dq2dt(n)));

% Energy balance
dTdt(1) = 1./(eplt.*rhog(1).*Cpg(1)+rhop*Cps) .* (epl.*Kl/h.*(((T(2)-T(1))/h)-2/h*(T(1)-Thalf(1))) - rhog(1).*Cpg(1).*epl.*u(1)/h.*(Thalf(2)-Thalf(1)) + rhop.*(deltaH(1)*dq1dt(1)+deltaH(2)*dq2dt(1)));
dTdt(2:n-1) = 1./(eplt.*rhog(2:n-1).*Cpg(2:n-1)+rhop*Cps) .* (epl.*Kl/h.*(((T(3:n)-T(2:n-1))/h)-2/h*(T(2:n-1)-T(1:n-2))) - rhog(2:n-1).*Cpg(2:n-1).*epl.*u(2:n-1)/h.*(Thalf(3:n)-Thalf(2:n-1)) + rhop.*(deltaH(1)*dq1dt(2:n-1)+deltaH(2)*dq2dt(2:n-1)));
dTdt(n) = 1./(eplt.*rhog(n).*Cpg(n)+rhop*Cps) .* (epl.*Kl/h.*(2/h*(Thalf(n+1)-T(n))-(T(n)-T(n-1))/h) - rhog(n).*Cpg(n).*epl.*u(n)/h.*(Thalf(n+1)-Thalf(n)) + rhop.*(deltaH(1)*dq1dt(n)+deltaH(2)*dq2dt(n)));

% total mass balance
dPdt(1) = D/h.*( ((P(2)-P(1))/h)-2/h*(P(1)-Phalf(1))) - P(1)./h.*(uhalf(2)-uhalf(1)) - u(1)/h.*(Phalf(2)-Phalf(1)) - P(1)./T(1).*( D/h.*(((T(2)-T(1))/h)-2/h*(T(1)-Thalf(1))) + dTdt(1) + u(1)/h.*(Thalf(2)-Thalf(1)) ) - ((1-epl)/epl)*rhop*R.*T(1).*(dq1dt(1)+dq2dt(1));
dPdt(2:n-1) = D/h.*( (P(3:n)-P(2:n-1))/h - (P(2:n-1)-P(1:n-2))/h) - P(2:n-1)./h.*(uhalf(3:n)-uhalf(2:n-1)) - u(2:n-1)/h.*(Phalf(3:n)-Phalf(2:n-1)) - P(2:n-1)./T(2:n-1).*( D/h.*(((T(3:n)-T(2:n-1))/h)-2/h*(T(2:n-1)-T(1:n-2))) + dTdt(2:n-1) + u(2:n-1)/h.*(Thalf(3:n)-Thalf(2:n-1)) ) - ((1-epl)/epl)*rhop*R.*T(2:n-1).*(dq1dt(2:n-1)+dq2dt(2:n-1));
dPdt(n) = D/h.*(2/h*(Phalf(n+1)-P(n))-(P(n)-P(n-1))/h) - P(n)./h.*(uhalf(n+1)-uhalf(n)) - u(n)/h.*(Phalf(n+1)-Phalf(n)) - P(n)./T(n).*( D/h.*(2/h*(Thalf(n+1)-T(n))-(T(n)-T(n-1))/h) + dTdt(n) + u(n)/h.*(Thalf(n+1)-Thalf(n)) ) - ((1-epl)/epl)*rhop*R.*T(n).*(dq1dt(n)+dq2dt(n));


Velocity=cat(1,uhalf,u);
Velocity([1:2:end,2:2:end])=Velocity;

%%%%%%%%%%%%%%%% Concatenate vectors of time derivatives
dydt = [dy1dt;dq1dt;dq2dt;dTdt;dPdt];
end

Hopefully this will show what im trying to get at with my orginial question

[Torsten]

Hopefully this will show what im trying to get at with my orginial question

No. I can only guess: You want to make a surface plot of "purityh" ?

[Thomas Stone-Wigg]

Sorry, i want to make a mesh plot of velocity, time and position in the bed.

Answer 1

[Torsten]

sol = adsorptionSolver(t,z,yF,TPF)

for i = 1:size(sol.y,2)
  [~,velocity(:,i)] = adsorptionModel(sol.x(i),sol.y(:,i),z(2)-z(1),numel(z),yF,TPF);  
end

Since "velocity" is an array of size 43x61 and not 21x61, I don't know how the 43 should be interpreted.

[Thomas Stone-Wigg]

Thank you Torsten and Dyuman thats what i was trying to do, sorry about my confusion.

Ill use the nodal values of velocity or manually add the z half points so i can include the velocity wall values in the plot.

Many thanks again

[Thomas Stone-Wigg]

Hi @Torsten,

Thanks for the help either, ive got the code working but have run into a problem i cant seem to fix. When when i run the simulation for depressurization there seems to be a large pressure gradient at the end of the bed and i cant work out why.

Running the code with with t being output i notice the solver gets to the finishing time t=2 and then goes back to 0 and runs a few more loops before finishing at 2 again; These points correspond to the large gradients. Would you have any suggestions as to why this is happening?

Ive checked the balance equtions and WENO quite a few times and they seem to be fine. Thanks again

clear,clc
%% MAIN
% Physical Parameters
L   = 1;            % Column length m 
r   = 0.25;         % Bed radius m 

t0  = 0;            % Initial Time
tf  = 2;           % Final time
%dt  = 0.01;          % Time step
t   = [t0 tf];     % Time vector
dz  = 0.05;          % Mesh size
z   = 0:dz:L;       % Mesh generation
n   = numel(z);

%%%%%%% blended gas parameters %%%%%%%%%%%
%%% calc u input

TPF = [1 300 8e5 2e5 1e5];   % Feed: Velocity (m/s), tempurature (K) and pressure (Bar) [Vfeed Tfeed PH PI PL] 

yF  = 0.7;              % Mole fraction for methane 

%%%%%%% Preliminary calculations %%%%%%%%%%%

A = pi*r^2;         % Bed cross-sectional area m2
V = A*L;            % Bed volume m^3
mBed = A*L*426.7;   % Mass of adsorbent (in bed) 426.7 = bed density

%%%%%%%%% loop around here to work out p and n and mole fraction for each
%%%%%%%%% loop

% y is an array size n*5 of y1 = 1:n, q1 = n+1:2*n,  
% q2 = 2*n+1:3*n, T = 3*n+1:4*n, P = 4*n+1:5*n
sol = adsorptionSolver(t,z,yF,TPF);

for i = 1:size(sol.y,2)
  [~,velocity(:,i),~,~] = adsorptionModel(sol.x(i),sol.y(:,i),z(2)-z(1),numel(z),yF,TPF); 
end

for i = 1:size(sol.y,2)
  [~,~,vhalf(:,i),~] = adsorptionModel(sol.x(i),sol.y(:,i),z(2)-z(1),numel(z),yF,TPF);  
end

for i = 1:size(sol.y,2)
  [~,~,~,pspan(:,i),] = adsorptionModel(sol.x(i),sol.y(:,i),z(2)-z(1),numel(z),yF,TPF);  
end

% sol.x is time steps
% sol.y is 
Moley2 = 1-sol.y(1:n,1:end);

format long

purityh = 100*Moley2(n,end) / (sol.y(n,end) + Moley2(n,end));
fprintf('Purity of hydrogen is %f%% after %4.2fs\n', purityh, tf)

% purityh = sol.y(3*n,:) ./(sol.y(n,:) + sol.y(3*n,:));
%purityh = (sol.y(n,:) + sol.y(3*n,:)); % mole fraction
% figure
% plot(sol.x, purityh)
% xlabel('time')
% ylabel('purity of hydrogen')

dzhalf = dz/2;
zhalf   = (-dzhalf:dz:L+dzhalf)';
znodal = z';
zspan=cat(1,zhalf,znodal)';
zspan([1:2:end,2:2:end])=zspan; 
 
 
    figure(1)
    
    subplot(1,2,1)
    mesh(sol.x,z,sol.y(1:n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('mole fraction y1')
    title('mole fraction of methane')
    
    subplot(1,2,2)
    mesh(sol.x,z,Moley2)
    xlabel('t')
    ylabel('bed length')
    zlabel('mole fraction y2')
    title('mole fraction of hydrogen')
    
     figure(2)
    
    subplot(1,2,1)
    mesh(sol.x,z,sol.y(n+1:2*n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('loading q1 mol/kg')
    title('loading of methane')
    
    subplot(1,2,2)
    mesh(sol.x,z,sol.y(2*n+1:3*n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('loading q2 mol/kg')
    title('loading of hydrogen')
    
    figure(3)
    
    subplot(1,2,1)
    mesh(sol.x,z,sol.y(3*n+1:4*n,1:end))
    xlabel('t')
    ylabel('bed length')
    zlabel('Temp')
    title('temp of system')
    subplot(1,2,2)
    mesh(sol.x,z,sol.y(4*n+1:5*n,1:end))
    xlabel('time')
    ylabel('bed length')
    zlabel('Pressure')
    title('Pressure of system')
    
    figure(4)
    subplot(1,2,1)
    mesh(sol.x,z,velocity)
    xlabel('t')
    ylabel('bed length')
    zlabel('Velocity')
    title('Intersistial velocity of system')
    subplot(1,2,2)
    mesh(sol.x,zspan,vhalf)
    xlabel('t')
    ylabel('bed length')
    zlabel('Velocity')
    title('Intersistial velocity of system')
    
    
    figure(5)
    mesh(sol.x,z,pspan)
    xlabel('t')
    ylabel('bed length')
    zlabel('dpdz')

    
%% Solver function
function out = adsorptionSolver(t,z,yFeed,TPFeed)
n = numel(z);      % Size of mesh grid
h = diff(z);
h = h(1);
% start up condidtions
Tw = 300;                   % Ambient Tempurature K
Pw = 8e5;                % Ambient Pressure Pa
y1w = 0.7;

% Start up conditions Conditions
y1_0       = zeros(n,1) + y1w;     

q1_0        = zeros(n,1);
q2_0        = zeros(n,1);

T_0         = zeros(n,1) + Tw;

P_0         = zeros(n,1) + Pw;

% y is an array size n*4 of c1 = 1:n, q1 = n+1:2*n, c2 = 2*n+1:3*n, q2 = 3*n+1:4*n
y0      = [y1_0; q1_0; q2_0; T_0; P_0];

% % % Creating mass matrix
% M               = eye(5*n);                         % Mass matrix for temporal domain integrator
% M(1,1)          = 0;                                % Algebraic equation for yi1 at z = 0;
% M(n,n)          = 0;                                % Algebraic equation for yi1 at z = L
% M(2*n+1,2*n+1)  = 0;                                % Algebraic equation for yi2 at z = 0;
% M(3*n,3*n)      = 0;                                % Algebraic equation for yi2 at z = L
% M(4*n+1,4*n+1)  = 0;
% opts = odeset('Mass',M,'MassSingular','yes');       % creates mass matrix for the function ode15s

% % Solving
% te=t(2)

out = ode15s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0);

%[Velocity, ~] = adsorptionModel(t,y,h,n,yFeed,TPFeed)
% Solving
% out = ode23s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0);
end
%% Adsorption model
function [dydt, Velocity, Vhalf, Pspan] = adsorptionModel(t,y,h,n,yiFeed,TPFeed)

% Variables allocation
% y is an array size n*5 of y1 = 1:n, q1 = n+1:2*n,  
% q2 = 2*n+1:3*n, T = 3*n+1:4*n, P = 4*n+1:5*n
y1 = y(1:n);
y2 = 1 - y1;

q1 = y(n+1:2*n);
q2 = y(2*n+1:3*n);

T  = y(3*n+1:4*n);
P  = y(4*n+1:5*n);

% Z half points
yhalf = zeros(n+1,1);
Thalf = zeros(n+1,1);
Phalf = zeros(n+1,1);
uhalf = zeros(n+1,1);
u = zeros(n,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Constant physical properties and basic simulation parameters
R = 8.3145;                 % ideal gas constant J/mol/K 
D   = 1.3e-5;               % Axial Dispersion coefficient m2/s
epl = 0.4043;               % Void fraction of bed
eplp = 0.546;               % Void fraction of particle
eplt = epl + (1-epl)*eplp;  % Total void fraction 
k   = [0.136;0.259];        % (lumped) Mass Transfer Coefficient 1/s
rhop = 716.3;               % Particle density kg/m3
% rhob = 426.7;               % Bed density kg/m3 = (1-epl)*rhop
Cps = 1046.7;               % heat capacity of solid J/kg/K 
Kl = 1.2e-6;                % Thermal diffusity J/m/s/K
Kz = 0.09;                   % Effective gas thermal conductivity 
deltaH = [24124; 8420];     % heat of adsorption J/mol
Vfeed = TPFeed(1);
Tfeed = TPFeed(2);
PH = TPFeed(3); 
PI = TPFeed(4);
PL = TPFeed(5);
lambda = 0.5;                % rate of pressure change (1/s)

% Ergun equation parameters
visc = 3.73e-5;     % gas viscosity kg/m/s
Rp = 5.41e-3;        % particle radius m

% Constants for CH4
c1a = -0.703029; c1b = 108.4773; c1c = -42.52157; c1d = 5.862788; c1e = 0.678565;
% Constants for H2
c2a = 33.066178; c2b = -11.363417; c2c = 11.432816; c2d = -2.772874; c2e = -0.158558;


% Cpg = heat capacity (J/mol*K)
Cpg1 = c1a + c1b*(T/1e3) + c1c*(T/1e3).^2 + c1d*(T/1e3).^3 + c1e./((T/1e3).^2);
Cpg2 = c2a + c2b*(T/1e3) + c2c*(T/1e3).^2 + c2d*(T/1e3).^3 + c2e./((T/1e3).^2);
Cpg = (y1.*Cpg1 + y2.*Cpg2);
% cool prop for viscosity eqs

%%%%%%%% Boundary conditions: j = 0.5 and n + 0.5 %%%%%%%%%%%%%%%%%%%%%%%%    
    %% Pressurization
    
% Phalf(1) = PH - (PH - PL)*exp(-lambda*t);
% Phalf(n+1) = P(n);
% 
% uhalf(1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((P(1)-Phalf(1)));
% %vhalf(n+1) = -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*(Phalf(n+1)-P(n));
% uhalf(n+1) = 0;
% 
% %yhalf(1) = 0.7;
% yhalf(1) = (y1(1)+yiFeed*uhalf(1)*h/D/2) / (1+uhalf(1)*h/D/2);
% %yhalf(1) = uhalf(1)*h/(2*D)*(yiFeed-y1(1)) + y1(1)
% 
% yhalf(n+1) = y1(n);
% 
%     % Inlet gas condtions
%     rhogInlet = ((Phalf(1)) .* (yhalf(1)*16.04e-3 + (1-yhalf(1))*2.02e-3)) ./ (R*Tfeed);
%     Cpg1in = c1a + c1b*(Tfeed/1e3) + c1c*(Tfeed/1e3).^2 + c1d*(Tfeed/1e3).^3 + c1e./((Tfeed/1e3).^2);
%     Cpg2in = c2a + c2b*(Tfeed/1e3) + c2c*(Tfeed/1e3).^2 + c2d*(Tfeed/1e3).^3 + c2e./((Tfeed/1e3).^2);
%     CpgInlet = (yhalf(1).*Cpg1in + (1-yhalf(1)).*Cpg2in);
%     %
% 
% Thalf(1) = (T(1) + Tfeed*uhalf(1)*rhogInlet*CpgInlet*h/Kl/2)/(1 + uhalf(1)*rhogInlet*CpgInlet*h/Kl/2);
% Thalf(n+1) = T(n);
    
    %% Adsorption

%     uhalf(1) = TPFeed(1);
%     %Phalf(1) = PH;
%     Phalf(1)  =  P(1) + (uhalf(1)*h*visc/2) / (4/150*(epl/(1-epl))^2*Rp^2);
%     Phalf(n+1) = PH;
%     uhalf(n+1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((Phalf(n+1)-P(n)));
%     yhalf(1) = (y1(1)+yiFeed*uhalf(1)*h/D/2) / (1+uhalf(1)*h/D/2);
%     yhalf(n+1) = y1(n);
%     
%      % Inlet gas condtions
%         rhogInlet = ((Phalf(1)) .* (yhalf(1)*16.04e-3 + (1-yhalf(1))*2.02e-3)) ./ (R*Tfeed);
%         Cpg1in = c1a + c1b*(Tfeed/1e3) + c1c*(Tfeed/1e3).^2 + c1d*(Tfeed/1e3).^3 + c1e./((Tfeed/1e3).^2);
%         Cpg2in = c2a + c2b*(Tfeed/1e3) + c2c*(Tfeed/1e3).^2 + c2d*(Tfeed/1e3).^3 + c2e./((Tfeed/1e3).^2);
%         CpgInlet = (yhalf(1).*Cpg1in + (1-yhalf(1)).*Cpg2in);
%         %
% 
%     Thalf(1) = (T(1) + Tfeed*uhalf(1)*rhogInlet*CpgInlet*h/Kl/2)/(1 + uhalf(1)*rhogInlet*CpgInlet*h/Kl/2);
%     Thalf(n+1) = T(n);

    %% Depres

    Phalf(1) = P(1);
    Phalf(n+1) = PI + (PH - PI)*exp(-lambda*t);
    
    uhalf(1) = 0;
    %uhalf(1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((P(1)-Phalf(1)));
    uhalf(n+1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((Phalf(n+1)-P(n)));
    
    yhalf(1) = y1(1);
    yhalf(n+1) = y1(n);
    
    Thalf(1) = T(1);
    Thalf(n+1) = T(n);

    %% Purge
%     Phalf(1) = PL + (PI - PL)*exp(-lambda*t);
%     Phalf(n+1) = P(n);
%     
%     uhalf(1) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((P(1)-Phalf(1)));
%     uhalf(n+1) = 0; 
%     %uhalf(n) =   -(2/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((Phalf(n+1)-P(n)));
%     
%     yhalf(1) = y1(1);
%     yhalf(n+1) = y1(n);
%     
%     Thalf(1) = T(1);
%     Thalf(n+1) = T(n);
  
%%%%%%%% Velocity wall values: j 1.5 to n - 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
uhalf(2:n) =   -(1/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((P(2:n) - P(1:n-1)));

%%%%%%%% Wall values: j = 1.5 %%%%%%%%%%%%%%%%%%%%%%%%
alpha0P = (2/3) / (P(2)-P(1) + 1e-10)^4;
alpha1P = (1/3) / (2*(P(1)-Phalf(1)) + 1e-10)^4;
Phalf(2) = alpha0P/(alpha0P+alpha1P)*(0.5*(P(1)+P(2))) + alpha1P/(alpha0P+alpha1P)*(2*P(1)-Phalf(1));

alpha0y = (2/3) / (y1(2)-y1(1) + 1e-10)^4;
alpha1y = (1/3) / (2*(y1(1)-yhalf(1)) + 1e-10)^4;
yhalf(2) = alpha0y/(alpha0y+alpha1y)*(0.5*(y1(1)+y1(2))) + alpha1y/(alpha0y+alpha1y)*(2*y1(1)-yhalf(1));

alpha0T = (2/3) / (T(2)-T(1) + 1e-10)^4;
alpha1T = (1/3) / (2*(T(1)-Thalf(1)) + 1e-10)^4;
Thalf(2) = alpha0T/(alpha0T+alpha1T)*(0.5*(T(1)+T(2))) + alpha1T/(alpha0T+alpha1T)*(2*T(1)-Thalf(1));

%%%%%%%% Wall values: j = 2.5 to n - 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
alpha0P = (2/3) ./ (P(3:n)-P(2:n-1) + 1e-10).^4;
alpha1P = (1/3) ./ (P(2:n-1)-P(1:n-2) + 1e-10).^4;
Phalf(3:n) = alpha0P./(alpha0P+alpha1P).*(0.5.*(P(2:n-1)+P(3:n))) + alpha1P./(alpha0P+alpha1P).*(1.5.*P(2:n-1)- 0.5* P(1:n-2));

alpha0y = (2/3) ./ (y1(3:n)-y1(2:n-1) + 1e-10).^4;
alpha1y = (1/3) ./ (y1(2:n-1)-y1(1:n-2) + 1e-10).^4;
yhalf(3:n) = alpha0y./(alpha0y+alpha1y).*(0.5*(y1(2:n-1)+y1(3:n))) + alpha1y./(alpha0y+alpha1y).*(1.5.*y1(2:n-1)- 0.5* y1(1:n-2));

alpha0T = (2/3) ./ (T(3:n)-T(2:n-1) + 1e-10).^4;
alpha1T = (1/3) ./ (T(2:n-1)-T(1:n-2)+ 1e-10).^4;
Thalf(3:n) = alpha0T./(alpha0T+alpha1T).*(0.5*(T(2:n-1)+T(3:n))) + alpha1T./(alpha0T+alpha1T).*(1.5.*T(2:n-1)- 0.5* T(1:n-2));
   
% rhog = (P .* (y1*16.04e-3 + y2*2.02e-3)) ./ (R*T);  % ρ = PM/RT, M: H2 = 2.02 g/mol CH4 = 16.04 g/mol Air= 28.97 g/mol 50/50 H2 and CH4 = 9.03
% % Central velocity eqautions
u(1:n) =   -(1/visc/h)*(4/150*(epl/(1-epl))^2)*(Rp^2)*((Phalf(2:n+1) - Phalf(1:n)));
rhog = (P .* (y1*16.04e-3 + y2*2.02e-3)) ./ (R*T);  % ρ = PM/RT, M: H2 = 2.02 g/mol CH4 = 16.04 g/mol Air= 28.97 g/mol 50/50 H2 and CH4 = 9.03
%%%%%%%%%%%%%%% Langmuir Equation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % CH4
a11 = -9.5592;	a21 = 4638.5;	b11 = 3.725e-4/1e5;	b21 = 1.443e4;
% % H2
a12 = -23.131;	a22 = 8069.1;	b12 = 2.248/1e5;	b22 = -1.435e4;
qs1 = a11 + (a21./T);
qs2 = a12 + (a22./T);
B1 = b11.*exp(b21./(8.3145*T));
B2 = b12.*exp(b22./(8.3145*T));
qstar1 = (qs1.*B1.*P.*y1) ./ (1+((P.*y1.*B1)+(P.*(1-y1).*B2))); %mols/kg
qstar2 = (qs2.*B2.*P.*(1-y1)) ./ (1+((P.*y1.*B1)+(P.*(1-y1).*B2)));

LT1 = qstar1-q1;
LT2 = qstar2-q2;

%%%%%%%%%%%%%%% Eqautions for loading of component i %%%%%%%%%%%%%%
dq1dt = k(1)*LT1;
dq2dt = k(2)*LT2;

%%%%%%%%%%%%%% Balance equations

dPdt = zeros(n,1);
dy1dt = zeros(n,1);
dTdt = zeros(n,1);

% component mass balance
dy1dt(1) = D/h.*(((y1(2)-y1(1))/h) - 2/h*(y1(1)-yhalf(1))) - u(1)/h.*(yhalf(2)-yhalf(1)) - R*T(1)./P(1).*((1-epl)/epl).*rhop.*(dq1dt(1) - y1(1).*(dq1dt(1)+dq2dt(1)));
dy1dt(2:n-1) = D/h.*((y1(3:n)-y1(2:n-1))/h - (y1(2:n-1)-y1(1:n-2))/h) - u(2:n-1)/h.*(yhalf(3:n)-yhalf(2:n-1)) - R*T(2:n-1)./P(2:n-1).*((1-epl)/epl).*rhop.*(dq1dt(2:n-1) - y1(2:n-1).*(dq1dt(2:n-1)+dq2dt(2:n-1)));
dy1dt(n) = D/h.*(2/h*(yhalf(n+1)-y1(n))-(y1(n)-y1(n-1))/h) - u(n)/h.*(yhalf(n+1)-yhalf(n)) - R*T(n)./P(n).*((1-epl)/epl).*rhop.*(dq1dt(n) - y1(n).*(dq1dt(n)+dq2dt(n)));

% Energy balance
dTdt(1) = 1./(eplt.*rhog(1).*Cpg(1)+rhop*Cps) .* (epl.*Kl/h.*(((T(2)-T(1))/h)-2/h*(T(1)-Thalf(1))) - rhog(1).*Cpg(1).*epl.*u(1)/h.*(Thalf(2)-Thalf(1)) + rhop.*(deltaH(1)*dq1dt(1)+deltaH(2)*dq2dt(1)));
dTdt(2:n-1) = 1./(eplt.*rhog(2:n-1).*Cpg(2:n-1)+rhop*Cps) .* (epl.*Kl/h.*(((T(3:n)-T(2:n-1))/h)-(T(2:n-1)-T(1:n-2))/h) - rhog(2:n-1).*Cpg(2:n-1).*epl.*u(2:n-1)/h.*(Thalf(3:n)-Thalf(2:n-1)) + rhop.*(deltaH(1)*dq1dt(2:n-1)+deltaH(2)*dq2dt(2:n-1)));
dTdt(n) = 1./(eplt.*rhog(n).*Cpg(n)+rhop*Cps) .* (epl.*Kl/h.*(2/h*(Thalf(n+1)-T(n))-(T(n)-T(n-1))/h) - rhog(n).*Cpg(n).*epl.*u(n)/h.*(Thalf(n+1)-Thalf(n)) + rhop.*(deltaH(1)*dq1dt(n)+deltaH(2)*dq2dt(n)));

% total mass balance
dPdt(1) = D/h.*( ((P(2)-P(1))/h)-2/h*(P(1)-Phalf(1))) - P(1)./h.*(uhalf(2)-uhalf(1)) - u(1)/h.*(Phalf(2)-Phalf(1)) + P(1)./T(1).*( -D/h.*(((T(2)-T(1))/h)-2/h*(T(1)-Thalf(1))) + dTdt(1) + u(1)/h.*(Thalf(2)-Thalf(1)) ) - ((1-epl)/epl)*rhop*R.*T(1).*(dq1dt(1)+dq2dt(1));
dPdt(2:n-1) = D/h.*( (P(3:n)-P(2:n-1))/h - (P(2:n-1)-P(1:n-2))/h) - P(2:n-1)./h.*(uhalf(3:n)-uhalf(2:n-1)) - u(2:n-1)/h.*(Phalf(3:n)-Phalf(2:n-1)) + P(2:n-1)./T(2:n-1).*( -D/h.*(((T(3:n)-T(2:n-1))/h)-(T(2:n-1)-T(1:n-2))/h) + dTdt(2:n-1) + u(2:n-1)/h.*(Thalf(3:n)-Thalf(2:n-1)) ) - ((1-epl)/epl)*rhop*R.*T(2:n-1).*(dq1dt(2:n-1)+dq2dt(2:n-1));
dPdt(n) = D/h.*(2/h*(Phalf(n+1)-P(n))-(P(n)-P(n-1))/h) - P(n)./h.*(uhalf(n+1)-uhalf(n)) - u(n)/h.*(Phalf(n+1)-Phalf(n)) + P(n)./T(n).*( -D/h.*(2/h*(Thalf(n+1)-T(n))-(T(n)-T(n-1))/h) + dTdt(n) + u(n)/h.*(Thalf(n+1)-Thalf(n)) ) - ((1-epl)/epl)*rhop*R.*T(n).*(dq1dt(n)+dq2dt(n));

Velocity = u;

Vhalf=cat(1,uhalf,u);
Vhalf([1:2:end,2:2:end])=Vhalf;

% Pspan=cat(1,Phalf,P);
% Pspan([1:2:end,2:2:end])=Pspan;
Pspan = (Phalf(2:n+1) - Phalf(1:n));

%%%%%%%%%%%%%%%% Concatenate vectors of time derivatives
dydt = [dy1dt;dq1dt;dq2dt;dTdt;dPdt];
end

[Torsten]

Seems it's somehow caused by the implementation of the boundary conditions in your discretization scheme.

Does the same happen if you use a simple first-order upwind scheme instead of WENO ?

[Thomas Stone-Wigg]

I’ll give it a go tomorrow, thanks for the suggestion

[Thomas Stone-Wigg]

No luck with first-oder upwind scheme. The spike is still there after trying it with just P, and then P, T and y.

This is the method i used to find the wall values:

Phalf(2:n) = P(1:n-1);

[Torsten]

What are the equations you tried to implement with initial and boundary conditions ?

[Thomas Stone-Wigg]

I'm trying to implement the balance equations from the attached document with a few changes:

I use Darcy's equation for velocity not the Ergun equation (questionably assuming laminar flow but im trying to change to the Ergun)

The cycle steps are sligthly different, at the moment i have feed pressurization, adsorption, forward depressurization and reverse evacuation (also about to change an open-open purge step).

They use Danckwert’s boundary conditions for a dispersed plug flow system. I also use half cell approximations to find some edge gradients so i dont have to use ghost cells, eg:

P(1)-P(0) = 2*(P(1)-Phalf(1))

The pressure changing steps has pressure changing as a fucntion of time with lambda controlling the rate:

Phalf(1) = PH - (PH - PL)*exp(-lambda*t);

As i test each step indiviually i set the pressure to what it would be at after the previous step and the mole fraction is always 0.7 initially.

[Torsten]

Sorry, but it means too much effort for me to check.

But it cannot be all you have to do to switch to first-order upwind by defining

Phalf(2:n) = P(1:n-1);

I always build up a code in steps and check whether it produces reasonable results from step to step. So I think it cannot be true that you started with the difficult WENO scheme, can it ?

[Thomas Stone-Wigg]

No your right i started with a simpler FDM and only two of the balance equations which is working well. Thanks for your help, ill try break the code down again.