# 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 visualización (últimos 30 días)
Kshitija Jadhav el 8 de Jun. de 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 comentariosMostrar 1 comentario más antiguoOcultar 1 comentario más antiguo
the cyclist el 8 de Jun. de 2021
Editada: the cyclist el 8 de Jun. de 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 el 9 de Jun. de 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

Iniciar sesión para comentar.

Joel Lynch el 9 de Jun. de 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.
##### 5 comentariosMostrar 3 comentarios más antiguosOcultar 3 comentarios más antiguos
Joel Lynch el 9 de Jun. de 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 el 10 de Jun. de 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....

Iniciar sesión para comentar.

### Categorías

Más información sobre Structural Mechanics en Help Center y File Exchange.

R2021a

### Community Treasure Hunt

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

Start Hunting!

Translated by