You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Passing out additions parameters after ODE solver.
4 views (last 30 days)
Show older comments
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
4 Comments
Dyuman Joshi
on 23 Jan 2024
Edited: Dyuman Joshi
on 23 Jan 2024
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
on 23 Jan 2024
Edited: Thomas Stone-Wigg
on 23 Jan 2024
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
on 23 Jan 2024
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
on 23 Jan 2024
Sorry, i want to make a mesh plot of velocity, time and position in the bed.
Accepted Answer
Torsten
on 23 Jan 2024
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.
9 Comments
Thomas Stone-Wigg
on 23 Jan 2024
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
on 23 Jan 2024
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
on 23 Jan 2024
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
on 24 Jan 2024
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
on 24 Jan 2024
What are the equations you tried to implement with initial and boundary conditions ?
Thomas Stone-Wigg
on 24 Jan 2024
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
on 24 Jan 2024
Edited: Torsten
on 24 Jan 2024
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
on 25 Jan 2024
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.
More Answers (0)
See Also
Categories
Find more on Gas Dynamics in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)