Plot PV diagram for 2 arrays

Question

Parth Bipin 2024-10-13

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2159315-plot-pv-diagram-for-2-arrays

编辑： Pavl M. 2024-10-13

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

2 个评论
显示无隐藏无

Shubham 2024-10-13

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

Parth 2024-10-13

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. 2024-10-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2159315-plot-pv-diagram-for-2-arrays#answer_1531365

编辑：Pavl M. 2024-10-13

在 MATLAB Online 中打开

Dear Friends,

Dear prospective employer,

Attn.: Important Priority Level 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.
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 (shirsch@mathworks.com) 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

% 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

%+380990535261

%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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Plot PV diagram for 2 arrays

2 个评论
显示无隐藏无

回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Plot PV diagram for 2 arrays

2 个评论 显示 无隐藏 无

回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论