Stability Analysis of Explicit and Implicit method used for solving Transient 2D heat conduction equation using MATLAB

Objective:

Compare the effect of time step on the stability of the solution of 2D Transient heat conduction using implicit and explicit methods. Also, explain

1. How does an unstable solution look like?

2. At what CFL number does the solution become unstable?

Introduction:

As discussed in the previous challenge, a transient heat conduction equation can be solved by using Explicit or Implicit methods, so for this study, we will be focusing on the Explicit method and the SOR on GS Implicit method.

To solve the 2D heat conduction equation by using the point iterative or time marching techniques, the Boundary conditions for the problem is specified as follows:

Top Boundary = 600 K

Bottom Boundary = 900 K

Left Boundary = 400 K

Right Boundary = 800 K

The General governing partial differential equation for heat transfer is:

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

where T is the temperature, `alpha` is thermal diffusivity and x,y,z are space variables and t is time variable.

The governing partial differential equation for the 2D unsteady-state case is :

`(dT)/(dt) = alpha*((del^2T)/(delx^2) + (del^2T)/(dely^2))`

Discretization of GDE:

The unsteady-state GDE is hyperbolic in nature as it depends on the information from left to the right direction in the mesh therefore, Forward Differencing Scheme is used to discretize this second-order partial differential equation.

Let `T_(i,j)^(n,m)text( be the temperature at ) i^(th) text( row & )j^(th) text( column in the assumed matrix field at )n^(th) text( time-step and )m^(th)text( iteration.)`

This equation can be discretized as:

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

Where:

`T_P= T_(i,j) text( (at point p))`

`T_R= T_(i+1,j)text( (at right point))`

`T_L= T_(i-1,j)text( (at left point))`

`T_T= T_(i,j+1)text( (at top point))`

`T_B= T_(i,j-1)text( (at bottom point))`

on rearranging the equation, we get:

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

Let us consider `K1=alpha*(Deltat)/(Deltax^2)` and `K2=alpha*(Deltat)/(Deltay^2)`

By further simplification, we get:

`T_P^(n+1) = (1-2*K1-2*K2)*T_P^(n) + K1*{T_L^n+T_R^n} + K2*{T_B^n+T_T^n}` -------------------------(1)

Here, K1 and K2 are X and Y components of CFL number. The CFL number stands for Courant–Friedrichs–Lewy number. For a stable solution, the sum of components of CFL number should be equal to or less than 0.5. Therefore,

`K1=CFL_x`

`K2=CFL_y`

`CFL_x+CFL_y <=0.5`

`:. K1+K2<=0.5` ------------------------------(2)

From this condition, the size of time step can be calculated by following formula:

`dt= (Deltay^2*Deltax^2)/(2*alpha*(Deltax^2+Deltay^2))`

Where:

`dx=X/(nx) text( ,X is size of body in x direction in meters & nx is no. of divisions in x direction)`

`dy=Y/(ny) text( ,Y is size of body in y direction in meters & ny is no. of divisions in y direction)`

This equation can be solved by using fixed-point iterative solving methods like Jacobi, Gauss-Siedal, Successive Over Relaxation on Gauss-Siedal. To study the comparison between these solvers, we first derive there solving equations and then compare the solution, the number of iterations taken to converge on MATLAB.

Explicit Method for 2D Unsteady-state heat conduction:

`T_P^(n+1) = (1-2*K1-2*K2)*T_P^(n) + K1*{T_L^n+T_R^n} + K2*{T_B^n+T_T^n}`

Implicit (SOR on Gauss-Seidel) Method for 2D Unsteady-state heat conduction:

The SOR method is an advanced method that is used to scale the solution of other iterative solvers. If the relaxation factor "Omega" `omega` is greater than 1, then the solution will converge fast, and if it's less than 1 then the solution will take more iterations to reach the desired solution. The formula of SOR applied on GS method is:

`T_P^(n+1,m) =(1-omega)* T_P^(n+1,m-1) + omega*(T_(GS))`

or

`T_P^(n+1,m) =(1-omega)* T_P^(n+1,m-1) + omega*(T_P^(n) /(K3)+ (K1)/(K3)*{T_L^(n+1,m)+T_R^(n+1,m-1)} + (K2)/(K3)*{T_B^(n+1,m)+T_T^(n+1,m-1)})`

`text(where, )omega text( is relaxation factor, if )omega>1 text( solution is over relaxed or if )omega<1 text( solution is under relaxed)`

Experiment details:

Initial conditions are visualized by creating a 10*10 matrix and assigning the wall temperatures as mentioned, all inside points are initially at 320K. A filled contour is drawn from this array to visualize temps at all points.

This matrix is visualized using "contourf" command in MATLAB, the rainbow color scheme is used to show the gradient of temperature:

For the comparison of effect of `dt` on solvers, three different `CFL` numbers are chosen as the time step depends on the CFL number.

The calculated `dt` for every `CFL` number is shown below:

Program explanation:

The main program is used to assign all initial, variable, and constant values. These values are then passed on to the respective functions to calculate and plot the results, final answers are retrieved from every function to draw a comparison. The comparison between initial conditions and other types of solvers, a subplot is used.

The bar graphs are used to compare the number of iterations and errors at the last iteration.

• figure(5) will show a 2*2 subplot of filled contours of the initial condition and final conditions of Explicit Unsteady-state solver.
• figure(6) will show a 2*2 subplot of filled contours of the initial condition and final conditions of Implicit SOR on GS Unsteady-state solver.
• figure(7) will show a bar plot of comparison of the number of iteration required for different CFL number to solve or fail un-steady state problem using explicit
• figure(8) will show a bar plot of comparison of the number of iteration required for different CFL number to solve or fail un-steady state problem using implicit SORGS

Global Variables specified in the main program:

• No. of divisions in X & Y direction respectively
  • nx=10
  • ny=10
• `Deltax" & "Deltay` (dx & dy) is calculated by taking a difference between two consecutive points on X & Y direction respectively
• Thermal diffusivity constant `(alpha)` is set to a=1.5
• `Deltat` (dt) is time step size, which is calculated using values of `Deltax" & "Deltay` according CFL stability criteria
• Max. no. of the time step is specified as nt=1400
• Value of CFL numbers are 0.3,0.5,1 which is specified in an array
• error tolerance value for space loop of Un-Steady state solvers, tol=1e-4
• defining SOR factor: sor_factor=1.2

Main Program:

clear all
close all
clc
%This program calls iterative solver functions to solve a given problem
%Function 1 is to solve unsteady state problem using Explicit method
%Function 2 is to solve unsteady state problem using Implicit-SOR-GS iterative method
nx=10; %no. of grid points in X-direction
ny=10; %no. of grid points in Y-direction
t_initial=320*ones(ny,nx); %initial condition for inner nodes
t_initial(2:ny-1,1)=400;%initial condition for boundary nodes
t_initial(1,2:nx-1)=600;%initial condition for boundary nodes
t_initial(2:ny-1,nx)=800;%initial condition for boundary nodes
t_initial(ny,2:nx-1)=900;%initial condition for boundary nodes
t_initial(1,1)=500;%initial condition for corner nodes
t_initial(1,nx)=700;%initial condition for corner nodes
t_initial(ny,1)=650;%initial condition for corner nodes
t_initial(ny,nx)=850;%initial condition for corner nodes

% creating grid
x=linspace(0,1,nx); % x grid points
dx=abs(x(2)-x(1)); %distance between 2 points in x direction
y=linspace(0,1,ny); % y grid points
dy=abs(y(2)-y(1)); %distance between 2 points in y direction

figure(1)
contourf(t_initial,20) % plotting initial conditions
colorbar; colormap(jet);
axis ij %to flip the Y-axis
axis equal
grid on
xlabel('X-axis')
ylabel('Y-axis')
title_text1=sprintf('Initial conditions');
title_text2=sprintf('Initial condition: Top 600K, Left 400K, Right 800K, Bottom 900K, Rest 320K');
title({[title_text1],[title_text2]}) 

a=1.5; %thermal diffusivity
cfl=[0.3,0.5,1] ;%array of CFL nos.
o=length(cfl);
nt=500;%number of time steps
%difining tolerence values
tol=1e-4; %error tolerance value for space loop of Un-Steady state solvers
%defining sor factor
sor_factor=1.2;
%calulating time step
for z=1:o
    q=cfl(z);
    dt(z)= (q*(dy^2)*(dx^2))/(a*((dx^2)+(dy^2)));
    k1(z)= a*dt(z)/dx^2;
    k2(z)= a*dt(z)/dy^2;
    k3(z)= 1/(1+2*k2(z)+2*k1(z));
end
% calling functions to solve the problem
[unsteady_exp,counter1,error1,ts1]= explicit_func_stability(t_initial,nx,ny,k1,k2,nt,dt,o); %Function 1
[unsteady_sorgs,counter4,error4,ts4]= sorgs_imp_stability(t_initial,nx,ny,k1,k2,k3,tol,nt,dt,sor_factor,o) ;%Function 2

% plotting a subplot to compare the results using a filled contour
figure(5) %Explicit method
subplot(2,2,1) %initial
contourf(t_initial) % plotting initial conditions
colorbar; colormap(jet);
axis ij %to flip the Y-axis
axis equal
grid on
xlabel('X-axis')
ylabel('Y-axis')
title_text1=sprintf('Initial conditions, Thermal Diffusivity= %.1f',a);
title_text2=sprintf('Initial condition: Top 600K, Left 400K, Right 800K, Bottom 900K, Rest 320K');
title({[title_text1],[title_text2]}) 

subplot(2,2,2)
contourf(unsteady_exp(:,:,1),20) %timestep 1
colorbar; colormap(jet);
title_text1=sprintf('Unsteady-Explicit Method, iteration number = %d, CFL= %.1f',counter1(1),cfl(1));
title_text2=sprintf('Timestep size = %.10f sec, timestep = %d',dt(1),ts1(1));
title_text3=sprintf('at time = %.4f sec, temp. diff. = %.10f K',(dt(1)*ts1(1)),error1(1));
title({[title_text1],[title_text2],[title_text3]}) 
axis ij  %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')

subplot(2,2,3)
contourf(unsteady_exp(:,:,2),20) %timestep 2
colorbar; colormap(jet);
title_text1=sprintf('Unsteady State-Explicit Method, iteration number = %d, CFL= %.1f',counter1(2),cfl(2));
title_text2=sprintf('Timestep size = %.10f sec, timestep = %d',dt(2),ts1(2));
title_text3=sprintf('at time = %.4f sec, temp. diff. = %.10f K',(dt(2)*ts1(2)),error1(2));
title({[title_text1],[title_text2],[title_text3]}) 
axis ij  %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')

subplot(2,2,4)
contourf(unsteady_exp(:,:,3)) %timestep 3
colorbar; colormap(jet);
title_text1=sprintf('Unsteady State-Explicit Method, iteration number = %d, CFL= %.1f',counter1(3),cfl(3));
title_text2=sprintf('Timestep size = %.10f sec, timestep = %d',dt(3),ts1(3));
title_text3=sprintf('at time = %.4f sec, temp. diff. = %.10f K',(dt(3)*ts1(3)),error1(3));
title({[title_text1],[title_text2],[title_text3]}) 
axis ij  %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')
 
figure(6) % SOR-GS implicit method
subplot(2,2,1) %initial
contourf(t_initial,20) % plotting initial conditions
colorbar; colormap(jet);
axis ij %to flip the Y-axis
axis equal
grid on
xlabel('X-axis')
ylabel('Y-axis')
title_text1=sprintf('Initial conditions, Thermal Diffusivity= %.1f',a);
title_text2=sprintf('Initial condition: Top 600K, Left 400K, Right 800K, Bottom 900K, Rest 320K');
title({[title_text1],[title_text2]}) 

subplot(2,2,2)
contourf(unsteady_sorgs(:,:,1),20) %timestep 1
colorbar; colormap(jet);
title_text1=sprintf('Unsteady-Implicit-SOR+GS, iteration number = %d, CFL= %.1f',counter4(1),cfl(1));
title_text2=sprintf('Timestep size = %.10f sec, timestep = %d',dt(1),ts4(1));
title_text3=sprintf('at time = %.4f sec, temp. diff. = %.10f K',(dt(1)*ts4(1)),error4(1));
title({[title_text1],[title_text2],[title_text3]}) 
axis ij  %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')

subplot(2,2,3)
contourf(unsteady_sorgs(:,:,2),20) %timestep 2
colorbar; colormap(jet);
title_text1=sprintf('Unsteady-Implicit-SOR+GS, iteration number = %d, CFL= %.1f',counter4(2),cfl(2));
title_text2=sprintf('Timestep size = %.10f sec, timestep = %d',dt(2),ts4(2));
title_text3=sprintf('at time = %.4f sec, temp. diff. = %.10f K',(dt(2)*ts4(2)),error4(2));
title({[title_text1],[title_text2],[title_text3]}) 
axis ij  %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')

subplot(2,2,4)
contourf(unsteady_sorgs(:,:,3),20) %timestep 3
colorbar; colormap(jet);
title_text1=sprintf('Unsteady-Implicit-SOR+GS, iteration number = %d, CFL= %.1f',counter4(3),cfl(3));
title_text2=sprintf('Timestep size = %.10f sec, timestep = %d',dt(3),ts4(3));
title_text3=sprintf('at time = %.4f sec, temp. diff. = %.10f K',(dt(3)*ts4(3)),error4(3));
title({[title_text1],[title_text2],[title_text3]}) 
axis ij  %to flip the Y-axis
axis square
grid on
xlabel('X-axis')
ylabel('Y-axis')

figure(7) % comparison of number of iteration and error in diffrent CFL number to solve un-steady state problem using explicit
X2 = categorical({num2str(dt(1)),num2str(dt(2)),num2str(dt(3))});
b=bar(X2,counter1,'g','g','r');
b.FaceColor='flat';
b.CData(3,:)=[1 0 0];
xlabel('Time Step sizes')
ylabel({['Number of iteration done till 10^6>error<10^-^6 or max Timestep=500'],['(whichever is achieved first)']})
title('Unsteady-State Explicit method, comparison of effect of Time Step')

figure(8) % comparison of number of iteration and error in diffrent CFL number to solve un-steady state problem using implicit SORGS
X2 = categorical({num2str(dt(1)),num2str(dt(2)),num2str(dt(3))});
b=bar(X2,counter4);
xlabel('Time Step sizes')
ylabel({['Number of iteration done till 10^6>error<10^-^6 or max Timestep=500'],['(whichever is achieved first)']})
legend('No.of iterations for stable solution')
title('Unsteady-State Implicit- SOR+GS method, comparison of effect of Time Step')

Unsteady-State Explicit Method:

This function solves the Unsteady-state problem using the Explicit method.

"t_un_ex" is initialized with a 3D matrix of zeros whose dimensions are "nx,ny,o", where "o" is the number of different CFL numbers

An outermost FOR loop is used to change the CFL values

"t_un_ex" is then initialized by copying initial temps, then a copy of "t_un_ex" is defined as "told_unex" to store the values at the previous iteration.

value of error is initialized as `9^9`.

The second most outer FOR loop is used to change the time-step of the solver, once it reaches maximum timestep defined in the main program or if it converges before that, it will stop the solver.

A Y-direction space loop is defined inside the X-direction space loop. New temperature values are calculated as per the equation derived above.

Error is calculated inside the Y-direction space loop as a maximum of difference in temperatures of all points between consecutive iterations.

"told_unex" is updated with "t_un_ex" values after the space loop so that the old values changes after each iteration.

A filled contour plot of "t_un_ex" is plotted within the while loop to visualize the convergence.

an IF loop is used to check for the convergence.

Function:

function [unsteady_exp,counter1,error1,ts1] = explicit_func_stability(t_initial,nx,ny,k1,k2,nt,dt,o)
    %unsteady state
    %trnsient state solution
    %explicit method
    % dT/dt = a*((d^2T/dx^2)+(d^2T/dy^2))
    % (tpnew - tpold)/dt = a*((tl+2*tp+tr)/dx^2 + (tt-2*tp+tb)/dy^2)
    t_un_ex=zeros(nx,ny,o); %initializing array to zero
    
    for c=1:o  %changing time step loop
         t_un_ex(:,:,c)=t_initial;
        %creating a copy of t
        told_unex=t_un_ex(:,:,c);
        %convergence loop
        counter_unex(c)=1; % loop counter
            for n=1:nt  %time loop
                error(c)=9e9; %error re-initialisation for every time step
                current_time = dt(c)*n;
                for i=2:nx-1 %spatial loop x direction
                    for j=2:ny-1 %spatial loop y direction           
                        term1 = k1(c)*(told_unex(i-1,j) + told_unex(i+1,j) - 2*told_unex(i,j));
                        term2 = k2(c)*(told_unex(i,j-1) + told_unex(i,j+1) - 2*told_unex(i,j));                
                        t_un_ex(i,j,c)= told_unex(i,j) + term1 + term2; % calculating new temperature
                        error(c)= max(max(abs(told_unex(i,j)-t_un_ex(i,j,c))));%calculating error                         
                    end       
                end                                       
                figure(c+1)
                contourf(t_un_ex(:,:,c))
                axis ij
                grid on;
                colorbar; colormap(jet);
                title_text1=sprintf('Unsteady State-Explicit Method');
                title_text2=sprintf('dt= %.4f Time Step = %d, current time = %.4f',dt(c),n,current_time);
                title_text3=sprintf('Temp. diff. b/w time steps = %.10f',error(c));
                title({[title_text1],[title_text2],[title_text3]});                            
                
                told_unex=t_un_ex(:,:,c); % updating old temperatures     
                counter_unex(c)=counter_unex(c)+1
                
                if(error(c)<1e-4) % error out of bounds checking break loop
                    break
                end
                if(error(c)>1e6)
                    break
                end                
            end            
            g(c)=n;  %no. of timesteps for each time step         
    end    
    unsteady_exp=t_un_ex;
    counter1=counter_unex;
    error1=error;
    ts1=g;
end

Unsteady-State Successive Over Relaxation on Gauss-Seidel Method:

This function solves the Unsteady-state problem using the Implicit method using SOR on the Gauss-Seidel iterative solver.

"t_un_ex" is initialized with a 3D matrix of zeros whose dimensions are "nx,ny,o", where "o" is the number of different CFL numbers

An outermost FOR loop is used to change the CFL values

"t_un_imp" is initialized by copying initial temps, then a copy of "t_un_imp" is defined as "told_unimp" to store the values at the previous iteration.

value of error is initialized as `9^9`.

An outer FOR loop is used to change the time-step of the solver, once it reaches maximum timestep defined in the main program or if it converges before that, it will stop the solver.

An outer while loop is used to check the tolerance criteria, if the value of error goes less than the tolerance specified in the main program, it will stop the space solver

A Y-direction space loop is defined inside the X-direction space loop. New temperature values are calculated as per the equation derived above.

Error is calculated inside the Y-direction space loop as a maximum of difference in temperatures of all points between consecutive iterations.

"told_unimp" is updated with "t_un_imp" values after the space loop so that the old values changes after each iteration.

an IF loop is used to check for the convergence.

Function:

function[unsteady_imp,counter4,error4,ts4]= sorgs_imp_stability(t_initial,nx,ny,k1,k2,k3,tol,nt,dt,sor_factor,o)
    %unsteady state
    %transient state solution
    %implicit method, SOR on G-S iterative solver
    % dT/dt = a*((d^2T/dx^2)+(d^2T/dy^2))
    t_un_imp=zeros(nx,ny,o);   %initializing array to zero
    
    for c=1:o  %changing time step loop
        %creating a copy of t
        t_un_imp(:,:,c)=t_initial;
        told_unimp=t_un_imp(:,:,c);
        tprev=t_un_imp(:,:,c);
        %convergence loop 
        count(c)=0; % loop counter
        for n=1:nt %time loop
            current_time = dt(c)*n;  
            error(c)=9e9; %error(c) re-initialisation for every time step
            while(error(c)>tol) %tolerance checking for every space loop
                for i=2:nx-1 %spatial loop x direction
                    for j=2:ny-1 %spatial loop y direction           
                        term1 = k1(c)*(t_un_imp(i-1,j) + told_unimp(i+1,j));
                        term2 = k2(c)*(t_un_imp(i,j-1) + told_unimp(i,j+1));
                        term3 = tprev(i,j);   %temperature at that point in previous timestep
                        t_un_imp(i,j,c)= (1-sor_factor)*told_unimp(i,j) + sor_factor*k3((c))*(term3 + term1 + term2); % calculating new temperature
                        error(c)= max(max(abs(told_unimp(i,j)-t_un_imp(i,j,c)))); %calculating error(c)                        
                    end       
                end   
            count(c)=count(c)+1;
            told_unimp=t_un_imp(:,:,c); % updating old temperatures  
            end
            tprev=t_un_imp(:,:,c); 
                if(error(c)<1e-7)% error out of bounds checking break loop
                    break
                end
                if(error(c)>1e6)
                    break
                end                
        end 
        l(c)=n; %no. of timesteps for each time step 
    end
    unsteady_imp=t_un_imp;
    counter4=count;
    error4=error;
    ts4=l;
end

Results:

Figure(5) shows the initial conditions and the final conditions for every CFL number, we can see that for CFL 0.3 or 0.5 the solution is stable and reaches a steady-state whereas for the CFL equal to 1, the solution is blown out and does not reaches convergence which proves the CFL stability criteria.

Screen capture of the unstable solution is shown below, it shows how the solution gets unstable and crosses the error limit set in the function:

Figure(6) shows a bar plot of numbers of iterations required for Explicit method solver to get a converged result or to become unstable is shown below:

Figure(7) shows the initial conditions and the final conditions for every CFL number, we can see that for every CFL number the solution is stable and reaches a steady-state and the required number of time steps to achieve convergence decreases as the size of time step increases.

Figure(8) shows a bar plot of numbers of iterations or time steps required for SOR on GS Implicit method solver to get a converged result or to become unstable is shown below:

Explanation of Results:

In an explicit numerical method, the rate of change equation would be evaluated in terms of known quantities at the previous time step `n`. An implicit method, in contrast, would evaluate some or all of the terms in the rate of change equation in terms of unknown quantities at the new time step `n+1`. Since new quantities appear on both the left and right sides of the quantity whose values we are calculating for time step `n+1`.

Numerical stability has to do with the behavior of the solution as the time-step dt is increased. If the solution remains well behaved for arbitrarily large values of the time step, the method is said to be unconditionally stable. This situation never occurs with explicit methods, which are always conditionally stable.

In an implicit formulation, a solution for the unknowns at the new time step n+1 may be obtained for any size time step

For an Explicit method to prevent the situation in which oscillations appear, a limitation on the size of the time step D t should be maintained below a certain limit. The limit depends on `Deltax` &`Deltay` and other parameters that form a stability criterion. A general condition for stability of explicit difference formulae in the two-dimensional system is the requirement that no coefficient in (1) is negative, which means the coefficient (1-2*K1-2*K2) is greater than or equal to zero. Otherwise, the oscillations in the solution occur.
This means that the stability criterion for a two-dimensional node is equation (2):

`K1 + K2<=0.5`

where,

`K1=alpha*(Deltat)/(Deltax^2)`

and

`K2=alpha*(Deltat)/(Deltay^2)`

Conclusion of Results:

1. From the above visual comparisons, we conclude that the Explicit method is very sensitive to the size of the time step and will be unstable if the inappropriate time step is chosen for the solution. Explicit methods are conditionally stable.
2. On the other hand, the Implicit method was stable for all CFL numbers which prove that Implicit methods are unconditionally stable.