Remove part of (fsurf) surface

23 views (last 30 days)
Luke G.
Luke G. on 20 May 2020
Edited: Ameer Hamza on 20 May 2020
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 Comments
darova
darova on 20 May 2020
Is it possible to see this?
Luke G.
Luke G. on 20 May 2020
@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]

Sign in to comment.

Accepted Answer

Ameer Hamza
Ameer Hamza on 20 May 2020
Edited: Ameer Hamza on 20 May 2020
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 Comment
Luke G.
Luke G. on 20 May 2020
This is brilliant! And very concise. Thanks for answering!

Sign in to comment.

More Answers (1)

darova
darova on 20 May 2020
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 Comments
darova
darova on 20 May 2020
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. on 20 May 2020
@darova, thanks for these suggestions!

Sign in to comment.

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!