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

1 view (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
##### 3 CommentsShow 1 older commentHide 1 older comment
the cyclist on 8 Jun 2021
Edited: the cyclist on 8 Jun 2021
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 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.
Joel Lynch on 9 Jun 2021
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 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....

### Categories

Find more on Structural Mechanics in Help Center and File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!