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

Solution

Equation of two-dimensional heat conduction equation is

`frac{delT}{delt}=alpha**(frac{del^2T}{delx^2}+frac{del^2T}{delx^2})`

where `alpha` is thermal diffusion coefficient.

The equation is best described as a transient state two-dimensional heat conduction equation whereas if the time rate of change of temperature is zero it will be called a steady-state equation.

Applying numerical discretization for the above-mentioned equation

Transient state two-dimensional heat conduction equation 

Nomenclature

Nodes representing X-axis = "i".

Nodes representing Y axis="j".

At any (i,j) temperature is represented by `T(i,j)`

The temperature at the top of (i,j)=`T(i,j+1)`

The temperature at the bottom of (i,j)=`T(i,j-1)`

The temperature at the left of (i,j)=`T(i-1,j)`

The temperature at the right of (i,j)=`T(i+1,j)`

(Please note that the central difference scheme is applicable for the temperature gradient (space loop) as boundary conditions are known to us, temperature for the inner nodes to be calculated. Forward difference scheme is applied for the time loop)

The equation for two-dimensional transient state heat conduction is described below

`frac{T(i)^(n+1)-T(i)^n}{dt}=alpha**(frac{(Ti,j+1)+(Ti,j-1)-2Ti,j}{dy^2}+frac{(Ti+1,j)+(Ti-1,j)-2Ti,j}{dx^2})`

`T(i)^(n+1)=T(i)^n+frac{alpha**dt**((Ti,j+1)+(Ti,j-1)-2Ti,j)){dy^2}+frac{alpha**dt**((Ti+1,j)+(Ti-1,j)-2Ti,j)}{dx^2}`

Explicit equation can be calculated for n time steps.

`T(i,j)^(n+1)=T(i,j)^n+frac{alpha**dt**((T(i,j+1)^n)+(T(i,j-1)^n)-2T(i,j)^n)){dy^2}+frac{alpha**dt**((T(i+1,j)^n)+(T(i-1,j)^n)-2T(i,j)^n)}{dx^2}`

This is one of the easiest way to calculate the temperature for (n+1) time steps and at (i,j) nodes since right hand side equation terms are already known for previous time steps i,e "n".

The implicit equation can be calculated for n+1 time steps.

`T(i,j)^(n+1)=T(i,j)^n+frac{alpha**dt**((T(i,j+1)^(n+1))+(T(i,j-1)^(n+1))-2T(i,j)^(n+1))){dy^2}+frac{alpha**dt**((T(i+1,j)^(n+1))+(T(i-1,j)^(n+1))-2T(i,j)^(n+1))}{dx^2}`

Calculating K1 and K2

`K1=frac{alpha**dt}{dx^2}`

`K2=frac{alpha*dt}{dy^2}`

`T(i,j)^(n+1)**(1+(2**K1)+(2**K2))=T(i,j)^n+K2**((T(i,j+1)^(n+1))+(T(i,j-1)^(n+1))+K1**((T(i+1,j)^(n+1))+(T(i-1,j)^(n+1))`

`T(i,j)^(n+1)=frac{T(i,j)^n+K2**((T(i,j+1)^(n+1))+(T(i,j-1)^(n+1)))+K1**((T(i+1,j)^(n+1))+(T(i-1,j)^(n+1)))}{1+(2**K1)+(2**K2)}`

The implicit equation consist of lot of unknowns apart from LHS because RHS consists of all the terms in the present time step.

Steady state two-dimensional heat conduction equation 

For steady state `frac{delT}{delt}``=0`

The equation for two-dimensional steady state heat conduction is described below

`frac{del^2T}{delx^2}+frac{del^2T}{dely^2}=0`

Applying numerical discretisation for the above equation

`frac{Ti+1,j+Ti-1,j}{dx^2}+frac{Ti,j-1+Ti,j+1}{dy^2}=Ti,j**frac{2**(dx^2+dy^2)}{dx^2**dy^2}`

`Ti,j=frac{1}{K}**(frac{Ti+1,j+Ti-1,j}{dx^2}+frac{Ti,j-1+Ti,j+1}{dy^2})`

where

`K=frac{2**(dx^2+dy^2)}{dx^2**dy^2}`

The implicit equation can be solved with the help of iterative solvers which are listed below

1. Jacobi method : As the iteration marches forward jacobi will use the value of temperature of previous node every time till the point of convergence. Eg if the temperature at node 3 is to be calculated jacobi will use `T2`and `T4`old values to calculate node at `T3`
2. Gauss-Seidel method: For calculating the value of temperature at any node it will use the new value of previous node which is updated after iteration and the old value of next node since gauss seidel has the value of that node from previous iteration.
3. Successive over relaxation method : It is a variant of method number 2. It uses a relaxation factor `omega`. When `omega``gt1` it is over relaxation. When `omega``lt1` it is under relaxation.

Transient state two-dimensional heat conduction  implicit equation 

Jacobi method

• `T(i,j)^(n+1)=frac{Told(i,j)^n+K2**((Told(i,j+1)^(n+1))+(Told(i,j-1)^(n+1)))+K1**((Told(i+1,j)^(n+1))+(Told(i-1,j)^(n+1)))}{1+(2**K1)+(2**K2)}`

Gauss-seidel method

• `T(i,j)^(n+1)=frac{Told(i,j)^n+K2**((Told(i,j+1)^(n+1))+(T(i,j-1)^(n+1)))+K1**((Told(i+1,j)^(n+1))+(T(i-1,j)^(n+1)))}{1+(2**K1)+(2**K2)}`

Successive over relaxation method

• `T(i,j)^(n+1)=omega**frac{Told(i,j)^n+K2**((Told(i,j+1)^(n+1))+(T(i,j-1)^(n+1)))+K1**((Told(i+1,j)^(n+1))+(T(i-1,j)^(n+1)))}{1+(2**K1)+(2**K2)}+(1-omega)**Told(i,j)`

 Transient state two-dimensional heat conduction explicit equation 

•  `T(i,j)^(n+1)=Told(i,j)^n+frac{alpha**dt**((Told(i,j+1)^n)+(Told(i,j-1)^n)-2Told(i,j)^n)){dy^2}+frac{alpha**dt**((Told(i+1,j)^n)+(Told(i-1,j)^n)-2Told(i,j)^n)}{dx^2}`

Steady state two-dimensional heat conduction implicit equation 

Jacobi method

• `Ti,j=frac{1}{K}**(frac{T(old)i+1,j+T(old)i-1,j}{dx^2}+frac{T(old)i,j-1+T(old)i,j+1}{dy^2})`

Gauss-seidel method

• `Ti,j=frac{1}{K}**(frac{T(old)i+1,j+Ti-1,j}{dx^2}+frac{Ti,j-1+T(old)i,j+1}{dy^2})`

Successive over relaxation method

• `Ti,j=omega**frac{1}{K}**(frac{T(old)i+1,j+Ti-1,j}{dx^2}+frac{Ti,j-1+T(old)i,j+1}{dy^2})+(1-omega)**Told(i,j)`

Courant-Fredreichs-Lewy condition

It is a necessary condition for convergence for solving partial difference equation.It is an important criteria determining stability of the profile.

In this case

`CFLx=frac{alpha*dt}{dx^2}`

`CFLy=frac{alpha*dt}{dy^2}`

For temperature profile to be stable

`CFLle0.5`

`CFLx+CFLyle0.5`

Gathering other input parameters

• Top Boundary = 600 K
• Bottom Boundary = 900 K
• Left Boundary = 400 K
• Right Boundary = 800 K
• Absolute error criteria = 1e-4

Matlab coding algorithm

• Defining input parameters
• Error is calculated which is `max(max(ab s(Told-T)))` after nodal loop.
• Told=T is created for steady and transient state two-dimensional heat conduction equation.
• Told=T is updated after the nodal loop.
• Tprevious=T is created and updated after the convergence loop.
• The graph is plotted after the convergence loop to get the best possible result.
• To get the computation time tic and toc should be used at the starting and ending of for loop.

Matlab coding for two-dimensional steady-state equation solved implicitly.

% solving 2D steady state temperature heat conduction equation 
% iterative solvers
% implicit equation
% Jacobi_method

clear 
close all
clc

% length of domain
L=1;
%Number of node points
n=10;
% defining x value
x=linspace(0,L,n);
dx=x(2)-x(1);
y=linspace(0,L,n);
dy=y(2)-y(1);
% other input parameters
T=100*ones(n);
T(1,1:n)=600;
T(n,1:n)=900;
T(1:n,n)=800;
T(1:n,1)=400;
T(1,1)=500;
T(1,n)=700;
T(n,n)=750;
T(n,1)=650;
Told=T;
% iterative solver to be used (jacobian_iter=1,GS=2,SOR=3)
iterative_solver=1;
jacobi_iter=0;
% error and tolerance
error=9e9;
tol=1e-4;

% applying convergence loop
tic
while (error>tol)
    for i=2:n-1
        for j=2:n-1
    if iterative_solver==1
        figure(1)
        % jacobian
        T(i,j)=0.25*((Told(i+1,j)+Told(i-1,j)+Told(i,j+1)+Told(i,j-1)));
    end
        end
    end
    jacobi_iter_time=toc;
    error=max(max(abs(Told-T)));
    Told=T;
    [c,h]=contourf(x,y,T);
    colorbar
    colormap(jet)
    xlabel('xaxis')
    ylabel('yaxis')
    clabel(c,h)
    title_text=sprintf('2D steady state temperature gradient jacobian iteration=%d,computation time=%f',jacobi_iter,jacobi_iter_time);
    title(title_text)
    jacobi_iter=jacobi_iter+1;
end

% solving 2D steady state temperature heat conduction equation 
% iterative solvers
% implicit equation
% Gauss_seidel_method

clear 
close all
clc

% length of domain
L=1;
%Number of node points
n=10;
% defining x value
x=linspace(0,L,n);
dx=x(2)-x(1);
y=linspace(0,L,n);
dy=y(2)-y(1);
% other input parameters
T=100*ones(n);
T(1,1:n)=600;
T(n,1:n)=900;
T(1:n,n)=800;
T(1:n,1)=400;
T(1,1)=500;
T(1,n)=700;
T(n,n)=750;
T(n,1)=650;
Told=T;
% iterative solver to be used 
iterative_solver=1;
GS=0;
% error and tolerance
error=9e9;
tol=1e-4;

% applying convergence loop


  tic;
   
while (error>tol)
   
   for i=2:n-1
        for j=2:n-1
    if iterative_solver==1
        figure(1)
        % GS
        T(i,j)=0.25*((Told(i+1,j)+T(i-1,j)+Told(i,j+1)+T(i,j-1)));
    end
        end
    
    end
   error=max(max(abs(Told-T)));
    Told=T;
    GS=GS+1;
    GS_time=toc;
   [c,h]=contourf(x,y,T);
    colorbar
    colormap(jet)
    xlabel('xaxis')
    ylabel('yaxis')
    clabel(c,h)
    title_text=sprintf('2D steady state temperature gradient gauss seidel iteration=%d,computation time=%f',GS,GS_time);
    title(title_text)
  
end
% solving 2D steady state temperature heat conduction equation 
% iterative solvers
% implicit equation
% Successive_over_relaxation

clear 
close all
clc

% length of domain
L=1;
%Number of node points
n=10;
% defining x value
x=linspace(0,L,n);
dx=x(2)-x(1);
y=linspace(0,L,n);
dy=y(2)-y(1);
% other input parameters
T=100*ones(n);
T(1,1:n)=600;
T(n,1:n)=900;
T(1:n,n)=800;
T(1:n,1)=400;
T(1,1)=500;
T(1,n)=700;
T(n,n)=750;
T(n,1)=650;
Told=T;
% iterative solver to be used 
iterative_solver=1;
SOR=0;
% error and tolerance
error=9e9;
tol=1e-4;
omega=1.8;
tic;
% applying convergence loop
while (error>tol)
    for i=2:n-1
        for j=2:n-1
    if iterative_solver==1
      
        % SOR
        T(i,j)=omega*(0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1)))+(1-omega)*Told(i,j);
    end
        end
    end
    SOR_time=toc;
    error=max(max(abs(Told-T)));
    Told=T;
    figure(2)
    [c,h]=contourf(x,y,T);
    colorbar
    colormap(jet)
    xlabel('xaxis')
    ylabel('yaxis')
    clabel(c,h)
    title_text=sprintf('2D steady state temperature gradient SOR iteration=%d,computation time=%f',SOR,SOR_time);
    title(title_text)
    SOR=SOR+1;
end

• Matlab coding for two-dimensional transient-state equation solved explicitly.

% solving 2D transient state temperature heat conduction equation 
% iterative solvers
% explicit equation

clear 
close all
clc

% length of domain
L=1;
%Number of node points
n=10;
% defining x value
x=linspace(0,L,n);
dx=x(2)-x(1);
y=linspace(0,L,n);
dy=y(2)-y(1);
% time marching
t=1000;
dt=10;
nt=t/dt;

% thermal diffusivity
alpha=1.1;
K1=(alpha*dt*(dx^-2));
K2=(alpha*dt*(dy^-2));
% other input parameters
T=300*ones(n);
T(1,1:n)=600;
T(n,1:n)=900;
T(1:n,n)=800;
T(1:n,1)=400;
T(1,1)=500;
T(1,n)=700;
T(n,1)=650;
T(n,n)=850;
Told=T;
Tprevious=T;
term1=(1+(2*K1)+(2*K2))^-1;
term2=K1*term1;
term3=K2*term1;
explicit_time=1;
tic;
for k=1:nt
    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)=(Tprevious(i,j)*term1)+(H*term2)+(V*term3);
        end
    end
    error=max(max(abs(Told-T)));
    Told=T;
    explicit_time=explicit_time+1;
    ET=toc;
    [c,h]=contourf(x,y,T);
    colorbar
    colormap(jet)
    clabel(c,h)
    xlabel('xaxis')
    ylabel('yaxis')
    title_text=sprintf('2D transient state conduction explicit iteration=%d,computation time=%d',explicit_time,ET);
    title(title_text)
    pause(0.03)
end

• Matlab coding for two-dimensional transient-state equation solved implicitly.

% solving 2D transient state temperature heat conduction equation 
% iterative solvers
% imlicit equation
% jacobi_method

clear 
close all
clc

% length of domain
L=1;
%Number of node points
n=10;
% defining x value
x=linspace(0,L,n);
dx=x(2)-x(1);
y=linspace(0,L,n);
dy=y(2)-y(1);
% time marching
t=1000;
dt=10;
nt=t/dt;

% thermal diffusivity
alpha=1.1;
% error and tolerance
error=9e9;
tol=1e-4;
% other input parameters
T=300*ones(n);
T(1,1:n)=600;
T(n,1:n)=900;
T(1:n,n)=800;
T(1:n,1)=400;
T(1,1)=500;
T(1,n)=700;
T(n,n)=750;
T(n,1)=650;
Told=T;
Tprevious=T;
K1=(alpha*dt)/dx^2;
K2=(alpha*dt)/dy^2;
term1=(1+(2*K1)+(2*K2))^-1;
term2=K1*term1;
term3=K2*term1;
%  iteration
iterative_method=1;


if iterative_method==1
    jacobi_iter=1;
    tic
    for k=1:nt
        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)=(Tprevious(i,j)*term1)+(H*term2)+(V*term3); 
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            jacobi_iter=jacobi_iter+1;
        end
        Tprevious=T;
        jacobian_time=toc;
        [c,h]=contourf(x,y,T);
        colorbar
        colormap(jet)
        xlabel('xaxis')
        ylabel('yaxis')
        clabel(c,h)
        title_text=sprintf('2D temperature transient state jacobian iteration=%d,computational time=%f',jacobi_iter,jacobian_time);
        title(title_text)
        pause(0.03)
        
    end
   
end

% solving 2D transient state temperature heat conduction equation 
% iterative solvers
% imlicit equation
% Gauss_Seidel_method

clear 
close all
clc

% length of domain
L=1;
%Number of node points
n=10;
% defining x value
x=linspace(0,L,n);
dx=x(2)-x(1);
y=linspace(0,L,n);
dy=y(2)-y(1);
% time marching
t=1000;
dt=10;
nt=t/dt;

% thermal diffusivity
alpha=1.1;
% error and tolerance
error=9e9;
tol=1e-4;
% other input parameters
T=1000*ones(n);
T(1,1:n)=600;
T(n,1:n)=900;
T(1:n,n)=800;
T(1:n,1)=400;
T(1,1)=500;
T(1,n)=700;
T(n,n)=750;
T(n,1)=650;
Told=T;
Tprevious=T;
% iteration
iterative_method=1;
K1=(alpha*dt)/dx^2;
K2=(alpha*dt)/dy^2;
term1=(1+(2*K1)+(2*K2))^-1;
term2=K1*term1;
term3=K2*term1;

% applying time marching loop
if iterative_method==1
    GS=1;
  
    tic;
    for k=1:nt
        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)+Told(i+1,j));
                  V=(T(i,j-1)+Told(i,j+1));
                  T(i,j)=(Tprevious(i,j)*term1)+(H*term2)+(V*term3);    
                end
            end
            error=max(max(abs(Told-T)));
            Told=T;
            GS=GS+1;
        end
        Tprevious=T;
        GS_time=toc;
        [c,h]=contourf(x,y,T);
        colorbar
        colormap(jet)
        disp('Temperature')
        xlabel('x axis')
        ylabel('y axis')
        clabel(c,h)
        title_text=sprintf('GS iteration for 2D transient heat conduction=%d,computational time=%f',GS,GS_time);
        title(title_text)
        pause(0.03)
    end
end

% solving 2D transient state temperature heat conduction equation 
% iterative solvers
% imlicit equation
% SOR

clear 
close all
clc

% length of domain
L=1;
%Number of node points
n=10;
% defining x value
x=linspace(0,L,n);
dx=x(2)-x(1);
y=linspace(0,L,n);
dy=y(2)-y(1);
% time marching
t=1000;
dt=10;
nt=t/dt;
omega=1.8;

% thermal diffusivity
alpha=4e-3;
% error and tolerance
error=9e9;
tol=1e-4;
% other input parameters
T=100*ones(n);
T(1,1:n)=600;
T(n,1:n)=900;
T(1:n,n)=800;
T(1:n,1)=400;
T(1,1)=500;
T(1,n)=700;
T(n,n)=750;
T(n,1)=650;
Told=T;
Tprevious=T;
Tend=T;
% iteration
iterative_method=1;
SOR=1;
K1=(alpha*dt)/dx^2;
K2=(alpha*dt)/dy^2;
term1=(1+(2*K1)+(2*K2))^-1;
term2=K1*term1;
term3=K2*term1;
% applying time marching loop
if iterative_method==1
    SOR=1;
    tic;
    for k=1:nt
    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)+Told(i+1,j));
              V=(T(i,j-1)+Told(i,j+1));
              T(i,j)=(Tprevious(i,j)*term1)+(H*term2)+(V*term3);  
              Tend(i,j)=(1-omega)*Told(i,j)+(T(i,j)*omega);
            end
        end
        error=max(max(abs(Told-T)));
        Told=T;
        SOR=SOR+1;
    end
    SOR_time=toc;
    Tprevious=T;
    Tend=T;
    [c,h]=contourf(x,y,T);
    colorbar
    colormap(jet)
    xlabel('x axis')
    ylabel('y axis')
    clabel(c,h)
    title_text=sprintf('SOR iteration for 2D transient heat conduction=%d,computational time=%f',SOR,SOR_time);
    title(title_text)
    pause(0.03)
    end
end

• Two-dimensional steady-state heat conduction graph (Jacobi method) graph.

• Two-dimensional steady-state heat conduction graph (Gauss-Seidel method) graph.

• Two-dimensional steady-state heat conduction graph (SOR method) graph.

• Two-dimensional transient-state heat conduction graph (explicit method).

• Two-dimensional transient-state heat conduction graph (Jacobi method).

• Two-dimensional transient-state heat conduction graph (Gauss-Seidel method).

• Two-dimensional transient-state heat conduction graph (SOR method).

Conclusion

The following information can be inferred from the plot

It can be concluded that the number of iteration for SOR is less compared to the other two (the highest being the Jacobian method).

`Jacobian-method gt Gauss-Seid el-method gtSOR`

• The explicit approach is by far the easiest way to find the temperature at the inner nodes. Hence for a given value of dx dt must be very less as compared to some stability constraints. Hence CFL along x and y should be less than 0.5 to get a stable profile.
• The implicit approach contains information obtained from solving the simultaneous equation for the full grid for each time step. This is computationally more demanding but allows for larger time steps and better stability.