Unable to get the following DAE system of equations solved
조회 수: 1 (최근 30일)
이전 댓글 표시
%% DAE of solution
tSpan = [0.000001,1];
Y0= [0.02646;0.97354;1;101325000;0.0001;0.0001;0.0001;0.0001];
MassM=[1 0 0 0 0 0 0 0;
0 1 0 0 0 0 0 0;
0 0 1 0 0 0 0 0;
0 0 0 1 0 0 0 0;
0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0];
opt=odeset('Mass', MassM);
[tSol,YSol]=ode15i(@(t,Y) YTHFcountercurrentflowFun(t,Y), tSpan, Y0, opt);
%% Plot solutions
plot(tSol, YSol);
function dY = YTHFcountercurrentflowFun (t, Y)
% CROSSFLOW cross flow model
global c1 c2 z ph pl permN2 d do Nf Uf mu1 R T di pi
c1=0.00000000004;
c2=0.0004696;
d=0.0023876;
ph= 101325;
pl= 0.001;
permN2= 0.00000000000000104836;
do= 0.0002;
z = 0.2;
Nf= 10000;
Uf= 0.1;
mu1= 0.00001805;
R= 8.314;
T= 303.15;
di= 0.0001;
pi= 3.14159265;
s= (pi*do*z*Nf*permN2*ph)/(d*Uf);
A= (256*mu1*R*T*Uf*Uf*d)/(pi*pi*permN2*ph*ph*ph*do*di*di*di*di*Nf*Nf);
% Functions
dY(1,1)= (1/Y(3))*(Y(8)/permN2)*(Y(1)-(sqrt(Y(4))*Y(5)))*((Y(1)/Y(5))-1);
dY(2,1)= (1/Y(3))*(Y(2)-(sqrt(Y(4))*Y(6)))*((Y(2)/Y(6))-1);
dY(3,1)= (-(1-Y(4))/Y(7));
dY(4,1)= A* (1- Y(3));
dY(5,1)= Y(5)-(((Y(8)/permN2))*Y(1)*Y(7))/(1-sqrt(Y(4))+(Y(7)*(Y(8)/permN2)*sqrt(Y(4))));
dY(6,1)= Y(6)-(Y(2)*Y(7))/(1-sqrt(Y(4))+(Y(7)*sqrt(Y(4))));
dY(7,1)= Y(7)- Y(1)+ ((Y(8)/permN2)*Y(2));
dY(8,1)= Y(8)-c1*(exp(c2*Y(1)*(sqrt(Y(4))/pl)));
end
댓글 수: 0
답변 (2개)
Walter Roberson
2024년 1월 20일
ode15i passes three parameters to the given function: (t,y,y′)
You are not required to do anything with the third parameter... but if you do not do anything with the third parameter then it is a waste to use ode15i.
댓글 수: 0
Torsten
2024년 1월 21일
편집: Torsten
2024년 1월 21일
The time derivative of the first equation is 1e45 at the beginning. Check your parameters and their units.
%% DAE of solution
tSpan = [0.000001,1];
Y0= [0.02646;0.97354;1;101325000];
[tSol,YSol]=ode15s(@(t,Y) YTHFcountercurrentflowFun(t,Y), tSpan, Y0)
%% Plot solutions
plot(tSol, YSol(:,2));
function dY = YTHFcountercurrentflowFun (t, Y)
% CROSSFLOW cross flow model
c1=0.00000000004;
c2=0.0004696;
d=0.0023876;
ph= 101325;
pl= 0.001;
permN2= 0.00000000000000104836;
do= 0.0002;
z = 0.2;
Nf= 10000;
Uf= 0.1;
mu1= 0.00001805;
R= 8.314;
T= 303.15;
di= 0.0001;
pi= 3.14159265;
s= (pi*do*z*Nf*permN2*ph)/(d*Uf);
A= (256*mu1*R*T*Uf*Uf*d)/(pi*pi*permN2*ph*ph*ph*do*di*di*di*di*Nf*Nf);
Y(8) = -(-c1*(exp(c2*Y(1)*(sqrt(Y(4))/pl))));
Y(7) = -(- Y(1)+ ((Y(8)/permN2)*Y(2)));
Y(6) = -(-(Y(2)*Y(7))/(1-sqrt(Y(4))+(Y(7)*sqrt(Y(4)))));
Y(5) = -(-(((Y(8)/permN2))*Y(1)*Y(7))/(1-sqrt(Y(4))+(Y(7)*(Y(8)/permN2)*sqrt(Y(4)))));
% Functions
dY(1,1)= (1/Y(3))*(Y(8)/permN2)*(Y(1)-(sqrt(Y(4))*Y(5)))*((Y(1)/Y(5))-1);
dY(2,1)= (1/Y(3))*(Y(2)-(sqrt(Y(4))*Y(6)))*((Y(2)/Y(6))-1);
dY(3,1)= (-(1-Y(4))/Y(7));
dY(4,1)= A* (1- Y(3));
%dY
end
댓글 수: 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!