zibian1.m
%Isobutene(1) Methanol(2) MTBE(3)
function xdot=zibian1(t,x,flag)
if isempty(flag),
T=x(4);
% Universal gas constant(cal/(K·mol))8.314472J/(K·mol)
R=1.987;
P=8*(1.013*1e+5);
% Wilson model binary interaction parameters Wilson模型二元相互作用引數
a(1,1) = 0;
a(1,2) = 169.9953;
a(1,3) = -30.2477;
a(2,1) = 2576.8532;
a(2,2) = 0;
a(2,3) = 1483.2478;
a(3,1) = 271.5669;
a(3,2) = -406.3902;
a(3,3) = 0;
% molar volumes(cm3/mol)
V(1) = 93.33;
V(2) = 44.44;
V(3) = 118.8;
A(1,1)=1;
A(2,1)=V(1)/V(2)*exp(-a(2,1)/(R*T));
A(3,1)=V(1)/V(3)*exp(-a(3,1)/(R*T));
A(2,2)=1;
A(1,2)=V(2)/V(1)*exp(-a(1,2)/(R*T));
A(3,2)=V(2)/V(3)*exp(-a(3,2)/(R*T));
A(3,3)=1;
A(2,3)=V(3)/V(2)*exp(-a(2,3)/(R*T));
A(1,3)=V(3)/V(1)*exp(-a(1,3)/(R*T));
% Activity coefficients 活度系數
G(1)=exp(1-x(1)*A(1,1)/(x(1)*A(1,1)+x(2)*A(1,2)+x(3)*A(1,3))...
-x(2)*A(2,1)/(x(1)*A(2,1)+x(2)*A(2,2)+x(3)*A(2,3))...
-x(3)*A(3,1)/(x(1)*A(3,1)+x(2)*A(3,2)+x(3)*A(3,3)))...
/(x(1)*A(1,1)+x(2)*A(1,2)+x(3)*A(1,3));
G(2)=exp(1-x(1)*A(1,2)/(x(1)*A(1,1)+x(2)*A(1,2)+x(3)*A(1,3))...
-x(2)*A(2,2)/(x(1)*A(2,1)+x(2)*A(2,2)+x(3)*A(2,3))...
-x(3)*A(3,2)/(x(1)*A(3,1)+x(2)*A(3,2)+x(3)*A(3,3)))...
/(x(1)*A(2,1)+x(2)*A(2,2)+x(3)*A(2,3));
G(3)=exp(1-x(1)*A(1,3)/(x(1)*A(1,1)+x(2)*A(1,2)+x(3)*A(1,3))...
-x(2)*A(2,3)/(x(1)*A(2,1)+x(2)*A(2,2)+x(3)*A(2,3))...
-x(3)*A(3,3)/(x(1)*A(3,1)+x(2)*A(3,2)+x(3)*A(3,3)))...
/(x(1)*A(3,1)+x(2)*A(3,2)+x(3)*A(3,3));
% Vapor pressure using Antoine's Equation
A(1) = 20.64559; A(2) = 23.49989; A(3) = 20.71616;
B(1) = -2125.74886; B(2) = -3643.31362; B(3) = -2571.58460;
C(1) = -33.160; C(2) = -33.434; C(3) = -48.406;
Psat(1)=exp(A(1)+B(1)/(T+C(1)));
Psat(2)=exp(A(2)+B(2)/(T+C(2)));
Psat(3)=exp(A(3)+B(3)/(T+C(3)));
% Equilibrium constant for the reaction isobutene+methanol=MTBE 異丁烯+甲醇=MTBE反應的平衡常數
Keq=exp(-10.0982 + 4254.05/T + 0.2667*log(T));
% Modified Raoults Law
y(1)=Psat(1)*G(1)*x(1)/P;
y(2)=Psat(2)*G(2)*x(2)/P;
y(3)=Psat(3)*G(3)*x(3)/P;
RR=x(1)*G(1)*x(2)*G(2)/(Keq*x(3)*G(3));
v=[-1,-1,1];
kf=74.4*exp(-3187/T);
kfref=74.4*exp(-3187/333.35);
Da=0;
xdot(1)=x(1)-y(1)+kf/kfref*(v(1)-sum(v)*x(1))*Da*RR;
xdot(2)=x(2)-y(2)+kf/kfref*(v(2)-sum(v)*x(2))*Da*RR;
xdot(3)=x(3)-y(3)+kf/kfref*(v(3)-sum(v)*x(3))*Da*RR;
xdot(4)=x(1)+x(2)+x(3)-1;
xdot=xdot' % xdot must be a column vector
else
% Return M the mass matirx 回傳質量矩陣
M = zeros(4,4);
M(1,1) = 1;
M(2,2) = 1;
M(3,3) = 1;
xdot = M;
end
zibian2.m
%Isobutene(1) Methanol(2) MTBE(3)
function xdot=zibian2(t,x,flag)
if isempty(flag),
T=x(4);
% Universal gas constant(cal/(K·mol))8.314472J/(K·mol)
R=1.987;
P=8*(1.013*1e+5);
% Wilson model binary interaction parameters Wilson模型二元相互作用引數
a(1,1) = 0;
a(1,2) = 169.9953;
a(1,3) = -30.2477;
a(2,1) = 2576.8532;
a(2,2) = 0;
a(2,3) = 1483.2478;
a(3,1) = 271.5669;
a(3,2) = -406.3902;
a(3,3) = 0;
% molar volumes(cm3/mol)
V(1) = 93.33;
V(2) = 44.44;
V(3) = 118.8;
A(1,1)=1;
A(2,1)=V(1)/V(2)*exp(-a(2,1)/(R*T));
A(3,1)=V(1)/V(3)*exp(-a(3,1)/(R*T));
A(2,2)=1;
A(1,2)=V(2)/V(1)*exp(-a(1,2)/(R*T));
A(3,2)=V(2)/V(3)*exp(-a(3,2)/(R*T));
A(3,3)=1;
A(2,3)=V(3)/V(2)*exp(-a(2,3)/(R*T));
A(1,3)=V(3)/V(1)*exp(-a(1,3)/(R*T));
% Activity coefficients 活度系數
G(1)=exp(1-x(1)*A(1,1)/(x(1)*A(1,1)+x(2)*A(1,2)+x(3)*A(1,3))...
-x(2)*A(2,1)/(x(1)*A(2,1)+x(2)*A(2,2)+x(3)*A(2,3))...
-x(3)*A(3,1)/(x(1)*A(3,1)+x(2)*A(3,2)+x(3)*A(3,3)))...
/(x(1)*A(1,1)+x(2)*A(1,2)+x(3)*A(1,3));
G(2)=exp(1-x(1)*A(1,2)/(x(1)*A(1,1)+x(2)*A(1,2)+x(3)*A(1,3))...
-x(2)*A(2,2)/(x(1)*A(2,1)+x(2)*A(2,2)+x(3)*A(2,3))...
-x(3)*A(3,2)/(x(1)*A(3,1)+x(2)*A(3,2)+x(3)*A(3,3)))...
/(x(1)*A(2,1)+x(2)*A(2,2)+x(3)*A(2,3));
G(3)=exp(1-x(1)*A(1,3)/(x(1)*A(1,1)+x(2)*A(1,2)+x(3)*A(1,3))...
-x(2)*A(2,3)/(x(1)*A(2,1)+x(2)*A(2,2)+x(3)*A(2,3))...
-x(3)*A(3,3)/(x(1)*A(3,1)+x(2)*A(3,2)+x(3)*A(3,3)))...
/(x(1)*A(3,1)+x(2)*A(3,2)+x(3)*A(3,3));
% Vapor pressure using Antoine's Equation
A(1) = 20.64559; A(2) = 23.49989; A(3) = 20.71616;
B(1) = -2125.74886; B(2) = -3643.31362; B(3) = -2571.58460;
C(1) = -33.160; C(2) = -33.434; C(3) = -48.406;
Psat(1)=exp(A(1)+B(1)/(T+C(1)));
Psat(2)=exp(A(2)+B(2)/(T+C(2)));
Psat(3)=exp(A(3)+B(3)/(T+C(3)));
% Equilibrium constant for the reaction isobutene+methanol=MTBE 異丁烯+甲醇=MTBE反應的平衡常數
Keq=exp(-10.0982 + 4254.05/T + 0.2667*log(T));
% Modified Raoults Law
y(1)=Psat(1)*G(1)*x(1)/P;
y(2)=Psat(2)*G(2)*x(2)/P;
y(3)=Psat(3)*G(3)*x(3)/P;
RR=x(1)*G(1)*x(2)*G(2)/(Keq*x(3)*G(3));
v=[-1,-1,1];
kf=74.4*exp(-3187/T);
kfref=74.4*exp(-3187/333.35);
Da=0;
xdot(1)=-x(1)+y(1)+kf/kfref*(v(1)-sum(v)*x(1))*Da*RR;
xdot(2)=-x(2)+y(2)+kf/kfref*(v(2)-sum(v)*x(2))*Da*RR;
xdot(3)=-x(3)+y(3)+kf/kfref*(v(3)-sum(v)*x(3))*Da*RR;
xdot(4)=x(1)+x(2)+x(3)-1;
xdot=xdot' % xdot must be a column vector
else
% Return M the mass matirx 回傳質量矩陣
M = zeros(4,4);
M(1,1) = 1;
M(2,2) = 1;
M(3,3) = 1;
xdot = M;
end
zibianrun.m
c=rand(50,3);
for i=1:7,
i
xO1=0.3;
x02=0.2;
x03=0.5;
T0=350;
% We choose various initial points to solve the system of 7
% Differential Algebraic Equations
X01=0.1*2;
X02=0.1*(i+0.5);
X03=1-X01-X02;
% we stop at time = 20
tf=20;
x0 = [X01 X02 X03 T0];
opts = odeset('Mass','M','MassSingular','yes');
% postive warped time
[t,x] = ode15s('zibian1',[0 tf],x0,opts);
% We plot the transformed compositions in the liquid phase
x1=x(:,1);
x2=x(:,2);
figure(1);
axis([0 1 0 1])
line([0 1],[1 0])
x1=x(:,1);
x2=x(:,2);
hold on
plot(x1,x2,'color',[c(i+1,1) c(i+1,2) c(i+1,3)],'LineWidth',2);
end
for i=1:7,
i
xO1=0.1;
x02=0.1;
x03=0.8;
T0=340;
X01=0.1*2;
X02=0.1*(i+0.5);
X03=1-X01-X02;
tf=20;
x0 = [X01 X02 X03 T0];
opts = odeset('Mass','M','MassSingular','yes');
% negative warped time
[t,x] = ode15s('zibian2',[0 tf],x0,opts);
x1=x(:,1);
x2=x(:,2);
X03=1-X01-X02;
figure(1);
x1=x(:,1);
x2=x(:,2);
hold on
plot(x1,x2,'color',[c(i+1,1) c(i+1,2) c(i+1,3)],'LineWidth',2);
end
% we repeat the same calculations as above several times using
% different starting compositions
for i=1:4,
i
xO1=0.75;
x02=0.1;
x03=0.15;
T0=340;
X02=0.1*1;
X01=0.1*(i+0.5);
X03=1-X01-X02;
tf=100;
x0 = [X01 X02 X03 T0];
opts = odeset('Mass','M','MassSingular','yes');
9
[t,x] = ode15s('zibian1',[0 tf],x0,opts);
x1=x(:,1);
x2=x(:,2);
figure(1);
x1=x(:,1);
x2=x(:,2);
hold on
plot(x1,x2,'color',[c(i+1,1) c(i+1,2) c(i+1,3)],'LineWidth',2);
end
for i=1:4,
i
xO1=0.7;
x02=0.1;
x03=0.2;
T0=300;
X02=0.1*1;
X01=0.1*(i+0.5);
X03=1-X01-X02;
tf=100;
x0 = [X01 X02 X03 T0];
opts = odeset('Mass','M','MassSingular','yes');
[t,x] = ode15s('zibina2',[0 tf],x0,opts);
x1=x(:,1);
x2=x(:,2);
figure(1);
x1=x(:,1);
x2=x(:,2);
hold on
plot(x1,x2,'color',[c(i+1,1) c(i+1,2) c(i+1,3)],'LineWidth',2);
end
錯誤提示:
錯誤使用 daeic12 (line 76)
此 DAE 的索引大于 1。
出錯 ode15s (line 310)
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
出錯 zibianrun2 (line 27)
[t,x] = ode15s('zibian1',[0 tf],x0,opts);
請各位老師幫我看看怎么修改一下!實在是matlab小白不太懂
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/36732.html
標籤:其他開發語言
上一篇:之前評價網路性能的時候用的是單目標,我想用多目標報錯了
下一篇:ts檔案下載后不能播放
