清空筛选条件
清空筛选条件

Im trying you solve the component mass balance and energy balance equations for an adsorption system. Im getting some exponential growth in the plots.

16 次查看（过去 30 天）
Thomas Stone-Wigg
Thomas Stone-Wigg 2023-11-1
评论： Thomas Stone-Wigg 2023-11-3
Ran in:
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);
Warning: Failure at t=1.893223e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.684342e-14) at time t.
% 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

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2023-11-1
编辑：Torsten 2023-11-1
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 个评论
显示 1更早的评论
Torsten
Torsten 2023-11-1
编辑：Torsten 2023-11-2
If there are only two species present in your mixture (methane and hydrogen), you don't need two molar balances since y1 + y2 = 1.
And I forgot to mention that you should use upwinding also for dy/dz as I suggested for dTdz above.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by