MATLAB Examples

PCYLZEROCORR performs ZERO level correction of cylinder pressure data

Crank angle (C.A.) reference for the calculation is from 0 to 720 deg. Top Dead center (TDC) occurs at 360 deg.

The thermodynamic concept behind the ZERO level correction of the cylinder pressure trace is that the beginning of the compression stroke can be approximated by a polytropic process that obeys:

P*V^k = Constant

Or doing a logarithmic transformation:

log(P) = -k*log(V)+Constant

Where the slope k for a Diesel engine is approximately 1.374

The idea of the ZERO level correction algorithm is then to find the absolute pressure level during the compression stroke that has a slope of -1.374 on the logP-logV domain, within the region where the polytropic hypothesys is valid (usually between 115 and 65 C.A. deg BTDC).

The algorithm consists of a few steps (namely steps 1, 3,4 and 5) that are explained in more detail below.

Contents

Step 1: Defining relevant parameters

The C.A. interval used for the ZERO level correction is defined as well as the gas polytropic coefficient. Engine geometric parameters used to calulate the displaced volume as a function of C.A. are also defined in this step.

% Defining crank angle interval where P*V^k = constant
% (from 115 to 65 deg CA BTDC)
ZeroCAs1 = 245;         % Interval start angle
ZeroCAs2 = 295;         % Interval end angle
ZeroK  = 1.374;         % polytropic coefficient cp/cv

% Approx. crank angle for intake valve closing
WCAivo = 180;

% Defining Engine parameters
EngCyl = 6;             % number of cylinders
EngD   = 0.126;         % bore diameter	(m)
Engs   = 0.166;         % stroke(m)
EngL   = 0.251;         % conrod length	(m)
EngCR  = 18.5;          % compression ratio

Step 2: Loading cylinder pressure data and plotting trace

The pressure data is stored in a CSV file for the purpose of this example but should be readily available inside the data processing software. Some phase-preserving filtering method is recommended.

% Loading CSV file containing cylinder pressure data
% (column 1 <- C.A. [deg] , column 2 <- Pressure [Pa])
Data = csvread('Cylinder Pressure - CA Based.csv');

% Assigning pressure data
CA = Data(:,1);         % Column 1 -> C.A. [deg]
P  = Data(:,2);         % Column 2 -> Cylinder Pressure [deg]

% Filtering data with simple 3-point moving average
P = [P(1); (P(1:end-2)+P(2:end-1)+P(3:end))/3; P(end)];

% Creating figure 1 (for illustration purposes)
Fig1 = figure;
FigPos = Fig1.Position;
Fig1.Position = [FigPos(1) FigPos(2) 600 400];
Fig1.Name = 'Input Pressure Trace';
Fig1.NumberTitle = 'off';
Fig1.MenuBar = 'none';
Fig1.Color = [0.97 0.99 0.99];

% Plotting pressure data from -90 to 90 deg C.A. (for illustration purposes)
Axes1 = axes('Parent',Fig1);
Axes1.Position = [0.08 0.12 0.86 0.8];
Axes1.NextPlot = 'add';
Axes1.FontName = 'Arial';
Axes1.XLim = [-90 90];
Axes1.XTick = -90:30:90;
Axes1.XGrid = 'on';
Axes1.YGrid = 'on';
Plot = plot(Axes1,CA-360,P);
Plot.LineWidth = 1;
xlabel('Crank Angle (deg.)');
ylabel('Pressure (Pa)');

Step 3: Performing LogV calculation

During the initial phase of the compression stroke, the gas state obeys a polytropic transformation of the type P*V^k = constant. Or using the log transformation: log(P) = -k*log(V) + Const . In this step, LogV is calculated as a function of the C.A. array to be used in the next steps.

% Total displaced volume (m3)
Vd = (pi*EngD^2/4)*(2*(Engs/2));

% Dead volume (m3)
Vo = Vd/(EngCR-1);

% Calculating piston axial position as a function of C.A.
Xp = (Engs/2)*((1+EngL/(Engs/2))-cos(pi/180*CA) - ...
    ((EngL/(Engs/2))^2-(sin(pi/180*CA)).^2).^0.5);

% Calculating displaced volume
V = Vo+(pi*EngD^2/4)*Xp;

% Calculating log of displaced volume
LogV = log10(V);

Step 4: Sweeping pressure offset values

Some indeces in the C.A. array are determined for further calculations. The pressure trace is then offset to avoid negative abs. values during the intake stroke. Finally, a sweep of pressure trace offsets is done in order to enable the search for the offset that will have the correct slope -k. For each pressure offset value, a linear fit using least-squares is done to find the best line fit within the C.A. region of interest.

% Finding relevant indeces
iw1 = find(CA>=WCAivo,1);
iz1 = find(CA>=ZeroCAs1,1);
iz2 = find(CA>=ZeroCAs2,1);

% Avoiding negative absolute pressures by shifting to 0.25 bar abs.
% (mainly caused by pressure transducer drift or noise)
P = P - mean(P(1:iw1)) + 0.25e5;

% Creating logarithmically spaced array of 31 pressure offsets
% between 0.5E+4 to 0.5E+6 Pa for search of best offset value
POffset = 0.5*logspace(4,6,31);

% Preallocating variable
C = zeros(length(POffset),2);

% Looping through pressure offset values
for j = 1:length(POffset)

    % Calculating least squares fit coefficients
    C11 = sum(LogV(iz1:iz2).^2);
    C12 = sum(LogV(iz1:iz2));
    C21 = C12;
    C22 = length(LogV(iz1:iz2));
    D1 = sum(LogV(iz1:iz2).*log10(P(iz1:iz2)+POffset(j)));
    D2 = sum(log10(P(iz1:iz2)+POffset(j)));

    % Calculating coefficients of linear fit (y=C1x+C2)
    % (or logP = -k*logV + Const. Where C(:,1) contains slope values)
    C(j,2) = (D2-C21/C11*D1)/(C22-C21/C11*C12);
    C(j,1) = (D1-C12*C(j,2))/C11;

end

% Calculating error function as the difference between
% the calculated values of slopes for different pressure
% offsets and the target value ZeroK defined in step 1
y = C(:,1)'+ZeroK;

Step 5: Determining pressure trace offset (ZERO level correction)

The zero crossing of the error function corresponds to a pressure offset that will have the correct value of the slope. Linear interpolation around the crossing point is used to determine a more accurate value.

% Finding zero crossing point
i0 = find(y<0);

% Calculating pressure offset using linear
% interpolation around zero crossing point
if all(y>0)  ||  isempty(i0)

    % Case where no crossing occurs (all positive values)
    P0 = 0;

else

    if i0 < length(POffset)

        % Case where crossing occurs
        i0 = i0(end);
        P0 = POffset(i0) - (POffset(i0+1)-POffset(i0))/(y(i0+1)-y(i0))*y(i0);

    else

        % Case where no crossing occurs (all negative values)
        P0 = 0;

    end

end

% Applying pressure offset
PCorrected = P + P0;

Step 6: Plotting error function and zero crossing point

The error function as well as the zero crossing point are shown on the figure below along with the segment used for interpolation.

% Creating figure 3 (for illustration purposes)
Fig3 = figure;
FigPos = Fig3.Position;
Fig3.Position = [FigPos(1) FigPos(2) 600 400];
Fig3.Name = 'Error Function Zero Crossing';
Fig3.NumberTitle = 'off';
Fig3.MenuBar = 'none';
Fig3.Color = [0.97 0.99 0.99];

% Plotting ZeroK error function (for illustration purposes)
Axes3 = axes('Parent',Fig3);
Axes3.Position = [0.1 0.12 0.84 0.8];
Axes3.NextPlot = 'add';
Axes3.FontName = 'Arial';
Axes3.XGrid = 'on';
Axes3.YGrid = 'on';
Plot = plot(Axes3,POffset,y);
Plot.LineWidth = 1;
Plot.Marker = '.';
Plot = plot(Axes3,POffset(i0:i0+1),y(i0:i0+1));
Plot.LineWidth = 2;
Plot.Marker = '.';
Plot.MarkerSize = 10;
Plot = plot(Axes3,P0,0,'k+');
Plot.LineWidth = 2;
Plot.MarkerSize = 8;
Text = text(Axes3,P0,0,['   P_0 = ' int2str(P0) ' Pa']);
Text.FontName = 'Arial';
Text.VerticalAlignment = 'top';
xlabel('Pressure Offset [Pa]');
ylabel('ZeroK Error [ - ]');

Step 7: Plotting logP-logV Diagram

Finally, the logP-logV diagram is a good way to observe the effect of the correct offset of the cylinder pressure trace on the compression stroke of the cycle. The average slope in the highlighted region should match the slope of the ideal fit line (black) with slope -k

% Generating fit line for polytropic range of P-V cycle
LogV1 = LogV(iz1);
LogV2 = LogV(iz2);
C0 = C(i0,2) - (C(i0+1,2)-C(i0,2))/(y(i0+1)-y(i0))*y(i0);
LogP1 = -ZeroK*LogV1 + C0;
LogP2 = -ZeroK*LogV2 + C0;

% Creating figure 4 (for illustration purposes)
Fig4 = figure;
FigPos = Fig4.Position;
Fig4.Position = [FigPos(1) FigPos(2) 600 400];
Fig4.Name = 'P-V diagram';
Fig4.NumberTitle = 'off';
Fig4.MenuBar = 'none';
Fig4.Color = [0.97 0.99 0.99];

% Plotting logP-logV diagram with fitted line for PV^k = const. region
Axes4 = axes('Parent',Fig4);
Axes4.Box = 'on';
Axes4.Position = [0.1 0.12 0.84 0.8];
Axes4.NextPlot = 'add';
Axes4.FontName = 'Arial';
Plot = plot(Axes4,LogV,log10(PCorrected));
Plot.LineWidth = 0.75;
Plot = plot(Axes4,LogV(iz1:iz2),log10(PCorrected(iz1:iz2)));
Plot.LineWidth = 1.5;
Plot = plot(Axes4,[LogV1 LogV2],[LogP1 LogP2],'k');
Plot.LineWidth = 2.0;
Text = text(Axes4,0.25*LogV1+0.75*LogV2,(LogP1+LogP2)/2,['Line fit with slope = -' num2str(ZeroK)]);
Text.FontName = 'Arial';
Text.HorizontalAlignment = 'right';
xlabel('Log(V) [ - ]');
ylabel('Log(P) [ - ]');