Need help solving 2d heat equation using adi method
7 views (last 30 days)
Show older comments
someone please help me correct this code
% 2D HEAT EQUATION USING ADI IMPLICIT SCHEME
clear all;
clc;
close all;
%%% DEFINING PARAMETERS GIVEN %%%
a = 0.65; % alpha in ft^2/hr
w = 1; % width of bar (feet)
dx = 0.05; % step size in x-direction (feet)
nx = 1+w/dx; % no of steps in x-direction
L = 1; % length of bar (feet)
dy = 0.05; % step size in y-direction (feet)
ny = 1+L/dy; % no of steps in y-direction
tf = 0.1; % Final time
dt = 0.002; % Time Step
Nt = tf/dt+1; % Total number of time steps
%%% DEFINING ADDITIONAL CONSTANTS TO SIMPLIFY EQUATION
Dx = a*dt/(dx^2); % Diffusion number(dx)
A = 2*(1+Dx);
Dy = a*dt/(dy^2); % Diffusion number (dy)
B = 2*(1-Dy);
C = 2*(1+Dy);
D = 2*(1-Dx);
% Define 'x'
for i=1:nx
x(i)=(i-1)*dx;
end
% Define 'y'
for j=1:ny
y(j)=(j-1)*dy;
end
% Define 't'
for k=1:Nt
t(k)=(k-1)*dt;
end
%%% INITIAL CONDITIONS %%%
for k = 1:Nt
for i = 2:nx
for j = 2:ny
T(i,j,k) = 0;
end
end
end
%%% BOUNDARY CONDITIONS %%%
for k=1:Nt
for i = 1:nx
for j = 1:ny
if j==1
T(i,j,k) = 200; % Lower wall
elseif i==nx
T(i,j,k) = 0; % Right Wall
elseif i==1
T(i,j,k) = 200; % Left Wall
elseif j==ny
T(i,j,k) = 0; % Top Wall
end
end
end
end
% DEFINING MATRIX ON RHS
Tx=zeros(nx-2,1);
Ty=zeros(ny-2,1);
for k=1:Nt-1
for j=2:ny-1
for i=2:nx-1
if i==1
Tx(i,j) = Dy*T(i,j+1) + B*T(i,j) + Dy*T(i,j-1) + Dx*T(i,j+1);
elseif i==nx-1
Tx(i,j) = Dy*T(i,j+1) + B*T(i,j) + Dy*T(i,j-1) + Dx*T(i,j+1);
else
Tx(i,j) = dy*T(i,j+1) + B*T(i,j) + dy*T(i,j-1);
end
end
end
% Tridiagonal matrix in x-sweep
% First we transform scalar constant B and Dx into column vectors
BB_M = B*ones(nx-1,1);
DDx_U = Dx*ones(nx-2,1);
DDx_L = DDx_U;
ma_x = diag(BB_M) + diag(-DDx_U,1) + diag(-DDx_L,-1);
% Find the solution at the interior nodes
T_x = (inv(ma_x))*Tx;
% Tridiagonal matrix in y-sweep
CC_M = C*ones(ny-1,1);
DDy_U = Dy*ones(ny-2,1);
DDy_L = DDy_U;
ma_y = diag(CC_M) + diag(-DDy_U,1) + diag(-DDy_L,-1);
% to solve y sweep
T_y = inv(ma_y)*T_x;
% To start over
T=T_y;
end
pcolor(T)
0 Comments
Answers (1)
itrat fatima
on 29 Dec 2019
clear all;
clc;
close all;
%%% DEFINING PARAMETERS GIVEN %%%
a = 0.65; % alpha in ft^2/hr
w = 1; % width of bar (feet)
dx = 0.05; % step size in x-direction (feet)
nx = 1+w/dx; % no of steps in x-direction
L = 1; % length of bar (feet)
dy = 0.05; % step size in y-direction (feet)
ny = 1+L/dy; % no of steps in y-direction
tf = 0.1; % Final time
dt = 0.002; % Time Step
Nt = tf/dt+1; % Total number of time steps
%%% DEFINING ADDITIONAL CONSTANTS TO SIMPLIFY EQUATION
Dx = a*dt/(dx^2); % Diffusion number(dx)
A = 2*(1+Dx);
Dy = a*dt/(dy^2); % Diffusion number (dy)
B = 2*(1-Dy);
C = 2*(1+Dy);
D = 2*(1-Dx);
% Define 'x'
for i=1:nx
x(i)=(i-1)*dx;
end
% Define 'y'
for j=1:ny
y(j)=(j-1)*dy;
end
% Define 't'
for k=1:Nt
t(k)=(k-1)*dt;
end
%%% INITIAL CONDITIONS %%%
for k = 1:Nt
for i = 2:nx
for j = 2:ny
T(i,j,k) = 0;
end
end
end
%%% BOUNDARY CONDITIONS %%%
for k=1:Nt
for i = 1:nx
for j = 1:ny
if j==1
T(i,j,k) = 200; % Lower wall
elseif i==nx
T(i,j,k) = 0; % Right Wall
elseif i==1
T(i,j,k) = 200; % Left Wall
elseif j==ny
T(i,j,k) = 0; % Top Wall
end
end
end
end
% DEFINING MATRIX ON RHS
Tx=zeros(nx-2,1);
Ty=zeros(ny-2,1);
for k=1:Nt-1
for j=2:ny-1
for i=2:nx-1
if i==1
Tx(i,j) = Dy*T(i,j+1) + B*T(i,j) + Dy*T(i,j-1) + Dx*T(i,j+1);
elseif i==nx-1
Tx(i,j) = Dy*T(i,j+1) + B*T(i,j) + Dy*T(i,j-1) + Dx*T(i,j+1);
else
Tx(i,j) = dy*T(i,j+1) + B*T(i,j) + dy*T(i,j-1);
end
end
end
% Tridiagonal matrix in x-sweep
% First we transform scalar constant B and Dx into column vectors
BB_M = B*ones(nx-1,1);
DDx_U = Dx*ones(nx-2,1);
DDx_L = DDx_U;
ma_x = diag(BB_M) + diag(-DDx_U,1) + diag(-DDx_L,-1);
% Find the solution at the interior nodes
T_x = (inv(ma_x))*Tx;
% Tridiagonal matrix in y-sweep
CC_M = C*ones(ny-1,1);
DDy_U = Dy*ones(ny-2,1);
DDy_L = DDy_U;
ma_y = diag(CC_M) + diag(-DDy_U,1) + diag(-DDy_L,-1);
% to solve y sweep
T_y = inv(ma_y)*T_x;
% To start over
T=T_y;
end
pcolor(T)
0 Comments
See Also
Categories
Find more on Thermal Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!