How can i run this script to make a plot for a lot of T
조회 수: 1 (최근 30일)
이전 댓글 표시
function [f] = acidpretreatment(t,x,T,xyl_zero)
%Components
xyl= x(1);
arab= x(2);
acet= x(3);
acet_dot = x(4);
cel= x(5);
Xylo= x(6);
Acetic= x(7);
Acetic_dot = x(8);
Arabi = x(9);
glu = x(10);
furf= x(11);
%---------------------------------------------------------------------
%Data
R=8.314;
HA=0.07; % m=948.320,38 kg
i=100:180
T=375+1i; % kelvin "Pressure 10 bar"
k10=2.37; k20=2.17; k1ac=2.37; k30=2.37; %min^-1%
n1=1.51; n2=0.29;n1ac=0.604;n3=1.359;
E1=83.3*1e3;E2=143.5*1e3;E1ac=83.3*1e3;E3=94.962*1e3; %J/mol
Bmax=120;
Xc=0.40680; Xh1=0.22136; Xh2=0.03786;Xh3=0.03330; %wt
xylmax=Xh1*Bmax; arabmax=Xh2*Bmax;acetmax=Xh3*Bmax;celmax=Xc*Bmax; %g/L
%---------------------------------------------------------------------
% Effective coefficients
ef1=1/(1+(xyl/xylmax)^20);
ef2=1/(1+(arab/arabmax)^20);
ef3=1/(1+(acet/acetmax)^20);
ef4=1/(1+(cel/celmax)^20);
% Effective concentrations
xylef=0.95*ef1*xyl+((1-ef1)*xylmax);
arabef=0.95*ef2*arab+((1-ef2)*arabmax);
acetef=0.95*ef3*acet+((1-ef3)*acetmax);
celef=0.95*ef4*cel+((1-ef4)*celmax);
% K min^(-1)
K1=(k10*1e10)*(HA^n1)*exp(-E1/(R*T));
K2=(k20*1e15)*(HA^n2)*exp(-E2/(R*T));
K3=(k30*1e10)*(HA^n3)*exp(-E3/(R*T)) ;
K1ac=k1ac*(10^10)*(HA^n1ac)*(xyl/xyl_zero)^2*exp(-E1ac/(R*T));
% Rates
rxyl=-K1*xylef;
rarab=-K1*arabef;
racet=-K1ac*acetef;
rcel=-K3*celef;
rXylo=(1/0.88)*K1*xylef-K2*Xylo;
rArabi=(1/0.88)*K1*arabef-K2*Arabi;
rAcetic=(1/0.7)*K1ac*acetef;
rglu=(1/0.9)*K3*celef;
% Mass balance
dxyldt=rxyl;
darabdt=rarab;
dacetdt = acet_dot;
d2acetdt2 = (racet - 10*acet_dot)/400;
daceticdt = Acetic_dot;
d2aceticdt2 = (rAcetic - 10*Acetic_dot)/400;
dceldt=rcel;
dXylodt=rXylo;
dArabidt=rArabi;
dgludt=rglu;
dfurfdt=0.64*K2*(Xylo+Arabi);
%-----------------------------------
f =[ dxyldt; darabdt; dacetdt; d2acetdt2; dceldt; dXylodt; daceticdt; d2aceticdt2; dArabidt; dgludt; dfurfdt];
end
The Xo is another file that i run first and then I run then on cmd : [t, x] = ode45(@(t,y) acidpretreatment(t,y,xyl_zero), [0 240], [x0]);
How it is possible to sovle this ode for T from 100 to 180.
댓글 수: 0
답변 (1개)
VBBV
2021년 12월 2일
편집: VBBV
2021년 12월 2일
T = 100:10:180;
x0 =100.8;
xyl_zero = 4;
% x0 = rand(1,length(T))
for i = 1:length(T)
[t,x] = ode45(@acidpretreatment, [0 24], x0);
plot(t,x)
hold on
end
function [dxyldt darabdt dacetdt d2acetdt2 daceticdt d2aceticdt2 dceldt dXylodt dArabidt dgludt dfurfdt] = acidpretreatment(T,xyl_zero)
%Components
X = rand(1,11);
xyl= X(1);
arab= X(2);
acet= X(3);
acet_dot = X(4);
cel= X(5);
Xylo= X(6);
Acetic= X(7);
Acetic_dot = X(8);
Arabi = X(9);
glu = X(10);
furf= X(11);
%---------------------------------------------------------------------
%Data
R=8.314;
HA=0.07; % m=948.320,38 kg
%i=100:180
T=375+T; % kelvin "Pressure 10 bar"
k10=2.37;
k20=2.17;
k1ac=2.37;
k30=2.37; %min^-1%
n1=1.51;
n2=0.29;
n1ac=0.604;
n3=1.359;
E1=83.3*1e3;E2=143.5*1e3;E1ac=83.3*1e3;E3=94.962*1e3; %J/mol
Bmax=120;
Xc=0.40680; Xh1=0.22136; Xh2=0.03786;Xh3=0.03330; %wt
xylmax=Xh1*Bmax; arabmax=Xh2*Bmax;acetmax=Xh3*Bmax;celmax=Xc*Bmax; %g/L
%---------------------------------------------------------------------
% Effective coefficients
ef1=1/(1+(xyl/xylmax)^20);
ef2=1/(1+(arab/arabmax)^20);
ef3=1/(1+(acet/acetmax)^20);
ef4=1/(1+(cel/celmax)^20);
% Effective concentrations
xylef=0.95*ef1*xyl+((1-ef1)*xylmax);
arabef=0.95*ef2*arab+((1-ef2)*arabmax);
acetef=0.95*ef3*acet+((1-ef3)*acetmax);
celef=0.95*ef4*cel+((1-ef4)*celmax);
% K min^(-1)
K1=(k10*1e10)*(HA^n1)*exp(-E1/(R*T));
K2=(k20*1e15)*(HA^n2)*exp(-E2/(R*T));
K3=(k30*1e10)*(HA^n3)*exp(-E3/(R*T)) ;
K1ac=k1ac*(10^10)*(HA^n1ac)*(xyl/xyl_zero)^2*exp(-E1ac./(R*T));
% Rates
rxyl=-K1*xylef;
rarab=-K1*arabef;
racet=-K1ac*acetef;
rcel=-K3*celef;
rXylo=(1/0.88)*K1*xylef-K2*Xylo;
rArabi=(1/0.88)*K1*arabef-K2*Arabi;
rAcetic=(1/0.7)*K1ac*acetef;
rglu=(1/0.9)*K3*celef;
% Mass balance
dxyldt=rxyl;
darabdt=rarab;
dacetdt = acet_dot;
d2acetdt2 = (racet - 10*acet_dot)/400;
daceticdt = Acetic_dot;
d2aceticdt2 = (rAcetic - 10*Acetic_dot)/400;
dceldt=rcel;
dXylodt=rXylo;
dArabidt=rArabi;
dgludt=rglu;
dfurfdt=0.64*K2*(Xylo+Arabi);
%-----------------------------------
%f = [dxyldt; darabdt; dacetdt; d2acetdt2; dceldt; dXylodt; daceticdt; d2aceticdt2; dArabidt; dgludt; dfurfdt]
end
Change the number of return arguments to function and try it.
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!