Stability analysis in unsteady state heat conduction equation solve using implicit and explicit method

Stability analysis in unsteady state heat conduction equation solve using implicit and explicit method.

CFL Number - Courant Friedrich Lewy,it is a dimensionless number.It must be less than 1 to maintain the solution stability.

A physical explanation, for a CFD simulation, of the Courant number could be that it tells you something about how fluid is moving through your computational cells. Meaning that If the Courant number is <=1 then the fluid particles move from one cell to another within one time step (at most). While if the Courant number is >1 a fluid particles moves through two or more cells at each time step and this can affect convergence negatively.

CFL Number is based on the thermal diffusivity,time step and cell grid size.

Stability analysis for unsteady state heat conduction -implicit method:

CFL number is mainly based on timestep(dt) and cell grid size(dx=dy).

Obersvation:

From this stability analysis,What i have observed is with the use of any timestep the implicit method gives us the converged result always.

Stability analysis for unsteady state heat conduction - Explicit method:

Obesrvation:

From this stability analysis,i have observed that With the incremental change of timestep will make my result unstable.

And CFL no <= 0.2025 ;dt =0.001;themal diffusivity(alpha) =1.25 at this point my result got unstable.

Matlab code for unsteady state heat conduction equation solving explicitly:

%Unsteady state 2D heat conduction equation using explicit method
%================================================================


clear all
close all
clc


%length 
L= 1;

%Defining square domain
nx=10;
ny=nx;

n=100;
%Grids and nodes
%---------------
%Grids & nodes in X-axis 
x=linspace(0,L,nx);
dx=x(2)-x(1);

%Grids & nodes in Y-axis
y=linspace(0,L,nx);
dy=y(2)-y(1);


%nx X ny Matrix with boundary condition
T=ones(nx,ny);
T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;


dt = 0.001;

%Previous temperature
Told = T;



%Assumed thermal diffusivity
alpha =1.25;


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

CFL = alpha *(1/dx^2 + 1/dy^2)*dt;

[X, Y] = meshgrid(x,y);

k = 1;
%%Time loop
tic;
for nt = 1:n  
            for i=2:nx-1
                for j=2:ny-1
                    term1 = Told(i,j);
                    term2 = k1*(T(i-1,j)-2*Told(i,j)+Told(i+1,j));
                    term3 = k2*(T(i,j-1)-2*Told(i,j)+Told(i,j+1));
                    
                    T(i,j) = term1+term2+term3;
                end
            end
                   
                     Told=T;
                    
                     colormap(jet);
                     [C,h]=contourf(x,y,T);
                     clabel(C,h);
                     colorbar
                     title(sprintf('No.of iteration in excplicit method is %d',k))
                     k=k+1;
                     xlabel('x axis');
                     ylabel('Spatial axis');
                     pause(0.03)
                     
   toc;                  
   fprintf('Time taken for solving is %f\n',toc)          
      end
    
                  

Conclusion:

In implicit method at any time step we will get a converged solution and it will maintain the stability.

In explicit method, if the timestep increases the solution will become unstable.

And the CFL Number plays a important role to determine the correct timestep to maintain the stability.