BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

VVarun Gupta27 Apr 2021

Aim : To solve 2D heat conduction equation using Jacobi, Gauss-Seidal, SOR, Explicit and implicit schemes.

Theory : for solving 2D heat conduction equation, iterative solvers are required.

1. Jacobi - It uses the old values to find the new values.

2. Gauss-Seidal - It uses latest values to find the new values.

3. SOR - It uses a factor in GS solver or Jacobi solver to converge more faster.

these are solvers are used for complex and coupled equations like in implicit scheme.

however for explicit scheme, equations are simple and the value of previous equations are used in the next equation. So such solvers are not required in explicit scheme.

Governing equations

Jacobi : Let `Ax = B`

`where, A = D+L+U`

`x^(k+1) = D^-1(b-(L+U)x^k)`

Gauss- Seidal

Let `Ax = B`

`where, A = D+L+U`

`L1 = D + L`

`x^(k+1) = L1^-1(b-Ux^k)`

SOR:

`T3^(k+1) = T3^k(1-alpha) +alpha(T(GS))`

Code is as follows:

clear all
close all
clc;

tol = 1e-4;
L = 1;
n = 10;
dx = L/(n-1);
dy = L/(n-1);
k = 2*(dx^2+dy^2)/(dx^2*dy^2);
T = zeros(10,10);
T(:,1) = 400;
T(1,:) = 600;
T(:,10) = 800;
T(10,:)= 900;
T_old = T;
err = 1;
counter = 0;
iter = 7;
w = 1.2;
alpha = 0.3;
dt = 0.01;
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
T_previous = T;



%jacobi 
if iter == 1
    while(err > tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
             T(i,j) = (1/k)*(T_old(i-1,j)+T_old(i+1,j))/dx^2 + (1/k)*(T_old(i,j-1)+T_old(i,j+1))/dy^2; 
            
            end
        
        end
    
        counter = counter + 1;
    
        err = max(max(abs(T-T_old)));
    
        T_old = T;
        
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
        
    end
end

%GS
if iter == 2
    
    while(err > tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = (1/k)*(T(i-1,j)+T(i+1,j))/dx^2 + (1/k)*(T(i,j-1)+T(i,j+1))/dy^2; 
            
            end
        
        end
    
        counter = counter + 1;
    
        err = max(max(abs(T-T_old)));
    
        T_old = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
    
    end
    
end

%SOR
if iter == 3
        
    while(err > tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = (1-w)*T_old(i,j) + w*((1/k)*(T(i-1,j)+T(i+1,j))/dx^2 + (1/k)*(T(i,j-1)+T(i,j+1))/dy^2); 
            
            end
        
        end
    
        counter = counter + 1;
    
        err = max(max(abs(T-T_old)));
    
        T_old = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
    
    end
    
end

%explicit
if iter == 4
        
    for nt = 1:200
        
        %while(err>tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = T_old(i,j) + (alpha*dt/dx^2)*(T_old(i-1,j)+T_old(i+1,j)-2*T_old(i,j)) + (alpha*dt/dy^2)*(T_old(i,j-1)+T_old(i,j+1)-2*T_old(i,j)); 
            
            end
        
        end
    
        counter = counter + 1;
    
        %err = max(max(abs(T-T_old)));
    
        T_old = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
    
    %end
    
    end

end

term1 = 1/(2*k1+2*k2+1);
term2 = k1*term1;
term3 = k2*term1;


%implicit (Jacobi)
if iter == 5
        
    for nt = 1:200
        err = 10;
        while(err>tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = T_previous(i,j)*term1 + term2*(T_old(i-1,j)+T_old(i+1,j)) + term3*(T_old(i,j-1)+T_old(i,j+1)); 
            
            end
        
        end
    
        counter = counter + 1;
    
        err = max(max(abs(T-T_old)));
    
        T_old = T;
        
    
        end
        T_previous = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
    
    end

end

%implicit (GS)
if iter == 6
        
    for nt = 1:200
        err = 10;
        while(err>tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = T_previous(i,j)*term1 + term2*(T(i-1,j)+T(i+1,j)) + term3*(T(i,j-1)+T(i,j+1)); 
            
            end
        
        end
    
        counter = counter + 1;
    
        err = max(max(abs(T-T_old)));
    
        T_old = T;
        
    
        end
        T_previous = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
    
    end

end

%implicit (SOR)
if iter == 7
        
    for nt = 1:200
        err = 10;
        while(err>tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = (1-w)*T_previous(i,j) + w*(T_previous(i,j)*term1 + term2*(T(i-1,j)+T(i+1,j)) + term3*(T(i,j-1)+T(i,j+1))); 
            
            end
        
        end
    
        counter = counter + 1;
    
        err = max(max(abs(T-T_old)));
    
        T_old = T;
        
    
        end
        T_previous = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
    
    end

end

The graph obtained is :

the graph is same for all the schemes.

Results and Discussion:

All the schemes converged the solution. The only difference is the time required for convergence or number of iterations were different. Explicit scheme found to be more faster than implicit scheme. But explicit scheme has to fulfill the CFL condition where as implicit scheme is unconditionally stable.

In implicit scheme - SOR found to be more faster than GS or Jacobi scheme.

Now we will discuss the convergence rate of all these schemes:

Steady state includes:

Jacobi

Gauss-Seidal

Successive over relaxation

Jacobi:

Gauss-Seidal:

SOR:

Transient simulation include:

Explicit

Implicit - Jacobi

Implicit - Gauss-Seidal

Implicit - SOR

here we will compare simulation time:

Explicit

Implicit - Jacobi

Implicit - Gauss Seidal

Implicit - SOR

For steady state : From the above table, it is clear that SOR scheme is the fastest when we select the weight properly.

For unsteady state: Explicit scheme is fastest then implicit schemes.

The reason for less iterations in SOR is because it is multiplied with weight which increases convergence rate then other two.

The reason for less simulation time in explicit scheme is the equations are relatively simple and not coupled. So it takes less time than Implicit scheme where the equations are coupled and complex.

In Last part we will discuss the stability criteria of explicit and implicit schemes.

Aim: To study the unstabilty criteria of explicit and implicit scheme.

Theory: When the sudden variations are seen in the values and it starts oscillating very fast, it means system has become unstable.

stability criteria or CFL condition is

`alpha*dt/dx^2 <0.5` for two dimension`

where `alpha*dt/dx^2` is the CFL number.

as CFL number becomes greater than 0.5, solution becomes unstable.

explicit scheme has to follow the CFL criteria where as implicit scheme is unconditionally stable. Means at any value of dt it remains stable.

first we will see stable solution:

clear all
close all
clc;

tol = 1e-4;
L = 1;
n = 10;
dx = L/(n-1);
dy = L/(n-1);
k = 2*(dx^2+dy^2)/(dx^2*dy^2);
T = zeros(10,10);
T(:,1) = 400;
T(1,:) = 600;
T(:,10) = 800;
T(10,:)= 900;
T_old = T;
err = 1;
counter = 0;
iter = 4;
w = 1.2;
alpha = 0.3;
dt = 0.01;
k1 = alpha*dt/dx^2;
k2 = alpha*dt/dy^2;
T_previous = T;


%explicit
if iter == 4
        
    for nt = 1:200
        
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = T_old(i,j) + (alpha*dt/dx^2)*(T_old(i-1,j)+T_old(i+1,j)-2*T_old(i,j)) + (alpha*dt/dy^2)*(T_old(i,j-1)+T_old(i,j+1)-2*T_old(i,j)); 
            
            end
        
        end
    
        counter = counter + 1;
    
        T_old = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
        
      
    end

end

term1 = 1/(2*k1+2*k2+1);
term2 = k1*term1;
term3 = k2*term1;


%implicit (Jacobi)
if iter == 5
        
    for nt = 1:200
        err = 10;
        while(err>tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = T_previous(i,j)*term1 + term2*(T_old(i-1,j)+T_old(i+1,j)) + term3*(T_old(i,j-1)+T_old(i,j+1)); 
            
            end
        
        end
    
        counter = counter + 1;
    
        err = max(max(abs(T-T_old)));
    
        T_old = T;
        
    
        end
        T_previous = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
    
    end

end

%implicit (GS)
if iter == 6
        
    for nt = 1:200
        err = 10;
        while(err>tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = T_previous(i,j)*term1 + term2*(T(i-1,j)+T(i+1,j)) + term3*(T(i,j-1)+T(i,j+1)); 
            
            end
        
        end
    
        counter = counter + 1;
    
        err = max(max(abs(T-T_old)));
    
        T_old = T;
        
    
        end
        T_previous = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
    
    end

end

%implicit (SOR)
if iter == 7
        
    for nt = 1:200
        err = 10;
        while(err>tol)
    
        for i = 2:n-1
        
            for j = 2:n-1
            
                T(i,j) = (1-w)*T_previous(i,j) + w*(T_previous(i,j)*term1 + term2*(T(i-1,j)+T(i+1,j)) + term3*(T(i,j-1)+T(i,j+1))); 
            
            end
        
        end
    
        counter = counter + 1;
    
        err = max(max(abs(T-T_old)));
    
        T_old = T;
        
    
        end
        T_previous = T;
        contourf(T,'ShowText','on')
        set(gca,'YDIR','reverse')
    
    end

end

the graph of this code is as follows:

this is stable solution for explicit scheme. The CFL number in this case was 0.243 which is less than 0.5

Now, we will increase the dt from 0.01 to 0.1 seconds. It means the CFL number will become 2.43 which is more than 0.5.

the results shows:

this is unstable solution as some values are in positive and some are in negative.

Now we will discuss about implicit scheme. It is unconditionally stable. There is no CFL criteria need to fulfill.

implicit - jacobi solver

dt in this case kept - 0.01 sec.

graph shows:

This is stable solution.

Now we change dt to 0.1 sec.

Now the graph is:

the same stable solution. We can try this for higher time step and also for other schemes like GS or SOR. the result will remain same.

Thus it proves Implicit scheme is unconditionally stable while explicit is stable with CFL criteria.

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.