Plotting/visualizing results of a 3D heat diffusion equation?

Hey guys; I have been playing with this problem for a couple of days now. And I think now it's about time to get some help. :-)

I have a 3D space (i.e. a room) with different dimensions as the width (W), depth (D) and height (H) of the room (Please refer to the attached image, named "geometry"). There is a heat source within the geometry somewhere near the right-back-floor intersection (the location of the heat source is NOT the focus of my question). I am solving the 3D heat diffusion equation to calculate the variation of the temperature within the room, due to the heat source, as the time progresses.

Number of discretization along W, D and H directions are NW, ND and NH, respectively. Therefore, I go ahead and create a 3D matrix for the temperature field with the dimension of Tn=(NH,NW,ND). Note that the first (NH), second (NW) and third (ND) arrays of this matrix are representing height (rows), width (column) and depth (page) of the room (considering the fact that our 3D matrix is starting at point O in the schematic, as per MATLAB convention for the 3D matrix).

I solve the equation using finite difference method (using some initial and boundary conditions) and the result of the temperature field is according to my expectation (based on the inspection of the temperature matrix, Tn). However, I am having difficulty plotting/visualizing the results. Can you help me plot the temperature field in a coordinate system with its origin at point M (refer to the attached image), its x axis lying in W direction, its y in D direction and its z in H direction??

You may recommend meshgrid. But I already tried that with no success. So, let's say I define the x and y and z as three vectors with dimension (1*NW), (1*ND) and (1*NH), respectively. Using [X,Y,Z]=meshgrid(x,y,z) outputs three 3D matrices X, Y and Z with dimensions that is not the same as the Tn. So, I am really stuck.

(My .m code is attached to this question for your reference)

Thanks in advace

Answer 1

[mohammad faghih]

Ok, Thanks KKSV for the help. However, here is a problem. The order of elements in the vectors X(:), Y(:) and Z(:) is not the same as the order of element when you vectorize Tn, that's why when you plot them, they do not match.

So, I had to write a nested loop to re-arrange the order of elements in Tn and store it in another vector called TT to make it match with the order of X(:), Y(:) and Z(:). (By 'order', I mean the way they are arranged)

Here is the nested loop attached to this post for your reference.

I am sure there is a better and more efficient way of doing this, but this is all I can do about it.

[mohammad faghih]

Here is a useful link in regards to this issue:

https://eli.thegreenplace.net/2014/meshgrids-and-disambiguating-rows-and-columns-from-cartesian-coordinates/

Answer 2

[KSSV]

rho=1.2;
cp=50;

%Length in thee three directions (i.e. width, depth and height)
W=100;
D=50;
H=30;

%Number of nodes in the three directions
NW=101;
ND=51;
NH=31;

%Distance between nodes in each of the three directions
w=W/(NW-1);
d=D/(ND-1);
h=H/(NH-1);

%Number of time steps and time interval
Nt=18000;
dt=1;
Total_t=Nt*dt;

%Creating matrix for the temperature
Tn=zeros(NH,NW,ND);

%Creating matrix for thermal conductivity
K=10*ones(NH,NW,ND);

K(:,1,:)=0.001;
K(:,end,:)=0.001;

K(1,:,:)=0.001;
K(end,:,:)=0.001;

K(:,:,1)=0.001;
K(:,:,end)=0.001;

%Initial conditions
Tn(:,:,:)=0;
t=0;
[X,Y,Z] = ndgrid(1:NH,1:NW,1:ND);

%March in time
for k=1:Nt
    k
    %Update time
    t=t+dt;

      %Update the temperature
      Tc=Tn;

      %Calculating new temperature in x and y and z locations
      for i=2:NW-1
          for j=2:ND-1
              for m=2:NH-1
                  Tn(m,i,j)=Tc(m,i,j)+...
                      dt*(K(m,i,j)/rho/cp)*...
                      (((Tc(m+1,i,j)-2*Tc(m,i,j)+Tc(m-1,i,j))/h/h)+...
                      ((Tc(m,i+1,j)-2*Tc(m,i,j)+Tc(m,i-1,j))/w/w)+...
                      ((Tc(m,i,j+1)-2*Tc(m,i,j)+Tc(m,i,j-1))/d/d));
              end
          end
      end
       scatter3(X(:),Y(:),Z(:),1,Tn(:),'filled')
       colorbar
       title(sprintf('Time Step:%d',k))
       drawnow
      %Applying the source term
      if (t<2*3600)
          Tn(20,80,40)= Tn(20,80,40)+(dt*100000/rho/cp);
      end

      %Boundary condition
      Tn(:,1,:)=Tn(:,2,:);
      Tn(:,end,:)=Tn(:,end-1,:);

      Tn(1,:,:)=Tn(2,:,:);
      Tn(end,:,:)=Tn(end-1,:,:);

      Tn(:,:,1)=Tn(:,:,2);
      Tn(:,:,end)=Tn(:,:,end-1);

end

[mohammad faghih]

Thanks for the quick reply; But, what is A in S=size(A)???

[KSSV]

Typo error..it is not required.....removed the line.

[mohammad faghih]

Thanks;

But looking at the figure attached to this comment and considering the location of the source, it does not produce my desired coordinate system. I updated my schematic to show the location of the source.

Based on your coordinate system, heat source is occurring at (x=20,y=80,z=40), however the actual location must be (x=80,y=40,z=10). Refer to the updated schematics and the Coordinate.png file attached to this comment.

[KSSV]

The procedure should remain same...check the meshgrid dimensions.

[laddawan buakaw]

What if there is a cylinder? What is a cylinder command?

Answer 3

[zero trent]

Dear @mohammad faghih, did you correct the problem? Please can you send me the corrected version? I need it. Thank you

[mohammad faghih]

Hi zero trent;

Sorry this post is too old I don't really remember about it.

But I do remember I could find a work around for it, but unfortunately I don't have the files anymore.

Sorry :-(

[zero trent]

Thank you for your reply.

Can you help me in resolving it? Even if you don't have the files, you can't resolve it given that you found around for it.

Because my professor ask me to solve it, and my Matlab knowledges are very poor.

Thank you a lot.