Plot PV diagram for 2 arrays

Load the data from Fall2024_PV_Diagram.mat and plot the PV diagram. Note: You must first calculate the inst. cylinder volume from the crank angle and engine geometry data provided. How do I plot PV diagram I have calculated inst cylinder volume. The pressure and inst cylinder volume are arrays

[Shubham]

Could you show us what you have done and what issues are you facing?

[Parth]

I am given a .mat file which has 2 arrays crank angle and pressure. So i found out inst volume from crank angle bore stroke and connecting rod length(as the assignment said me to do). Now the assignment says plot PV diagram and when I am using plot it is coming as a graph.

Answer 1

[Pavl M.]

Dear Friends,

Dear prospective employer,

Attn.: Important Priority Level 1:

1. This is quite very valuable and interesting project of Engines Cycles quantitative visualization and may have many applications in HEV. I deal with Advanced Propulsion and HEV.
2. Which cycles would you like to dissect: Otto, Miller, Atkinson, which?

In order it comes as a diagram and more than a curve (like you see) it is needed to devide the curve into phases and invoke multiple matlab plot functions per intervals.

Please accept next my transcoded and modified with additions solution.

% Program to calculate different state variables in an Otto Cycle, calculate its thermal efficiency and plot its p-v diagram
clc;
clear;
close all;

%Simplified arrow addition function for plot (auxilliary help function):
function [xfigure, yfigure]=axescoord2figurecoord(varargin)
% AXESCOORD2FIGURECOORD Transform axes coordinates in current
% figure units coordinate to the figure for annotation location
% [xfigure, yfigure]=axescoord2figurecoord(xaxes,yaxes)
% [xfigure, yfigure]=axescoord2figurecoord(xaxes,yaxes,handle_axes)
%
% Ex.
%       % Create some data
% 		t = 0:.1:4*pi;
% 		s = sin(t);
%       % Add an annotation requiring (x,y) coordinate vectors
% 		plot(t,s);ylim([-1.2 1.2])
%       set(gcf,'Units','normalized');
% 		xa = [1.6 2]*pi;
% 		ya = [0 0];
% 		[xaf,yaf] = axescoord2figurecoord(xa,ya);
% 		annotation('arrow',xaf,yaf)
% Acknowledgments are due to Scott Hirsch ([e-mail]) for is
% function ds2nfu. Some part of the present function derived from ds2nfu.
%
% Valley Benoît / Jan 2007

% Process inputs
narginchk(2,3);
if nargin==2
    xaxes=varargin{1};
    yaxes=varargin{2};
    h_axes = get(gcf,'CurrentAxes');
else
    xaxes=varargin{1};
    yaxes=varargin{2};
    h_axes = varargin{3};
end
% get axes properties
funit=get(get(h_axes,'Parent'),'Units');
% get axes properties
aunit=get(h_axes,'Units');
darm=get(h_axes,'DataAspectRatioMode');
pbarm=get(h_axes,'PlotBoxAspectRatioMode');
dar=get(h_axes,'DataAspectRatio');
pbar=get(h_axes,'PlotBoxAspectRatio');
xlm=get(h_axes,'XLimMode');
ylm=get(h_axes,'YLimMode');
xd=get(h_axes,'XDir');
yd=get(h_axes,'YDir');
% set the right units for h_axes
set(h_axes,'Units',funit);
axesoffsets = get(h_axes,'Position');
x_axislimits = get(h_axes, 'xlim');     %get axes extremeties.
y_axislimits = get(h_axes, 'ylim');     %get axes extremeties.
x_axislength = x_axislimits(2) - x_axislimits(1); %get axes length
y_axislength = y_axislimits(2) - y_axislimits(1); %get axes length
% mananged the aspect ratio problems
set(h_axes,'units','centimeters');
asc=get(h_axes,'Position');
rasc=asc(4)/asc(3);
rpb=pbar(2)/pbar(1);
if rasc<rpb
    xwb=axesoffsets(3)/rpb*rasc;
    xab=axesoffsets(1)+axesoffsets(3)/2-xwb/2;
    yab=axesoffsets(2);
    ywb=axesoffsets(4);
elseif rasc==rpb
    xab=axesoffsets(1);
    yab=axesoffsets(2);
    xwb=axesoffsets(3);
    ywb=axesoffsets(4);
else
    ywb=axesoffsets(4)*rpb/rasc;
    yab=axesoffsets(2)+axesoffsets(4)/2-ywb/2;
    xab=axesoffsets(1);
    xwb=axesoffsets(3);
end

if strcmp(darm,'auto') & strcmp(pbarm,'auto')
    xab=axesoffsets(1);
    yab=axesoffsets(2);
    xwb=axesoffsets(3);
    ywb=axesoffsets(4);
end
% compute coordinate taking in account for axes directions
if strcmp(xd , 'normal')==1
    xfigure = xab+xwb*(xaxes-x_axislimits(1))/x_axislength;
else
    xfigure = xab+xwb*(x_axislimits(2)-xaxes)/x_axislength;
end
if strcmp(funit,'normalized');
    xfigure(find(xfigure>1))=1;
    xfigure(find(xfigure<0))=0;
end
if strcmp(yd , 'normal')==1
    yfigure = yab+ywb*(yaxes-y_axislimits(1))/y_axislength;
else
    yfigure = yab+ywb*(y_axislimits(2)-yaxes)/y_axislength;
end
if strcmp(funit,'normalized');
    yfigure(find(yfigure>1))=1;
    yfigure(find(yfigure<0))=0;
end
set(h_axes,'Units',aunit); % put axes units back to original state
end

%Visualization(Interface) function:
function addarrow(h, ev, varargin)

%if verLessThan('matlab', '8.4.0')
    if ishandle(h) && isempty(ev) && strcmp(get(h,'type'), 'line')
        % initial call or button-down, line handle is already h
    elseif strcmp(class(h), 'schema.prop') %
        h = ancestor(get(ev, 'AffectedObject'), 'line');
    end
%else
%    if ishandle(h) && strcmp(h.Type, 'line')
        % initial call and button-down
%    else
 %       h = ev.AffectedObject;
%    end
%%end

npv = length(varargin);

% Normalize figure

hfig = ancestor(h, 'figure');
figunit = get(hfig, 'units');
set(hfig, 'units', 'normalized');

% Create/update arrow
x = get(h, 'xdata');
y = get(h, 'ydata');

[xa, ya] = axescoord2figurecoord(x(end-1:end),y(end-1:end));

arrowexists = isappdata(h, 'arrow');
if arrowexists
    harrow = getappdata(h, 'arrow');
end
% Default options
p = inputParser;
if arrowexists
    p.addParameter('HeadLength', get(harrow, 'headlength'));
    p.addParameter('HeadWidth' , get(harrow, 'headwidth'));
    p.addParameter('HeadStyle' , get(harrow, 'headstyle'));
else
    p.addParameter('HeadLength', 10);
    p.addParameter('HeadWidth' , 10);
    p.addParameter('HeadStyle' , 'vback2');
end
p.addParameter('Color'     , get(h, 'color'));
p.addParameter('LineStyle' , get(h, 'linestyle'));
p.addParameter('LineWidth' , get(h, 'linewidth'));

% Override with user options
if npv > 0
    p.parse(varargin{:});
    Opt = p.Results;
end

if arrowexists
    set(harrow, 'x', xa, 'y', ya, ...
             'color', Opt.Color, ...
             'headlength', Opt.HeadLength, ...
             'headwidth', Opt.HeadWidth, ...
             'headstyle', Opt.HeadStyle, ...
             'LineStyle', Opt.LineStyle, ...
             'LineWidth', Opt.LineWidth);

else
    harrow = annotation('arrow', xa, ya, ...
             'color', Opt.Color, ...
             'headlength', Opt.HeadLength, ...
             'headwidth', Opt.HeadWidth, ...
             'headstyle', Opt.HeadStyle, ...
             'LineStyle', Opt.LineStyle, ...
             'LineWidth', Opt.LineWidth);
    setappdata(h, 'arrow', harrow);
end
% Return figure units to their starting value
set(hfig, 'units', figunit);

    addlistener(h, 'Color',     'PostSet', @addarrow);
    addlistener(h, 'LineStyle', 'PostSet', @addarrow);
    addlistener(h, 'LineWidth', 'PostSet', @addarrow);
    set(h, 'buttondownfcn', @addarrow)

end

% Function for determining the behaviour of volume during isentropic
% process (Main Logic Function)
function V = engine_kinematics(bore, stroke, cr_length, comp_ratio, start_crank_angle, end_crank_angle)
    % Engine geometry properties
    a = stroke / 2;
    R = cr_length / a;
    % Volume parameters
    v_swept = (pi / 4) * (bore^2) * stroke;
    v_clearance = v_swept / (comp_ratio - 1);
    theta_startcrank = deg2rad(start_crank_angle);
    theta_endcrank = deg2rad(end_crank_angle);
    num_val = 100;
    dtheta = (theta_endcrank - theta_startcrank) / (num_val - 1);
    V = zeros(1, num_val);

    for i = 0:num_val-1
        theta = theta_startcrank + i * dtheta;
        term_1 = 0.5 * (comp_ratio - 1);
        term_2 = R + 1 - cos(theta);
        term_3 = sqrt(R^2 - sin(theta)^2);
        V(i + 1) = (1 + term_1 * (term_2 - term_3)) * v_clearance;
    end
end

% Taking Inputs variable
% Engine geometry properties
% load('EngGeom.mat')
bore = 0.18;
stroke = 0.2;
cr_length = 0.8;
comp_ratio = 12;
gamma = 1.4;                 % for air gamma= 1.4
c_v = 0.714;                    % specfic heat of air at constant volume KJ/Kg K
% State variables assumed
p1 = 101325;
T1 = 350;
T3 = 2000;
% Calculating swept volume and clearance volume
v_swept = (pi/4)*(bore^2)*stroke;
v_clearance = v_swept/(comp_ratio-1);
v1 = v_swept+v_clearance;
v2 = v_clearance;
% Calculating state variables at different points
% State variables at point 2
% Process 1-2 is an isentropic comp process (pv^gamma=constant)
p2 = p1*(comp_ratio^gamma);
% (p1*v1)/T1 = (p2*v2)/T2 | T2 = p2*v2*T1/(p1*v1)
T2 = p2 * v2 * T1 / (p1 * v1);
V_comp = engine_kinematics(bore, stroke, cr_length, comp_ratio, 180, 0);
C_const = p1 * (v1^gamma);
P_comp = zeros(1, length(V_comp));
P_comp = C_const./(V_comp.^gamma);
% State variables at point 3
% Process 2-3 is an isochoric heat addition process (volume=constant)
v3 = v2;
% Isochoric process    ''' (p3*v3)/T3 = (p2*v2)/T2 | p3 = p2*T3/(T2) '''
p3 = p2 * T3 / T2;
V_expan = engine_kinematics(bore, stroke, cr_length, comp_ratio, 0, 180);
C_const = p3*(v3^gamma);
P_expan = zeros(1, length(V_expan));
P_expan = C_const./(V_expan.^gamma);
% State variables at point 4
% Process 3-4 is an isentropic expan process (pv^gamma=constant)
v4 = v1;
p4 = p3 * (v3^gamma) / (v4^gamma);
% (p3*v3)/T3 = (p4*v4)/T4 | T4 = p4*v4*T3/(p3*v3)
T4 = p4 * v4 * T3 / (p3 * v3);
% Calculating thermal efficiency
%thermal_efficiency = 1 - (1 / (comp_ratio^(gamma - 1)));
%fprintf('Thermal efficiency of given Otto Cycle is : %.2f%%\n', thermal_efficiency * 100);
% Plotting
figure;
h1 = plot(V_comp, P_comp,'DisplayName', '1-2 Isentropic comp');
%addarrow(h1,[],'color', 'r', 'headwidth', 4)
hold on
h2 = plot([v2, v3], [p2, p3],'DisplayName','2-3 Isochoric heat addition');
addarrow(h2,[],'headwidth', 4)
h3 = plot(V_expan, P_expan,'DisplayName','3-4 Isentropic expan');
addarrow(h3,[],'headwidth', 4)
h4 = plot([v4, v1],[p4, p1],'DisplayName', '4-1 Isochoric heat rejection');
addarrow(h4,[],'headwidth', 4)
addarrow(h1,[],'headwidth', 4) %substantial workaround
%ylim([-900000, 9500000]);
%xlim([0, 0.02]);
xlabel('Volume (m^3)');
ylabel('Pressure (Pa)');
title('P-V diagram of an Otto Cycle');
%To set as per engine's geometry
text(V_comp(1),P_comp(1), '1');
text(v2,p2, '2');
text(v3,p3, '3');
text(v4, p4, '4');
legend('Location', 'best');
hold off;

%Constructed from needing help code by
    %https://independent.academia.edu/PMazniker
    %[phone]
    %https://join.skype.com/invite/oXnJhbgys7oW
    %https://diag.net/u/u6r3ondjie0w0l8138bafm095b
    %https://github.com/goodengineer
    %https://orcid.org/0000-0001-8184-8166
    %https://willwork781147312.wordpress.com/portfolio/cp/
    %https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w