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

Stability Analysis in Unsteady(Transient) State problem

MMohammed Mansoor Ali Khan22 Sept 2020

Stability Analysis in the unsteady(Transient) problem:

To study: Effect of time step dt on stability of the solution of 2D heat conduction equation for both implict and explicit schemes

In solving a unsteady state,the problem is that the solution changes with time that means the temp values change with time i.e `frac{delT}{delt}` and become stable after some particular time`

Transient 2D heat conduction equation

`frac{delT}{delt}=alpha[frac{del^2T}{delx^2}+frac{del^2T}{dely^2}]`

CFL (courant-Friedrichs-lewy) is an essential condition in solving the partial differential equation. CFL number is controlled by altering the space interval and time interval.

CFL for Transient problem CFL<=0.5 solution is stable.

CFL>0.5 solution is unstable only for explicit scheme

so in explicit scheme

`k1=alphafrac{Deltat}{Deltax^2} ;k2=alphafrac{Deltat}{Deltay^2}`

Transient problem:

we have two schemes to solve and check for stability for time steps:

  • Implicit scheme- in this the values are not knowns so we use iterative techniques to find out the values.

Jacobi method

Gauss siedel method

Successive over-relaxation method

  • Explicit scheme- In this we do not use any of the above methods as we know the values.

stable solution:

for dt=0.0008 and `alpha=1.4`the CFL number is 0.224 ;since its less tha 0.5 the solutuion is stable for 900 iterations.

code for explicit transient problem

clear all 
close all
clc
%inputs
n=10;
x=linspace(0,1,n);
y=linspace(0,1,n);
dx = 0.1;
dy = 0.1;
dt=0.0008;
alpha=1.4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
CFL=k1+k2;
T = ones(n);
nt=900;

%boundary conditions
T(1,1:n) = 600;%top bc
T(n,1:n) = 900; %bottom bc
T(1:n,1) = 400;%left bc
T(1:n,n) = 800;%right bc
Told=T;
iterative_solver = 1;
tol=1e-4;

tic;
    
 for v=1:nt
     for i = 2:n-1
        for j = 2:n-1
         
     term1=Told(i,j);
     term2=k1*(Told(i-1,j) - 2*Told(i,j)+ Told(i+1,j));
     term3=k2*(Told(i,j-1) - 2*Told(i,j)+ Told(i,j+1));
     T(i,j)=(term1+term2+term3);
    end
  end
   Told=T;
    end
    
 c= toc;
  
figure(1)
[P,Q]=contourf(x,y,T);
clabel(P,Q)
colormap('jet');
colorbar;
axis ij;
xlabel('x axis')
ylabel('y axis')
title_text1=sprintf('Transient explicit');
title_text2=sprintf('iteration number=%d dt = %s ',nt,dt);
title_text3=sprintf('CFL number=%s',CFL);
title({title_text1,title_text2,title_text3});

Output plot

Unstable solution :

The solution becomes unstable for the explicit transient problem after certain time steps:

parameters: for `dt=0.002`,`alpha=1.4` the CFL comes out to be `CFL=0.56` the solution becomes unstable.

code :

clear all 
close all
clc
%inputs
n=10;
x=linspace(0,1,n);
y=linspace(0,1,n);
dx = 0.1;
dy = 0.1;
dt=0.002;
alpha=1.4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
CFL=k1+k2;
T = ones(n);
nt=900;

%boundary conditions
T(1,1:n) = 600;%top bc
T(n,1:n) = 900; %bottom bc
T(1:n,1) = 400;%left bc
T(1:n,n) = 800;%right bc
Told=T;
iterative_solver = 1;
tol=1e-4;

tic;
    
 for v=1:nt
     for i = 2:n-1
        for j = 2:n-1
         
     term1=Told(i,j);
     term2=k1*(Told(i-1,j) - 2*Told(i,j)+ Told(i+1,j));
     term3=k2*(Told(i,j-1) - 2*Told(i,j)+ Told(i,j+1));
     T(i,j)=(term1+term2+term3);
    end
  end
   Told=T;
    end
    
 c= toc;
  
figure(1)
[P,Q]=contourf(x,y,T);
clabel(P,Q)
colormap('jet');
colorbar;
axis ij;
xlabel('x axis')
ylabel('y axis')
title_text1=sprintf('Transient explicit');
title_text2=sprintf('iteration number=%d dt = %s ',nt,dt);
title_text3=sprintf('CFL number=%s',CFL);
title({title_text1,title_text2,title_text3});

Plot for Unstable solution:

The plot shows that for time step of 0.002 , alpha=1.4 and CFL of 0.56 the solution becomes unstable

Implicit method:

these methods are exceptionally stable than that of explicit . even for step size of 0.3 and alpha of 1.4 these have a stable solution

Stability in Jacobi:

step size dt=0.02,alpha=1.4 the CFL comes out to be CFL=0.56 for nt=900 iteration.

code:

clear all 
close all
clc
%inputs
n=10;
x=linspace(0,1,n);
y=linspace(0,1,n);
dx = 0.1;
dy = 0.1;
dt=0.02;
alpha=1.4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=(1/(2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
CFL=k1+k2;
T = ones(n);
nt=900;

%boundary conditions
T(1,1:n) = 600;%top bc
T(n,1:n) = 900; %bottom bc
T(1:n,1) = 400;%left bc
T(1:n,n) = 800;%right bc
Told=T;
T_prev=Told;
tol=1e-4;

tic;
    
 for v=1:nt
     for i = 2:n-1
        for j = 2:n-1
         
     H=(Told(i-1,j)+Told(i+1,j));
     V=(Told(i,j-1)+Told(i,j+1));
     T(i,j)=(term1*T_prev(i,j)) + (term2*H) + (term3*V);
    end
  end
   Told=T;
    end
    
 c= toc;
  
figure(1)
[P,Q]=contourf(x,y,T);
clabel(P,Q)
colormap('jet');
colorbar;
axis ij;
xlabel('x axis')
ylabel('y axis')
title_text1=sprintf('Transient implicit jacobi');
title_text2=sprintf('iteration number=%d dt = %s ',nt,dt);
title_text3=sprintf('CFL number=%s',CFL);
title({title_text1,title_text2,title_text3});
pause(0.3)

Stability in Gauss siedel:

step size dt=0.02,alpha=1.4 the CFL comes out to be CFL=0.56 for nt=900 iteration.

code:

clear all 
close all
clc
%inputs
n=10;
x=linspace(0,1,n);
y=linspace(0,1,n);
dx = 0.1;
dy = 0.1;
dt=0.02;
alpha=1.4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=(1/(2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
CFL=k1+k2;
T = ones(n);
nt=900;

%boundary conditions
T(1,1:n) = 600;%top bc
T(n,1:n) = 900; %bottom bc
T(1:n,1) = 400;%left bc
T(1:n,n) = 800;%right bc
Told=T;
T_prev=Told;
tol=1e-4;

tic;
    
 for v=1:nt
     for i = 2:n-1
        for j = 2:n-1
         
     H=(T(i-1,j)+Told(i+1,j));
     V=(T(i,j-1)+Told(i,j+1));
     T(i,j)=(term1*T_prev(i,j)) + (term2*H) + (term3*V);
    end
  end
   Told=T;
    end
    
 c= toc;
  
figure(1)
[P,Q]=contourf(x,y,T);
clabel(P,Q)
colormap('jet');
colorbar;
axis ij;
xlabel('x axis')
ylabel('y axis')
title_text1=sprintf('Transient implicit gs');
title_text2=sprintf('iteration number=%d dt = %s ',nt,dt);
title_text3=sprintf('CFL number=%s',CFL);
title({title_text1,title_text2,title_text3});
pause(0.3)

Stability in SOR:

step size dt=0.02,alpha=1.4 the CFL comes out to be CFL=0.56 for nt=900 iteration.

code:

clear all 
close all
clc
%inputs
n=10;
w=1.3;
x=linspace(0,1,n);
y=linspace(0,1,n);
dx = 0.1;
dy = 0.1;
dt=0.3;
alpha=1.4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=(1/(2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
CFL=k1+k2;
T = ones(n);
nt=900;

%boundary conditions
T(1,1:n) = 600;%top bc
T(n,1:n) = 900; %bottom bc
T(1:n,1) = 400;%left bc
T(1:n,n) = 800;%right bc
Told=T;
T_prev=Told;
tol=1e-4;

tic;
    
 for v=1:nt
     for i = 2:n-1
        for j = 2:n-1
         
     H=(T(i-1,j)+Told(i+1,j));
         V=(T(i,j-1)+Told(i,j+1));
         T(i,j)=(1-w)*Told(i,j) + w*((term1*T_prev(i,j)) + (term2*H) + (term3*V));
       end
  end
   Told=T;
    end
    
 c= toc;
  
figure(1)
[P,Q]=contourf(x,y,T);
clabel(P,Q)
colormap('jet');
colorbar;
axis ij;
xlabel('x axis')
ylabel('y axis')
title_text1=sprintf('Transient implicit SOR');
title_text2=sprintf('iteration number=%d dt = %s ',nt,dt);
title_text3=sprintf('CFL number=%s',CFL);
title({title_text1,title_text2,title_text3});
pause(0.3)

output plot for jaccobi,Gauss siedel and SOR for nt=900,alpha=1.4:

code for different implicit scheme for nt=900,dt=0.1; the CFL=28

jacobi code:

clear all 
close all
clc
%inputs
n=10;
x=linspace(0,1,n);
y=linspace(0,1,n);
dx = 0.1;
dy = 0.1;
dt=0.1;
alpha=1.4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=(1/(2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
CFL=k1+k2;
T = ones(n);
nt=900;
%boundary conditions
T(1,1:n) = 600;%top bc
T(n,1:n) = 900; %bottom bc
T(1:n,1) = 400;%left bc
T(1:n,n) = 800;%right bc
Told=T;
T_prev=Told;
tol=1e-4;
tic;
  
 for v=1:nt
     for i = 2:n-1
        for j = 2:n-1
         
     H=(Told(i-1,j)+Told(i+1,j));
     V=(Told(i,j-1)+Told(i,j+1));
     T(i,j)=(term1*T_prev(i,j)) + (term2*H) + (term3*V);
    end
  end
   Told=T;
    end
c= toc;
 figure(1)
[P,Q]=contourf(x,y,T);
clabel(P,Q)
colormap('jet');
colorbar;
axis ij;
xlabel('x axis')
ylabel('y axis')
title_text1=sprintf('Transient implicit jacobi');
title_text2=sprintf('iteration number=%d dt = %s ',nt,dt);
title_text3=sprintf('CFL number=%s',CFL);
title({title_text1,title_text2,title_text3});

code for Gauss sidel:

clear all 
close all
clc
%inputs
n=10;
x=linspace(0,1,n);
y=linspace(0,1,n);
dx = 0.1;
dy = 0.1;
dt=0.1;
alpha=1.4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=(1/(2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
CFL=k1+k2;
T = ones(n);
nt=900;

%boundary conditions
T(1,1:n) = 600;%top bc
T(n,1:n) = 900; %bottom bc
T(1:n,1) = 400;%left bc
T(1:n,n) = 800;%right bc
Told=T;
T_prev=Told;
tol=1e-4;

tic;
    
 for v=1:nt
     for i = 2:n-1
        for j = 2:n-1
         
     H=(T(i-1,j)+Told(i+1,j));
     V=(T(i,j-1)+Told(i,j+1));
     T(i,j)=(term1*T_prev(i,j)) + (term2*H) + (term3*V);
    end
  end
   Told=T;
    end
    
 c= toc;
  
figure(1)
[P,Q]=contourf(x,y,T);
clabel(P,Q)
colormap('jet');
colorbar;
axis ij;
xlabel('x axis')
ylabel('y axis')
title_text1=sprintf('Transient implicit gs');
title_text2=sprintf('iteration number=%d dt = %s ',nt,dt);
title_text3=sprintf('CFL number=%s',CFL);
title({title_text1,title_text2,title_text3});

code for SOR:

clear all 
close all
clc
%inputs
n=10;
w=1.3;
x=linspace(0,1,n);
y=linspace(0,1,n);
dx = 0.1;
dy = 0.1;
dt=0.1;
alpha=1.4;
k1=(alpha*dt)/dx^2;
k2=(alpha*dt)/dy^2;
term1=(1/(2*k1+2*k2));
term2=k1*term1;
term3=k2*term1;
CFL=k1+k2;
T = ones(n);
nt=900;

%boundary conditions
T(1,1:n) = 600;%top bc
T(n,1:n) = 900; %bottom bc
T(1:n,1) = 400;%left bc
T(1:n,n) = 800;%right bc
Told=T;
T_prev=Told;
tol=1e-4;

tic;
    
 for v=1:nt
     for i = 2:n-1
        for j = 2:n-1
         
     H=(T(i-1,j)+Told(i+1,j));
         V=(T(i,j-1)+Told(i,j+1));
         T(i,j)=(1-w)*Told(i,j) + w*((term1*T_prev(i,j)) + (term2*H) + (term3*V));
       end
  end
   Told=T;
    end
    
 c= toc;
  
figure(1)
[P,Q]=contourf(x,y,T);
clabel(P,Q)
colormap('jet');
colorbar;
axis ij;
xlabel('x axis')
ylabel('y axis')
title_text1=sprintf('Transient implicit SOR');
title_text2=sprintf('iteration number=%d dt = %s ',nt,dt);
title_text3=sprintf('CFL number=%s',CFL);
title({title_text1,title_text2,title_text3});

output plot for nt=900,dt=0.1; the CFL=28:

Explanation for code:

  • declare nodes n,x axis and yaxis,dx,dy,alpha,omega in case of SOR.
  • write the equation for k1,k2 and CFL
  • declare the boundary conditions
  • set Told=T and Tprev =Told.
  • initiate a for loop for nt ,another two for i,j
  • write the descritized eqaution for individual methods.
  • close the for loop and increment the counter.
  • set again Told=T
  • plot the coutour for each indivual schemes and methods.

Explanation for results/plots:

The solution become unstable for explicit at dt=0.002 ,alpha=1.4 as CFL is 0.56.

As the CFL>0.5 the solution becomes unstable for transient explicit scheme.

Implicit schemes are exceptionally stable for larger time steps.

for dt=0.02,alpha=1.4, CFL=0.56 the solution is stable for implicit scheme in transient state.

Even for a larger time step dt=0.1, alpha =1.4 ,CFL=28 the solution is exceptionally stable for implicit scheme in transient state.

conclusion:

So CFL number plays a vital role in determining the stability of solution depending on the methods and schemes used.

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.