Week 5.1 - Mid term project - Solving the steady and unsteady 2D heat conduction problem

Challenge – Solving the steady and unsteady 2D heat conduction problem.

The challenge is to solve 2D heat conduction equation.

In physics and mathematics, the heat equation is a partial differential equation that describe how the distribution of some quantity (such as heat) evolves over time in a solid medium, as it spontaneously flows from places where it is higher towards where it is lower. It is a special case of the diffusion equation.

This challenge will focus on solving the equation with three solvers, Jacobian, gauss-seidel, and SOR iteration for two case steady and transient.

While we use iterative solvers to solve steady state PDE and both explicit and implicit methods are used to solve transient PDE.

Given,

Generating a 2D Profile.

The Temperature at 4 boundaries.

Top – 400 K,

Bottom – 900 K,

Left – 600 K,

Right – 800 K,

Solving 2D heat conduction equation of steady state:

Since the temperature doesn’t change with time the above equation reduces to zero.

01. Jacobian Iterative Solver – In Jacobian method, each iteration is done using a set of previous value, even if new values are available in the middle of an iteration, like when you are solving for n+1 iteration and for I th node, you will use values of previous (n) iteration even if you have the updated values of the nodes in current iteration. So, this method is called method of simulations displacement.

The 2D steady state heat conduction equation for this method is.

02. Gauss-Seidel Iterative solver - In Gauss-seidel, as soon as you find a new value in the current iteration you will use it for finding other values in the current iteration.

Like when you are solving for n+1 iteration and for I th node you will use all the available values up to i-1 node of current iteration (n+1) and use I + 1 node value and above values of previous iteration(n). So, we get more converged values, so this method is called method of spontaneous displacements.

Equation for this method is.

03. Successive over – Relaxation (SOR) Method – It is a variant of the gauss seidel method for solving a linear system of equations, resulting in faster convergence. You can also use jacobi method. Here an extra term called relaxation factor (Omega) is multiplied to gauss seidel generation term to scale the convergence rate and solve the equation faster. To calculate new value, we take old value and take a parameter omega and take rate of change.

Equation for this method is.

When omega is 0, we get new value is equal to old value which gives wrong result.

When omega goes from 0 to 1, we are switching from no update to an update using gauss siedal or jacobi method.

When omega >1, then we are taking the update fro gauss siedal method and scaling up the convergence rate. The value we give is called over relaxation factor.

When omega <1, the value we give is called under relaxation factor. This is used when the solution given by gauss seidel is over relaxed we give omega value less than 1 to get the converged solution.

Input Variable:

1. Nodes; we are going to take no. of nodes in x – direction, nx and in Y – direction, ny = 20.

2. Length of the domain: This is taken as 1 along x – axis and y – axis.

3. Given, Tolerance assumed to be 1e-4.

4. And we also assumed error as 999. We can assume any value but it must be greater than the required tolerance.

By taking these as input we create a function for each solver and use those function in the main program and solve the 2D steady state equation on a square domain.

Main Program:

clear all

close all

clc

L = 1;

nx = 20;

ny = nx;

tolerance = 1e-4; %temperature matrix

error = 999;

solver_number = input('enter_the_solver_number: ')

if solver_number==1

%function for steady_state_jacobian_iteration_method

steady_jacobian_iteration_method(L,nx,ny,error,tolerance)

end

if solver_number==2

% function for steady_state_gauss_seidel_iteration_method

steady_gauss_seidel_iteration_method(L,nx,ny,error,tolerance)

end

if solver_number==3

% function for steady_state_SOR_method applied on gauss seidel method

steady_gauss_SOR_method(L,nx,ny,error,tolerance)

end

Function Program:

% Function to solve 2D steady state heat conduction equation using jacobi iterative solver.

function steady_jacobian_iteration_method(L,nx,ny,error,tolerance)

tic

x = linspace(0,L,nx);

y = linspace(0,L,ny);

dx = L/(nx-1);

dy = L/(ny-1);

T = 300*ones(nx,ny); %Temperature matrix

T_L = 400;

T_R = 800;

T_T = 600;

T_B = 900;

T(:,1) = T_L;

T(:,end) = T_R;

T(1,:)=T_T;

T(end,:)=T_B;

T(1,1) = (T_L+T_T)/2;

T(1,end) = (T_T+T_R)/2;

T(end,1) = (T_L+T_B)/2;

T(end,end) = (T_B+T_R)/2;

k=(2*(dx^2+dy^2))/(dx^2*dy^2)

Told=T;

j_iteration_number=1;

%convergence loop

while (error>tolerance)

%nodal loop

for i= 2: (nx-1)

for j= 2:(ny-1)

T(i,j)=(1/k)*((Told(i-1,j)+Told(i+1,j))/dx^2 + (Told(i,j-1)+Told(i,j+1))/(dy^2));

end

end

error=max(max(abs(Told-T)));

%updating old values

Told=T;

j_iteration_number= j_iteration_number+1;

%creating a plot

figure(1)

[C,h]=contourf(x,fliplr(y),T);

clabel(C,h);

colormap(jet)

colorbar;

time_required=toc

title_text = sprintf('solving 2D steady state heat equation using jacobi iterative solver \n j_iteration number=%d ; computation time: %0.3f s', j_iteration_number,time_required),

title(title_text);

xlabel('X axis')

ylabel('Y axis')

end

end

% Function to solve 2D steady state heat conduction equation using gauss seidel iterative solver.

function steady_gauss_seidel_iteration_method(L,nx,ny,error,tolerance)

tic

x=linspace(0,L,nx);

y=linspace(0,L,ny);

dx=L/(nx-1);

dy=L/(ny-1);

T= 300*ones(nx,ny);

T_L = 400;

T_R = 800;

T_T = 600;

T_B = 900;

T(:,1)=T_L;

T(:,end)=T_R;

T(1,:)=T_T;

T(end,:)=T_B;

%edge refining

T(1,1) = (T_L+T_T)/2;

T(1,end) = (T_T+T_R)/2;

T(end,1) = (T_L+T_B)/2;

T(end,end) = (T_B+T_R)/2;

k=(2*(dx^2+dy^2))/(dx^2*dy^2);

Told=T;

gs_iteration_number=1;

while (error>tolerance)

for i= 2: (nx-1)

for j= 2:(ny-1)

T(i,j)=(1/k)*((T(i-1,j)+Told(i+1,j))/dx^2 + (T(i,j-1)+Told(i,j+1))/(dy^2));

end

end

error=max(max(abs(Told-T)));

Told=T;

gs_iteration_number= gs_iteration_number+1;

figure(1)

[C,h]=contourf(x,fliplr(y),T);

clabel(C,h)

colormap(jet)

colorbar;

time_required=toc;

title_text = sprintf('solving 2D steady state heat equation using gauss seidel iterative solver \n iteration number=%d ; computation time: %0.1f s', gs_iteration_number,time_required)

title(title_text);

xlabel('X axis')

ylabel('Y axis')

end

end

% Function to solve 2D steady state heat conduction equation using SOR iterative solver.

function steady_gauss_SOR_method(L,nx,ny,error,tolerance)

tic

x=linspace(0,L,nx);

y=linspace(0,L,ny);

dx=L/(nx-1);

dy=L/(ny-1);

T= 300*ones(nx,ny);

T_L = 400;

T_R = 800;

T_T = 600;

T_B = 900;

T(:,1)=T_L;

T(:,end)=T_R;

T(1,:)=T_T;

T(end,:)=T_B;

%edge refining

T(1,1) = (T_L+T_T)/2;

T(1,end) = (T_T+T_R)/2;

T(end,1) = (T_L+T_B)/2;

T(end,end) = (T_B+T_R)/2;

k=(2*(dx^2+dy^2))/(dx^2*dy^2);

Told=T;

omega=1.72;

sor_iteration_number=1;

while (error>tolerance)

for i= 2: (nx-1)

for j= 2:(ny-1)

Tgs(i,j)=(1/k)*((T(i-1,j)+Told(i+1,j))/dx^2 + (T(i,j-1)+Told(i,j+1))/(dy^2));

T(i,j)=(1-omega)*(Told(i,j))+omega*(Tgs(i,j));

end

end

error=max(max(abs(Told-T)));

Told=T;

sor_iteration_number= sor_iteration_number+1;

figure(1)

[C,h]=contourf(x,fliplr(y),T);

clabel(C,h)

colormap(jet)

colorbar;

time_required=toc;

title_text = sprintf('solving 2D steady state heat equation using SOR iterative solver \niteration number=%d ; computation time: %0.1f s', sor_iteration_number,time_required);

title(title_text);

xlabel('X axis')

ylabel('Y axis')

end

end

In the above three plot results you can see that the number of iterations and the computation time for jacobi iteration method is high i.e, 826 iteration and 320.655 seconds time and then it is foe gauss seidel method with 442 iteration and 158.1 seconds time and then for SOR. On addition of over relaxation factor to the gauss seidel method. Ex SOR method you can see that the solution is calculated in just 59 iteration and 18.3 second’s time. The SOR method depends on relaxation factor and the solution is optimum with a smaller number of iterations and computational time when the omega value is 1.72 and the value above or below this the number iterations gets increased and when the value is above 1.99 the solution is getting blow up. So, the maximum value of omega must be less than 2. You can see the blown-up result if it is given greater than 2 in the below image.

Solving 2D heat conduction equation of transient state:

The equation for 2D heat conduction equation of transient state is –

Where alpha is the thermal diffusivity.

Then discretize the equation by using forward difference scheme for time and central difference for space.

Here we have a classification.

01. If we consider RHS terms of nth time step while calculating for (n+1) time step then we call it as explicit form.

Program:

% Solving the transient state 2D heat conduction problem using explicit method over a square domain.

clear all

close all

clc

tic

%intial Conditions

L=1;

nx=20;

ny=nx;

% TEMPERATURE matrix

x=linspace(0,L,nx);

y=linspace(0,L,ny);

dx=L/(nx-1);

dy=L/(ny-1);

% TEMPERATURE matrix

T= 300*ones(nx,ny);

T_L = 400;

T_R = 800;

T_T = 600;

T_B = 900;

T(:,1)=T_L;

T(:,end)=T_R;

T(1,:)=T_T;

T(end,:)=T_B;

%edge refining

T(1,1) = (T_L+T_T)/2;

T(1,end) = (T_T+T_R)/2;

T(end,1) = (T_L+T_B)/2;

T(end,end) = (T_B+T_R)/2;

Tprev=T;

dt=1e-4;

t=0.15;

nt=t/dt;

alpha=1;

k1=alpha*(dt/dx^2);

k2=alpha*(dt/dy^2);

for k=1:nt

for i=2:nx-1

for j=2:ny-1

T(i,j)=Tprev(i,j)*(1-2*k1-2*k2)+k1*(Tprev(i-1,j)+Tprev(i+1,j))+k2*(Tprev(i,j-1)+Tprev(i,j+1));

end

end

Tprev=T;

figure(1)

[C,h]=contourf(x,fliplr(y),T);

clabel(C,h)

colormap(jet)

colorbar;

time_of_simulation=toc;

title_text = sprintf('solving 2D transient state heat equation explicitly \n time of simulation=%0.3f',time_of_simulation);

title(title_text)

xlabel('X axis')

ylabel('Y axis')

end

Here we get the solution when it satisfies the CFL number condition i.e.

Since in this case dx = dy the condition modifies to

If the condition is not satisfied, then you will see that the solution is blown up.

When we consider RHS terms of (n + 1)th time step while calculating for (n + 1) time step then we call it as implicit form.

Main Program:

% Solving the transient state 2D heat conduction problem using iterative solvers over a square domain.

clear all

close all

clc

tic

%intial Conditions

L=1;

nx=20;

ny=nx;

tolerance=1e-4;

solver_number= input('enter_the_solver_number: ')

if solver_number==1

%function for transient_state_jacobian_iteration_method

transient_jacobian_iteration_method(L,nx,ny,tolerance)

end

if solver_number==2

% function for transient_state_gauss_seidel_iteration_method

transient_gauss_seidel_iteration_method(L,nx,ny,tolerance)

end

if solver_number==3

% function for transient SOR_iteration_method applied on gauss seidel method

transient_gauss_SOR_iteration_method(L,nx,ny,tolerance)

end

Function Program :

% Function to solve 2D transient state heat conduction equation using jacobi iterative solver

function transient_jacobian_iteration_method(L,nx,ny,tolerance)

tic

x=linspace(0,L,nx);

y=linspace(0,L,ny);

dx=L/(nx-1);

dy=L/(ny-1);

% TEMPERATURE matrix

T= 300*ones(nx,ny);

T_L = 400;

T_R = 800;

T_T = 600;

T_B = 900;

T(:,1)=T_L;

T(:,end)=T_R;

T(1,:)=T_T;

T(end,:)=T_B;

%edge refining

T(1,1) = (T_L+T_T)/2;

T(1,end) = (T_T+T_R)/2;

T(end,1) = (T_L+T_B)/2;

T(end,end) = (T_B+T_R)/2;

Tprev=T;

Told =T;

dt=1e-4;

t=0.15;

nt=t/dt;

alpha=1;

k1=alpha*(dt/dx^2);

k2=alpha*(dt/dy^2);

j_iteration_number=1

%time loop

for k=1:nt

error =999

% convergence loop

while(error>tolerance)

% nodal loop

for i=2:nx-1

for j=2:ny-1

T(i,j)=(Tprev(i,j)+k1*(Told(i-1,j)+Told(i+1,j))+k2*(Told(i,j-1)+Told(i,j+1)))/(1+2*k1+2*k2);

end

end

error=max(max(abs(T-Told)));

%updating old values

Told=T;

j_iteration_number= j_iteration_number+1

end

%updating previous values

Tprev=T;

figure(1)

[C,h]=contourf(x,fliplr(y),T);

clabel(C,h);

colormap(jet)

colorbar;

time_required=toc

title_text = sprintf('solving 2D transient state heat equation using jacobi iterative solver \n j_iteration number=%d ; computation time: %0.3f s', j_iteration_number,time_required),

title(title_text);

xlabel('X axis')

ylabel('Y axis')

end

end

% Function to solve 2D transient state heat conduction equation using gauss seidel iterative solver

function transient_gauss_seidel_iteration_method(L,nx,ny,tolerance)

tic

x=linspace(0,L,nx);

y=linspace(0,L,ny);

dx=L/(nx-1);

dy=L/(ny-1);

% TEMPERATURE matrix

T= 300*ones(nx,ny);

T_L = 400;

T_R = 800;

T_T = 600;

T_B = 900;

T(:,1)=T_L;

T(:,end)=T_R;

T(1,:)=T_T;

T(end,:)=T_B;

%edge refining

T(1,1) = (T_L+T_T)/2;

T(1,end) = (T_T+T_R)/2;

T(end,1) = (T_L+T_B)/2;

T(end,end) = (T_B+T_R)/2;

Tprev=T;

Told =T;

dt=1e-4;

t=0.25;

nt=t/dt;

alpha=1;

k1=alpha*(dt/dx^2);

k2=alpha*(dt/dy^2);

gs_iteration_number=1

% timeloop

for k=1:nt

error =999

% convergence loop

while(error>tolerance)

% nodal loop

for i=2:nx-1

for j=2:ny-1

T(i,j)=(Tprev(i,j)+k1*(T(i-1,j)+Told(i+1,j))+k2*(T(i,j-1)+Told(i,j+1)))/(1+2*k1+2*k2);

end

end

error=max(max(abs(T-Told)));

% upgating old values

Told=T;

gs_iteration_number= gs_iteration_number+1;

end

% updating previous values

Tprev=T;

figure(1)

% creating the plot

[C,h]=contourf(x,fliplr(y),T);

clabel(C,h);

colormap(jet)

colorbar;

time_required=toc

title_text = sprintf('solving 2D transient state heat equation using gauss seidel solver \n gs_iteration number=%d ; computation time: %0.3f s', gs_iteration_number,time_required),

title(title_text);

xlabel('X axis')

ylabel('Y axis')

end

end

% Function to solve 2D transient state heat conduction equation using SOR solver

function transient_gauss_SOR_iteration_method(L,nx,ny,tolerance)

tic

x=linspace(0,L,nx);

y=linspace(0,L,ny);

dx=L/(nx-1);

dy=L/(ny-1);

% TEMPERATURE matrix

T= 300*ones(nx,ny);

T_L = 400;

T_R = 800;

T_T = 600;

T_B = 900;

T(:,1)=T_L;

T(:,end)=T_R;

T(1,:)=T_T;

T(end,:)=T_B;

%edge refining

T(1,1) = (T_L+T_T)/2;

T(1,end) = (T_T+T_R)/2;

T(end,1) = (T_L+T_B)/2;

T(end,end) = (T_B+T_R)/2;

Tprev=T;

Told =T;

dt=1e-4;

t=0.15;

nt=t/dt;

alpha=1;

% Tgs=ones(nx,ny)

k1=alpha*(dt/dx^2);

k2=alpha*(dt/dy^2);

SOR_iteration_number=1

omega=1.06

for k=1:nt

error =999;

while(error>tolerance)

for i=2:nx-1

for j=2:ny-1

Tgs(i,j)=(Tprev(i,j)+k1*(T(i-1,j)+Told(i+1,j))+k2*(T(i,j-1)+Told(i,j+1)))/(1+2*k1+2*k2);

T(i,j)=(1-omega)*Told(i,j)+omega*Tgs(i,j);

end

end

error=max(max(abs(T-Told)));

Told=T;

SOR_iteration_number= SOR_iteration_number+1;

end

Tprev=T;

figure(1)

[C,h]=contourf(x,fliplr(y),T);

clabel(C,h);

colormap(jet)

colorbar;

time_required=toc;

title_text = sprintf('solving 2D transient state heat equation using SOR iterative solver \n SOR_iteration number=%d ; computation time: %0.3f s', SOR_iteration_number,time_required);

title(title_text);

xlabel('X axis')

ylabel('Y axis')

end

end

From the above plots you can conclude that explicit method takes less time to execute but when it comes to stability it mainly depends on CFL condition if this condition is not satisfied the solution gets blow up. Where as in the implicit methods Jacobian method gives the solution in 7970 iterations whereas gauss method takes 10115 iterations but using SOR method you get the solution in 5098 iterations but the computational time high for gauss and then Jacobian but for the SOR method is faster than other two. Here in SOR method if you use omega value more or less than 1.06 you can see the iteration count increases when it crosses 1.9 then you see that the solution got blew up.so its value should be given less than 2.you can also observe that for any values of alpha dx and dt the solution will be stable, like explicit it is not dependent of CFL condition.