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 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
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 件のコメント
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Particle & Nuclear Physics についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!