Week 5 - Mid term project - Solving the steady and unsteady 2D heat conduction problem
RESULTS AND DISCUSSION:
The given problem is solved using 'Jacobi ,gauss siedal and successive over relaxation methods',for steady and transient states.
For discretizing the 2D heat conduction, we are employing the 'Central differencing scheme' as it is more accurate than the skewed schemes.
In the 'steady state' temperature remains constant over all the simulation time and hence temperature changes only in the spacial coordinates.
In 'Tansient state' temperature changes with time along with spacial cordinates.Two schemes are used to solve 1)Explicit scheme 2) Implicit scheme
Out of the above two methods explicit scheme gives the converges faster and uses no iterative solvers as we are not using forward time step temperature values to find the current time step temperatures.
In case of implicit scheme, Three solvers are used out of which SOR solver converges faster as it acts as a enhancer for the other two methods along with a value 'omega' which is a scaling factor .Here omega is greater than 1 for which it is known as over relaxation method and also should not be more than an optimal value so that it doesnot over shoot in reaching the required temperature.
After SOR solver, GS solver is faster because it uses the previously solved temp
As we are solving for 10X10 grid points(nx=ny=10) for L=1m ,number of curves in the plots are accordingly.As always heat conducts from high temperature region to lower temperature region unless and until there is no external means of work done on it.Initially the right handside of the mesh is at higher temperature according to given boundary conditions and the core of the mesh is at low temperature when heat conduction starts from the right side of mesh ,temperature increases and cooler core moves to the left at a pace 'alpha' which is known as Thermal diffusivity and the solution becomes stable when error loop criteria is satisfied,where convergence of solution takes place and loop runs till error value is greater than the given tolerance value.
For getting a stable solution total CFL number should be <0.5 otherwise the solution becomes unstale and result in an inappropriate solution or no solution.
CFL=CFL1+CFL2;
CFL1=k1;
CFL2=k2;
where k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
Transient-Implicit scheme iterations table
|
Jacobi-iterations |
Guass -Siedal iterations | SOR- iterations |
| 5676 | 5027 | 4333 |
Steady state iterations table.
| Jacobi-iterations | Guass -Siedal iterations | SOR-iterations |
| 208 | 112 | 93 |
clear all
close all
clc
%solving 2D transient heat conducton equation using implicit scheme
% number of grid points
nx=10;
ny=10;
% defining mesh size
x=linspace(0,1,nx);
y=linspace(0,1,ny);
%defining variables
alpha=1.3 %thermal diffusivity constant
dt=3e-4; %time step
t=0.5 % simulation time
dx=x(2)-x(1);
dy=y(2)-y(1);
% pre defining T
T=300*ones(nx,ny);
%applying boundary conditions
T(:,1)=400;
T(:,end)=800;
T(1,:)=600;
T(end,:)=900;
T(1,1)=500;
T(1,end)=700;
T(end,1)=650;
T(end,end)=850;
% copy of T which changes with in the time step
Told=T;
% copy of T which does not changes with in the time step(to be updated for every next time step)
T_prev_dt=T;
k1=((alpha*dt)/(dx^2));
k2=((alpha*dt)/(dy^2));
CFL=k1+k2
error=9e9;
tol=1e-4;
%number of time steps
nt=t/dt
iterative_solver=1
%jacobi method
if iterative_solver==1
jacobi_iter=1
for nts=1:nt;
error=1;
while error>tol
for i=2:nx-1;
for j=2:ny-1;
term1=(1+(2*k1)+(2*k2))^-1;
term2=k1*term1;
term3=k2*term1;
H=(Told(i-1,j)+Told(i+1,j));
V=(Told(i,j-1)+Told(i,j+1));
T(i,j)=(T_prev_dt(i,j)*term1)+(H*term2)+(V*term3);
end
end
error=max(max(abs(Told-T)));
Told=T;
jacobi_iter=jacobi_iter+1;
end
T_prev_dt=T;
end
end
figure(1)
contourf(x,y,T)
colormap(jet)
colorbar
title_text=sprintf("transient state heat conduction,Jacobi iter=%d",jacobi_iter);
title(title_text);
xlabel("X_axis");
ylabel("Y_axis");
%Gauss siedel method
if iterative_solver==2
gs_iter=1;
for nts=1:nt
error=1;
while error>tol
for i=2:nx-1;
for j=2:ny-1;
term1=(1+(2*k1)+(2*k2))^-1;
term2=(T(i-1,j)+T(i+1,j))*(term1);
term3=(T(i,j-1)+T(i,j+1))*(term1);
T(i,j)=(T_prev_dt(i,j)*term1)+(k1*(term2))+(k2*(term3));
end
end
error=max(max(abs(Told-T)));
Told=T;
gs_iter=gs_iter+1;
end
T_prev_dt=T;
end
end
figure(2);
contourf(x,y,T);
colormap(jet);
colorbar;
title_text=sprintf("transient state heat conduction,GS iterations=%d",gs_iter);
title(title_text);
xlabel("X_axis");
ylabel("Y_axis");
%SOR method
omega=1.1
if iterative_solver==3
sor_iter=1
for nts=1:nt
error=1;
while error>tol
for i=2:nx-1;
for j=2:ny-1;
term1=((1+(2*k1)+(2*k2))^-1);
term2=(Told(i-1,j)+T(i+1,j))*(term1);
term3=(Told(i,j-1)+T(i,j+1))*(term1);
T(i,j)=(1-omega)*Told(i,j)+omega*((T_prev_dt(i,j)*term1)+(k1*(term2))+(k2*(term3)));
end
end
error=max(max(abs(Told-T)));
Told=T;
sor_iter=sor_iter+1;
end
T_prev_dt=T;
end
end
figure(3);
contourf(x,y,T);
colormap(jet);
colorbar;
title_text=sprintf("transient state heat conduction,SOR iter=%d",sor_iter);
title(title_text);
xlabel("X_axis");
ylabel("Y_axis");
%%%% END OF PROGRAM %%%%
%%%%%solving steady state 2d heat conduction equation%%%
clear all
close all
clc
%solving steady state 2d heat conduction equation
% number of grid points
nx=10;
ny=10;
%pre defining T
T=300*ones(nx,ny);
% defining boundary conditions
T(:,1)=400;
T(:,end)=800;
T(1,:)=600;
T(end,:)=900;
T(1,1)=500;
T(1,end)=700;
T(end,1)=650;
T(end,end)=850;
%creating a copy of T
Told=T;
%mesh size
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
k=2*(dx^2+dy^2)/(dx^2*dy^2)
error=10;
tol=1e-4;
% jacobi
iterative_solver=1
if iterative_solver==1
jacobi_iter=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;
jacobi_iter=jacobi_iter+1;
end
end
%plotting
figure(1)
contourf(x,y,T)
colormap(jet);
colorbar;
title_text=sprintf("steady state conduction,jacobi iter=%d",jacobi_iter);
title(title_text);
xlabel("X_axis");
ylabel("Y_axis");
%GS_method
if iterative_solver==2
gs_iter=1
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;
gs_iter=gs_iter+1;
end
end
figure(2)
contourf(x,y,T);
colormap(jet);
colorbar;
title_text=sprintf("steady state conduction,GS iter=%d",gs_iter);
title(title_text);
xlabel("X_axis");
ylabel("Y_axis");
%SOR_method
omega=1.1
if iterative_solver==3
sor_iter=1
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)= (1-omega)*Told(i,j)+omega*((term1+term2)/k);
end
end
error=max(max(abs(Told-T)));
Told=T;
sor_iter=sor_iter+1;
end
end
figure(3)
contourf(x,y,T);
colormap(jet);
colorbar;
title_text=sprintf("steady state conduction,SOR iter=%d",sor_iter);
title(title_text);
xlabel("X_axis");
ylabel("Y_axis");
%%%% END OF PROGRAM %%%%%%%
%%%% %solving 2D transient heat conducton equation using explicit scheme %%%%%
clear all
close all
clc
%solving 2D transient heat conducton equation using explicit scheme
nx=10;
ny=10;
x=linspace(0,1,nx);
y=linspace(0,1,ny);
alpha=1.3
dt=3e-4;
t=0.5
dx=x(2)-x(1);
dy=y(2)-y(1);
%defining T
T=300*ones(nx,ny);
%applying boundary conditions
T(:,1)=400;
T(:,end)=800;
T(1,:)=600;
T(end,:)=900;
T(1,1)=500;
T(1,end)=700;
T(end,1)=650;
T(end,end)=850;
% copy of T
Told=T;
T_prev_dt=T;
k1=((alpha*dt)/(dx^2));
k2=((alpha*dt)/(dy^2));
CFL=k1+k2
error=9e9;
tol=1e-4;
nt=t/dt;
if iter=1
iter=1;
for nts=1:nt;
error=1;
while error>tol;
for i=2:nx-1;
for j=2:ny-1;
term1=Told(i-1,j)-2*T(i,j)+T(i+1,j);
term2=Told(i,j-1)-2*T(i,j)+T(i,j+1);
T(i,j)=T_prev_dt(i,j)+(k1*term1)+(k2*term2);
end
end
error=max(max(abs(Told-T)));
Told=T;
iter=iter+1;
end
T_prev_dt=T;
end
end
%plotting
contourf(x,y,T)
colormap(jet)
colorbar
title_text=sprintf("Transiet state heat conduction Explicit scheme,iter=%d,nts=%d",iter,nts);
title(title_text);
xlabel("X_axis");
ylabel("Y_axis");
%%%%% END OF PROGRAM %%%%%