MATLAB Finite Difference Method using Matlab

AI Thread Summary
A user seeks assistance with a Matlab project involving 2-D transient finite difference methods for a heat transfer class. They express frustration with their professor's teaching quality and inquire about the feasibility of using separation of variables in Matlab, questioning if it aligns with finite difference methods. Participants in the discussion suggest providing a clearer problem description and share links to relevant resources, including example codes and assignments. They emphasize the importance of adjusting boundary conditions and simplifying the code for better clarity and efficiency. Suggestions include using matrix representations to streamline calculations and ensuring stability by adhering to specific conditions for time and space discretization. The user shares their modified code, seeking feedback on errors and boundary conditions, while others offer guidance on necessary adjustments to achieve accurate results. Overall, the conversation focuses on troubleshooting Matlab code for heat transfer simulations and improving understanding of finite difference methods.
mpm
Messages
82
Reaction score
0
I am curious to know if anyone has a program that will solve for 2-D Transient finite difference.

I have a project in a heat transfer class and I am supposed to use Matlab to solve for this. However, when I took the class to learn Matlab, the professor was terrible and didnt teach much at all.

Any help would be appreciated.

Also, can I possibly use like a separation of variables process in Matlab and that count as finite difference?

I know there is always more than one way to do things in Matlab and still get the same thing but I don't know that those two things can qualify as the same thing.
 
Physics news on Phys.org
mpm said:
I am curious to know if anyone has a program that will solve for 2-D Transient finite difference.

I have a project in a heat transfer class and I am supposed to use Matlab to solve for this. However, when I took the class to learn Matlab, the professor was terrible and didnt teach much at all.

Any help would be appreciated.

Also, can I possibly use like a separation of variables process in Matlab and that count as finite difference?

I know there is always more than one way to do things in Matlab and still get the same thing but I don't know that those two things can qualify as the same thing.

I think you should describe how is your problem and how you plan to solve it, and maybe we could advice you.
 
I guess a more specific problem would help. I don't want to type it all out so here is a link.

http://www.me.ttu.edu/files/pantoya/S05%20-%20Matlab%20Finite%20Diff%20DP2%20asgnmnt.pdf

The professor gave us a simliar example and he wrote some code for it. I think I might be able to manipulate it into fitting my problem but I am not exactly sure because the conditions are different on each side of the block.

Here is a link to his code.

http://www.me.ttu.edu/files/pantoya/ftcs_block_cavity.m

Also here is a link to his problem he solved.

http://www.me.ttu.edu/files/pantoya/Block%20with%20cavity%20example.doc

Can maybe someone tell me what part needs to be changed to alter the conditions I need? If I know which part needs to be changed maybe I can figure it out. I think I know which parts to remove for the hole in the block of his problem.


Any help is appreciated.
 
Holy cow that's some messy code. Do yourself a favor and represent the entire surface with a single vector. Then your lapacian becomes a tridiagonal matrix, and the boundary conditions are much easier to enforce that way (you don't need a million if elseif statements).

But if you don't know how to do that, it looks like you just pretty much use the code that was given, but you need to alter the boundary conditions to match your problem. The problem you are given has much simpler conditions than the one that the code is given for, so it shouldn't be that hard.

Then it looks likes he's using a Forward Euler scheme for the time integration, and Forward Euler often blows up when it's not supposed to. There's a parameter involving delta T and the coefficients of the equation that controls the stability of Forward Euler. You want to place this value in the unit disc on the left half of the real axis to make Forward Euler stable.
 
LeBrad said:
Holy cow that's some messy code.

The fact is that professor has shown not to have any idea about elegance at programming. That problem can be solved in 60 lines or so.
 
Last edited:
I don't see the need for matrix inversion or Gauss-Siedel because he's already solving explicitly for the values at the new time. The heat equation explicitly in time is u^{n+1} = u^n + Au^n, and he's given u(x,y,0). I was just suggesting that by making A a matrix and u a vector he could reduce all the stepping around the surface to one matrix-vector product.

Now if he were solving Poisson's equation, Au=f, then he would need matrix inversion or Gauss-Siedel. But even then, I would recommend Krylov subspace methods over both and still go with matrix-vector form.
 
LeBrad said:
I don't see the need for matrix inversion or Gauss-Siedel because he's already solving explicitly for the values at the new time.

What inversion are you talking about?

Take a look at my post, Have I ever mentioned inverting anywhere?? :-p :biggrin: :smile:

:biggrin: :biggrin: :biggrin:
 
Ok. I have edited my professor's code and taken out the cavity. Or at least I think I did. I am sure I need to change some boundaries but I am not quite sure what I need to change yet. Since you say it is messy, maybe there is an easier way for me to do this.

I noticed you mention Laplace and triadiagonal. I know of a code that does this I think but I am not sure how I need to manipulate it.

Here is the link.
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=5894&objectType=file

I apologize for having to click the download button. This was just on the Matlab file sharing thing and I think this may be what you speak of.

If this is what you speak of, maybe you can help me manipulate it into working.

I appreciate everyones help here. I feel a little resentful in the fact that I was jipped of learning Matlab because my professor was angry about getting a sophomore level class when he usually teaches senior level classes.

Thanks to all.
 
Last edited by a moderator:
I'm afraid I have only created more confusion here. Forget everything I said about matrices. You can do it the way your professor did and it will work fine. Why don't you show us what you have changed his code to and we can tell you if you're on the right track.
 
  • #10
Here is my code that I have. Sorry I cut and pasted it on here, but I don't have a server to put the file on. I am getting it to run, but there are problems somewhere because my plot doesn't come out right.

Any help with my errors would be appreciated. I probably need to change a coordinate with an i or j somewhere and maybe some other things.

Thanks.

mpm


clear
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;

dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes

Cx = round(0.4/dx); Cy = round(0.25/dy);

tf = 60; %units of seconds
dt = 10; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny);
T(1,:,1)=100;

time(1) = 0;
time(2) = dt;
t=1;
while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction
if i==1; %identify left boundary
T(i,j,t+1)=100;
elseif i>=2 & i<=(Cx-1);
if j==1; %lower boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny; %upper boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
else %all other vertical intermediate nodes on left boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
T(i,1,t+1)=T(i,2,t+1);
if j==1; %lower boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny; %top boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
if j==1; %bottom
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
elseif i==Nx;
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + -10)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx - 10 + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + -10 + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);
end
end
end
end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end

figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e)
T4(e) = T(end,center_y,e)
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')
 
  • #11
Ok, first of all it looks like Cx and Cy were coordinates for the cavity in the other example, so those shouldn't show up anywhere in your code. You erased most of them but still have a crucial elseif depending on one. In fact, if you just change that Cx to Nx, I think you will find that you get a very nice looking solution. Unfortunately it's not the solution to the problem you want to be solving.

Second, when you cut out one elseif statement about i==Cx you also need to kill the stuff below it. As it is now your first elseif statement contains repeated equations for the top and bottom boundaries.

Most importantly though, rather than just changing the geometry, you also need to change the actual boundary conditions. Your professer's problem held the left side at a constant temperature and had a flux leaving the right side. Your problem appears to be heating the left side and holding the right side constant. So your staments regarding the right side, i==Nx, should just set T = 100. Now since you are heating the left side of the box, i == 1, you want to be increasing the temperature at that boundary by some amount at each time step. You can figure out exactly what to put by thinking about what 2000 W means.

Hopefully this points you in the right direction.
 
  • #12
I just want to make sure I get a couple of things straight.

In the Original code, CLx and Cly were used for the cavity dimensions. Are you sure Cx and Cy are for the cavity as well?

If that is the case, do I need to remove any line of code that involves Cx and Cy or do I need to change that to Nx and Ny.

As for the elseif statements, do I need to remove all of that stuff involved with elseif and Cx? Or do I need to remove all elseif statements? Just want to be clear on this.

Once I get to this part, I may can handle the rest about the boundary conditions. I will just have to see. I just want to make sure all of my geometry coding is correct before I mess with boundary conditions.

One last question, if I mess dx and dt, should I keep those relatively small (ie.e < 1cm) or larger (>1cm)? It seems like there might be a stability problem if I go too small. Is that the problem or is it just the coding causing the problem?

Thanks for all your help LeBrad.
 
  • #13
I think Clx and Cly are the size of the cavity, but Cx and Cy are where it starts or stops on the domain. I admit that I didn't look at it long enough to figure out what they were referring to, but since your domain is just a big rectangle, the only boundaries you should need are Nx and Ny.

You are working your way from bottom to top and then left to right (assuming (1,1) is the bottom left corner), so I think the general outline should be something like this.

if x == left wall,
if y = bottom, do something - bottom left corner​
elseif y = middle, do something - left side​
elseif y = top, do something - top left corner​
end​
elseif x is in the middle
if y == bottom, do something - bottom​
elseif y = middle, do something -- This will handle most of the points​
elseif y = top, do something - top​
end​
elseif x is on the right boundary
if y = bottom do something - bottom right corner​
elseif y = middle do something - right side​
elseif y = top do something - top right corner​
end​
end

Those statements should only depend on Nx, Ny, 1, and 1.

As for timesteps, this is just a wild guess from my limited knowledge of semi-discrete Fourier analysis, but the fact that you have a Laplacian makes me say set \Delta T &lt; \Delta X^2 to keep it stable. You'll figure it out when you do the stability analysis.
 
Last edited:
  • #14
LeBrad,


Thanks for your time and help on this. Hopefully I can figure it out from what you have told me. At least I have a direction to go towards. Hopefully I am not annoying you with my questions.

I have until Thursday to finish this so hopefully I can figure it out.
 
  • #15
Also if anyone else has any suggestions that will make things clearer, feel free to add comments.
 
  • #16
%ftcs_block_cavity.m
%
%forward in time centered in space
%based on heat transfer energy balance at all nodes
%
%two dimensional finite difference solution to a rectangular block with
%a rectangular crack or cavity

clear
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;


dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes


tf = 60; %units of seconds
dt = 1; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny);
T(1,:,1)=100;

time(1) = 0;
time(2) = dt;
t=1;
while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction
if i==1; %identify left boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + 2000)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);
end

if j==1; %lower boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny; %upper boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
else %all other vertical intermediate nodes on left boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

T(i,1,t+1)=T(i,2,t+1);
if j==1; %lower boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny; %top boundary
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

if j==1; %bottom
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end

if j==1;
elseif i==Nx;
T(i,j,t+1)=100;
end
end
end
end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end

figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e);
T4(e) = T(end,center_y,e);
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')
 
  • #17
Here is my latest code

while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction

if i==1; %identify left boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + 2000)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);

elseif i>1 & i<Nx; %identify middle boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);

elseif i ==Nx; %identify right boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j>1 & j<Ny;
T(i,j,t+1) = [100];
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);

end
end
end
end
end

end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end

figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e);
T4(e) = T(end,center_y,e);
end
 
  • #18
Here is my last post of code.

%ftcs_block_cavity.m
%
%forward in time centered in space
%based on heat transfer energy balance at all nodes
%
%two dimensional finite difference solution to a rectangular block with
%a rectangular crack or cavity

clear
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;


dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes


tf = 60; %units of seconds
dt = 1; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny);
T(1,:,1)=100;

time(1) = 0;
time(2) = dt;
t=1;
while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction

if i==1; end %identify left boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/(dx) + dx*(T(i,j+1,t)-T(i,j,t))/dy + 2000)+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + 2000 + dx*(T(i,j-1,t)-T(i,j,t))/dy)+T(i,j,t);
end
end

if i>1 & i<Nx; end %identify middle boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j>1&j<Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
end

if i==Nx; end %identify right boundary
if j==1;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
elseif j>1 & j<Ny;
T(i,j,t+1) = [100];
elseif j==Ny;
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*(dy*(T(i-1,j,t)-T(i,j,t))/dx + dy*(T(i+1,j,t)-T(i,j,t))/dx + dx*(T(i,j-1,t)-T(i,j,t))/(dy))+T(i,j,t);
end
end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end

figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e);
T4(e) = T(end,center_y,e);
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')
 
  • #19
LeBrad, Here is my coding as of 3/23/05 1:06 AM Central Time.

LeBrad,

I took your adivce and copied and pasted my professors T equations. I manipulated them into what I thought was right. 5 of the 9 equations work I think. They don't give me problems at least when I run them.

The other 4 I guess I didnt manipulate in the write way. They are marked as comments. if you have any suggestions feel free to give them.

Here is my code.

Thanks for everything.


%ftcs_block_cavity.m
%
%forward in time centered in space
%based on heat transfer energy balance at all nodes
%
%two dimensional finite difference solution to a rectangular block with
%a rectangular crack or cavity

clear
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;


dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes


tf = 60; %units of seconds
dt = 1; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny);
T(1,:,1)=100;

time(1) = 0;
time(2) = dt;
t=1;
%while time(t+1)<=tf %index time step increments
for i=1:Nx; %index in x direction
for j=1:Ny; %index in y direction

if i==1; end %identify left boundary
if j==1; %bottom left corner
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t) + 2000 +T(i,j,t);
elseif j>1&j<Ny; %left side
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + 2000 + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
elseif j==Ny; %top left corner
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + 2000 + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
end
end

if i>1 & i<Nx; end %identify middle boundary
if j==1; %bottom
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t) + T(i+1,j,t)-T(i,j,t) + T(i,j,t);
elseif j>1&j<Ny; %middle
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t) + T(i-1,j,t) + T(m,j+1,t) + T(i,j-1,t) + T(i,j,t);
elseif j==Ny; %top
% T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i+1,j,t)-T(i,j,t) + T(i,j-1,t)-T(i,j,t) + T(i,j,t);
end
end

if i==Nx; end %identify right boundary
if j==1; %bottom right corner
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t)+T(i,j,t);
elseif j>1 & j<Ny; % right side
T(i,j,t+1) = [100];
elseif j==Ny; %top right corner
T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
end
end
t = t+1;
time(t+1) = time(t)+dt;
end

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
x(i) = x(i-1)+dx;
end
for j=2:Ny;
y(j) = y(j-1)+dy;
end
time = [0:dt:tf];
figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
T1(e) = T(2,center_y,e);
T4(e) = T(end,center_y,e);
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')
 
  • #20
You need to delete the END statements that appear like this
if i ==1; end;
You have 3 of those that shouldn't be there.

I fixed your time counting, and I changed the line near the top that says
T(Nx,:,1) = 100; That sets the right side to 100 to start. Before it was setting the left side.

You need another set of parentheses in your equations so your coefficients multiply all of the differences.

This is basically the heat equation
a*\frac{\partial T}{\partial t} = \nabla^2 T= \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}

It is saying that the change in temperature with respect to time depends on second spatial derivatives of the temperature. Since we don't know the function T, we can't take its derivatives. But we can approximate them as differences. The derivative of a line just tells you its slope. And the slope is just the rise over the run. So we can approximate a first derivative by saying \frac{\partial T}{\partial x} = \frac{T_{i+1} - T_i}{\Delta x} where T(i+1) and T(i) are the temperatures at two nodes adjacent in the x direction.

To get a second derivative, you just take the derivative of two first derivative points, \frac{\partial^2 T}{\partialx x^2} = \frac{\frac{T_{i+1} - T_i}{\Delta x} - \frac{T_{i}-T_{i-1}}{\Delta x}}{\Delta x}

To get the formulas for the derivatives with respect to y, just replace the i's with j's.

So when you put it all together, your equation for the temperature at a node T(i,j) will be a function of T_{i+1,j}, T_{i,j}, T_{i-1,j}, T_{i,j+1}, T_{i,j-1}

This means that the demperature at a point depends on the temperatures of the 4 points around it at the previous step. Now, the heat equation is only definded in the middle of your rectange and ALONG the edges, but not normal to the edges. So any point in the middle of the rectangle will have a formula involving those 5 T() values.

When you're on the left edge, you can't use T(i-1,j) because there is no node there. You can still use j+1 and j-1 though, which means you're still using the second derivative along the boundary. So if you keep those and ditch the T(i-1) term, you're left with T(i,j), T(i,j+1), T(i,j-1), T(i+1,j). This means you have a second derivative along the boundary and a first derivative normal to the boundary, which is what the heat equation and your boundary conditions say.

So, when putting the T equations in the right spot, look at where you are on the boundary and figure out which derivatives you need to worry about there and then that will tell you which indexed terms need to be in that particular equation. You can also think about how a point on the top of the rectangle cannot access T(i,j+1), so that term shouldn't be there. That's another good way to check if you have the right equation in the right spot.

The full equations you will need are mostly the same as the ones from your professor's code. So use this to find which ones go where and understand why.

Here is what you had with a few small things fixed. (If you poste code again, put it between the code tags to save space and preserve indentations.)

Code: 
%ftcs_block_cavity.m
%
%forward in time centered in space
%based on heat transfer energy balance at all nodes
%
%two dimensional finite difference solution to a rectangular block with
%a rectangular crack or cavity

clear all
clc

%define material properties
k = 0.401; %conduction coefficient units = W/cm*K
Cp = 385; %specific heat units = J/kg*K
rho = 0.008933; %density units = kg/cm^3

%define geometry
Lx = 50; Ly = 20;


dx = 1; %differential x length of cells
Nx = Lx/dx; %number of x nodes
dy = 1; %differential y length of cells
Ny = Ly/dy; %Number of Y nodes


tf = 60; %units of seconds
dt = 1; %seconds

Box_x = 0; %define lower left corner of box
Box_y = 0;

%Initialize mesh
T(:,:,1)=25*ones(Nx,Ny); % Set T initial
T(Nx,:,1)=100; % Set right side to initial value


time = [dt:dt:tf]; % make a vector of time values
for t = 1:length(time) % loop through the time values
   % time(count) = count*dt;
    for i=1:Nx %index in x direction
        for j=1:Ny %index in y direction

            if i==1 %identify left boundary
                if j==1 %bottom left corner
                    % T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t) + 2000 +T(i,j,t);
                elseif j>1&j<Ny %left side
                    % T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + 2000 + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
                elseif j==Ny %top left corner
                    % T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t)-T(i,j,t) + 2000 + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
                end
            end

            if i>1 & i<Nx; %identify middle boundary
                if j==1; %bottom
                 %   T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t) + T(i+1,j,t)-T(i,j,t) + T(i,j,t);
                elseif j>1&j<Ny %middle
                    T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i+1,j,t) + T(i-1,j,t) + T(i,j+1,t) + T(i,j-1,t) + T(i,j,t);
                elseif j==Ny %top
                    % T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i+1,j,t)-T(i,j,t) + T(i,j-1,t)-T(i,j,t) + T(i,j,t);
                end
            end

            if i==Nx %identify right boundary
                if j==1 %bottom right corner
                    T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j+1,t)-T(i,j,t)+T(i,j,t);
                elseif j>1 & j<Ny % right side
                    T(i,j,t+1) = [100];
                elseif j==Ny %top right corner
                    T(i,j,t+1) = (dt*k/(rho*Cp*dx*dy))*T(i-1,j,t)-T(i,j,t) + T(i,j-1,t)-T(i,j,t)+T(i,j,t);
                end
            end % end right boundary if
        end % end y loop
    end % end x loop
end % end time loop

x(1) = Box_x;
y(1) = Box_y;
for i=2:Nx;
    x(i) = x(i-1)+dx;
end
for j=2:Ny;
    y(j) = y(j-1)+dy;
end
time = [0:dt:tf];
figure(1); clf;
Tfinal = T(:,:,end);
mesh(y,x,Tfinal)
title(strcat('Solution at time = ',num2str(time(end))));
xlabel('y')
ylabel('x')

true_time = time(1:end-1);
clear time;
time = true_time;

figure(2); clf;
center_y = round(Ly/2);
for e = 1:length(time)
    T1(e) = T(2,center_y,e);
    T4(e) = T(end,center_y,e);
end
plot(time,T1,'k',time,T4,'b')
title('Temp history of specific center line nodes')
xlabel('time')
ylabel('T ^oC')
 
  • #21
Thanks LeBrad for keeping your word.

My code runs now, I just need to mess with my boundary conditions a little more. I am getting a line similar to what I need but I don't think its the right one.

I will try and look at your comments more in depth and see if I can figure it out.

Thanks again, you have been a great help.
 
  • #22
Oh, one more thing. I noticed whenever I change my dt and my dx and dy, my program starts running nuts. Do I need to keep one or all of these in a certain range?

For my assignment I have to change these but I am not sure if there is a boundary I need to remember.

Im guessing this hasto do with the stability criterion.

Is this correct?
 
  • #23
I think I have figured out all of my code. My only concern now is the stability criteria. I am not sure what values for dt I need to use. Also I am not sure if dx effects this as well.

I believe the solution is supposed to converge but its not right now. I havea feeling this has somethign to do with the stability criteria.

Any thoughts?
 
  • #24
Is there any method that you were told to use for finding the stability criteria? I will show you what I remember of one way, but there may be an easier way to do it.

Ok, suppose you're working with the equation \frac{\partial u}{\partial t} = C* \frac{\partial^2 u}{\partial x^2}

Assume a solution of the form \hat{U}(t) = e^{ikx}

Now if you plug it into the discrete equation

\frac{u^{n+1}-u^n}{\Delta t} = \frac{C}{\Delta x^2} (u^n_{i+1} - 2u^n_i + u^n_{i-1})

you will get

\frac{\hat{U}^{n+1}-\hat{U}^n}{\Delta t} e^{ikx} = \frac{C}{\Delta x^2} \hat{U}^n ( e^{ik(x+\Delta x)} - 2e^{ikx} + e^{ik(x-\Delta x)})
Then pull out exp(ikx) on the right and move some stuff from the left to the right to get

\hat{U}^{n+1} e^{ikx}= \hat{U}^n e^{ikx}+ \frac{C \Delta t}{\Delta x^2}\hat{U}^n e^{ikx} (e^{ik\Delta x} -2 + e^{-ik\Delta x})

Now cancel the exponential and pull out U to get

\hat{U}^{n+1} = \alpha \hat{U}^n

where alpha is all the other constant terms gathered up, something like (1+C*delt/delx^2*stuff).

Now look at that equation. It says that the U for the next time step is just U from the previous time step multiplied by a constant. So if you want your method to be stable (meaning as time gets large the U values don't go to infinity), you want alpha to be less than one. So alpha has a dependence on delt/delx^2 in it, so if you write it out and set it less than one, you will find a value that the ratio must be less than. So based on a deltaX, you can pick a deltat to keep it stable.

If you play with it for a while you should be able to find alpha. I don't want to take all the fun out of it for you. A hint is that it involes a trig function (sin^2 I think).
 
  • #25
May i ask a question here?
If I have a 3-D wave equation, I need to solve it by finite difference method...
How can I do?
Thanks
 
  • #26
You are trying to solve \nabla^2 u = c\frac{\partial^2u}{\partial t^2} ?

Over what domain? What are the boundary / initial conditions? Do you have a specific question about it, or are you not sure where to start?
 
  • #27
Can ADI Methods be applied in Cylindrical coordinate system?
Thank you
 
  • #28
Would "any body" like to give me some references about numerical method using MATLAB (journal or e-book)? I have had some references, but very limited edition.
Thanks
 
  • #29
Dear All,

I am also trying to solve, non-linear partial differential equations using finite difference method. They are parabolic partial differential equations but both coupled. Can anybody help me write coding in matlab. I did for single variable (temperature, that change with time and depth) nonlinear equations. But now I have two variable (Stress and temperature). They are interdependent too. Please anybody help me to write the coding

Thanks
Best Regards
Sujeevan
 

Similar threads

MATLAB Can You Solve This Finite Difference Equation Using Matlab?
Replies
1
Views
1K
MATLAB Finite difference numerical integration or ode45?
Replies
9
Views
3K
MATLAB 3D Finite different method using matlab
Replies
3
Views
3K
Matlab report help please (Finite Difference Method)
Replies
5
Views
2K
MATLAB Solving Poisson's Equation in Finite Difference Method with Matlab
Replies
3
Views
2K
MATLAB Writing a matlab program to solve the Kepler Problem (Universal Vars)
Replies
1
Views
3K
MATLAB How Does hObject Work in MATLAB GUI Callbacks?
Replies
6
Views
2K
MATLAB Finite element skeleton code for matlab
Replies
3
Views
2K
MATLAB Using MATLAB outside of the lab
Replies
4
Views
1K
MATLAB Solving transportation problems (LPP) using linprog in MATLAB
Replies
3
Views
4K

Hot Threads

Recent Insights

Back
Top