Stability analysis in the unsteady problem
The objective of this challenge is to study the stability criteria for the transient case, i.e Explicit and Implicit functions. The stability depends on the CFL number. Here, we can analyse the stability for various CFL number to see how the solution changes/appears.
What is a CFL number?
When timesteps are involved in the problem, the solution might become unstable for a smaller time set. To overcome this instability in the solution Courant-Friedrichs-Levy (CFL) is introduced.
If velocity is involved,
CFL=velocity⋅ΔtΔxCFL=velocity⋅ΔtΔx < 1
If diffusivity is involved, like here,
CFL=α⋅ΔtΔx2CFL=α⋅ΔtΔx2 < 0.5
Here k1 is the CFL along x
k2 is the CFL along y
CFL_x + CFL_y < 0.5
CODE FOR EXPLICIT FUNCTION:
clear all
close all
clc
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = dx;
tol = 1e-4;
dt = 1e-3;
t = 300*ones(nx,ny);
t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;
t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;
told = t;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
tol = 1e-4;
iter = 1;
if iter == 1
for k = 1:1400
error = 1;
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
A = told(i,j);
B = k1*(told(i+1,j) -2*(told(i,j)) + told(i-1,j));
C = k2*(told(i,j+1) -2*(told(i,j)) + told(i,j-1));
t(i,j) = A + B + C;
end
end
error = max(max(abs(told-t)));
told = t;
iter = iter + 1;
figure(1)
contourf(x,y,t);
clabel(contourf(x,y,t));
set(gca, 'ydir', 'reverse')
colormap(jet)
colorbar
title_text = sprintf('No of iter = %d' , iter);
title(title_text)
pause(0.03)
end
end
end
AT CFL <0.5:
Here the dt = 1e-3; CFL = 0.089; Alpha = 1.1;
AT CFL>0.5
Here the dt = 1e-2; CFL along x & y = 0.89; Alpha = 1.1;
This is an unstable solution.
Here the dt = 1e-1; CFL along x & y = 8.910; Alpha = 1.1;
CODE FOR IMPLICIT JACOBI:
clear all
close all
clc
nx = 10;
ny = 10;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
t = 300*ones(nx,ny);
t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;
t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;
told = t;
t_prev = t;
dt = 1e-3;
tol = 1e-4;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
jacobi_iter = 1;
for k = 1:1400
error = 9e9;
while(error > tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1/(1+ (2*k1) + (2*k2));
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) = (t_prev(i,j)*term1) + (H*term2) + (V*term3);
end
end
error = max(max(abs(told-t)));
told = t;
jacobi_iter = jacobi_iter + 1;
end
t_prev = t;
figure(1)
contourf(x,y,t)
clabel(contourf(x,y,t))
colorbar
set(gca, 'ydir', 'reverse')
title_text = sprintf('Jacobi iterations = %d' , jacobi_iter);
title(title_text)
end
AT CFL <0.5
Here the dt = 1e-3; CFL along x & y = 0.089; Alpha = 1.1;
AT CFL>0.5
Here the dt = 1e-1; CFL along x & y = 8.910; Alpha = 1.1;
CODE FOR IMPLICIT GAUSS SEIDEL:
clear all
close all
clc
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = dx;
t = 300*ones(nx,ny);
t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;
t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;
told = t;
t_prev = t;
dt = 1e-3;
k1 = 1.1*(dt/(dx^2));
k2 = 1.1*(dt/(dy^2));
error = 9e9;
tol = 1e-4;
GS_iter = 1;
for k = 1:1400
error = 1;
while(error > tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1/(1+ (2*k1) + (2*k2));
term2 = k1*term1;
term3 = k2*term1;
H = t(i-1,j) + told(i+1,j);
V = t(i,j+1) + told(i,j-1);
t(i,j) = (t_prev(i,j)*term1) + (H*term2) + (V*term3);
end
end
error = max(max(abs(told-t)));
told = t;
GS_iter = GS_iter + 1;
end
t_prev = t;
figure(1)
contourf(x,y,t);
set(gca, 'ydir', 'reverse')
clabel(contourf(x,y,t));
colorbar
title_text = sprintf('Gauss seidal iterations = %d' , GS_iter);
title(title_text)
end
AT CFL <0.5
Here the dt = 1e-3; CFL along x & y = 0.089; Alpha = 1.1;
AT CFL>0.5
Here the dt = 1e-1; CFL along x & y = 8.910; Alpha = 1.1;
CODE FOR IMPLICIT SOR:
clear all
close all
clc
nx = 10;
ny = nx;
x = linspace(0,1,nx);
y = linspace(0,1,ny);
dx = x(2) - x(1);
dy = dx;
t = 300*ones(nx,ny);
t(2:ny-1,1) = 400;
t(2:ny-1,nx) = 800;
t(1,2:nx-1) = 600;
t(ny,2:nx-1) = 900;
t(1,1) = 500;
t(1, ny) = 700;
t(nx,1) = 650;
t(nx, ny) = 850;
told = t;
t_prev = t;
t_sor = t;
dt = 1e-3;
k1 = 1.1*(dt)/(dx^2);
k2 = 1.1*(dt)/(dy^2);
omega = 1.1;
tol = 1e-4;
SOR_iter = 0;
for k = 1:1400
error = 1;
while(error > tol)
for i = 2:nx-1
for j = 2:ny-1
term1 = 1/(1+ (2*k1) + (2*k2));
term2 = k1*term1;
term3 = k2*term1;
H = t(i-1,j) + told(i+1,j);
V = t(i,j+1) + told(i,j-1);
t(i,j) = (t_prev(i,j)*term1) + (H*term2) + (V*term3);
t_sor(i,j) = ((1-omega)*told(i,j)) + (t(i,j)*omega);
end
end
error = max(max(abs(told-t)));
told = t_sor;
SOR_iter = SOR_iter + 1;
end
t_prev = t_sor;
figure(1)
contourf(x,y,t);
set(gca, 'ydir', 'reverse')
clabel(contourf(x,y,t));
colorbar
title_text = sprintf('SOR iterations = %d' , SOR_iter);
title(title_text)
end
AT CFL <0.5
Here the dt = 1e-3; CFL along x & y = 0.089; Alpha = 1.1;
AT CFL>0.5
Here the dt = 1e-1; CFL along x & y= 8.910; Alpha = 1.1;
CONCLUSIONS:
- The explicit function strictly obeys the CFL criteria which are to be less than 0.5.
- In order to keep the solution stable large time step is preferred in an explicit method.
- The implicit function doesn't care about the CFL criteria it looks stable for any given timestep