Stability analysis in the unsteady problem
Objective: To observe the effect of time step on transient state methods.
Description:
`CFL=α(Δt)/(Δx^2)`
i.e.. `CFL∝Δt`
If Δt increases then CFL increases
Case1 - Explicit Method:
close all
clear all
clc
nx = 20; %nx & ny will be nodes along the length of x & y respectively
ny = nx; %let us consider nx = ny
x = linspace(0,1,nx); %to create nx number of points from origin to length of domain =1
y = linspace(0,1,ny); %to create ny number of points from origin to length of domain =1
dx=(x(2)-x(1));
dy=(y(2)-y(1));
thermal_diffusivity = 1.4;
cfl = [0.4:0.001:0.7];
time=5;
for a=1:length(cfl)
dt = (cfl(a)*dx^2)/(2*thermal_diffusivity);
k1 = thermal_diffusivity*dt/dx^2;
k2 = thermal_diffusivity*dt/dy^2;
nt = time/dt;
%let us consider the boundary conditions
top=600;
bottom=300;
left=400;
right=800;
%let us start by taking all values as ones - nx rows and ny columns
Tfirst = ones(nx,ny);
%assiging those boundary conditions to the T matrix excepting the corners
Tfirst(1,2:nx-1) = top; %to assign the values for first row of T matrix
Tfirst(ny,2:nx-1) = bottom; %to assign the values for last row of T matrix
Tfirst(2:ny-1,1) = left; %to assign the values for first column of T matrix
Tfirst(2:ny-1,nx) = right; %to assign the values for last column of T matrix
%assigning the values at the corner positions of the T matrix by taking average of adjacent values
Tfirst(1,1)=(top+left)/2;
Tfirst(ny,1)=(bottom+left)/2;
Tfirst(1,nx)=(right+top)/2;
Tfirst(ny,nx)=(right+bottom)/2;
%dummy variables are created for updating the values
Tprev = Tfirst;
Tnew = Tfirst;
%every method should be strted with iteration number 1
plane_iter = 1;
term1 = (1-(2*k1)-(2*k2));
tolerance=1e-7;
%transient state - explicit scheme===================================
for t=1:nt
error = 9e9;%time loop
for i = 2:nx-1 %space loop
for j = 2:ny-1
H1=k1*(Tprev(i+1,j)+Tprev(i-1,j));
V1=k2*(Tprev(i,j+1)+Tprev(i,j-1));
Tnew(i,j) = (term1*Tprev(i,j))+(H1)+(V1);
end
end
error = max(max(abs(Tprev-Tnew)));
Tprev = Tnew;
plane_iter = 1+plane_iter;
if error==0; % to stop the loop after reaching convergence
break
end
end
% fprintf(' nComputational time using Explicit method = %g n',a);
% fprintf('No. of time levels approached in explicit method = %d n',plane_iter);
figure(1)
contourf(Tnew); %to get output as contour
colormap(jet); %color of contour to get rainbow colors
colorbar; %side bar representing the temperature with color
title_text = sprintf('CFL number = %g & No. of time levels = %g',cfl(a),plane_iter ); %to show the time loop number
title(title_text) %title of plot
set(gca, 'YDir','reverse')
M(a)=getframe(gcf); %to fetch frame for each plot at a
end
%to open and create a video using video writer
video=VideoWriter('CFL analysis','MPEG-4');
video.FrameRate=2; %to set the frame rate for video
video.Quality=100; %to get higher quality video
open(video); %to open the video
writeVideo(video,M); %to attach fetched frames together
close(video); % to close the process
Output:
Result:
From the above video, we can observe that from CFL number upto 0.504 the solution is showing similar pattern which is stable and from 0.504 we can see some random patterns which shows its instability.
Justification:
As we know that the selection of dt value in explicit method is based on CFL number. We must make sure that the CFL number should be less than or equal to 0.5 (for diffusion term) to get stable solution if it is more than 0.5 the solution will blow up.
Case2 - Implicit Method:
close all
clear all
clc
nx = 20; %nx & ny will be nodes along the length of x & y respectively
ny = nx; %let us consider nx = ny
x = linspace(0,1,nx); %to create nx number of points from origin to length of domain =1
y = linspace(0,1,ny); %to create ny number of points from origin to length of domain =1
dx=(x(2)-x(1));
dy=(y(2)-y(1));
thermal_diffusivity = 1.4;
cfl = [0.4:0.01:0.7];
time=5;
for a=1:length(cfl)
dt = (cfl(a)*dx^2)/(2*thermal_diffusivity);
k1 = thermal_diffusivity*dt/dx^2;
k2 = thermal_diffusivity*dt/dy^2;
nt = time/dt;
%let us consider the boundary conditions
top=600;
bottom=300;
left=400;
right=800;
%let us start by taking all values as ones - nx rows and ny columns
T = ones(nx,ny);
%assiging those boundary conditions to the T matrix excepting the corners
T(1,2:nx-1) = top; %to assign the values for first row of T matrix
T(ny,2:nx-1) = bottom; %to assign the values for last row of T matrix
T(2:ny-1,1) = left; %to assign the values for first column of T matrix
T(2:ny-1,nx) = right; %to assign the values for last column of T matrix
%assigning the values at the corner positions of the T matrix by taking average of adjacent values
T(1,1)=(top+left)/2;
T(ny,1)=(bottom+left)/2;
T(1,nx)=(right+top)/2;
T(ny,nx)=(right+bottom)/2;
term1 = 1/(1+2*k1+2*k2);
term2 = k1* term1;
term3 = k2* term1;
%dummy variables are created for updating the values of next time step
Tprev1=T; %for jacobi
Tprev2=T; %for Gauss-seidel
Tprev3=T; %for SOR
%dummy variables are created for updating the values
%T & Told are used for jacobi method
%Told1 & Told2 are used for Gauss-seidel method
%Told3 & Told4 are used for SOR method
Told=T;
Told1=T;
Told2=T;
Told3=T;
Told4=T;
tolerance = 1e-7;
%every method should be strted with iteration number 1
j_iter=1;
GS_iter=1;
SOR_iter=1;
%every method should be strted with 1st time loop
j_time=1;
GS_time=1;
SOR_time=1;
%for SOR method
alpha=1.2; %alpha should be between 1 & 2
k=2* (dx^2 + dy^2)/(dx^2 * dy^2);
%jacobi method - transient state - implicit scheme=====================================
for t=1:nt %time loop
error_jacobi = 9e9;
while (error_jacobi > tolerance)
for i = 2:nx-1 %space loop
for j = 2:ny-1
H1 = Told(i+1,j)+Told(i-1,j);
V1 = Told(i,j+1)+Told(i,j-1);
T(i,j) = (term1*Tprev1(i,j))+(term2*H1)+(term3*V1);
end
end
error_jacobi= max(max(abs(Told-T)));
Told= T;
j_iter = 1+j_iter;
end
Tprev1=T;
j_time=1+j_time;
if (error_jacobi==0)
break
end
end
% fprintf(' nComputational time using implicit Jacobi method = %g n',a);
% fprintf('No. of iterations in implicit Jacobi method = %d n',j_iter);
% fprintf('No. of time levels approached in implicit Jacobi method = %d n',j_time);
subplot(3,1,1)
contourf(T); %to get output as contour
colormap(jet); %color of contour to get rainbow colors
colorbar; %side bar representing the temperature with color
title_text = sprintf('CFL number=%g, No. of time levels(Jacobi)=%g',cfl(a),j_time); %to show the iteration number
title(title_text) %title of plot
set(gca, 'YDir','reverse')
M(a)=getframe(gcf); %to fetch frame for each plot at a
%Gauss_seidel method - transient state - implicit scheme===============================
for t=1:nt %time loop
error_gauss = 9e9;
while (error_gauss > tolerance)
for i = 2:nx-1 %space loop
for j = 2:ny-1
H2 = Told2(i+1,j)+Told1(i-1,j);
V2 = Told2(i,j+1)+Told1(i,j-1);
Told1(i,j) = (term1*Tprev2(i,j))+(term2*H2)+(term3*V2);
end
end
error_gauss= max(max(abs(Told2-Told1)));
Told2= Told1;
GS_iter = 1+GS_iter;
end
Tprev2=Told2;
GS_time=1+GS_time;
if (error_gauss==0)
break
end
end
% fprintf(' nComputational time using implicit GS method = %g n',b);
% fprintf('No. of iterations in implicit GS method = %d n',GS_iter);
% fprintf('No. of time levels approached in implicit GS method = %d n',GS_time);
subplot(3,1,2)
contourf(Told1); %to get output as contour
colormap(jet); %color of contour to get rainbow colors
colorbar; %side bar representing the temperature with color
title_text = sprintf('CFL number=%g, No. of time levels(Gauss-seidel)=%g',cfl(a),GS_time); %to show the iteration number
title(title_text) %title of plot
set(gca, 'YDir','reverse')
%SOR method - transient state - implicit scheme===============================
for t=1:nt %time loop
error_SOR = 9e9;
while (error_SOR > tolerance)
for i = 2:nx-1 %space loop
for j = 2:ny-1
H3 = Told4(i+1,j)+Told3(i-1,j);
V3 = Told4(i,j+1)+Told3(i,j-1);
Told3(i,j) = (1-alpha)*Told4(i,j)+(alpha)*((term1*Tprev3(i,j))+(term2*H3)+(term3*V3));
end
end
error_SOR= max(max(abs(Told4-Told3)));
Told4=Told3;
SOR_iter = 1+SOR_iter;
end
Tprev3=Told4;
SOR_time=1+SOR_time;
if (error_SOR<1e-10)
break
end
end
% fprintf(' nComputational time using implicit SOR method = %g n',c);
% fprintf('No. of iterations in implicit SOR method = %d n',SOR_iter);
% fprintf('No. of time levels approached in implicit SOR method = %d n',SOR_time);
subplot(3,1,3)
contourf(Told3); %to get output as contour
colormap(jet); %color of contour to get rainbow colors
colorbar; %side bar representing the temperature with color
title_text = sprintf('CFL number=%g, No. of time levels(SOR)=%g',cfl(a),SOR_time); %to show the iteration number
title(title_text) %title of plot
set(gca, 'YDir','reverse')
end
%to open and create a video using video writer for Jacobi method
video=VideoWriter('CFL analysis-1','uncompressed AVI');
video.FrameRate=2; %to set the frame rate for video
open(video); %to open the video
writeVideo(video,M); %to attach fetched frames together
close(video); % to close the process
Output:
Result:
From the above video, we can observe that the stability of the solution is not changing along with CFL number.
Justification:
As we know that the selection of dt value in implicit method doesnt depend on CFL number. So for larger values of dt also the solution will be unconditionally stable.
=======================================================================
link for Matlab code of contours in the above video
https://skill-lync.com/projects/solving-the-steady-and-unsteady-2d-heat-conduction-problem-62
link for Comparision of above steady and transient state results along with detailed description about steady and transient states and methods of solving steady(itterative solvers) and transient states(explicit and implicit)
https://skill-lync.com/projects/comparison-of-steady-state-and-transient-state-results