Heat diffusion problem. Why my program is showing error when I am trying to change 'nx' value i.e. no of rows beyond 13 ? Is there any problem with time or timesteps??

% clearing the workspace screen and variables
clear all;
close all;
clc; % to clean the command window

% Stainless steel 304 grade material is considered for a rectangular domain 
% Assumed parameter values for SS material 
k= 15;                  % Thermal conductivity in W/mK
rho = 8000;             % Density in kg/m^3
cp = 502;               % Specific heat in J/kg K
h = 500;                % Heat transfer coefficient

% formation of 1D matrix of size 5mm by 5mm
L = 0.05;
H = 0.05;
W = 0.001;              %thickness of plate is 1 mm

%number of nodes
nx=13; % Rows
ny=20; % Columns

% x and y vectors
x=linspace(0,L,nx);
y=linspace(0,H,ny);

dx=L/nx;
dy=H/ny;

% Initialization of temperature also called empty matrix
T=zeros(nx,ny);

% Initial Boundary conditions (All temperatures in Kelvin)
Tb1 = 1000;   % Left boundary Temperature 
Tb2 = 100;    % Right boundary Temperature 

%parameters to solve the equation
tmax=60;           %Total time in s
nt = 200;          %Total no of time steps
dt = tmax/nt;      %Each time step in s,It is an incremental change in time for which the governing equations are being solved.
T(:)=Tb2; 

for time=1:nt
    for i=1:nx
        for j=1:ny
                  if(i==1)&&(j==1) % 1st Corner
                  delta_T(i,j)= T(i,j)+(((k*(T(i+1,j)-2*T(i,j)+T(i,j+1))/dx)+(h*(Tb1-T(i,j))))/(dx*rho*cp));
                  
                  elseif(i==1)&&(j==ny) % 2nd Corner
                  delta_T(i,j)= T(i,j)+(((k*(T(i+1,j)-2*T(i,j)+T(i,j-1))/dx)+(h*(T(i,j)-Tb2)))/(dx*rho*cp));
                                   
                  elseif(i==nx)&&(j==1) % 3rd Corner
                  delta_T(i,j)= T(i,j)+(((k*(T(i-1,j)-2*T(i,j)+T(i,j+1))/dx)+(h*(Tb1-T(i,j))))/(dx*rho*cp));
                  
                  elseif(i==nx)&&(j==ny) %4th Corner
                  delta_T(i,j)= T(i,j)+(((k*(T(i-1,j)-2*T(i,j)+T(i,j-1))/dx)+(h*(T(i,j)-Tb2)))/(dx*rho*cp));                      
                                     
                  elseif(j==1)
                  % Left Edge
                  delta_T(i,j)=T(i,j)+(((k*(T(i-1,j)+T(i+1,j)-3*T(i,j)+T(i,j+1))/dx)+(h*(Tb1-T(i,j))))/(dx*rho*cp));
                  
                  elseif(j==ny)
                  % Right Egde
                  delta_T(i,j)=T(i,j)+(((k*(T(i-1,j)+T(i+1,j)-3*T(i,j)+T(i,j-1))/dx)+(h*(T(i,j)-Tb2)))/(dx*rho*cp));
                  
                  elseif(i==1) 
                  % Top Edge
                  delta_T(i,j)=T(i,j)+((k*(T(i+1,j)-3*T(i,j)+T(i,j-1)+T(i,j+1)))/(dx*dx*rho*cp));
                                   
                  elseif(i==nx)
                  % Bottom Edge
                  delta_T(i,j)=T(i,j)+((k*(T(i-1,j)-3*T(i,j)+T(i,j-1)+T(i,j+1)))/(dx*dx*rho*cp));
              
                  else
                  %Middle control volumes
                  delta_T(i,j)=T(i,j)+(k*(T(i-1,j)+T(i+1,j)-4*T(i,j)+T(i,j-1)+T(i,j+1))/(dx*dx*rho*cp));
                  
                  end 
        end 
    end 
    
    T=delta_T; %for updating a temperature each time
  
end

[SALAH ALRABEEI]

what is the error msg!

[the cyclist]

When you say "showing error", what do you mean? I am able to run your code for larger values of nx, and I don't get any MATLAB errors. (The code runs to completion.)

Do you mean that it gives you a result you are not expecting?

[Kshitija Jadhav]

yes its not showing the result which I want.

It should show maximum temperature from left surface then it should reduce until the last coulmn. Why the values are coming nagative on alternate rows once I change 'nx' beyond 13 ? I am unable to find the exact error of my program

Answer 1

[Joel Lynch]

It looks like you are not including the dt value in your stencil. You also need to be mindful of stability when setting the time step size.

Your dx/dy calculation has a fence post error, there are nx-1 increments in an array of nx elements.

You should also eliminate the i/j for loops, they will run very slowly. Best MATLAB practice is to vectorize.

[Kshitija Jadhav]

I am new to this programming things. Never did this before. So I am unable to understand few things that you have explained.

[Joel Lynch]

Sure, it looks like an euler explicit 2-dimensional heat diffusion problem. This is a common problem in numerical methods, so there are many different solutions widely available to reference.

The basic idea is that at every time step, you are updating the nodes of the T matrix, , using the derivative from the previous time step multiplied by change in time . Your code does not include that time-increment on the right-hand side of your equations, so it cannot be integrating correctly. Check your units, they should be equal on both sides of the equation.

A fence-post error is a type of "off-by-one" error. This is an easy mistake to make, as

x=linspace(0,L,nx); 

Has nx points (or fence-posts), but nx-1 increments (fence sections). You should be calculating dx using

dx = L/(nx-1); % or alternativley dx = x(2)-x(1);

Your approach is a good one for learning to program, as it shows the stencils and numerical methods in an intuitive way. However, it will be very slow as the problem scales up, because it updates each i/j value in sequence. The "right" way to implement this scheme in MATLAB (so it is simpler and faster) is to update as many values in a single line of code as possible. For your interior points, this can be done outside of the i/j for loops by creating arrays of indicies:

idx = 2:nx-1; % Array of interior i values
idy = 2:ny-1; % Array of interior j values

Which can then be used to access all interior stencil points at once by shifting them +/- 1 increment.

for time=1:nt
    delta_T(idx,idy) = T(idx,idy)+ (k*dt/(dx*dx*rho*cp)) .* ( ...
                                   T(idx-1,idy) + T(idx+1,idy) -4.*T(idx,idy) + ...
                                   T(idx,idy-1) + T(idx,idy+1) );
end

The boundary conditions can also be done in a similar way, without using two nested for loops.

Note, you are now performing element-wise multiplcation and divison, and MATLAB needs to know you are not trying to do matrix mulitplication. This means we have to use element-wise operators (.* and ./) for any non-scalar math. Technically, you are still doing scalar calculations, but its good practice to learn to use these operators when working with arrays and matricies.

[Kshitija Jadhav]

Thank you very much Joel for your clear explaination. I will try to modify the program.

[Joel Lynch]

No problem! Also, in some of the examples online, they discuss stability conditions for euler explicit.

Right now your time step stize is arbitrary, but there are very specific requirements (CFL conditions) on the time step, depending on the dx/dy grid spacing and the thermal diffusivity (or charactiertic speed) of the problem.

Essentially, if you step too quickly, gradients propogate faster than your numerical methods can account for.

[Kshitija Jadhav]

Okay. All right. I will try to understand the different examples more clearly then I will implement it in my program. Thanks a lot again....