CFD || Solving Steady State 2D Heat Conduction Equation

OBJECTIVE

To solve a Steady state 2D heat conduction equation using different Iterative techniques (Jacobi,Gauss Seidal,SOR)

ASSUMPTIONS

• Assume that the domain is a unit square.
• Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]
• Boundary conditions for steady case - Left:400K, Top:600K, Right:800K, Bottom:900K

OVERVIEW

We know 2D Steady State Heat Conduction Equation as;

`(partial^2T)/(partialx^2)+(partial^2T)/(partialy^2)=0`

Let,

• `T_p = T(i,j);`
• `T_L = T(i-1,j);`
• `T_R = T(i+1,j);`
• `T_T = T(i,j+1);`
• `T_B = T(i,j-1);`
• `k = [(2(Deltax^2+Deltay^2))/(Deltax^2+Deltay^2)]`

Discretizing the above equation using central difference;

`(T_R-2T_P+T_L)/(Deltax^2)+(T_T-2T_P+T_B)/(Deltay^2)=0`

`implies T_P = ([(T_R+T_L)/(Deltax^2)+(T_T+T_B)/(Deltay^2)])/k` ...(eqn)`

To the above eqn different iterative techniques are apllied

• Jacobi

`T_P = ([(T_(oldR)+T_(oldL))/(Deltax^2)+(T_(oldT)+T_(oldB))/(Deltay^2)])/k`

• Gauss-siedal

`T_P = ([(T_(oldR)+T_(L))/(Deltax^2)+(T_(oldT)+T_(B))/(Deltay^2)])/k`

• SOR

`T_P = omega*([(T_(oldR)+T_(L))/(Deltax^2)+(T_(oldT)+T_(B))/(Deltay^2)])/k+(1-w)T_(oldP)`

MATLAB FUNCTION

%Solving for steady state
%dt^2/dx^2 + dt^2/dy^2=0

clear all
close all
clc

nx=10;
ny=10;
T=ones(nx,ny);

%BC
T(1,:)=900;    %Temp bottom
T(:,nx)=800;  %Temp right
T(ny,:)=600;   %Temp top
T(:,1)=400;   %Temp Left

%Intialisation
Told=T;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);

tol=1e-4;
error=10;
k=(2*(dx^2+dy^2)/(dx^2*dy^2));

iter_solve=3; %%%%(Choose the solver)

%__________________________________________________________________________

%_________Jacobi solver______________%
if iter_solve==1
    tic
    j_iter=1;
    
    while(error>tol)
        for i=2:nx-1
            for j=2:ny-1
               T(i,j)=(((Told(i+1,j)+Told(i-1,j))/(dx^2))+((Told(i,j+1)+Told(i,j-1))/(dy^2)))/k;
            end
        end
        error=max(max(abs(T-Told)));
        Told=T;
        j_iter=j_iter+1;
        
        %plotting
        [x1,y1]=contourf(x,y,T);
        colormap(jet);
        clabel(x1,y1);
        xlabel('xaxis')
        ylabel('yaxis')
        text_title=sprintf('Temperature-Jacobi Iteration=%d',j_iter); 
        title(text_title)
        pause(0.1)
    end
    time_jacobi = toc;
end

%__________________________________________________________________________


    %_________GS solver__________%
if iter_solve==2
    tic
    gs_iter=1;
    
    while(error>tol)
        for i=2:nx-1
            for j=2:ny-1
             T(i,j)=(((Told(i+1,j)+T(i-1,j))/(dx^2))+((Told(i,j+1)+T(i,j-1))/(dy^2)))/k;   
            end
        end
        error=max(max(abs(T-Told)));
        Told=T;
        gs_iter=gs_iter+1;
        
        %plotting
        
        [x1,y1]=contourf(x,y,T);
        colormap(jet);
        clabel(x1,y1);
        xlabel('xaxis')
        ylabel('yaxis')
        text_title=sprintf('Temperature-GS Iteration=%d',gs_iter); 
        title(text_title)
        pause(0.01)
    end
    time_gs = toc;
end

%__________________________________________________________________________


        %____________SOR______________%       
        if iter_solve==3
            tic
            sor_iter=1;
            w=1.4;
            while(error>tol)
                for i=2:nx-1
                    for j=2:ny-1
                         T(i,j)=(w*(((Told(i+1,j)+T(i-1,j))/(dx^2))+((Told(i,j+1)+T(i,j-1))/(dy^2)))/k)+((1-w)*Told(i,j));
                    end
                end
                error=max(max(abs(T-Told)));
                Told=T;
                sor_iter=sor_iter+1;
                
                %plotting
               
                [x1,y1]=contourf(x,y,T);
                colormap(jet);
                clabel(x1,y1);
                xlabel('xaxis')
                ylabel('yaxis')
                text_title=sprintf('Temperature-SOR Iteration=%d',sor_iter); 
                title(text_title)
                pause(0.001)
            end
         time_sor = toc;    
        end

RESULTS

• Time taken by Jacobi = 44 sec.
• Time taken by Gauss-siedal = 20 sec.
• Time taken by SOR = 12 sec.

CONCLUSION

Its clear from the above data that the no.of iteration and time taken to reach steady state by different iterative techniques is in the order Jacobi>Gauss-siedal>SOR.