Stability analysis in the unsteady problem
Aim:- The aim of this project is to observe stability analysis in the unsteady problem.
Theory:-The heat conduction equation is a partial differential equation that describes the distribution of heat (or the temperature field) in a given body over time.
Assumption & Boundary conditions:-Following assumption are made to solve the steady and unsteady equation which as follows:-
1. Assuming that the domain is a unit square.
2. Assuming that nx = ny [Number of points along the x direction is equal to the number of points along the y direction]
3. Boundary conditions for steady and unsteady case
- Left:400K
- Top:600K
- Right:800K
- Bottom:900K
4. Tolerance=1e-4
Before proceeding further lets see some of the following important points which are as follows:-
1.How does solution become unstable?
Actually,solution becomes unstable when there is error in it.The error occur in the solution is because of following reasons.
- When we discretize the equation we approximate the term,and in the approximation we ignore the terms which is called as truncation error.
2.What does CFL number means and its criteria for implicit and explicit scheme?
Definition:-Courant–Friedrichs–Lewy (CFL) condition is a necessary condition for convergence while solving certain partial differential equation (usually hyperbolic PDE) numerically. It arises in the numerical analysis of explicit time integration schemes, when these are used for the numerical solution. As a consequence, the time step must be less than a certain time in many explicit time-marching simulation, otherwise the simulation produces incorrect results.
- Following is the formula for CFl number for 2D case is,
C=`alpha*(frac{1}{dx^2}+frac{1}{dy^2})dt`
here, `alpha`is thermal diffusivity
dx=dy=number of grids in x and y direction
dt=Time-step
CFl number criteria:-
- The condition for CFL number is,
`Cle0.25`
Explicit scheme:-If CFL number is more then the 0.25 for the solution become unstable.
Implicit scheme:-The implicit scheme is inherently stable which means if CFL number is more then its condition value it can be tolerated.
Now, to understand above concept in detail the we have solve 2D heat conduction equation for unsteady state:-
Govering equation for 2D heat conduction equation for unsteady state is
....................(2)
We will solve above equation by using two scheme,
- Explict scheme
- Implicit scheme
a.Explict scheme :-
After applying first order forward differencing scheme on L.H.S and central differencing scheme on R.H.S on eqaution(2),we get
Let says,
And,on simplification we get ...(3)
Where k1 & k2 is,
The equation (2) will be solve by using point iterative method.The point iterative method which will be use to solve the equation(2) are as follows:-
- Jacobi
- Gauss seidal
- Successive over-relaxation
Matlab code:-
- In the following program we have solve eqaution(3) by using diiferent iterative method which is mentioned above.
- First,we have assume that number of grid in x direction is equal to number of grid in y direction
- And after that we have initialize the array and assigining the boundary condition.
- Then time loop is set in in which equation (3) will be solve by using different iterative method until the time loop get finish and convergence result is obtain.
- After,then result is plotted using contourf command.
clear all
close all
clc
%Solving transient 2D heat conduction equation
%Input parameter
nx=40; %no of grids in x-direction
ny=40; %no of grids in y_direction
rf=1.5; %relaxation factor
dt=change;
alpha=1.4;
nt=400;
%Calculation
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1); %Time step in x-direction
dy=y(2)-y(1); %Time step in y-direction
%Initializing array and assigning BC
T=ones(nx,ny);
T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;
Told=T;
%Constant
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
cfl_num=alpha*(dt/dx^2)
iterative_solver=input
%Jacobi_iteration
if iterative_solver==1
tic
for k=1:nt % time loop
for i=2:nx-1
for j=2:ny-1
Term1=(1-2*k1-2*k2);
Term2=k1;
Term3=k2;
H=(Told(i-1,j)+Told(i+1,j));
V=(Told(i,j-1)+Told(i,j+1));
T(i,j)=(Told(i,j)*Term1)+(Term2*H)+(Term3*V);
end
end
Told=T;
end
time1=toc;
figure(1)
colormap(jet)
contourf(x,y,T,'showText','on')
colorbar
xlabel('x axis')
ylabel('Y axis')
text1=sprintf('Simulation time =%d s',time1)
text2=sprintf('CFL number=%d',cfl_num)
text3=sprintf('Time step=%d',dt)
title(['Transient state explicit using jacobi';text1;text2;text3])
colorbar
pause(0.03)
end
if iterative_solver==2
%gs_iteration
tic
for k=1:nt
for i=2:nx-1
for j=2:ny-1
Term1=(1-2*k1-2*k2);
Term2=k1;
Term3=k2;
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j)=(Told(i,j)*Term1)+(Term2*H)+(Term3*V);
end
end
Told = T;
end
time2=toc;
figure(2)
colormap(jet)
contourf(x,y,T,'showText','on')
xlabel('X axis')
ylabel('Y axis')
text1=sprintf('Simulation time =%d s',time2)
text2=sprintf('CFL number=%d',cfl_num)
text3=sprintf('Time step=%d',dt)
title(['Transient state using gs_solver';text1;text2;text3])
colorbar
pause(0.3)
end
%sor_iteration
if iterative_solver==3
for k=1:nt
for i=2:nx-1
for j=2:ny-1
Term1=(1-2*k1-2*k2);
Term2=k1;
Term3=k2;
H=(T(i-1,j)+Told(i+1,j));
V=(T(i,j-1)+Told(i,j+1));
T(i,j)=Told(i,j)*(1-rf)+rf*((Told(i,j)*Term1)+(Term2*H)+(Term3*V));
end
end
Told = T;1
end
time3=toc;
figure(3)
colormap(jet)
contourf(x,y,T,'showText','on')
xlabel('X axis')
ylabel('Y axis')
text1=sprintf('Simulation time =%d s',time3)
text2=sprintf('CFL number=%d',cfl_num)
text3=sprintf('Time step=%d',dt)
title(['Transient state using sor_solver';text1;text2;text3])
colorbar
pause(0.3)
endResult:-Following plots shows the results for 2D heat conduction for Transient-state explicit using different iterative methods.
Observation:-
- By observing,above plots its quite clear that jacobi_iteration takes more number of iteration and computation time as compared to other two iterative solver this is because jacobi method takes the old values.
- While the other two iterative method which is Gauss seidel & Successive over relaxiation iterative method shows most convergence result for same value of number of time steps as compared to jacobi method this is beacuase Gauss seidel & Successive over relaxiation takes the updated values.Hence the new values is passed in new iteration as the solution proceeds.Thus in this type of iterative solver we get convergence result in less number of iteration as well as in less computation time.
b.Implicit scheme:-
Govering equation for 2D heat conduction equation for unsteady state is,
.....................(3)
After applying first order forward differencing scheme on L.H.S and central differencing scheme on R.H.S on eqaution(3),we get
Let says,
And,on simplification we get
where k1 & k2 are,
..................(4)
The equation (4) will be solve by using point iterative method.The point iterative method which will be use to solve the equation(4) are as follows:-
- Jacobi
- Gauss seidal
- Successive over-relaxation
Matlab code:-
- In the following program we have solve eqaution(4) by using diiferent iterative method which is mentioned above.
- First,we have assume that number of grid in x direction is equal to number of grid in y direction
- Then error and tolerance is set in the time loop and then equation (4) is solve one by one by diiferent iterative method in a loop.And each iterative method solve equation(4) till convergence result is obtain.
- After,then result is plotted using contourf command.
clear all
close all
clc
%Solving transient 2D heat conduction equation
%Input parameter
nx=40; %no of grids in x-direction
ny=40; %no of grids in y_direction
rf=1.5; %relaxation factor
dt=change;
alpha=1.4;
nt=400;
%Calculation
x=linspace(0,1,nx);
y=linspace(0,1,ny);
dx=x(2)-x(1); %Time step in x-direction
dy=y(2)-y(1); %Time step in y-direction
%Initializing array and assigning BC
T=ones(nx,ny);
T(:,1)=400;
T(1,:)=600;
T(:,end)=800;
T(end,:)=900;
Told=T;
T_prev=T;
%Constant
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
%CFL_number
cfl_num=alpha*(dt/dx^2)
error=9e9; %error
tol=1e-4; %tolerance
iterative_solver=input
tic
jacobi_iter=1
if iterative_solver==1
for k=1:nt
error=9e9;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
Term1=(1+2*k1+2*k2)^-1;
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)+(Term2*H)+(Term3*V);
end
end
error=max(max(abs(Told-T)));
Told = T;
jacobi_iter=jacobi_iter+1
end
T_prev=T;
end
time1=toc;
figure(1)
colormap(jet)
contourf(x,y,T,'showText','on')
xlabel('X axis')
ylabel('Y axis')
text1=sprintf('Number of iteration =%d',jacobi_iter)
text2=sprintf('Simulation time =%d s',time1)
text3=sprintf('CFL number=%d',cfl_num)
text4=sprintf('Time step=%d',dt)
title(['Transient state using jacobi_solver';text1;text2;text3;text4])
colorbar
pause(0.3)
end
tic
gs_iter=1
if iterative_solver==2
for k=1:nt
error=9e9;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
Term1=(1+2*k1+2*k2)^-1;
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)+(Term2*H)+(Term3*V);
end
end
error=max(max(abs(Told-T)));
Told = T;
gs_iter=gs_iter+1
end
T_prev=T;
end
time2=toc;
figure(2)
colormap(jet)
contourf(x,y,T,'showText','on')
xlabel('X axis')
ylabel('Y axis')
text1=sprintf('Number of iteration =%d',gs_iter)
text2=sprintf('Simulation time =%d s',time2)
text3=sprintf('CFL number=%d',cfl_num)
text4=sprintf('Time step=%d',dt)
title(['Transient state using gs_solver';text1;text2;text3;text4])
colorbar
pause(0.3)
end
tic
sor_iter=1
if iterative_solver==3
for k=1:nt
error=9e9;
while(error>tol)
for i=2:nx-1
for j=2:ny-1
Term1=(1+2*k1+2*k2)^-1;
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)=rf*((T_prev(i,j)*Term1)+(Term2*H)+(Term3*V))+(1-rf)*Told(i,j);
end
end
error=max(max(abs(Told-T)));
Told = T;
sor_iter=sor_iter+1
end
T_prev=T;
end
time3=toc;
figure(3)
colormap(jet)
contourf(x,y,T,'showText','on')
xlabel('X axis')
ylabel('Y axis')
text1=sprintf('Number of iteration =%d',sor_iter)
text2=sprintf('Simulation time =%d s',time3)
text3=sprintf('CFL number=%d',cfl_num)
text4=sprintf('Time step=%d',dt)
title(['Transient state using sor_solver';text1;text2;text3;text4])
colorbar
pause(0.3)
endResult:-Following plots shows the results for 2D heat conduction for Transient-state explicit using different iterative methods.
Observation:-
- Tough we are getting same result for different iterative method so, the different between them is number of iteration and computation time taken by iterative solver to get convergence result.
- By observing,above plots its quite clear that jacobi_iteration takes more number of iteration and computation time as compared to other two iterative solver this is because jacobi_iteration takes the old values.
- While other two iterative solver(Gauss seidel & Successive over relaxiation) takes the updated values.Hence the new values is passed in new iteration as the solution proceeds.Thus in this type of iterative solver we get convergence result in less number of iteration as well as in less computation time.
Now,observing stability analysis of above programs :-
1.Explicit scheme:-Observing effect of CFL number and dt on 2D heat conduction of Transient-state explicit scheme using different iterative methods.
a.When dt=1e-4 and CFL=0.21294
- The following result looks stable and converge,since CFL is less then its stability criteria which 0.25.
b.When dt=1.2e-4 and CFl=0.259
- When we slightly increase the value of dt,the jacobi started to become unstable,since CFL number become little bit above its stability criteria,however other two method which is gauss seidel and successive over relaxation has successfully achieve their result.
c.When dt=4.7e-4 and CFL=1.00082
- When dt is little more increase,the result for all solver for explicit become completely unstable.
2.Implicit scheme:-Since,implicit scheme is inherently stable which means high value CFL number cannot effect the solution which is shown in below plots.