Im trying you solve the component mass balance and energy balance equations for an adsorption system. Im getting some exponential growth in the plots.
15 views (last 30 days)
Show older comments
Thomas Stone-Wigg
on 1 Nov 2023
Commented: Thomas Stone-Wigg
on 3 Nov 2023
clear,clc
%% MAIN
% Parameters
L = 1; % Column length
t0 = 0; % Initial Time
tf = 20; % Final time
dt = 20; % Time step
t = t0:dt:tf; % Time vector
dz = 0.25; % Mesh size
z = 0:dz:L; % Mesh generation
n = numel(z);
% % Parameters
% L = 1; % Column length
% t0 = 0; % Initial Time
% tf = 20; % 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 %%%%%%%%%%%
P = 6e5; % feed pressure in pascals
TFeed = 310; %tempurature of feed (degrees C)
yiF_1 = 0.95; % Mole fraction for methane
yiF_2 = 1-yiF_1; % Mole fraction for hydrogen
yiF = [yiF_1;yiF_2];
% 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 T = 4*n+1:5*n
sol = adsorptionSolver(t,z,yiF,TFeed);
% sol.x is time steps
% sol.y is
return
figure
subplot(4,1,1)
plot(sol.x,sol.y(n,:))
xlabel('t')
ylabel('c1')
title('yi1(t) at L')
subplot(4,1,2)
plot(sol.x,sol.y(3*n,:))
xlabel('t')
ylabel('c2')
title('yi2(t) at L')
subplot(4,1,3)
plot(sol.x,sol.y(4*n,:))
xlabel('t')
ylabel('q2')
title('q2(t) at L')
subplot(4,1,4)
plot(sol.x,sol.y(5*n,:))
xlabel('t')
ylabel('T')
title('T(t) at L')
%% Solver function
function out = adsorptionSolver(t,z,yiFeed,TFeed)
n = numel(z); % Size of mesh grid
h = diff(z);
h = h(1);
Tw = 300; % Ambient Tempurature K
% Initial Conditions
yi1_0 = zeros(n,1);
yi1_0(1) = yiFeed(1);
q1_0 = zeros(n,1);
yi2_0 = zeros(n,1);
yi2_0(1) = yiFeed(2);
q2_0 = zeros(n,1);
T_0 = zeros(n,1) + Tw;
T_0(1) = TFeed;
% 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 = [yi1_0; q1_0; yi2_0; q2_0; T_0];
% Creating mass matrix
M = eye(5*n); % Mass matrix
M(1,1) = 0; % Algebraic equation for c1 at z = 0;
M(n,n) = 0; % Algebraic equation for c1 at z = L
M(2*n+1,2*n+1) = 0; % Algebraic equation for c2 at z = 0;
M(3*n,3*n) = 0; % Algebraic equation for c2 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,yiFeed,TFeed),t,y0,opts);
end
%% Adsorption model
function dydt = adsorptionModel(~,y,h,n,yiFeed,TFeed)
% Variables allocation
% 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, T = 4*n+1:5*n
yi1 = y(1:n);
qi1 = y(n+1:2*n);
yi2 = y(2*n+1:3*n);
qi2 = y(3*n+1:4*n);
T = y(4*n+1:5*n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Physical properties and basic simulation parameters
P = 6e5; % feed pressure in pascals
R = 8.3145; % ideal gas constant J/mol/K
D = 1.3e-5; % Axial Dispersion coefficient m2/s
u = 3e-1; % Superficial velocity m/s
epl = 0.4043; % Void fraction of bed
% eplp = 0.546; % Void fraction of particle
% eplt = epl + (1-epl)*eplp; % Toatl void fraction
k = [0.136;0.259]; % (lumped) Mass Transfer Coefficient 1/s
rhop = 716.3; % Particle density kg/m3
rhog = (P * (yi1*16.04e-3 + yi2*2.02e-3)) ./ (R*T);
rhob = 426.7; % Bed density kg/m3
Cps = 1046.7; % heat capacity of solid J/kg/K
Cpg = yi1*2226 + yi2*14310; % specific heat capacity of gas % https://www.engineeringtoolbox.com/hydrogen-d_976.html
Kl = 1.2e-6; % Thermal diffusity J/m/s/K
deltaH = [24124; 8420]; % heat of adsorption J/mol
%%%%%%%%%%%%%%% Langmuir Equation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CH4
alp11 = 0.0086; alp12 = -0.2155; bet11 = 0.0004066; bet12 = -0.010604;
% H2
alp21 = -0.0000379; alp22 = -0.01815; bet21 = 2.2998; bet22 = -0.05993;
a1 = alp11.*exp(alp12*T);
a2 = alp21*exp(alp22*T);
b1 = bet11*exp(bet12*T);
b2 = bet21*exp(bet22*T);
qstar1 = (P*a1.*yi1) ./ (1+((P*b1.*yi1)+(P*b2.*yi2))); %mols/kg
qstar2 = (P*a2.*yi2) ./ (1+((P*b1.*yi1)+(P*b2.*yi2))); %mols/kg
LT1 = qstar1-qi1;
LT2 = qstar2-qi2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Eqautions for loading of component i %%%%%%%%%%%%%%
dqi1dt = k(1)*LT1;
dqi2dt = k(2)*LT2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Temperature Equations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dTdz = zeros(n,1);
dTdz(2:n-1) = (T(3:n)-T(2:n-1))/h;
d2Tdz2 = zeros(n,1);
d2Tdz2(2:n-1) = (T(3:n)-2*T(2:n-1)+T(1:n-2))/h^2;
dTdt = (ones(n,1)./ (epl*rhog.*Cpg+rhob*Cps)) .* ( (-rhog.*Cpg*epl*u.*dTdz) + (Kl*d2Tdz2) + (rhob*(dqi1dt*deltaH(1)+dqi2dt*deltaH(2))) );
dTdt(1) = T(1)-TFeed(1);
dTdt(n) = T(n)-T(n-1);
% change in 1/t with z
J = 1./T;
dJdz = zeros(n,1);
dJdz(2:n-1) = (J(3:n)-J(2:n-1))/h;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Concentration equations %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change concentration of 1 (methane) with change of z
dyi1dz = zeros(n,1);
dyi1dz(2:n-1) = (yi1(3:n)-yi1(2:n-1))/h;
% Rate of change concentration of 1 (methane) with change of z
d2yi1dz2 = zeros(n,1);
d2yi1dz2(2:n-1) = (yi1(3:n)-2*yi1(2:n-1)+yi1(1:n-2))/h^2;
% Change concentration of 1 (methane) with change of time
dyi1dt = D * (d2yi1dz2 + 2*T.*dyi1dz.*dJdz) - u*dyi1dz - (R*T/P)*((1-epl)/epl)*rhop.* (dqi1dt - yi1.*(dqi1dt+dqi2dt));
dyi1dt(1) = yi1(1)-yiFeed(1);
dyi1dt(n) = yi1(n)-yi1(n-1);
% Change concentration of 2 (hydrogen) with change of z
dyi2dz = zeros(n,1);
dyi2dz(2:n-1) = (yi2(3:n)-yi2(2:n-1))/h;
% Rate of change concentration of 2 (hydrogen) with change of z
d2yi2dz2 = zeros(n,1);
d2yi2dz2(2:n-1) = (yi2(3:n)-2*yi2(2:n-1)+yi2(1:n-2))/h^2;
% Change concentration of 2 (hydrogen) with change of time
dyi2dt = D * (d2yi2dz2 + 2*T.*dyi2dz.*dJdz) - u*dyi2dz - (R*T/P)*((1-epl)/epl)*rhop.* (dqi2dt - yi2.*(dqi1dt+dqi2dt));
dyi2dt(1) = yi2(1)-yiFeed(2);
dyi2dt(n) = yi2(n)-yi2(n-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Concatenate vectors of time derivatives
dydt = [dyi1dt;dqi1dt;dyi2dt;dqi2dt;dTdt];
end
0 Comments
Accepted Answer
Torsten
on 1 Nov 2023
Edited: Torsten
on 1 Nov 2023
Use Upwind Differencing for the convection terms:
dTdz(2:n) = (T(2:n)-T(1:n-1))/h;
dJdz(2:n) = (J(2:n)-J(1:n-1))/h;
The term
-D_x*2*T*dy_i/dz*d(1/T)/dz
in your component equation looks rather strange. Where does it stem from and what does it model ?
3 Comments
More Answers (0)
See Also
Categories
Find more on Particle & Nuclear Physics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!