# 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??

4 views (last 30 days)
Kshitija Jadhav on 8 Jun 2021
Commented: Kshitija Jadhav on 10 Jun 2021
% 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
Kshitija Jadhav on 9 Jun 2021
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

Joel Lynch on 9 Jun 2021
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 on 10 Jun 2021
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....

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!