Solving the steady and unsteady 2D heat conduction problem in MATLAB
STEADY STATE HEAT CONDUCTION EQUATION
A function solver2d is made to solve the steady state heat conduction equation implicitly. The initial values are initialised i.e. number of nodes along x and y axis, temperature matrix T, x, y, dx, dy, error and tolerance. The default value of k is also assigned. The iteration counter, count starts from 0. Boundary conditions are assigned of the left, right, top and bottom of the matrix. The corner element values are determined by taking an average of the adjacent values.
If the solver variable called with the function is 1, then Jacobi method is used. The while loop is used for when the error value is greater than tolerance. For loop is created for x and y term, and the solution is divided into 2 terms, which are calculated separately. The 2 terms are added and then divided by the constant k. The error is calculated, the temperature matrix is updated and count is increemented.
If the solver variable called with the function is 2, then Gauss Siedal method is used. It is similar to Jacobi method but in Jacobi method, the old terms of temperature matrix are used to determine updated value but in Gauss Siedal method, the updated terms are used to get the updated value of temperature.
If the solver variable called with the function is 3, then SOR method is used. The method involves used the Jacobi or GS final value multiplied with a constant,a and then added with the product of old temperature term and (1 - a).
Plotting is done used the contourf function along with colorbar, colormap function. Based on the type of solver, the plot is made. The time is calculated at the end of the loop.
function solver2d(solver)
%code to solve steady state heat conduction implicitly
%starting values
ny=10; nx=10;
T=298*ones(nx,ny);
Told=T;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=1; tol=1e-4;
k=(2*(dx^2+dy^2)/(dx^2*dy^2));
[xx,yy]=meshgrid(x,y);
count=0;
%Boundary Conditions
%Tl = 400; Tr=800; Tt=600; Tb=900;
T(1,:)=900; T(end,:)=600; T(:,1)=400; T(:,end)=800;
T(1,1)=(400+900)/2; T(1,end)=(900+800)/2; T(end,1)=(400+600)/2; T(end,end)=(800+600)/2;
Told=T;
tic
%jacobi method
if solver==1
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=((Told(i-1,j)+Told(i+1,j))/(dx^2));
term2=((Told(i,j-1)+Told(i,j+1))/(dy^2));
T(i,j)=(term1+term2)/k;
end
end
error=max(max(abs(Told-T)));
Told=T;
count=count+1;
end
end
%gauss siedal method
if solver==2
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=((T(i-1,j)+T(i+1,j))/(dx^2));
term2=((T(i,j-1)+T(i,j+1))/(dy^2));
T(i,j)=(term1+term2)/k;
end
end
error=max(max(abs(Told-T)));
Told=T;
count=count+1;
end
end
%SOR method
if solver==3
a=1.2;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=((T(i-1,j)+T(i+1,j))/(dx^2));
term2=((T(i,j-1)+T(i,j+1))/(dy^2));
T(i,j)=a*((term1+term2)/k)+ Told(i,j)*(1-a);
end
end
error=max(max(abs(Told-T)));
Told=T;
count=count+1;
end
end
total_time=toc;
%Contour plot
figure(1)
[C,h]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
xlabel(\'X Axis\');
ylabel(\'Y Axis\');
if solver==1
title(sprintf(\'Contour plot of steady state heat conduction using GAUSS JACOBI METHOD. \\n No. of Iterations: %d. \\n Computation Time: %0.4f s.\',count,total_time));
end
if solver==2
title(sprintf(\'Contour plot of steady state heat conduction using GAUSS SEIDAL METHOD. \\n No. of Iterations: %d. \\n Computation Time: %0.4f s.\',count,total_time));
end
if solver==3
title(sprintf(\'Contour plot of steady state heat conduction using SOR METHOD. \\n No. of Iterations: %d. \\n Computation Time: %0.4f s.\',count,total_time));
end
end
TRANSIENT IMPLICIT HEAT CONDUCTION EQUATION
A function tisolver is made to solve the transient heat conduction equation implicitly. The initial values are initialised i.e. number of nodes along x and y axis, temperature matrix T, x, y, dx, dy, error and tolerance. The default value of k1, k2, k are also assigned. The iteration counter, count starts from 0. Boundary conditions are assigned of the left, right, top and bottom of the matrix. The corner element values are determined by taking an average of the adjacent values. The method is same as the steady state solver but the equation is different for transient and steady state equation.
The terms required for time loop are defined, i.e. alpha, dt, t_start, t_stop and number of time steps. A for loop is created for required time steps.
If the solver variable called with the function is 1, then Jacobi method is used. The while loop is used for when the error value is greater than tolerance. For loop is created for x and y term, and the solution is divided into 2 terms, which are calculated separately. The 2 terms are added along with the previous temperature term and then divided by the constant k. The error is calculated, the temperature matrix is updated and count is increemented.
If the solver variable called with the function is 2, then Gauss Siedal method is used. It is similar to Jacobi method but in Jacobi method, the old terms of temperature matrix are used to determine updated value but in Gauss Siedal method, the updated terms are used to get the updated value of temperature.
If the solver variable called with the function is 3, then SOR method is used. The method involves used the Jacobi or GS final value multiplied with a constant,a and then added with the product of old temperature term and (1 - a).
Plotting is done used the contourf function along with colorbar, colormap function. Based on the type of solver, the plot is made. The time is calculated at the end of the loop.
function tisolver(solver)
%code to solve transient heat conduction equation implicitly
%starting values
ny=10; nx=10;
T=298*ones(nx,ny);
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=1; tol=1e-4;
count=0;
%time loop values
alpha=1e-3; dt=0.01;
t_start=0; t_stop=1/alpha;
nt=(t_stop - t_start)/dt;
k1=(alpha*dt)/dx^2; k2=(alpha*dt)/dy^2;
k=1+(2*k1)+(2*k2);
[xx,yy]=meshgrid(x,y);
%Boundary Conditions
%Tl = 400; Tr=800; Tt=600; Tb=900;
T(1,:)=900; T(end,:)=600; T(:,1)=400; T(:,end)=800;
T(1,1)=(400+900)/2; T(1,end)=(900+800)/2; T(end,1)=(400+600)/2; T(end,end)=(800+600)/2;
Told=T; Tprev=Told;
tic;
for l=1:nt
%jacobi
if solver==1
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=k1*(Told(i-1,j)+Told(i+1,j));
term2=k2*(Told(i,j-1)+Told(i,j+1));
T(i,j)=(Tprev(i,j)+term1+term2)/k;
end
end
error=max(max(abs(Told-T)));
Told=T;
count=count+1;
end
end
%gauss siedal
if solver==2
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=k1*(T(i-1,j)+Told(i+1,j));
term2=k2*(T(i,j-1)+Told(i,j+1));
T(i,j)=(Tprev(i,j)+term1+term2)/k;
end
end
error=max(max(abs(Told-T)));
Told=T;
count=count+1;
end
end
%SOR
if solver==3
while(error>tol)
for i=2:nx-1
for j=2:ny-1
term1=((1-alpha)*T(i-1,j)+alpha*(T(i-1,j)+Told(i+1,j)));
term2=((1-alpha)*T(i,j-1)+alpha*(T(i,j-1)+Told(i,j+1)));
T(i,j)=((Tprev(i,j))+(k1*term1)+(k2*term2))/k;
end
end
error=max(max(abs(Told-T)));
Told=T;
count=count+1;
end
end
Tprev=T;
end
total_time=toc;
%plotting
[C,h]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
xlabel(\'X Axis\');
ylabel(\'Y Axis\');
if solver==1
title(sprintf(\'Contour plot of implicit transient heat conduction using GAUSS JACOBI METHOD. \\n Computation Time: %0.4f s.\',total_time));
end
if solver==2
title(sprintf(\'Contour plot of implicit transient heat conduction using GAUSS SEIDAL METHOD. \\n Computation Time: %0.4f s.\',total_time));
end
if solver==3
title(sprintf(\'Contour plot of implicit transient heat conduction using SOR METHOD. \\n Computation Time: %0.4f s.\',total_time));
end
end
TRANSIENT EXPLICIT HEAT CONDUCTION EQUATION
A function tosolver is made to solve the transient heat conduction equation explicitly. The initial values are initialised i.e. number of nodes along x and y axis, temperature matrix T, x, y, dx, dy, error and tolerance. The default value of k1, k2, k are also assigned. The iteration counter, count starts from 0. Boundary conditions are assigned of the left, right, top and bottom of the matrix. The corner element values are determined by taking an average of the adjacent values. The method is same as the steady state solver but the equation is different for transient and steady state equation.
The terms required for time loop are defined, i.e. alpha, dt, t_start, t_stop and number of time steps. A for loop is created for required time steps.
There is no variable called with the function. The value is determined by obtaining 2 terms,i.e. term1 is the sum of the left and right term subtracted by twice of current term multiplied with constant k1, term2 is the sum of top and bottom term subtracted by twice of current term multiplied with constant k2. term1 and term2 are added with old temperature term to obtain new temperature term.
The temperature matrix and count is updated. Plotting is done used the contourf function along with colorbar, colormap function. Based on the type of solver, the plot is made. The time is calculated at the end of the loop.
function tosolver
%code to solve transient heat conduction equation implicitly
%starting values
ny=10; nx=10;
T=298*ones(nx,ny);
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
error=1; tol=1e-4;
count=0;
%time loop values
alpha=1e-3; dt=0.01;
t_start=0; t_stop=1/alpha;
nt=(t_stop - t_start)/dt;
k1=(alpha*dt)/dx^2; k2=(alpha*dt)/dy^2;
[xx,yy]=meshgrid(x,y);
%Boundary Conditions
%Tl = 400; Tr=800; Tt=600; Tb=900;
T(1,:)=900; T(end,:)=600; T(:,1)=400; T(:,end)=800;
T(1,1)=(400+900)/2; T(1,end)=(900+800)/2; T(end,1)=(400+600)/2; T(end,end)=(800+600)/2;
Told=T;
tic;
for l = 1:nt
for i = 2:nx-1
for j = 2:ny-1
term1 = k1*(Told(i-1,j)-2*Told(i,j)+Told(i+1,j));
term2 = k2*(Told(i,j-1)-2*Told(i,j)+Told(i,j+1));
T(i,j) = Told(i,j)+term1+term2;
end
end
Told = T;
count = count+1;
total_time=toc;
end
%plotting
[C,h]=contourf(xx,yy,T);
colorbar;
colormap(jet);
clabel(C,h);
xlabel(\'X Axis\');
ylabel(\'Y Axis\');
title(sprintf(\'Contour plot of explicit transient heat conduction using GAUSS JACOBI METHOD. \\n Computation Time: %0.4f s.\',total_time));
end