Solution & Analysis of Steady and Unsteady 2D Heat Conduction using MATLAB

Abstract : 

The heat conduction equation is a partial differential equation that describes the distribution of heat (or the temperature field) in a given body over time. Detailed knowledge of the temperature field is very important in thermal conduction through materials. The equation have been solved using finite difference method, using this discretized equations and boundary conditions a matlab program is created, where both steady and unsteady states of this equation are solved using iterative solvers like Jacobi, Gauss seidel and Successive over relaxation methods. 

Study has been conducted to understand the Computational time taken to solve each of these solvers and comparison is done. 

Study on the effect of time step on the stability of the Transient state of the heat conduction equation of both Implicit and explicit methods is done and understanding how and why an unstable solution forms, how CFL condition is developed using Von neumann analysis and how solution behaves using it.

                                       

I . Objectives :- 

1. Solution and Analysis of Steady and Unsteady 2D Heat Conduction Problem
2. In which, Solving the 2D heat conduction equation by using point iterative methods, such as : 
3.   a.  Jacobi method  b. Gauss-Seidelmethod  c. Successive Over-relaxation method
4. For transient case, Solving by both Implicit & Explicit Schemes.
5. Comparing Fastest steady state simulation vs fastest transient state simulation and analyzing results. 
6. Effect of time step on the stability of the Transient solution
7. Evaluating unstable solution and why it forms. 
8. Evaluating Courant-Friedrichs-Lewy condition for the solution.  

 

II. 2D Heat Conduction (Governing equation) :    `(delT)/(delt) = alpha((del^2T)/(delx^2))+((del^2T)/(dely^2))`

where,    

T : Temperature

t : time

`alpha` = Thermal expansion coefficient

 

II. A) Solving the equation for Steady state using Finite difference method :

Equation for Steady state : `((del^2T)/(delx^2))+((del^2T)/(dely^2)) = 0`

Now, using central difference scheme, eqn becomes, 

`(T_(i+1,j)^n-2T_(i,j)^n+T_(i-1,j)^n)/ (Deltax^2) + (T_(i,j+1)^n-2T_(i,j)^n+T_(i,j-1)^n)/ (Deltay^2) = 0`

on simplification, we get, 

`T_(i,j)^n = 1/k*[(T_(i+1,j)^n + T_(i-1,j)^n)/(Deltax^2) + (T_(i,j+1)^n + T_(i,j-1)^n)/(Deltay^2)]`    --(1)

    where, `k = 2*((Deltax^2+Deltay^2)/(Deltax^2*Deltay^2))` 

 

II. B) Solving the equation for Unsteady (Transient) state using Finite difference method : 

• Explicit Scheme : 

Equation for Transient state : `(delT)/(delt) + alpha((del^2T)/(delx^2)+(del^2T)/(dely^2)) = 0`

For ease, lets say  -

 `T_(i+1,j)^n = R`

`T_(i-1,j)^n = L`

`T_(i,j+1)^n = T`

`T_(i,j-1)^n = B`

`T_(i,j)^n = P`

 

using first order forward difference in time derivative and central difference in x,y derivative, we get,

`T_P^(n+1) = T_P^n + alpha*Deltat[(T_R^n - 2T_P^n + T_L^n)/(Deltax^2) + (T_T^n - 2T_P^n + T_B^n)/(Deltay^2)]` `

On Further Simplification, 

`T_P^(n+1) = T_P^n(1-2k_1-2k_2) + k_1T_L^n+ k_1T_R^n + k_2T_B^n + k_2T_T^n`      --(2) 

 where, `k_1 =( alpha*Deltat)/(Deltax^2)` & `k_2 = ( alpha*Deltat)/(Deltay^2)` 

 Similarly for, 

 Implicit Scheme : 

From the transient state equation, solving using FDM, we get, 

`T_P^(n+1)*(1+2k_1 + 2k_2) = T_P^n + k_1(T_L^(n+1)+ T_R^(n+1)) + k_2(T_T^(n+1) + T_B^(n+1))`  -- (3)

 

II. C) Point Iterative Methods : 

These equations obtained are solved using iterative methods such as Jacobi, Gauss-Seideland S.O.R, in matlab. 

 

III.  Assumptions & Boundary Condition for the problem : 

1. Assume that the domain is a unit square. 

2. Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]

3. Boundary conditions for steady and transient case

    1. Left:400K

    2. Top:600K

    3. Right:800K

    4. Bottom:900K

4. Tolerance = 1e-4 (for all solvers and schemes)

 

IV. Matlab Program : 

IV A.) Main Program : 

% Program to solve a 2D Heat Conduction equation
% dT/dt = alpha*(dT^2/dx^2 + dT^2/dy^2
% alpha = Thermal conductivity
% rf = relaxation factor (SOR method)
% dx = dy = grid space (domain = unit square)
% nx = ny [Number of points along the x direction is equal to the number of points along the y direction]

% Boundary conditions for steady and transient case (both schemes)
%
%     1. Left:400K
%     2. Top:600K
%     3. Right:800K
%     4. Bottom:900K
%
%                                      Code by AADIL SHAIKH
%                                          DATE : 7/29/2019

close all
clear all
clc

disp('1. Steady state solver')
disp('2. Explicit Scheme - Transient solver')
disp('3. Implicit Scheme - Transient solver')

% Choose a solver from the above 3 displayed solvers .
Solver = input('Choose a solver from above : ');

% After choosing solvers , choose the Iterative method for which the
% solution is needed .
% 1. Jacobi method
% 2. Gauss Seidel method
% 3. Successive over relaxation method

if Solver == 1
    [x,y,T] = Steady_state ;
    
elseif Solver == 2
    [x,y,T] = Explicit_transient ;
    
elseif Solver == 3
    [x,y,T] = Implicit_transient ;
end

This is the Main program which basically calls anyone of the three solvers, according to the users choice. 

Namely,

Steady state , Unsteady - Explicit scheme & Unsteady state - Implicit scheme .

 

IV. B) Steady state solver : Iterative methods - Jacobi, Gauss Seidel , S.O.R

 

function  [x,y,T] = Steady_state

%Domain  inputs
nx = 40;
ny = 40;

rf = 1.5 ;   % Relaxation factor

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2)-x(1);
dy = y(2)-y(1);

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

% assigining bc
T = ones(nx,ny);
T(:,1) = 400;
T(1,:) = 600;
T(:,end) = 800;
T(end,:) = 900;

Told = T;

disp('1. Jacobi method')
disp('2  Gauss Seidel method')
disp('3. S.O.R method')

iterative_solver = input('Choose which Iterative method to solve with : ');
error = 9e9;
tol = 1e-4;

%jacobi method
tic
iter_ctr = 1;
if iterative_solver ==1
    while (error > tol)
        
        for i = 2:nx-1
            for j = 2:ny-1
                H = (Told(i+1,j)+Told(i-1,j))/(dx^2);
                V = (Told(i,j+1)+Told(i,j-1))/dy^2;
                T(i,j) = (H + V)/k;
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        iter_ctr = iter_ctr + 1;
        
    end
    Simulation_time = toc;
    figure(1)
    colormap(jet)
    contourf(x,y,T,':','showText','on')
    set(gca, 'YDir','reverse')
    xlabel('x axis')
    ylabel('y axis')
    title({['steady state - Jacobi'];['Iteration number : ',num2str(iter_ctr)];['Simulation time : ',num2str(Simulation_time)]})
    colorbar
    % pause(0.0001)
    saveas(1,'ss_jacobi','png')
end



%Gauss Seidel method
tic
iter_ctr = 1;
if iterative_solver == 2
    while (error > tol)
        
        for i = 2:nx-1
            for j = 2:ny-1
                H = (Told(i+1,j)+T(i-1,j))/(dx^2);
                V = (Told(i,j+1)+T(i,j-1))/dy^2;
                T(i,j) = (H + V)/k;
            end
        end
        error = max(max(abs(Told-T)));
        
        Told = T;
        iter_ctr = iter_ctr + 1;
    end
    Simulation_time = toc;
    
    figure(2)
    colormap(jet)
    contourf(x,y,T,':','showText','on')
    set(gca, 'YDir','reverse')
    xlabel('x axis')
    ylabel('y axis')
    title({['steady state  - Gauss Seidel'];['Iteration number : ',num2str(iter_ctr)];['Simulation time : ',num2str(Simulation_time)]})
    colorbar
    % pause(0.0001)
    saveas(2,'ss_gs','png')
end



% Successive Over Relaxation method
tic
iter_ctr = 1;
if iterative_solver == 3
    while (error > tol)
        
        for i = 2:nx-1
            for j = 2:ny-1
                H = (Told(i+1,j)+T(i-1,j))/(dx^2);
                V = (Told(i,j+1)+T(i,j-1))/dy^2;
                T(i,j) = Told(i,j)*(1-rf) + rf*((H + V)/k);
                
            end
        end
        error = max(max(abs(Told-T)));
        Told = T;
        iter_ctr = iter_ctr + 1;
    end
    Simulation_time = toc;
    
    figure(3)
    colormap(jet)
    contourf(x,y,T,':','showText','on')
    set(gca, 'YDir','reverse')
    xlabel('x axis')
    ylabel('y axis')
    title({['Steady state - Successive Over Relaxation'];['Iteration number : ',num2str(iter_ctr)];['Simulation time : ',num2str(Simulation_time)]})
    colorbar
    % pause(0.0001)
    saveas(3,'ss_sor','png')
end
end

1.  Once the user chooses Steady state solver, it calls this function. It solves the steady state 2D heat conduction equation (1) derived using finite difference method as shown in the start of the report. 
2. The number of points chosen along x and y direction are 40 and dx,dy the space difference for each grid is determined from this parameter. 
3. Error and Tolerance are set and the equation (1) is placed inside the Convergence loop where the solution solves until convergence takes place.  
4. Iteration counter counts the number of iterations taken for the solver to reach convergence.
5. Simulation / Computational Time is calculated. 

 Results :                                                       

                                                         

                                                   

                   

1. As observed, The jacobi method takes most computational time, as well as number of iterations. 
2. This happens because in jacobi method, The old values are taken for convergence and from that the new value is derived. Then that becomes the old value and new one is derived. 
3. And in the Gauss Seidel method, The values we solve for are updated in both directions of the matrix (in the equation).  Hence the new value is passed on each iterations as the solution proceeds .  Thus convergenece is obtained faster and computational time is less as well as the iterations are lesser than jacobi. 

IV. B-1) SOR :

1. Successive over relaxation method can be used on either gauss Seidel or jacobi . 
2. Here we have used SOR on gauss Seidel , as its a faster method and more accurate than jacobi, with sor it will converge more accurately and computational time can be faster. If it takes more computational time in case of explicit scheme, it will give higher convergence.
3. As can be observed in the sor plot above . 
4. Sor method is like a performance enhancer. 
5. The term rf, used in the equation stands for Relaxation factor. 
6. rf > 1 - aggressive (sor)   ;  rf < 1 - not agressive (successive under relaxation)
7. The relaxation factor enhances the solution to converge faster. 
8. The Formula used for sor , as generalized for this project is  : 
9. `T_P^(n+1) = T_P^n*(1-rf) + rf*(T_(gauss seidel) or T_(jacobi))` 
10. The optimum value of rf can also be found using : `rf = 2/(1+sin(pi*dx)`
11. The sor converges for rf values 0<rf<2 . so simply 1.5 factor has been chosen in the problems. 
12. For Other theorems and for models, The rf formula is slightly moded and can be found in the references at the end of the report. 
13. where an rf of more than 1 enhances the convergence rate.  

 

IV. C) Transient state solver - Explicit Scheme : Iterative methods - Jacobi, Gauss seidel, S.O.R

function [x,y,T] = Explicit_transient
%Domain  inputs
nx = 40;
ny = 40;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2)-x(1);
dy = y(2)-y(1);

dt = 0.0001;
alpha = 1.4;
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);

nt = 400;

rf = 1.5 ;
CFL_Num = alpha*dt/(dx^2);

% assigining bc
T = ones(nx,ny);
T(:,1) = 400;
T(1,:) = 600;
T(:,end) = 800;
T(end,:) = 900;

Told = T;
Tprev_dt = T;

disp('1. Jacobi method')
disp('2. Gauss Seidel method')
disp('3. S.O.R method')
iterative_solver = input('Choose which Iterative method to solve with : ');
error = 9e9;
tol = 1e-4;

%jacobi method
tic
if iterative_solver == 1
    for k = 1:nt
        for i = 2:nx-1
            for j = 2:ny-1
                Term1 = (1 -2*k1 - 2*k2);
                Term2 = k1;
                Term3 = k2;
                H = (Told(i+1,j)+Told(i-1,j));
                V = (Told(i,j+1)+Told(i,j-1));
                T(i,j) = (Told(i,j)*Term1) + (Term2*H) + (Term3*V);
            end
        end
        Told = T;
        
    end
    Simulation_time = toc;
    
    
    figure(1)
    colormap(jet)
    contourf(x,y,T,':','showText','on')
    set(gca, 'YDir','reverse')
    xlabel('x axis')
    ylabel('y axis')
    title({['Explicit - Transient solver using Jacobi '];['Simulation time : ',num2str(Simulation_time)]})
    colorbar
    saveas(1,'EP_jacobi','png')
    
end



%   Gauss Seidel method
tic
if iterative_solver == 2
    for k = 1:nt
        for i = 2:nx-1
            for j = 2:ny-1
                Term1 = (1 -2*k1 - 2*k2);
                Term2 = k1;
                Term3 = k2;
                H = (Told(i+1,j)+T(i-1,j));
                V = (Told(i,j+1)+T(i,j-1));
                T(i,j) = (Told(i,j)*Term1) + (Term2*H) + (Term3*V);
            end
        end
        Told = T;
    end
    Simulation_time = toc;
    
    
    figure(2)
    colormap(jet)
    contourf(x,y,T,':','showText','on')
    set(gca, 'YDir','reverse')
    xlabel('x axis')
    ylabel('y axis')
    title({['Explicit - Transient solver using Gauss Seidel '];['Simulation time : ',num2str(Simulation_time)]})
    colorbar
    saveas(2,'EP_GS','png')
end


%   S.O.R method
tic
if iterative_solver == 3
    for k = 1:nt
        for i = 2:nx-1
            for j = 2:ny-1
                Term1 = (1 -2*k1 - 2*k2);
                Term2 = k1;
                Term3 = k2;
                H = (Told(i+1,j)+T(i-1,j));
                V = (Told(i,j+1)+T(i,j-1));
                T(i,j) = Told(i,j)*(1-rf) + rf*((Told(i,j)*Term1) + (Term2*H) + (Term3*V));
            end
        end
        Told = T;
        
    end
    Simulation_time = toc;
    
    
    figure(3)
    colormap(jet)
    contourf(x,y,T,':','showText','on')
    set(gca, 'YDir','reverse')
    xlabel('x axis')
    ylabel('y axis')
    title({['Explicit - Transient solver using Successive Over Relaxation '];['Simulation time : ',num2str(Simulation_time)]})
    colorbar
    saveas(3,'EP_SOR','png')
end
end

1.  This function is called when solver 2 is chosen by the user. It solves the 2D Transient state heat conduction equation (2), by point iterative methods mentioned above. 
2. Since there is a time derivative, this equation is solved using Time marching method in FDM. 
3. And as its an explicit scheme, from eqn (2) all the values in the rhs are from same time step & known and we solve for the Temperature in next time step on the lhs, the unknown.  
4. The explicit scheme has no convergence loop because of this. 
5. The solution runs iterations for the number of time loop. 400 in this case.
6. It is solved for all the three iterative methods . 

Results : 

                                  Transient State - Explicit scheme : Jacobi (nt = 400)

               

 

                               Transient State - Explicit scheme : Gauss Seidel (nt = 400)

                

 

                    Transient State - Explicit scheme : Successive Over Relaxation (nt = 400)             

               

1.  As there is no convergence loop for explicit scheme, it runs for the number of time loop (nt = 400)
2. The solution tries to approach steady state with jacobi but hasnt fully reached for nt = 400. but when this is increased, The observation is that it shows the plot close to steady state of the equation. 
3. But for the same nt = 400, Gauss Seidel method shows more accurate distribution throughout the plot. It takes more simulation time than jacobi but gives better result. 
4. Successive over relaxation plot is the most enhanced and closer out of 3 results for explicit scheme. 
5. rf is set to 1.5. and Sor is performed over Gauss Seidel . sor takes the most computational time in explicit case.  The formula used is modified here, for explicit scheme. Details given in : (IV. B-1) section.

 

IV. D) Transient state solver - Implicit Scheme : Iterative methods - Jacobi, Gauss seidel, S.O.R

function  [x,y,T] = Implicit_transient ;

%Domain  inputs
nx = 40;
ny = 40;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2)-x(1);
dy = y(2)-y(1);

dt = 0.001;
alpha = 1.4;
k1 = alpha*dt/(dx^2);
k2 = alpha*dt/(dy^2);

CFL_Num = alpha*dt/(dx^2)

nt = 400;
rf = 1.5 ;               % relaxation factor

% assigining bc
T = ones(nx,ny);
T(:,1) = 400;
T(1,:) = 600;
T(:,end) = 800;
T(end,:) = 900;

% creating a copy of u
%uold = u ;


Told = T;
Tprev_dt = T;

disp('1. Jacobi method')
disp('2 Gauss Seidel method')
disp('3. S.O.R method')
iterative_solver = input('Choose which method to solve with : ');
error = 9e9;
tol = 1e-4;

%jacobi method
tic
jacobi_iter = 1;
if iterative_solver ==1
    for k = 1:nt
        error = 9e9;
        while(error > tol)
            for i = 2:nx-1
                for j = 2:ny-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_dt(i,j)*Term1) + (Term2*H) + (Term3*V);
                end
            end
            
            error = (max(abs(T-Told)));
            error = max(error);
            Told = T;
            jacobi_iter = jacobi_iter + 1;
            
        end
        Tprev_dt = T;
        
        
    end
    Simulation_time = toc;
    
    figure(1)
    colormap(jet)
    contourf(x,y,T,':','showText','on')
    set(gca, 'YDir','reverse')
    xlabel('x axis')
    ylabel('y axis')
    title({['Implicit - Transient solver using Jacobi '];['Iteration number : ',num2str(jacobi_iter)];['Simulation time : ',num2str(Simulation_time)]})
    colorbar
    saveas(1,'im_jacobi','png')
end


% Gauss Seidel method
tic
gs_iter = 1;
if iterative_solver == 2
    for k = 1:nt
        error = 9e9;
        while(error > tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    Term1 = (1 +2*k1 + 2*k2)^-1;
                    Term2 = k1*Term1;
                    Term3 = k2*Term1;
                    H = (Told(i+1,j)+T(i-1,j));
                    V = (Told(i,j+1)+T(i,j-1));
                    T(i,j) = (Tprev_dt(i,j)*Term1) + (Term2*H) + (Term3*V);
                end
            end
            
            error = (max(abs(Told-T)));
            error = max(error);
            Told = T;
            gs_iter = gs_iter + 1;
            
        end
        Tprev_dt = T;
    end
    Simulation_time = toc;
    
    figure(2)
    colormap(jet)
    contourf(x,y,T,':','showText','on')
    set(gca, 'YDir','reverse')
    xlabel('x axis')
    ylabel('y axis')
    title({['Implicit - Transient solver using Gauss Seidel '];['Iteration number : ',num2str(gs_iter)];['Simulation time : ',num2str(Simulation_time)]})
    colorbar
    saveas(2,'im_gs','png')
    
end


% Successive over relaxation method
tic
sor_iter = 1;
if iterative_solver == 3
    for k = 1:nt
        error = 9e9;
        while(error > tol)
            for i = 2:nx-1
                for j = 2:ny-1
                    Term1 = (1 +2*k1 + 2*k2)^-1;
                    Term2 = k1*Term1;
                    Term3 = k2*Term1;
                    H = (Told(i+1,j)+T(i-1,j));
                    V = (Told(i,j+1)+T(i,j-1));
                    T(i,j) = Told(i,j)*(1 - rf) + rf*((Tprev_dt(i,j)*Term1) + (Term2*H) + (Term3*V));
                end
            end
            
            error = (max(abs(Told-T)));
            error = max(error);
            Told = T;
            sor_iter = sor_iter + 1;
            
        end
        Tprev_dt = T;
    end
    Simulation_time = toc;
    
    figure(3)
    colormap(jet)
    contourf(x,y,T,':','showText','on')
    set(gca, 'YDir','reverse')
    xlabel('x axis')
    ylabel('y axis')
    title({['Implicit - Transient solver using S.O.R '];['Iteration number : ',num2str(sor_iter)];['Simulation time : ',num2str(Simulation_time)]})
    colorbar
    saveas(3,'im_sor','png')
    
end

end

1. This function is called to solve Transient state - Implicit scheme using iterative solvers .
2. From eqn (3), This function is solved, This being an implicit scheme, we have unkown values for the next time step and one value with previous time step using which we solve for the next 'P' time step on the lhs. 
3. This scheme having time derivative is solved using time marching method in FDM. 
4. Being a transient implicit case, there is a time loop and a convergence loop, The convergence loop updates the values of the new time step values i.e at (n+1) with each iteration, while the Time loop updates the value of the previous time step.
5. The solution runs till convergence in reached with regards to the tolerance set.
6. Contour is plotted. 

Results :                                        

                                                     

                                  

                

1. From Observation of the results, all the iterative solvers reach convergence. 
2. However, Jacobi method takes most computational time and iterations whereas Sor takes the least time and No. of iterations is almost less than half of gauss seidel's iterations. 
3. We can conclude the sor enchances the solution fastest and in least number of iterations.
4. The implicit scheme performs most iterations amongst the steady state and explicit cases.
5. Implicit scheme gives better results as compared to the explicit scheme in any of the iterative solvers.

 

v) Comparing the Fastest Steady state simulation vs the Fastest Transient state simulation 

 

 V-1 ) Steady state results : 

1.  Program details and insights on section (IV.B) regarding steady state. 
2. As observed the fastest Simulation for steady state solver is from Iterative technique - Successive over relaxation . (more details on sor - section (IV .B-1))
3. The sor enhances the solution (time and convergence in this case) for the iterative techniques and thats what we observe here. 
4. Steady state converges the fastest in it , as compared to Gauss seidel and then jacobi.
5. As we can see the number of iterations lowers drastically from jacobi to sor method. 
6. The Sor form used here is updated from the Gauss seidel model. 

 

V-2) Unsteady state Results : 

Explicit scheme : 

1. Program details and results on section (IV-C)
2. As observed in explicit scheme the jacobi has given the fastest simulation as compared to the sor and Gauss seidel. 
3. As the explicit scheme does not have convergence loop, it runs for iterations of the time loop (400 in our case)
4. So we can see here for 400 iterations, the better solution which is nearest to the steady state is actually given by the sor method. so practically speaking the sor method will be time saving when highly convergent results are desired because for 400 iterations both the other methods were still in approach to reach full steady state convergence. 
5. But for comparison sake of simulation, the jacobi comes out to be fastest. 

 

Implicit Scheme : 

1.  Program details and results on section (IV-D)
2. As observed in implicit scheme, The sor method gives the fastest simulation.
3. Although implicit scheme takes most simulation time out of all the three solvers.
4. The equation and loops are complex in the implicit scheme and more time demanding.
5. The number of iterations are also more than other two. And iterations drop down from jacobi to sor with gauss seidel in between.
6. All the iterative solvers in implicit scheme reach high convergence.

 

V-3 ) Conclusion for Fastest Steady state simulation vs Transient state : 

1. Explicit scheme - Jacobi is the fastest in the transient state and S.O.R method is the fastest in steady state.
2. Above section explains individually how these are the fastest simulation. 
3. In comparison Explicit jacobi is faster than Steady state sor. 
4. But the transient solver for explicit and implicit scheme has been calculated at nt = 400 (time loop) to keep common terms when comparing. so the explicit scheme's solution is only for 400 iterations . 
5. As the iterations are increased the solution approaches very closely to steady state S.O.R in convergence and simulation time . 
6.                                      Example :    for ( nt = 1000)
7.                            

                                             Results In a tabular form 

                     

 

VI ) Stability Analysis in the Unsteady Problem

A solution can go unstable if there are errors in it . hence by studying those errors we can understand what makes it go unstable & hence its stability criteria

 The numerical solution for such an equation is influenced by two sources of error. 

1. Discretization error - difference between exact analytical solution of the pde and exact solution of the corresponding difference equation .  ( the truncation error) plus any errors introduced by numerical treatment of boundacy conditions.

2. Round off error : Computer rounding off repetitive numbers of calculation.

From equation 

`T_P^(n+1) = T_P^n + alpha*Deltat[(T_R^n - 2T_P^n + T_L^n)/(Deltax^2) + (T_T^n - 2T_P^n + T_B^n)/(Deltay^2)]` 

using a discrete Fourier decomposition , where range of x and y is defined seperately for each direction . Defining the amplification matrix G as in one dimensional case - we move forward, obtaining : 

`alpha*(1/dx^2 + 1/dy^2)*dt <= 1/2`

This stability condition is sufficient and analogous to normal condition. 

but when dx = dy , it puts a severe requirement on the time step , it is reduced further by a factor of 2 compared with one dimensional case, we get,

 `(alpha*dt)/dx^2 <= 1/4`  or    `dt <= dx^2/(4*alpha)`  

where, dx = dy : no. grid points in x,y direction

This equation gives the stability requirement for the solution of finite difference equation above to be stable. . 

Clearly for the value dx the value of dt must be small enough to satisfy the condition. 

As long as this condition is met, error will not grow with subsequent steps and numerical solution will proceed in a stable manner otherwise the error will grow larger and cause the solution to blow up . 

This analysis is called Von Neumann stability method. 

From using the Von neumann stability analysis we can determine Courant-Friedrichs-Levy condition for the stability of an explicit solution of a PDE.  

Thus that condition from the above analysis is :

`C <= 0.25`            , (C/CFL number/ courant number)

If in explicit scheme, CFL number is more than 0.25 the solution will blow up and become unstable. 

However Implicit solvers are unconditionally stable, which means the they can tolerate CFL numbers more larger than CFL condition. whereas the explicit solvers are conditionally stable and can never tolerate values more than CFL condition. 

 

VI-1 ) Observing Stability analysis of the above program : 

A)  Explicit solver: 

 The changing cfl number and time steps are mentioned in the plots itself and as we can see the cfl number is under 0.25 for nt = 400 iterations . The results are stable and approaching the steady state.

Thermal diffusivity used is 1.4 for all conditions. 

 

• From here on the time steps are varied to see how solution acts when the cfl condition is crossed.

 For just a little increase over the cfl , we see that th jacobi has become unstable. it tried to solve a we can see from the initial plot shows the temperature 500 but failed further and blew up. however gauss seidel and SOR managed to solve. 

 

 

 For a little more increase in the CFL number , jacobi has become fully unstable, however Gauss seidel method was successful in achieving the result , however the SOR updated from Gauss seidel became some what unstable and blew up .

 An increase in approximately 8 times the cfl condition totally blew up all the solvers for explicit. and became totally unstable that it couldnt produce it in the plot. 

 

B ) Implicit solver : 

Since implicit solvers are unconditionally stable, increasing the cfl number did not affect the stability of the solution . Experiments with really high courant number was done but it still did not become unstable. 

 VII ) References : 

1. CFD textbook - John D anderson
2. http://pubs.sciepub.com/ajams/4/4/3/index.html
3. https://www.flow3d.com/resources/cfd-101/numerical-issues/implicit-versus-explicit-numerical-methods/
4. https://www.simscale.com/blog/2017/08/cfl-condition/
5. https://pdfs.semanticscholar.org/a5cb/6355d1b1b01569f11214328441e5728984a7.pdf
6. https://nptel.ac.in/courses/111107063/module10/lecture2/lecture2.pdf
7. https://en.m.wikipedia.org/wiki/Von_Neumann_stability_analysis

keywords - NUMERICAL-ANALYSIS, MATLAB, MATLAB-BASICS, FINITE DIFFERENCE, 2D HEAT CONDUCTION