Search code examples
matlabplotfigurepde

MatLab PDE plotting issue


I am pretty new to Matlab enviroment and I am working on a heat transfer simulation in Matlab (2014b). I intend to have a multilayered wall of different materials (for now, there is only a single material - copper) and display the results in a single plot. Everything is working fine until I try to append the third layer of the wall to the plot. Below is the definition of basic variables and geometry for PDE solver (r3-r0):

k = 400; 
rho = 8960; 
specificHeat = 386; 
thick = .01; 
stefanBoltz = 5.670373e-8; 
hCoeff = 1;
ta = 300;
emiss = .5;

c = thick*k;
a = sprintf('2*%g + 2*%g*%g*u.^3', hCoeff, emiss, stefanBoltz);
f = 2*hCoeff*ta + 2*emiss*stefanBoltz*ta^4;
d = thick*rho*specificHeat;

r0 = [3 4 0 1 1 0 1 1 1.3 1.3];
r1 = [3 4 0 1 1 0 0.6 0.6 1 1];
r2 = [3 4 0 1 1 0 0.3 0.3 0.6 0.6];
r3 = [3 4 0 1 1 0 0 0 0.3 0.3];

Now, following piece of code calculates the heat transfer through the bottom layer (r3) of the wall, with the input temperature of 1000 K:

gdm = r3';
g = decsg(gdm, 'R3', ('R3')');
hmax = .1; % element size
[p, e, t] = initmesh(g, 'Hmax', hmax);    

numberOfPDE = 1;
pb = pde(numberOfPDE);
pg = pdeGeometryFromEdges(g);
uBottom = pdeBoundaryConditions(pg.Edges(1),'u',1000);
pb.BoundaryConditions = uBottom;

u = pdenonlin(pb,p,e,t,c,a,f, 'jacobian', 'lumped');
fprintf('Temperature at the top edge of the plate = %5.1f degrees-K\n', ...
 u(4));

figure
pdeplot(p, e, t, 'xydata', u, 'contour', 'on', 'colormap', 'jet')

hold on

After "hold on", I basically repeat the previous code for the "r2" rectangle, with the input temperature of "u(4)" (output of the previous layer), followed by this final bit of code:

hold off
axis([0,1,0,2])
caxis manual
caxis([u(4) 1000]);
colorbar; 

As I said, this is all working, and results for the both first and second layer are in the same plot. But after I repeat the process for the third layer (r1), and plot the results into the original figure (with the "hold off" bit being at the very end of the code, of course), the plot displays only the result for the third layer. I'm not sure if this is some limitation of Matlab or if my solution is wrong, so I would like to ask for a bit of help or some direction. Thanks in advance

Below is the full code for better understanding:

k = 400;
rho = 8960;
specificHeat = 386;
thick = .01;
stefanBoltz = 5.670373e-8;
hCoeff = 1; 
ta = 300;
emiss = .5;

c = thick*k;
a = sprintf('2*%g + 2*%g*%g*u.^3', hCoeff, emiss, stefanBoltz);
f = 2*hCoeff*ta + 2*emiss*stefanBoltz*ta^4;
d = thick*rho*specificHeat;

r0 = [3 4 0 1 1 0 1 1 1.3 1.3];
r1 = [3 4 0 1 1 0 0.6 0.6 1 1];
r2 = [3 4 0 1 1 0 0.3 0.3 0.6 0.6];
r3 = [3 4 0 1 1 0 0 0 0.3 0.3];
%---------------------------------------------------------

gdm = r3';
g = decsg(gdm, 'R3', ('R3')');
hmax = .1; % element size
[p, e, t] = initmesh(g, 'Hmax', hmax);

numberOfPDE = 1;
pb = pde(numberOfPDE);
pg = pdeGeometryFromEdges(g);
uBottom = pdeBoundaryConditions(pg.Edges(1),'u',1000);
pb.BoundaryConditions = uBottom;

u = pdenonlin(pb,p,e,t,c,a,f, 'jacobian', 'lumped');
fprintf('Temperature at the top edge of the plate = %5.1f degrees-K\n', ...
 u(4));

figure
pdeplot(p, e, t, 'xydata', u, 'contour', 'on', 'colormap', 'jet')

hold on

%----------------------------------------------------------------------
gdm = r2';
g = decsg(gdm, 'R2', ('R2')');
hmax = .1; % element size
[p, e, t] = initmesh(g, 'Hmax', hmax);

numberOfPDE = 1;
pb = pde(numberOfPDE);
pg = pdeGeometryFromEdges(g);
uBottom = pdeBoundaryConditions(pg.Edges(1),'u',u(4));
pb.BoundaryConditions = uBottom;

u = pdenonlin(pb,p,e,t,c,a,f, 'jacobian', 'lumped');
fprintf('Temperature at the top edge of the plate = %5.1f degrees-K\n', ...
 u(4));

pdeplot(p, e, t, 'xydata', u, 'contour', 'on', 'colormap', 'jet')

%----------------------------------------------------------------------------
gdm = r1';
g = decsg(gdm, 'R1', ('R1')');
hmax = .1; % element size
[p, e, t] = initmesh(g, 'Hmax', hmax);

numberOfPDE = 1;
pb = pde(numberOfPDE);
pg = pdeGeometryFromEdges(g);
uBottom = pdeBoundaryConditions(pg.Edges(1),'u',u(4));
pb.BoundaryConditions = uBottom;

u = pdenonlin(pb,p,e,t,c,a,f, 'jacobian', 'lumped');
fprintf('Temperature at the top edge of the plate = %5.1f degrees-K\n', ...
 u(4));

pdeplot(p, e, t, 'xydata', u, 'contour', 'on', 'colormap', 'jet')

%----------------------------------------------------------------------------
gdm = r0';
g = decsg(gdm, 'R0', ('R0')');
hmax = .1; % element size
[p, e, t] = initmesh(g, 'Hmax', hmax);

numberOfPDE = 1;
pb = pde(numberOfPDE);
pg = pdeGeometryFromEdges(g);
uBottom = pdeBoundaryConditions(pg.Edges(1),'u',u(4));
pb.BoundaryConditions = uBottom;

u = pdenonlin(pb,p,e,t,c,a,f, 'jacobian', 'lumped');
fprintf('Temperature at the top edge of the plate = %5.1f degrees-K\n', ...
 u(4));

pdeplot(p, e, t, 'xydata', u, 'contour', 'on', 'colormap', 'jet')

%hold off

axis([0,1,0,2])
%caxis manual
caxis([u(4) 1000]);

Solution

  • For some reason (which are not completely clear to me) the pdeplot function internally calls hold off.

    So, to get the desired result you need to add a hold on after each call to pdeplot.