MATLAB Answers

ode45 does not return the right answer for internal combusion engine pressure profile

조회 수: 2(최근 30일)
Kostas
Kostas 13 Jul 2019
Commented: Stephan 13 Jul 2019
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):
Thanks for your help in advance,

  댓글 수: 4

표시 이전 댓글 수: 1
Kostas
Kostas 13 Jul 2019
Good point, hadn't thought about that to be honest. I am taking it down as it is not very kind to the author.

로그인 to comment.

답변 수 (0)

이 질문에 답변하려면 로그인을(를) 수행하십시오.

제품


Translated by