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

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,

[Stephan]

Are you sure, that you are permitted to publish a full ebook here?

[KostasK]

Good point, hadn't thought about that to be honest. I am taking it down as it is not very kind to the author.

[Star Strider]

If it is open-source, that should not be a problem.

[Stephan]

unfortunally it is not open source