Remove part of (fsurf) surface

조회 수: 19 (최근 30일)
Luke G.
Luke G. 2020년 5월 20일
편집: Ameer Hamza 2020년 5월 20일
Hi All, I'm plotting a symbolic function using fsurf, and it yields the result I want. However, roughly half of the graph is irrelevant and I would like to not display that half of the results:
I started to use fill3 in a crude attempt (below) to hide the information I'm not interested in, but then it also obscures the data I am interested in and doesn't look very clean.
Does anyone have any clever solutions to this problem? The code used is below:
clc, clear all, close all
% DEFINE GEOMETRY CONSTs.
r = 20E-3; % absorber radius, [m]
Ah = pi*r^2; % heating surface area, [m^2]
Ac = Ah; % cooling surface area, [m^2]
l_stroke = 20E-3; % piston stroke, [m]
V1 = 250/1E6; % volume @ position1, [m^3]
% ^^ let's say V1 = 250 mL or one Red Bull can
V2 = Ah*l_stroke+V1; % volume @ position 2 (end of stroke)
V3 = V2; % for now, assuming very little expansion here
V4 = V1;
% AIR CONSTs.
R = 8.314; % IG const., [J/mol-K]
T_amb = 20+273.15; % ambient temp., [K]
P_atm = 101325; % atm. pressure, [Pa]
n = (P_atm*V1)/(R*T_amb); % moles of air, arbitrarily chosen
% THERMAL BOUNDARY CONDITIONS
q = 1E3; % input energy, [W/m^2]
h = 10; % convection coeff., [W/m^2-K]
k = 0.8; % condunction coeff. glass cylinder, [W/m-K]
t = 3E-3; % piston wall thickness, [m]
% DEFINING WORK FUNCTIONS
syms T12 T34
W12(T12) = n*R*T12*log(V2/V1); % Work, [J]
dW12(T12) = q*Ah - (T12-T_amb)*(t/(k*Ac) + 1/(h*Ac))^-1; % Work rate, [J/s]
t12 = W12/dW12; % Time, [s]
W34(T34) = n*R*T34*log(V4/V3); % Work, [J]
dW34(T34) = -(T34-T_amb)*(t/(k*Ac) + 1/(h*Ac))^-1; % Work rate, [J/s]
t34 = W34/dW34; % Time, [s]
W_tot(T12,T34) = W12(T12) + W34(T34);
t_stroke(T12,T34) = t12(T12) + t34(T34);
pwr = W_tot/t_stroke; % our final expression of power, [J/s]
%% 3D PLOT: Power Surface Plot
figure
xylims = [290 1E3 290 1E3];
ph = fsurf(pwr,xylims,'ShowContours','on');
ylabel('T_1_-_2 (K)'); xlabel('T_3_-_4 (K)'); zlabel('Power (J s^-^1)')
zlim([0 q*Ah].*3)
% Image Quality Stuff
brighten(0.3)
ph.EdgeColor = [1 1 1];
ph.AmbientStrength = 0.8;
ax = gca;
cb_ctrl = colorbar;
cb_ctrl.Limits = [0 q*Ah];
caxis([0 q*Ah]);
ax.Box = 'on';
Many thanks in advance,
Luke G.
  댓글 수: 2
darova
darova 2020년 5월 20일
Is it possible to see this?
Luke G.
Luke G. 2020년 5월 20일
@darova, definitely! See below:
clc, clear all, close all
% DEFINE GEOMETRY CONSTs.
r = 20E-3; % absorber radius, [m]
Ah = pi*r^2; % heating surface area, [m^2]
Ac = Ah; % cooling surface area, [m^2]
l_stroke = 20E-3; % piston stroke, [m]
V1 = 250/1E6; % volume @ position1, [m^3]
% ^^ let's say V1 = 250 mL or one Red Bull can
V2 = Ah*l_stroke+V1; % volume @ position 2 (end of stroke)
V3 = V2; % for now, assuming very little expansion here
V4 = V1;
% AIR CONSTs.
R = 8.314; % IG const., [J/mol-K]
T_amb = 20+273.15; % ambient temp., [K]
P_atm = 101325; % atm. pressure, [Pa]
n = (P_atm*V1)/(R*T_amb); % moles of air, arbitrarily chosen
% THERMAL BOUNDARY CONDITIONS
q = 1E3; % input energy, [W/m^2]
h = 10; % convection coeff., [W/m^2-K]
k = 0.8; % condunction coeff. glass cylinder, [W/m-K]
t = 3E-3; % piston wall thickness, [m]
% DEFINING WORK FUNCTIONS
syms T12 T34
W12(T12) = n*R*T12*log(V2/V1); % Work, [J]
dW12(T12) = q*Ah - (T12-T_amb)*(t/(k*Ac) + 1/(h*Ac))^-1; % Work rate, [J/s]
t12 = W12/dW12; % Time, [s]
W34(T34) = n*R*T34*log(V4/V3); % Work, [J]
dW34(T34) = -(T34-T_amb)*(t/(k*Ac) + 1/(h*Ac))^-1; % Work rate, [J/s]
t34 = W34/dW34; % Time, [s]
W_tot(T12,T34) = W12(T12) + W34(T34);
t_stroke(T12,T34) = t12(T12) + t34(T34);
pwr = W_tot/t_stroke; % our final expression of power, [J/s]

댓글을 달려면 로그인하십시오.

채택된 답변

Ameer Hamza
Ameer Hamza 2020년 5월 20일
편집: Ameer Hamza 2020년 5월 20일
The data for the surfaces plotted using fsurf is Read-only, so it is not possible to modified. However, you can use surf() like this and mask the unwanted region using nan
%% 3D PLOT: Power Surface Plot
figure
xylims = [290 1E3 290 1E3];
[T34g, T12g] = meshgrid(linspace(xylims(1), xylims(2), 100), linspace(xylims(3), xylims(4), 100));
PHg = double(subs(pwr, [T12, T34], {T12g, T34g}));
mask = T34g < T12g;
PHg(mask) = nan;
ph = surf(T12g, T34g, PHg);
ylabel('T_1_-_2 (K)'); xlabel('T_3_-_4 (K)'); zlabel('Power (J s^-^1)')
zlim([0 q*Ah].*3)
% Image Quality Stuff
brighten(0.3)
ph.EdgeColor = [1 1 1];
ph.AmbientStrength = 0.8;
ax = gca;
cb_ctrl = colorbar;
cb_ctrl.Limits = [0 q*Ah];
caxis([0 q*Ah]);
ax.Box = 'on';
  댓글 수: 1
Luke G.
Luke G. 2020년 5월 20일
This is brilliant! And very concise. Thanks for answering!

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

darova
darova 2020년 5월 20일
Try numerical solution
F = matlabFunction(pwr);
[X,Y] = meshgrid(300:20:1000);
Z = F(X,Y);
Z1 = (Y-3/4*X)/10; % line equations y = 3/4*x
Z(abs(Z) > 100) = nan; % remove high values
Z(Z1<0) = nan; % remove values above the line
contour(X,Y,Z1,[0 0],'b','linew',3)
surface(X,Y,Z,'facecolor','none')
view(3)
axis vis3d
  댓글 수: 2
darova
darova 2020년 5월 20일
Or bettet using initmesh
F = matlabFunction(pwr);
%%
x = [300 400 1000 1000 300]; % contour
y = [300 300 750 1000 1000];
gd = [2; length(x); x(:); y(:)]; % geometry
dl = decsg(gd); % decomposition
[p,e,t] = initmesh(dl,'hmax',20); % triangulate inside contour
Z = F(p(1,:),p(2,:)); % calculate Z value
ff.vertices = [p' Z(:)];
ff.faces = t(1:3,:)';
ff.facecolor= 'y';
plot(x,y,'linew',3)
patch(ff)
view(3)
axis vis3d
75
Luke G.
Luke G. 2020년 5월 20일
@darova, thanks for these suggestions!

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Surface and Mesh Plots에 대해 자세히 알아보기

태그

제품


릴리스

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by