Mathematical engine modeling (FHM) on Matlab

[Altairs]

I am trying to replicate a JAVA based applet on MATLAB. Applet can be seen here (http://www.engr.colostate.edu/~allan/thermo/page6/EngineParm1/engine.html).

The theory behind the functions of the applet is here (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)

[x,P] = ode45 ('final', [-180 180], 1)
and then
plot(x,P)

Note the amplitude

final.m
function dydx = final(x,P)
dydx = -1.4*(P/v1(x))*dvol(x)+((1.4-1)/v1(x))*HRF(x);

v1.m
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
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
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 :(

[laraque2003]

Hi,

I get the same problem...

Did you get the solution of this problem ?

[Altairs]

I do not remember axactly, as this is quite an old post, but I think the problem had something to do with angle conversion from degree to radian or vice versa...that is why the amplitude was not coming out correct...but once correct, the results were beautiful :D
Try playing around with angle..if you are using deg convert to radians...if you are using radians try degree...

[laraque2003]

Thanks !

Do you still have your code in order to compare with mine ?

[gmosbo01]

Must work with radians and then convert at the end back to degrees for the plot of pressure vs crank angle. Also, the burn fraction (xb) is not a constant but rather a function of crank angle as well. See the website below for the calculation.
http://www.engr.colostate.edu/~allan/thermo/page6/page6.html

In order for units to match, keep in mind that 1 bar is 10x5 Pa or 100 kPa.