Solve a system of coupled differential equations with ODE23s
16 views (last 30 days)
Show older comments
Thomas Stone-Wigg
on 11 Jan 2024
Commented: William Rose
on 12 Jan 2024
Hi,
Im trying to solve the system of coupled differential equations (TMB, CMB and EB at the bottom of the code) with ODE23s but I'm unsure of the correct method to set up the equations to pass them to the solver.
This snippet is just to show how they are coupled:
dPdt(1:n) = ...* dTdt(1:n);
dy1dt(1:n) = ...*dPdt(1:n) ... *dTdt(1:n);
dTdt(1:n) = ...*dPdt(1:n);
I've attached the three equations and their non-dimensional and discretized forms for reference
At the moment im not too worried about the efficiency of the code but any extra advice would be much appreciated
Many thanks
Tom
clear,clc
%% MAIN
% Physical Parameters
L = 1; % Column length m
r = 0.25; % Bed radius m
t0 = 0; % Initial Time
tf = 500; % Final time
dt = 0.1; % Time step
t = t0:dt:tf; % Time vector
dz = 0.1; % Mesh size
z = 0:dz:L; % Mesh generation
n = numel(z);
%%%%%%% blended gas parameters %%%%%%%%%%%
%%% calc u input
TPF = [1 300 6e5 1.75e5 1e5]; % Feed: Velocity (m/s), tempurature (K) and pressure (Pa) [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*4 of y1 = 1:n, q1 = n+1:2*n, y2 = 2*n+1:3*n,
% q2 =3*n+1:4*n T = 4*n+1:5*n
sol = adsorptionSolver(t,z,yF,TPF);
% sol.x is time steps
% sol.y is
format long
purityh = 100*sol.y(3*n,end) /(sol.y(n,end) + sol.y(3*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,sol.y(2*n+1:3*n,1:end))
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(3*n+1:4*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(4*n+1:5*n,1:end))
xlabel('t')
ylabel('bed length')
zlabel('Temp')
title('temp of system')
subplot(1,2,2)
mesh(sol.x,z,sol.y(5*n+1:6*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 = 5.95e5; % Ambient Pressure Pa
y1w = 0.7;
% Initial 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 = ode23s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0);
end
%% Adsorption model
function dydt = 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);
vhalf = zeros(n+1,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);
%%%%%%%% Boundary conditions: j = 0.5 and n + 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
%%%% Pressurization
Phalf(1) = PH - (PH - PL)*exp(-lambda*t);
Phalf(n+1) = P(n);
vhalf(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));
yhalf(1) = (y1(1)+yiFeed*vhalf(1)*h) / (2*D+vhalf(1)*h);
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*vhalf(1)*epl*rhogInlet*CpgInlet*h)/(2*Kz + vhalf(1)*epl*rhogInlet*CpgInlet*h);
Thalf(n+1) = T(n);
%%%%%%%% Velocity wall values: j 1.5 to n - 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
vhalf(2:n) = -(2/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)-y(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
%%%%%%%%%%%%%%% Langmuir Equation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % CH4
a11 = -9.5592; a21 = 4638.5; b11 = 3.725e-4/100000; b21 = 1.443e4;
% % H2
a12 = -23.131; a22 = 8069.1; b12 = 2.248/100000; 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);
% TMB
dPdt(1:n) = (T(1:n)./h).*(vhalf(2:n+1).*Phalf(2:n+1)./Thalf(2:n+1)-vhalf(1:n).*Phalf(1:n)./Thalf(1:n)) - ((R*T(1:n)*(1-epl))/epl).*(dq1dt(1:n)+dq2dt(1:n)) + P(1:n)/T(1:n) * dTdt(1:n);
% CMB
dy1dt(1) = (D.*T(1)./(h.*P(1))) * ( (Phalf(2)./Thalf(2)).*((y(2)-y(1))/h)-(Phalf(1)./Thalf(1)).*(2/h*(y(1)-yhalf(1)))) - (T(1)./P(1)/h).*(vhalf(2).*Phalf(2).*yhalf(2)./Thalf(2) - vhalf(1).*Phalf(1).*yhalf(1)./Thalf(1)) - ((R*T(1)*(1-epl))/epl./P(1)).*dq1dt(1) - (y1(1)/P(1))*dPdt(1) + (y1(1)/T(1))*dTdt(1);
dy1dt(2:n-1) = (D.*T(2:n-1)./(h.*P(2:n-1))) .* ( (Phalf(3:n)./Thalf(3:n)).*((y(3:n)-y(2:n-1))/h)-(Phalf(2:n-1)./Thalf(2:n-1)).*((y(2:n-1)-y(1:n-2))/h)) - (T(2:n-1)./P(2:n-1)/h).*(vhalf(3:n).*Phalf(3:n).*yhalf(3:n)./Thalf(3:n) - vhalf(2:n-1).*Phalf(2:n-1).*yhalf(2:n-1)./Thalf(2:n-1)) - ((R.*T(2:n-1).*(1-epl))./epl./P(2:n-1)).*dq1dt(2:n-1) - (y1(2:n-1)./P(2:n-1)).*dPdt(2:n-1) + (y1(2:n-1)./T(2:n-1)).*dTdt(2:n-1);
dy1dt(n) = (D.*T(n)./(h.*P(n))) * ( (Phalf(n+1)./Thalf(n+1)).*(2/h*(yhalf(n+1)-y(n)))-(Phalf(n)./Thalf(n)).*((y(n)-y(n-1))/h)) - (T(n)./P(n)/h).*(vhalf(n+1).*Phalf(n+1).*yhalf(n+1)./Thalf(n+1) - vhalf(n).*Phalf(n).*yhalf(n)./Thalf(n)) - ((R*T(n)*(1-epl))/epl./P(n)).*dq1dt(n) - (y1(n)/P(n))*dPdt(n) + (y1(n)/T(n))*dTdt(n);
% EB
Cpa = (q1.*Cpg1)./(q1+q2+1e-10) + (q2.*Cpg2)./(q1+q2+1e-10); % 1e-10 is there to stop nan
ebcoeff = 1./( ((1-epl)/epl).*(rhop.*Cps + Cpa.*(dq1dt+dq2dt)) );
dTdt(1) = ebcoeff(1)*(Kl/epl/h*((T(2)-T(1))/h-2*(T(1)-Thalf(1))/h)-Cpg(1)./(R*h).*(Phalf(2).*vhalf(2)-Phalf(1).*vhalf(1))-((1-epl)/epl).*Cpa(1).*T(1).*(dq1dt(1)+dq2dt(1)) + ((1-epl)/epl)*(-deltaH(1)*dq1dt(1)-deltaH(2)*dq2dt(1)) - Cpg(1)/R.*dPdt(1));
dTdt(2:n-1) = ebcoeff(2:n-1).*(Kl/epl/h*((T(3:n)-2.*T(2:n-1)+T(1:n-2))./h)-Cpg(2:n-1)./(R*h).*(Phalf(3:n).*vhalf(3:n)-Phalf(2:n-1).*vhalf(2:n-1))-((1-epl)/epl).*Cpa(2:n-1).*T(2:n-1).*(dq1dt(2:n-1)+dq2dt(2:n-1)) + ((1-epl)/epl).*(-deltaH(1).*dq1dt(2:n-1)-deltaH(2).*dq2dt(2:n-1)) - Cpg(2:n-1)./R.*dPdt(2:n-1));
dTdt(n) = ebcoeff(n)*(Kl/epl/h*(2*(Thalf(n+1)-T(n))/h-(T(n)+T(n-1))/h)-Cpg(n)./(R*h).*(Phalf(n+1).*vhalf(n+1)-Phalf(n).*vhalf(n))-((1-epl)/epl).*Cpa(n).*T(n).*(dq1dt(n)+dq2dt(n)) + ((1-epl)/epl)*(-deltaH(1)*dq1dt(n)-deltaH(2)*dq2dt(n)) - Cpg(n)/R.*dPdt(n));
%%%%%%%%%%%%%%%% Concatenate vectors of time derivatives
dydt = [dy1dt;dq1dt;dq2dt;dTdt;dPdt];
end
0 Comments
Accepted Answer
William Rose
on 12 Jan 2024
Edited: William Rose
on 12 Jan 2024
Does your code run? Do you want help making it run at all, or help making it run better? If better, better how? You said speed is not an issue. If it doesn't run, what is the error message?
I notice that you use t in multiple roles on the right hand side below. That may not turn out well.
t = t0:dt:tf; % Time vector
% Solving
out = ode23s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),t,y0);
I would do the following instead
tspan = t0:dt:tf; % Time vector
% Solving
out = ode23s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),tspan,y0);
Good luck.
2 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!