Bad result in 2D Transient Heat Conduction Problem Using BTCS Finite Difference Method implicitly

조회 수: 20 (최근 30일)
Everything seem's ok, but my solution's is wrong.
(Thank's to youtuber Sam R. & Bhartendu for Gauss-Seidel Method)
%% 2D HEAT EQUATION WITH CONSTANT TEMP. AT BC'S
%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////
%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////
%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////
%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////
%\\\\\\\\\\\\\\\\\\\\\\\///////////////////////
%% INITIALIZING
clc
clear
close all
%% PARAMETERIZATION
a = 1e-4; % THERMAL DIFFUSIVITY
t = 200; % TOTAL TIME
nt = 2; % TOTAL NUMBER OF TIME STEPS
dt = t/nt; % TIMESTEP
L = 1; % X = Y
nx = 4; % TOTAL NUMBER OF SPATIAL GRIDS
dx = L/nx; % dX = dY
%% BC'S
T1 = 100; % BCS OF DOMAIN
T2 = 100;
T3 = 100;
T4 = 100;
T5 = 200; % INITIAL TEMP.
Tmax = max([T1,T2,T3,T4,T5]);
%% DEFINITION
m = (L/dx)+1; % NO. OF X STEP (X = Y)
r = (t/dt)+1; % NO. OF TIME STEP
d = a*dt/(dx)^2; % DIFF. NO.
%% CREATING BC's & IC'S
B = zeros(m-2,m-2);
[nr,nc] = size(B);
nT = nr * nc;
B = zeros(m,m);
B(:,end) = T1;
B(end,:) = T2;
B(1,:) = T3;
B(:,1) = T4;
B(2:end-1,2:end-1) = T5; %% EDGES MUST BE REFINED
B1 = B(2:end-1,2:end-1)';
BC_v = B1(:); %IC's
%% CREATE TDMA LHS
AA2(1:nT) = 1+4*d;
AA3(1:nT-1) = -d;
AA5(1:nT-3) = -d;
AA = diag (AA2,0) + diag (AA3, -1) + diag (AA3,1) + ...
diag(AA5,-3)+diag(AA5,3); %% LHS HAS BEEN CREATEAD
%% CREATE RHS
T13(:,:,:) = zeros(size(B,1),size(B,1),r); % Define a cubic matrix
n = size(AA,1);
T13(:,:,1) = B;
BB(n,1) = 0;
for i = 1:n
if mod(i,2) == 1
BB(i,1) = T1+T2;
end
if mod(i,2) ~= 1
BB(i,1) = T1;
end
if mod(n,2)~=1
BB((n+1)/2,1) = 0;
elseif mod(n,2)==1
BB((n+1)/2,1) = 0;
end
end
RHS = BC_v + d*BB;
for k = 1:r
%\\\\\Solution of x in Ax=b using Gauss Seidel Method/////
%\\\\\Solution of x in Ax=b using Gauss Seidel Method/////
C = RHS; % constants vector
n = length(C);
X = zeros(n,1);
Error_eval = ones(n,1);
% Start the Iterative method
iteration = 0;
while max(Error_eval) > 1e-7
iteration = iteration + 1;
Z = X; % save current values to calculate error later
for i = 1:n
j = 1:n; % define an array of the coefficients' elements
j(i) = []; % eliminate the unknow's coefficient from the remaining coefficients
Xtemp = X; % copy the unknows to a new variable
Xtemp(i) = []; % eliminate the unknown under question from the set of values
X(i) = (C(i) - sum(AA(i,j) * Xtemp)) / AA(i,i);
end
Error_eval = sqrt((X - Z).^2);
end
%\\\\\Solution of x in Ax=b using Gauss Seidel Method/////
%\\\\\Solution of x in Ax=b using Gauss Seidel Method/////
RHS = X + d*BB;
T13(2:end-1,2:end-1,k+1) = (vec2mat(RHS',m-2)); % VEC. TO MAT.
T13(:,end,k+1) = T1;
T13(end,:,k+1) = T2;
T13 (1,:, k + 1) = T3;
T13 (:, 1, k + 1) = T4;
end
%% Animation
x = 1:m;
y = 1:m;
T21 = T13 (:,:, 1);
for k = 1:r-1
figure(1)
h = imagesc(x,y,T21);
shading interp
axis([0 m 0 m 0 Tmax])
title({['Transient Heat Conduction'];['time = ',...
num2str((k)*dt),' s']})
colorbar;
drawnow;
xlabel(['T = ',num2str(T2),'C'])
yyaxis left
ylabel(['T = ',num2str(T4),'C'])
yyaxis right
str = 'T0 = 20C' ;
dim = [.6 0.55 .3 .3];
pause(0.1);
refreshdata(h);
if k~=r
T21 = T13(:,:,k+1);
else
break;
end
end
  댓글 수: 8
Leroy Coelho
Leroy Coelho 2020년 7월 7일
I have written the code for this problem
I am more then happy to send the code your way for your refrence
Ragul Kumar
Ragul Kumar 2020년 11월 6일
Hello experts,
I am trying to solve the finite difference methof for crank nicolson scheme to 2d heat equation. please let me know if you have any MATLAB CODE for this
Boundary condition are
If you can kindly send me the matlab code, it will be very useful for my research work . thank you very much.

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

채택된 답변

Leroy Coelho
Leroy Coelho 2020년 7월 8일
편집: Leroy Coelho 2020년 7월 8일
%2-D Transient Heat Conduction With No Heat Generation [FDM][BTCS]
clc;
clear all;
%% Variable Declaration
n = 21; %number of nodes
L = 1; %length of domain
W = 1; %width of domain
alpha = 1e-4; %thermal diffusivity (m^2/s)
m =(n-2)*(n-2); %variable used to construct penta diagonal matrix
sm = sqrt(m);
dx = L/(n-1); %delta x domain
dy = W/(n-1); %delta y domain
x = linspace(0,L,n); %linearly spaced vectors x direction
y = linspace(0,W,n); %linearly spaced vectors y direction
[X,Y]=meshgrid(x,y);
Tin = 200; %internal temperature
T = ones(m,1)* Tin; %initilizing Space
Ta = zeros(n,n);
A = zeros(m,m);
Ax = zeros(m,1);
B = zeros(sm,sm);
dt = 0.1; %time step
tmax = 500; %total Time steps (s)
t = 0 : dt : tmax;
r = alpha * dt /(dx^2); %for stability, must be 0.5 or less
Tt = 100; %Top Wall
Tb = 100; %Bottom Wall
Tl = 100; %Left Wall
Tr = 100; %Right Wall
%% Boundry Conditions
Ta(1,1:n) = Tt; %Top Wall
Ta(n,1:n) = Tb; %Bottom Wall
Ta(1:n,1) = Tl; %Left Wall
Ta(1:n,n) = Tr; %Right Wall
%% Setup Matrix
for ix = 1 : 1 : m
for jx =1 : 1 : m
if (ix == jx)
A(ix,jx) = (1 + 4*r);
elseif ( (ix == jx + 1) && ( (ix - 1) ~= sm * round( (ix-1)/sm) ) ) %RHS
A(ix,jx) = -r;
elseif ( (ix == jx - 1) && ( ix ~= sm*round(ix/sm) ) )
A(ix,jx) = -r;
elseif (ix == jx + sm)
A(ix,jx) = -r;
elseif (jx == ix + sm)
A(ix,jx) = -r;
else
A(ix,jx) = 0;
end
end
end
for iy = 1 : 1 : sm
for jy = 1 : 1 : sm
if (iy == 1) && (jy == 1)
B(iy,jy) = r * (Tl + Tt);
elseif (iy == 1) && (jy == sm)
B(iy,jy) = r * (Tt + Tr); %LHS
elseif (iy == sm) && (jy == sm)
B(iy,jy) = r * (Tb + Tr);
elseif (iy == sm) && (jy == 1)
B(iy,jy) = r * (Tb + Tl);
elseif (iy == 1) && (jy == sm)
B(iy,jy) = r * (Tt + Tr);
elseif (iy == 1)&&(jy > 1 || jy < sm)
B(iy,jy) = r * Tt;
elseif (jy == sm) && (iy > 1 || iy < sm)
B(iy,jy) = r * Tr;
elseif (iy == sm) && ( jy > 1 || jy < sm)
B(iy,jy) = r * Tb;
elseif (jy == 1) && ( iy > 1 || jy < sm)
B(iy,jy) = r * Tl;
else
B(iy,jy) = 0;
end
end
end
Bx = reshape(B,[],1); %Convert matrix to vector
%% Solution
for l = 2 : length(t) %time steps
Xx = ( T + Bx );
Ax = A \ Xx;
T( 1 : m ) = Ax( 1 : m );
fprintf('Time t=%d\n',l-1);
end
Tx = reshape( Ax , sm , sm); %convert vector to matrix
for i = 2 : 1 : n-1
for j = 2 : 1 :n-1
Ta(i,j) = Tx ( i-1 , j-1 );
end
end
%% Plot
contourf(X,Y,Ta,50,'edgecolor','none');
h = colorbar;
ylabel(h, 'Temperature °C')
colormap jet
axis equal
title(['Top (Tt)= ',num2str(Tt),'°C']);
xlabel(['Bottom (Tb)= ',num2str(Tb),'°C'])
yyaxis left
ylabel(['Left (Tl)= ',num2str(Tl),'°C'])
yyaxis right
ylabel(['Right (Tr)= ',num2str(Tr),'°C'])
  댓글 수: 4
Héctor Antonio González Pomares
¿Cómo podría simular esto? Cuando el flujo es bidimensional de calor por conducción en condiciones transitorias para esta placa plana por medio de la aplicacióndel Método de Diferencias Finitas.
temperaturas en °C donde t1=250, t2=270, t3=230, t4=150, t5=210, t6=210, t7=270; t8=110; t9=300; t10=100; t11=290; t12=130 Cada cuadricula con dimensiones a x a donde a=0.25 m

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

추가 답변 (0개)

Community Treasure Hunt

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

Start Hunting!

Translated by