2D Heat conduction equation using a Transient solver and a Steady state solver using Iterative techniques Jacobi Gauss Seidel SOR

JORGE LOZANO

Introduction to CFD using Matlab and OpenFOAM

Midterm project - Solving the steady and unsteady 2D heat conduction problem

OBJECTIVE

Solve the conduction equation using a Transient solver and a Steady state solver using Iterative techniques (Jacobi, Gauss Seidel, SOR)

1. Assume that the domain is a unit square.
2. Assume 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 transient case.
4. Left:400K
5. Top:600K
6. Right:800K
7. Bottom:900K

Absolute error criteria = 1e-4

1. You will write code for both implicit and explicit schemes (This is for the Transient Analysis only)
2. Compare your fastest steady state simulation vs your fastest transient simulation and explain your results
3. For both implicit and explicit methods look at the effect of time step on the stability of the solution. You need to answer the following questions
4. How does an unstable solution look like?
5. At what CFL number does your solution become unstable?
6. Write a detailed explanation for your results, plots or codes without explanation will awarded zero marks.

THEORY

The conduction equations that solve these types of problems are specifically for a 2D heat conduction steady and unsteady case.

The following equation represents a two-dimensional heat conduction equation:

Also viewed along the alpha variable, being an important value of the material property in the geomerty:

When applying this equation to a steady state case the time component is eliminated as seen next.

Implicit vs. Explicit Numerical Methods

Numerical solution schemes are often referred to as being explicit or implicit. When a direct computation of the dependent variables can be made in terms of known quantities, the computation is said to be explicit. When the dependent variables are defined by coupled sets of equations, and either a matrix or iterative technique is needed to obtain the solution, the numerical method is said to be implicit.

Iterations are used to advance a solution through a sequence of steps from a starting state to a final, converged state. This is true whether the solution sought is either one step in a transient problem or a final steady-state result. In either case, the iteration steps resemble a time-like process. Of course, the iteration steps usually do not correspond to a realistic time-dependent behavior. In fact, it is this aspect of an implicit method that makes it attractive for steady-state computations, because the number of iterations required for a solution is often much smaller than the number of time steps needed for an accurate transient that asymptotically approaches steady conditions.

The Stability Issue

Numerical stability has to do with the behavior of the solution as the time-step dt is increased. If the solution remains well behaved for arbitrarily large values of the time step, the method is said to be unconditionally stable. This situation never occurs with explicit methods, which are always conditionally stable. It is easy to see this by dividing the Q-equation by dt and then letting dt approach infinity.

CFL: COURANT-FRIEDRICHS-LEWY condition

The convergence condition by Courant–Friedrichs–Lewy is a necessary condition for convergence while solving certain partial differential equations (usually hyperbolic PDEs) numerically. It arises in the numerical analysis of explicit time integration schemes, when these are used for the numerical solution. Therefore, the time step must be less than a certain time in many explicit time-marching computers simulations, otherwise the simulation produces incorrect results. The condition is named after Richard Courant, Kurt Friedrichs, and Hans Lewy who described it in their 1928 paper.

The equation used for obtaining the time step and the CFL value are the following equations, the first one being a general 2D specific equation and the second with the condition that the distance between each node are the same.

SIMULATIONS

STEADY-STATE THERMAL ANALYSIS

The code presented is specifically for a steady state implicit 2D case where there are 30 nodes for dx and dy, having a square with determined temperature boundary conditions. The implicit methods applied are the Jacobian method, Gauss-Seidel method, and the Successive over relaxation (SOR) method.

clear all
close all
clc
x=input(\'x and y nodes number =\');
y=x;
m=linspace(0,1,x);
n=linspace(0,1,y);
dx=1/(x-1);
dy=1/(y-1);
T_prev= 500*ones(x,y);
T_old=T_prev;
T_prev(:,1)=400;
T_prev(:,y)=800;
T_prev(1,:)=600;
T_prev(x,:)=900;
error=9e9;
tol=1e-4;
solver=input(\'number for solver =\');
alpha=1.5;
k=2*(dx^2 + dy^2)/((dx^2)*(dy^2));
tic;
%JACOBI METHOD
if solver==1
  j_iter=1;
  while(error>tol)
  for i=2 :x-1 
    for j=2 :y-1
      T_prev(i,j)=(k^-1)*(((T_old(i-1,j)+T_old(i+1,j))/(dx^2))+((T_old(i,j+1))+T_old(i,j-1))/(dy^2));
    endfor
  endfor
  error=max(max(abs(T_prev-T_old)));
  T_old=T_prev;
  j_iter=j_iter+1;
end
toc
sst=toc;
figure(1)
[X,Y]=contourf(m,n,T_prev,25);
colormap(bone)
colorbar
xlabel(\'X AXIS TEMPERATURE\');
ylabel(\'Y AXIS TEMPERATURE\');
title_text=sprintf(\'STEADY STATE JACOBIAN \\nITERATION NUMBER=%d \\n SIMULATION TIME=%d SECONDS\',j_iter,sst);
title(title_text)
end
%Gauss seidel method new values are used
if solver==2
  g_iter=1;
  while(error>tol)
  for i=2 : x-1
    for j=2 : y-1
      T_prev(i,j)=(k^-1)*((T_prev(i-1,j)+T_old(i+1,j))/(dx^2)+(T_old(i,j+1)+T_prev(i,j-1))/(dy^2));
    endfor
  endfor
  error=max(max(abs(T_prev-T_old)));
  T_old=T_prev;
  g_iter=g_iter+1;
end
toc
sst=toc;
figure(2)
[X,Y]=contourf(m,n,T_prev,25);
colormap(bone)
colorbar
xlabel(\'X AXIS TEMPERATURE\')
ylabel(\'Y AXIS TEMPERATURE\')
title_text=sprintf(\'STEADY STATE GAUSS SEIDEL \\nITERATION NUMBER=%d \\n SIMULATION TIME=%d SECONDS\',g_iter,sst);
title(title_text)
end
%Successive over relaxation
if solver==3
  sor_iter=1;
  while(error>tol)
  for i=2 : x-1
    for j=2 : y-1
      k1=(T_prev(i-1,j)+T_old(i+1,j))/(dx^2);
      k2=(T_old(i,j+1)+T_prev(i,j-1))/(dy^2);  
      T_prev_g(i,j)=(k1+k2)/k;
      T_prev(i,j)=T_old(i,j)*(1-alpha)+alpha*(T_prev_g(i,j));
    endfor
  endfor
  error=max(max(abs(T_prev-T_old)));
  T_old=T_prev;
  sor_iter=sor_iter+1;
end
toc
sst=toc;
figure(3)
[X,Y]=contourf(m,n,T_prev,25);
colormap(bone)
colorbar
xlabel(\'X AXIS TEMPERATURE\');
ylabel(\'Y AXIS TEMPERATURE\');
title_text=sprintf(\'STEADY STATE SUCCESSIVE OVER RELAXATION \\nITERATION NUMBER=%d \\n SIMULATION TIME=%d SECONDS \',sor_iter,sst);
title(title_text)
end

As seen in the postprocess data of the simulation the cases are similar in the results of the temperature which represents having a good accuracy and application for each method, although the Jacobian method takes the longest time to obtain a convergence in the results, the second fastest is the Gauss-Seidel method and the fastest overall is the Successive over relaxation method, that though has more mathematic equations to solver for each iterations, due to the approximation over using the newest values is faster to obtain the final result.

UNSTEADY-STATE THERMAL ANALYSIS

The code presented is specifically for an unsteady state EXPLICIT 2D case where there are 30 nodes for dx and dy, having a square with determined temperature boundary conditions, having as analyzed times: 50sec. and 500sec.

EXPLICIT METHOD

clear all
close all
clc


nx=15;
ny=15;
dt=1;
t=500;
nt=t/dt;
alpha=1e-4;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);
CFL=alpha*dt*((1/dx^2)+(1/dy^2));
T=600*ones(ny,nx);
T(:,1)=400;
T(:,ny)=800;
T(nx,:)=600;
T(1,:)=900;
T_old=T;
k1=alpha*dt/dx^2;
k2=alpha*dt/dy^2;

iter=0;
tic

for k=1:nt
  for j=2:ny-1
    for i=2:nx-1
    H=(T_old(i-1,j)-2*T_old(i,j)+T_old(i+1,j));
    V=(T_old(i,j-1)-2*T_old(i,j)+T_old(i,j+1));
    T(i,j)=T_old(i,j)+(H*k1)+(V*k2);
    end
  end
  T_old=T;
  iter=iter+1;
end

time=toc;
figure(1)
[c,h]=contourf(x,y,T,25);
colorbar
colormap(bone)
xlabel(\'X AXIS TEMPERATURE\');
ylabel(\'Y AXIS TEMPERATURE\');
formatSpec=\"UNSTEADY STATE-EXPLICIT \\n ITERATION NUMBER= %d \\n SIMULATION TIME= %f SECONDS \\n TIME ANALYZED= %d SECONDS\\nCFL= %d\";
title(sprintf(formatSpec,iter,time,t,CFL));

The following two images represent the postprocess for the explicit method by using a CFL number of 0.01682 for both.

The following two images represents a less defined CFL number as 0.841 being a really wrong number to use for these simulatio, and a calculated CFL number as 0.5.

IMPLICIT METHOD

The code presented is specifically for an unsteady state IMPLICIT 2D case where there are 30 nodes (ofr the 500s 20 nodes are used) for dx and dy, having a square with determined temperature boundary conditions, having as analyzed times: 50sec. and 500sec.

clear all
close all
clc

nx=input(\'x and y nodes number =\');
ny=nx;
solver=input(\'number for solver =\');
dt=.2;
t=200;
nt=t/dt;
alpha=1e-4;
m=linspace(0,1,nx);
dx=m(2)-m(1);
n=linspace(0,1,ny);
dy=n(2)-n(1);
CFL=(alpha*dt)/((1/dx^2)+(1/dy^2));
T=500*ones(nx,ny);
T(:,1)=400;
T(nx,:)=600;
T(:,ny)=800;
T(1,:)=900;
T_old=T;
T_prev=T;
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
term1=1/(1+2*k1+2*k2);
term2=k1/term1;
term3=k2/term1;
tol=1e-4;
tic

%JACOBI
if solver==1
M=(\"GAUSS JACOBI\");
for iter=0:nt-2
  error=1;
 while error>tol
    for j=2:ny-1
     for i=2:nx-1        
          H=(T_old(i-1,j)+T_old(i+1,j));
          V=(T_old(i,j-1)+T_old(i,j+1));
          T(i,j)=(T_prev(i,j)*term1)+(term2*H)+(term3*V);
  end
  end
  error=max(max(abs(T-T_old)));
  T_old=T;
  iter=iter+1;
end
T_prev=T;
end

%GS
elseif solver==2
M=(\"GAUSS SIEDEL\");
for iter=0:nt-2
  error=1;
    while error>tol
     for j=2:ny-1
       for i=2:nx-1    
          H=(T(i-1,j)+T_old(i+1,j));
          V=(T(i,j-1)+T_old(i,j+1));
          T(i,j)=(T_prev(i,j)*term1)+(term2*H)+(term3*V);
  end
  end
  error=max(max(abs(T-T_old)));
  T_old=T;
iter=iter+1;
end
T_prev=T;
end

%SOR
elseif solver==3
M=(\"SOR\");
omega=1.1;
for iter=0:nt-2
  error=1;
    while error>tol
     for j=2:ny-1
        for i=2:nx-1   
          H=(T(i-1,j)+T_old(i+1,j));
          V=(T(i,j-1)+T_old(i,j+1));
          T(i,j)=(T_prev(i,j)*term1)+(H*term2)+(V*term3);
          T(i,j)=T_old(i,j)*(1-omega)+omega*T(i,j);
  end
  end
  error=max(max(abs(T-T_old)));
  T_old=T;
  iter=iter+1;
end
T_prev=T;
end
end


time=toc;
figure(1);
[c,h]=contourf(m,n,T,25);
colorbar
colormap(bone)
xlabel(\'X AXIS TEMPERATURE\');
ylabel(\'Y AXIS TEMPERATURE\');
formatSpec=\"UNSTEADY STATE-IMPLICIT %s\\n ITERATION NUMBER= %d \\n SIMULATION TIME=%f SECONDS\";
title(sprintf(formatSpec,M,iter,time));

IMPLICIT JACOBIAN METHOD

IMPLICIT GAUSS-SEIDEL METHOD

IMPLICIT SUCCESSIVE OVER RELAXATION

1. Compare your fastest steady state simulation vs your fastest transient simulation and explain your results

As expected, the steady state simulations compared to the transient ones are faster. According to steady state simulations, the Jacobian method is the lowest and performs more iterations, the second slowest is the Gauss-Seidel method, having the fastest method as the SOR, with less iterations performed and less computational time. In the transient simulations, the Gauss-Seidel method was the fastest for both times analyzed, the second was the Jacobian method and the slowest by much was the SOR method for both times.

By comparing the fastest for steady and unsteady simulations having a computational time of 11 seconds (steady-state), and 388 seconds (unsteady-state), taking more computational time as transient methods.

1. How does an unstable solution look like?

An unstable solution, where the CFL number and therefore the time step is greater than the stability number results in greater and unknown values as follows:

1. At what CFL number does your solution become unstable?

For an explicit unsteady method:

By applying the next formula, for a 2D thermal case, and CFL as 0.5 when dx equals dy:

• Alpha:1e-4
• CFL=0.5
• dx = dy = (1/29)

Dt=2.97265

For an implicit unsteady method:

By using the Jacobian method, by using dt as 0.2, where dx and dy the same, by using the same formula:

• Alpha:1e-4
• dx = dy = (1/19)

CFL=0.01444

As researched, an implicit unsteady method doesn’t apply the rule for the CFL condition being less or equal than 0.25, in this case for an even smaller value, the results show a non-convergence along the temperatures, as the next image.

CONCLUSIONS

Several 2D heat conduction steady and unsteady numeric simulations were performed by varying several important variables and coefficients affecting directly the postprocessing results.

The CFL stability coefficient to determine the time step only applies for explicit unsteady simulations, whereas seen, it must be below or equal to 0.5 for this geometry. For these cases analyzed, the fastest numeric unsteady simulation obtained is the explicit solver.

Dt=2.97265 is the value obtained by using the CFL equation for 0.5, which is the greatest time step value available to apply without obtaining unstable results, only being able to apply for explicit unsteady simulations, having the fastest unsteady computational time. Although the SOR method results as being the fastest steady-state computational time.

REFERENCES

• ‘Solution of the Two-Dimensional Steady State Heat Conduction using the Finite Volume Method”, Rico Morasata 1, Kürsad Melih Güleren 2.