Event function ode45
Afficher commentaires plus anciens
I have this set of differential equations like this.This has to be integrated till y(2) is greater than a certain value
function dy = sys1(t,y)
dy=zeros(10,1);
ms=2.840414792361269e-04;
phi=0.8;
l=8;
m=18;
n=0;
soc=330;
ah=0;
a=(l+m/4-n/2)/phi ;
x=ah*(1+4.76*a)/(1-ah*1.5) ;
nr=1+1.5*x+4.76*a;
xfr=1/nr;
xh2r=x/nr;
xn2r=3.76*a/nr;
xo2r=(a+x/2)/nr;
rpm = 1400;
omega = (2*pi*rpm/60);
t1f=(soc-180)/(6*rpm);
t2f=(340-180)/(6*rpm);
t3f=(540-180)/(6*rpm);
Q=[ - 49.8508 2.9486 -0.3781;
133.9972 -0.7494 0.2096;
- 55.4971 -0.0264 0.0048;
-157.4988 -20.1404 3.3960;
1516.7802 13.3876 1.7898;
545.4127 33.0797 -3.9679;
-241.5668 -15.6721 1.3055;
-646.0548 -42.9936 -1.3144;
10.4278 28.0143 0.1638] ;
coffal=Q.';
als=[1 phi phi^2 ah ah^2 phi*ah phi^2*ah phi*ah^2 phi^2*ah^2];
% engine configurations
Rad = 42.5*10^-3 ; % radius of crank
L = 141*10^-3; % length of connecting rod
f =Rad/L; % crank radius to connecting rod length ratio
B = 77.4*10^-3; % bore of cylinder
r = 10; % compression ratio
Vd = (pi*(B^2)*Rad)/2;
V0 = Vd/(r-1); % clearance volume
r_cyl = B/2 ; % radius of cylinder
Lv = 0.12*B ; % maximum valve lift
yr = 0 ; % exhaust mole fraction
% air fuel ratio
% Sl0 = 0.358 % constant for finding lam v
Tref = 298 ; % reference temperature
Pref = 101325 ; % reference pressure
Mu=(xfr*114+xo2r*32+xn2r*28+xh2r*2)/1000;
ru = 8.314/Mu;
Tu(180) = 340 ;
di = 101325/(ru*Tu(180)) ; % average density at intake
vp = (rpm*2*Rad)/30; % mean piston velocity
Ui = 0.6*vp*(4); % average inlet gas speed .6 vol effcny
% constants
alst=als' ;
Twg = 1000;
sab=coffal*(alst);
sl0=sab(1,1) ;
alpha=sab(2,1);
beta=sab(3,1);
Ru=8.314 ;
fw=1+0.41737*ah ;
Mu=(xfr*114+xo2r*32+xn2r*28+xh2r*2)/1000;
ru = 8.314/Mu;
[cpu ,hu ,cpb ,hb,rb]=cph(y(8),y(9),y(7));
CA=180+6*t*rpm;
dy(6)=(V0/2) *(r-1) * (sind(CA) + (sind(CA)*cosd(CA) / sqrt((1/(f^2))-(sind(CA)^2)) ))*omega;
tu=y(8)
tb=y(9)
Auh=2*y(6)/r_cyl;
dy(4)=Auh*7.67*10^-3* (vp^(1/3)) *(y(7) * y(8))^0.5 * -(y(8) - Twg );
gu=1/(1-(8.314/(Mu*cpu)));
dy(7)=(1/y(6))*((gu-1)*dy(4)-gu*y(7)*dy(6));
dy(8)=((dy(4)+y(6)*dy(7)) /(ms*cpu) );
end
I tried doing this
function [value, isterminal, direction] = event1(t, y)
ms=2.840414792361269e-04;
value = ((y(2)/ms)<0.7);
isterminal = 1; % Stop the integration
direction = 0;
end
opt1 = odeset('Events', @event1);
[b B ]=ode45(@sys2,[t1f t2f],A(end,:),opt1)
but it is not stopping.Please help....
Réponses (0)
Catégories
En savoir plus sur Ordinary Differential Equations dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!