Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Finite Difference Method using Matlab

  1. Mar 20, 2005 #1

    mpm

    User Avatar

    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 dont know that those two things can qualify as the same thing.
     
  2. jcsd
  3. Mar 21, 2005 #2

    Clausius2

    User Avatar
    Science Advisor
    Gold Member

    I think you should describe how is your problem and how you plan to solve it, and maybe we could advice you.
     
  4. Mar 21, 2005 #3

    mpm

    User Avatar

    I guess a more specific problem would help. I dont want to type it all out so here is a link.

    http://www.me.ttu.edu/files/pantoya/S05 - Matlab Finite Diff DP2 asgnmnt.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 Im 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 with cavity example.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.
     
  5. Mar 21, 2005 #4
    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.
     
  6. Mar 22, 2005 #5

    Clausius2

    User Avatar
    Science Advisor
    Gold Member

    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: Mar 22, 2005
  7. Mar 22, 2005 #6
    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 [tex] u^{n+1} = u^n + Au^n[/tex], 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, [tex] Au=f [/tex], 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.
     
  8. Mar 22, 2005 #7

    Clausius2

    User Avatar
    Science Advisor
    Gold Member

    What inversion are you talking about???

    Take a look at my post, Have I ever mentioned inverting anywhere?? :tongue2: :biggrin: :rofl:

    :biggrin: :biggrin: :biggrin:
     
  9. Mar 22, 2005 #8

    mpm

    User Avatar

    Ok. I have edited my professor's code and taken out the cavity. Or at least I think I did. Im 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.
     
  10. Mar 22, 2005 #9
    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.
     
  11. Mar 22, 2005 #10

    mpm

    User Avatar

    Here is my code that I have. Sorry I cut and pasted it on here, but I dont have a server to put the file on. Im getting it to run, but there are problems somewhere because my plot doesnt 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')
     
  12. Mar 22, 2005 #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.
     
  13. Mar 22, 2005 #12

    mpm

    User Avatar

    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.
     
  14. Mar 22, 2005 #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 [tex] \Delta T < \Delta X^2[/tex] to keep it stable. You'll figure it out when you do the stability analysis.
     
    Last edited: Mar 22, 2005
  15. Mar 22, 2005 #14

    mpm

    User Avatar

    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.
     
  16. Mar 22, 2005 #15

    mpm

    User Avatar

    Also if anyone else has any suggestions that will make things clearer, feel free to add comments.
     
  17. Mar 22, 2005 #16

    mpm

    User Avatar

    %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')
     
  18. Mar 22, 2005 #17

    mpm

    User Avatar

    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
     
  19. Mar 22, 2005 #18

    mpm

    User Avatar

    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')
     
  20. Mar 23, 2005 #19

    mpm

    User Avatar

    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 dont 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')
     
  21. Mar 23, 2005 #20
    You need to delete the END statments 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
    [tex]a*\frac{\partial T}{\partial t} = \nabla^2 T= \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} [/tex]

    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 [tex] \frac{\partial T}{\partial x} = \frac{T_{i+1} - T_i}{\Delta x}[/tex] 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, [tex]\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}[/tex]

    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 [tex]T_{i+1,j}, T_{i,j}, T_{i-1,j}, T_{i,j+1}, T_{i,j-1}[/tex]

    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 (Text):
     
    %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')

     
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Finite Difference Method using Matlab
Loading...