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

Objective: Solve the discretized 2-dimentional Heat conduction equation and obtain the temperature field through time marching code in both implicit and explicit forms in MATLAB

Theory:

The heat conduction equation is a partial differential equation that describes the distribution of heat (or temperature field) in a given body over time.

It describes the flow of heat from areas of high temperature to areas of low temperature in a solid medium. The extent of heat transfer can be visualized through the temperature field. It is a function of the Cartesian (x,y,z) coordinates and the time elapsed (t) .

The boundary conditions define an initial state of temperatures and the heat flows over the solid with the passages of time until a steady state is reached at which point temperature variation due to heat transfer is negligibly low.

Implict , Explicit schemes and Iterative solvers

When equations containing derivatives with respect to time , a so called – ‘time marching’ solution scheme needs to be employed in which the time is increased from a reference zero( initial conditions) in additions of a given time step .

In the discretized form of the equation we solve for the variable in this case the temperature in terms of its own values in the previous time step. This is the basis of the ‘Explicit Solver’

In a different case, the value of the variable is expressed in terms of its values in different parts of the space in the present time step. So here an iterative solver is used to arrive at a convergent result. This is the 'implicit solver'

An iterative solver is used when the given matrix and the variables are so large that conventional techniques like gauss-elimination would be unrealistic. If in a given equation of..

`Ax=B`

An initial guess value of x is used, and the value of B is calculated and a residual is obtained.

A linear transformation is applied to the guess values (x1, x2 ,x3…) of x and the new values are obtained . If ‘T’ is the transformation matrix and ‘C’ a constant

`overset(n+1)X = T*overset(n)X + C`

(n- time generation)

This way a new generation of value of X is obtained and this process continues until a convergence is obtained wherein the difference of the new and old values tends to zero.

`overset(n+1)X - overset(n)X rarr 0`

A convergent solution is obtained.

The two approaches employed to solve the 2D equation are steady state and transient.

1. Steady State Solution schemes (Only Implicit):

• Jacobi
• Gauss Siedel
• Successive Over-relaxation

2. Transient Solution schemes:

• Matrix Inversion (Implicit)
• Explicit method

Jacobi , Gauss Siedal and Successive over-relaxation

These are the three solving techniques used. The update transformation matrix ‘T’ mentioned earlier contained values of X matrix. Therefore when solving for one variable other variables are also involved, their values given by the latest guess value(previous time step) .

In the Jacobi solver all equations use the previous time step values of the variables to calculate the next time step value of one variable.

In Gauss Seidel, the variables whose values in the next iteration is already calculated are used instead of the guess value. This way a more aggressive update is given to the variable values and the number of iterations required for convergence is lower.

In successive over-relaxation, the value being calculated from the R.H.S of the matric transformation is multiplied by a factor ‘w’. For the correct values of ‘w’ the iterations required for convergence is lowered. This is because a more aggressive convergence route is taken.

If ‘w’ is >1 over-relaxed.

If ‘w’ is <1 under-relaxed.

The heat conduction equation

`alpha*(del^2T/(del^2x) +del^2T/(del^2y)) = delT/delt`

in 2D

for the steady state condition,

`delT/(delt)=0`

and through central discretization and using the fact that

`Deltax=Deltay`

we get,

`underset(i,j)T=(1/4)*(underset(i-1,j)T + underset(i+1,j)T + underset(i,j-1)T + underset(i,j+1)T)`

for the transient case

`delT/(delt)!=0`

In this case both explicit and implicit solvers are used ,

Explicit Transient -

`alpha*(((underset(i-1,j)T - 2*underset(i,j)T + underset(i+1,j)T)^n/(Deltax^2)) +((underset(i,j-1)T - 2*underset(i,j)T + underset(i,j+1)T)^n/(Deltay^2)))=((underset(i,j)T^(n+1) - underset(i,j)T^(n))/(Deltat))`

Implicit Transient-

`alpha*(((underset(i-1,j)T - 2*underset(i,j)T + underset(i+1,j)T)^(n+1)/(Deltax^2)) +((underset(i,j-1)T - 2*underset(i,j)T + underset(i,j+1)T)^(n+1)/(Deltay^2)))=((underset(i,j)T^(n+1) - underset(i,j)T^(n))/(Deltat))`

Results plotted in contor graphs

The boundary conditions are 400,800,600,900K .

Steady State

Initial state

Final State

Transient Implicit

Initial state

Final state

Transient Explicit

Initial

Final

Conclustion

A comparision of the number of iterations taken for convergence in both steady state and transient state implicit solvers using Jacobi , Gauss Siedal and Successive over-relaxation(SOR).

Steady State

Jacobi -2667

Gauss Siedal-1396

SOR with 'w'=1.1 -1156

Transient implicit solver

Jacobi -11763

Gauss Seidal -7544

SOR with 'w'=1.7 -3564

Therefore with the right value of relaxation factor the number of iterations required for convergence is in the order of

JACOBI>GAUSS SEIDAL>SUCCESSIVE OVER-RELAXATION

1.Steady state
clc
clear all
close all
%The number of elemets in the square
n=30;
number=1;
t=ones(n,n);
number=1;

x=linspace(0,1,n);
dx=x(2)-x(1);
dy=dx;
k= 2*((dx*dx) + (dy*dy))/(dy*dy*dx*dx);
error=9e9;
method =1;
% 1-Jacobi , 2- GS , 3-SOR
t(2:n-1,1)=400;
t(2:n-1,n)=800;
t(1,2:n-1)=600;
t(n,2:n-1)=900;
t(1,1)=500;
t(1,n)=700;
t(n,n)=850;
t(n,1)=650;
told=t;
tol = (1e-6);
%Jacobi method
if(method ==1)
    while(error>tol)
    for i=2:n-1
        for j=2:n-1
            term1=(told(i,j-1) + told(i,j+1))*(1/k)*(1/(dx*dx));
            term2=(told(i-1,j) + told(i+1,j))*(1/k)*(1/(dx*dx));
            t(i,j)=term1+term2;
        end
    end
    error1= max(abs(told - t))
    error= max(error1)
    told=t;
    number=number+1;
    end
end
%GS
if(method==2)
    while(error>tol)
    for i=2:n-1
        for j=2:n-1
            term1=(t(i,j-1) + t(i,j+1))*(1/k)*(1/(dx*dx));
            term2=(t(i-1,j) + t(i+1,j))*(1/k)*(1/(dx*dx));
            t(i,j)=term1+term2;
        end
    end
    error1= max(abs(told - t))
    error= max(error1)
    told=t;
    number=number+1;
    end
end
%SOR
if(method==3)
    w=1.1;
    while(error>tol)
    for i=2:n-1
        for j=2:n-1
            term1=(t(i,j-1) + t(i,j+1))*(1/k)*(1/(dx*dx));
            term2=(t(i-1,j) + t(i+1,j))*(1/k)*(1/(dx*dx));
            t(i,j)= ((1-w)*told(i,j)) + w*(term1+term2);
        end
    end
    error1= max(abs(told - t))
    error= max(error1)
    told=t;
    number=number+1;
    end
end


2.Explicit solver for transient
close all
clc
clear all
n=30;
num=1;
t=ones(n,n);
t(2:n-1,1)=400;
t(2:n-1,n)=800;
t(1,2:n-1)=600;
t(n,2:n-1)=900;
t(1,1)=500;
t(1,n)=700;
t(n,n)=850;
t(n,1)=650;
told=t;
dt=0.001;
nt=2000;
dx=1/(n-1);
dy=dx;
for r=1:nt
    for i =2:n-1
        for j=2:n-1
            term1 = (told(i-1,j) - (2*told(i,j)) + (told(i+1,j)));
            term2 = ((told(i,j+1) - (2*told(i,j)) + told(i,j-1)));
            t(i,j) = told(i,j) + (1.4*term1*dt) + (1.4*term2*dt);
        end
    end
    told=t;
end


3. Implicit solver for transient case
clc
clear all
close all
n=30;
t=ones(n,n);
t(2:n-1,1)=400;
t(2:n-1,n)=800;
t(1,2:n-1)=600;
t(n,2:n-1)=900;
t(1,1)=500;
t(1,n)=700;
t(n,n)=850;
t(n,1)=650;
num=1;
method=3;
told=t;
tprev=t;
tol = (1e-4);
error1=ones(1,n);
dx=1/(n-1);
dy=1/(n-1);
dt=0.01;

k1=dt*1.4/dx/dy;
k2=1.4*dt/dy/dx;
%Jacobi
if(method==1)
for nt =1:2000
    error=9e9;
    while(error>tol)
        for i=2:n-1
            for j=2:n-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)=(tprev(i,j)*term1)+(H*term2)+(V*term3);
            end
        end
        error1=max(abs(told-t));
        error=max(error1);
        num=num+1;
       
        told=t;
    end
    tprev=t;
end
end
%GS
if(method==2)
    for nt =1:2000
    error=9e9;
    while(error>tol)
        for i=2:n-1
            for j=2:n-1
                term1=(1+(2*k1)+(2*k2))^-1;
                term2=k1*term1;
                term3=k2*term1;
                H=(t(i-1,j)+t(i+1,j));
                V=(t(i,j-1)+t(i,j+1));
                t(i,j)=(tprev(i,j)*term1)+(H*term2)+(V*term3);
            end
        end
        error1=max(abs(told-t));
        error=max(error1);
        num=num+1;
       
        told=t;
    end
    tprev=t;
end
end
%SOR
if(method==3)
   for nt =1:2000
    error=9e9;
    w=1.7;
    while(error>tol)
        for i=2:n-1
            for j=2:n-1
                term1=(1+(2*k1)+(2*k2))^-1;
                term2=k1*term1;
                term3=k2*term1;
                H=(t(i-1,j)+t(i+1,j));
                V=(t(i,j-1)+t(i,j+1));
                t(i,j)=((1-w)*told(i,j)) + (w*((tprev(i,j)*term1)+(H*term2)+(V*term3)));
            end
        end
        error1=max(abs(told-t));
        error=max(error1);
        num=num+1;
       
        told=t;
    end
    tprev=t;
end
end