ode45 does not return the right answer for internal combusion engine pressure profile
1 次查看(过去 30 天)
显示 更早的评论
Hi all,
I have been getting plenty of practice the last few days with ODE45 and there is one problem I am trying to solve that it does not return the right answer. So in order to rule out things, the aim of this question is to see if my code's sytanx in the way that I use ODE45 is responsible for any mistakes.
So here is the physical problem (Below is a brief explanation, but the detailed ODE is eqution 2.34 found on page 46 of Colin R Ferguson & Allan Kirkpatric's Internal Combustion Engines (Applied Thermosciences) 3rd Ed):
I have a 1st order ODE which is differentiating pressure (P) with respect to angle (th), inclusive also of two variable coefficients (c1 and c2) that are also functions of the angle (th). The equation is: dP/dth = P*c1(th) + c2(th). I know for a fact that I want to have a tspan from -pi to pi and that the initial condition should be y0=1. The last piece of the puzzle is that the coefficient c2(th) should be set to zero for an angle (th) less than a prespecified angle called thetas (which indicate the start of combusion). So here is my code:
% Inputs
Q_in = 1764 ; % Input heat (J)
P1 = 1e5 ; % Pressure ad BDC (Pa)
b = 0.1 ; % Bore (m)
s = 0.1 ; % Stroke (m)
r = 10 ; % Compression ratio
thetas = -20 ; % Start of combustion
thetad = 40 ; % Duration of combusion
gamma = 1.4 ; % Ratio of specific heats
a = 5 ; % Weibe function energy release constant a
n = 3 ; % Weibe function energy release constant n
% Conversion to rad
thetas = deg2rad(thetas) ;
thetad = deg2rad(thetad) ;
% Non-dimensionalise values
% Heat
v_d = s * pi * b^2 / 4 ;
q_in = Q_in ./ (P1 * v_d / ( 1 - 1/r)) ;
% Solve ODE
tspan = linspace(-pi, pi) ;
y0 = 1 ;
[th, P] = ode45(@(th, P) odefun(th, P, gamma, r, q_in, n, a, thetas, thetad), tspan, y0) ;
figure
plot(rad2deg(th),P )
title('My Pressure Profile')
ylabel('Pressure (non-dimensional)')
xlabel('Crank Angle (\theta)')
grid on
% ODE function during heat release
function dPdth = odefun(th, P, gamma, r, q_in, n, a, thetas, thetad)
%C1 Variable Coefficient
c1num = -gamma * (r - 1 ) * sin(th) / (2 * r) ;
c1den = 1 + (r - 1) * (1 - cos(th)) / (2 * r) ;
c1 = c1num ./ c1den ;
%C2 Variable Coefficient
if thetas < th
x_b = 1 - exp( -a * ((th - thetas) / thetad)^n) ;
c2num = (gamma - 1) * q_in * n * a * ( 1 - x_b) * ((th - thetas) / thetad)^(n - 1) ;
c2den = c1den ;
c2 = c2num ./ c2den ;
else
c2 = 0 ;
end
dPdth = P * c1 + c2 ;
end
So what I end up getting is a 'pressure' profile graph which is drastically different from the solution to the problem (the solution pretains to engine1 on the correct pressure profile graph):
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/229673/image.png)
Thanks for your help in advance,
4 个评论
回答(0 个)
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!