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(최근 30일)
Kshitija Jadhav
Kshitija Jadhav 2021년 6월 8일
댓글: Kshitija Jadhav 2021년 6월 10일
% 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
Kshitija Jadhav
Kshitija Jadhav 2021년 6월 9일
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
Joel Lynch 2021년 6월 9일
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
Kshitija Jadhav
Kshitija Jadhav 2021년 6월 10일
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....

댓글을 달려면 로그인하십시오.

추가 답변(0개)

제품


릴리스

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by