Passing out additions parameters after ODE solver.
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
Dyuman Joshi
2024-1-23
编辑:Dyuman Joshi
2024-1-23
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
2024-1-23
编辑:Thomas Stone-Wigg
2024-1-23
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
2024-1-23
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
2024-1-23
Sorry, i want to make a mesh plot of velocity, time and position in the bed.
采纳的回答
Torsten
2024-1-23
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 个评论
Thomas Stone-Wigg
2024-1-23
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
2024-1-23
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
2024-1-23
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
2024-1-24
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
2024-1-24
What are the equations you tried to implement with initial and boundary conditions ?
Thomas Stone-Wigg
2024-1-24
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
2024-1-24
编辑:Torsten
2024-1-24
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
2024-1-25
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Applications 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)