- #1
Altairs
- 127
- 0
I am trying to replicate a JAVA based applet on MATLAB. Applet can be seen http://www.engr.colostate.edu/~allan/thermo/page6/EngineParm1/engine.html" .
The theory behind the functions of the applet is http://www.engr.colostate.edu/~allan/thermo/page6/page6.html" .
It is actually a finite heat release model to give pressure, volume and work curves with respect to crank angle from -180 to 180.
It is quite simple actually. Simple integration using Runge-Kutta Order 4 method and this ODE is also of the first order. I am going to post the code below.
Problem is that the shape of the curve is alright but the amplitude (pressure) is totally unrealistic...it shoots to X10^76 or something where as originally in the applet its just 8800 KPa.
Can someone please check what I am doing wrong ? I have spent about a week and still no luck and something tells me that whatever it is it is very small (thats why I can see it :D).
Call the function by writing following on main
Initial conditions used are the same as used in the applet. (theta(x) from -180 to 180 and initial pressure is 1 bar)
and then
Note the amplitude
final.m
v1.m
dvol.m
HRF.m
Please check. I am so stuck at this :(
The theory behind the functions of the applet is http://www.engr.colostate.edu/~allan/thermo/page6/page6.html" .
It is actually a finite heat release model to give pressure, volume and work curves with respect to crank angle from -180 to 180.
It is quite simple actually. Simple integration using Runge-Kutta Order 4 method and this ODE is also of the first order. I am going to post the code below.
Problem is that the shape of the curve is alright but the amplitude (pressure) is totally unrealistic...it shoots to X10^76 or something where as originally in the applet its just 8800 KPa.
Can someone please check what I am doing wrong ? I have spent about a week and still no luck and something tells me that whatever it is it is very small (thats why I can see it :D).
Call the function by writing following on main
Initial conditions used are the same as used in the applet. (theta(x) from -180 to 180 and initial pressure is 1 bar)
Code:
[x,P] = ode45 ('final', [-180 180], 1)
Code:
plot(x,P)
Note the amplitude
final.m
Code:
function dydx = final(x,P)
dydx = -1.4*(P/v1(x))*dvol(x)+((1.4-1)/v1(x))*HRF(x);
v1.m
Code:
function vol1=v1(x)
R=3;
Vd=(pi/4)*(0.1^2)*(0.1);
vol1=(Vd/9)+(Vd/2)*(R+1-cosd(x)-(R^2-sind(x)^2)^(1/2));
dvol.m
Code:
function dvolume = dvol(x)
R=3;
Vd=(pi/4)*(0.1^2)*(0.1);
dvolume=(Vd/2)*sind(x)*[1+cosd(x)*(R^2-sind(x)^2)^(-1/2)];
HRF.m
Code:
function hrf = HRF(x)
n=3;
a=5;
thetaD=40;
thetaS=0;
Qin=1800;
xb=0.99;
if(x<thetaS || x>(thetaS+thetaD))
hrf=0;
else
hrf = (n*a*(Qin/thetaD)*(1-xb)*((x-thetaS)/thetaD)^(n-1));
end;
Please check. I am so stuck at this :(
Last edited by a moderator: