AMIRKABIR UNIVERSITY OF TECHNOLOGY AEROSAPCE ENGINEERING DEPARTMENT
Numerical Solution of 2D Heat Equation by ADI and SLOR methods Computational Fluid Dynamics Course Assignment Instructor: Dr. Jahangirian Ata Ghasemi Esfahani student ID: 92129071
WINTER
2014
Contents 1. Introduction ........................................................................................................................ 3 2. Problem definition ............................................................................................................... 3 3. Discretization of 2D heat equation and MATLAB codes ............................................................. 4 3.1 ADI method ................................................................................................................... 4 3.2 SLOR method ................................................................................................................. 8 4. Results and discussion .........................................................................................................12 4.1 ADI method ..................................................................................................................12 4.2 SLOR method ................................................................................................................21 4.3 Discussion.....................................................................................................................24 5. References .........................................................................................................................25
2
1. Introduction This report summarizes the results obtained from implementation of various implicit methods for solving two-dimensional heat conduction equation with a specific set of boundary and initial conditions. The most crucial advantage associated with explicit method is its relative ease of implementation. On the other hand, in order for the solution to be stable, for a fixed value of ∆x, one must consider limited values of ∆t such that the stability criteria are met. This usually results in long computer run times. Implicit methods often maintain stability even for much larger values of ∆t, but the massive matrix manipulations needed at each time step which results in an increase in computation time. Also the implicit methods might not be as accurate in yielding transient solutions since truncation error is large. The methods investigated in this assignment, namely ADI and SLOR, have been devised to take advantage of the so-called tri-diagonal matrix solution method that can be coded and implemented easily. Instead of attempting to manipulate massive matrices, the ADI method resorts to a method called two-level solution. This approach ensures that at each time step, there are no more than three unknowns to solve for thus the tri-diagonal solution routine, also known as Thomas algorithm, can be employed. The SLOR (SOR by line, where SOR stands for successive over-relaxation) is a line-iterative method. Although relaxation techniques can be formulated in explicit or implicit forms, the specific discretization form presented here is implicit and again, relies on Thomas algorithm for matrix manipulations. Within this framework, some quantities are assumed to be known at initial iterative step and some other are updated as the iterative solution progresses. Although both the ADI and SLOR methods are classified as implicit methods, care must be taken in selecting grid size and time steps. As will be seen in what follows, careful selection of the mentioned parameters will result in a physically sound solution and improper selection of those yields inappropriate results.
2. Problem definition
Constant Temperature T=100 °C
Constant Temperature T=0 °C
Constant Temperature T=300 °C
Constant Temperature T=100 °C
3
3. Discretization of 2D heat equation and MATLAB codes 3.1 ADI method The unsteady two-dimensional heat conduction equation (parabolic form) has the following form:
A forward time, central space scheme is employed to discretize the governing equation as described in the next page. The idea of the ADI-method (alternating direction implicit) is to alternate direction and thus solve two one-dimensional problem at each time step. The first step keeps y-direction fixed:
In the second step we keep x-direction fixed:
The ADI calculation procedure is depicted in the figure below:
Both equations can be written in a triadiagonal form. If we define: 4
Then we get:
Instead of five-band matrix in BTCS method, here each time step can be obtained in two sweeps. Each sweep can be done by solving a tri-diagonal system of equations. The ADI method is second order accurate with a T.E. of the order O[ Δt2 , Δx2, Δy2]. The amplification factor is as follows:
Where
It is obvious that this method is unconditionally stable, although the 3D ADI scheme is only conditionally stable. The first step of the ADI scheme is reduced to the following form that can be implemented as a MATLAB code: ATI-1,jn+1/2 – BTI,jn+1/2 +ATi+1,jn+1/2 = Ki Where: A= (α∆t/2∆x2) B= 1+ (α∆t/∆x2) Ki = -Ti,jn - (α∆t/2∆x2) [ Ti,j+1n -2Ti,jn +Ti,j-1n ]
The second step of the ADI scheme is implemented as follows: CTI-1,jn+1/2 – DTI,jn+1/2 +CTi+1,jn+1/2 = Ki
5
Where: C= (α∆t/2∆y2) D= 1+ (α∆t/∆y2) Lj = -Tin+1/2 - (α∆t/2∆x2) [ Ti+1,j+n+1/2 -2Ti,jn+1/2 +Ti-1,jn+1/2 ]
The matrix formed has the following form:
The MATLAB code written to solve the 2D diffusion equation by ADI method can be seen below: %Heat equation solution by ADI method % by Ata Ghasemi Esfahani 92129071 clear all clc % constants alfa=0.02; dt=0.5; dx=0.2; dy=0.2; a=(alfa*dt)/(2*dx^2); b=1+2*a; al=(alfa*dt)/(2*dy^2); bl=1+2*a; m=10; n=10000; %initial conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=2:m for i=2:m te(j,i)=0; end end
Initial conditions: zero temperature throughout the domain
6
%boundary conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:m+1 te(1,i)=300; te(m+1,i)=0; end
Boundary conditions
for j=2:m+1 te(j,1)=100; te(j,m+1)=100; end te(:,:) % formation of matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=1:m for i=1:m if (i==j) ma(j,i)=-b; end if (i==j+1) ma(j,i)=a; end if (i==j-1) ma(j,i)=a; end end end
Formation of coefficients matrix
for t=1:n for j1=2:m for i=2:m k(i-1)=-te(j-1,i)-((alfa*dt)/(2*dy^2))*(te(j-1,i+1)2*te(j-1,i)+te(j-1,i-1)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Thomas Algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% md(1,1)=-b; for j=1:m-1 for i=1:m-1 if (j~=1) if (i==j) md(j,i)=ma(j,i)-md(j-1,i)*a/md(j-1,j-1); end end if (i==j+1) md(j,i)=a; end end end 7
for i=2:9 k(i)=k(i)-k(i-1)*a/md(i-1,i-1); end te(j1,m)=k(m-1)/md(m-1,m-1); for i=1:m-2 i=m-i; te(j1,i)=(k(i-1)-a*te(j1,i+1))/md(i-1,i-1); end end for i=2:m for j=2:m l(j-1)=-te(j,i-1)-((alfa*dt)/(2*dx^2))*(te(j+1,i-1)2*te(j,i-1)+te(j-1,i-1)); end for i=2:9 l(i)=l(i)-l(i-1)*a/md(i-1,i-1); end te(m,i)=l(m-1)/md(m-1,m-1); for j=1:m-2 j=m-j; te(j,i)=(l(j-1)-a*te(j+1,i))/md(j-1,j-1); end end end
3.2 SLOR method Assuming the 2D steady heat equation has the following form (Laplace's equation):
8
We now have to describe the SLOR method. The general idea behind successive overrelaxation method, also classified as an iterative method, is to introduce an arbitrary correction to the intermediate values of the unknown considered in the following way:
where k denotes iteration level, ui,jk+1 is the most recent value of ui,,j calculated from the Gauss-Seidel procedure, ui,jk' is the value from the previous iteration as adjusted by previous application of this formula if the over-relaxation is being applied successively (at each iteration), and ui,jk+1' is the newly adjusted or "better guess" for ui,,j, at the k + 1 iteration level. That is, we expect ui,jk+1' to be closer to the final solution than the unaltered value ui,jk+1 from the Gauss-Seidel calculation. Here, w is the relaxation parameter, and when 1 < w < 2, overrelaxation is being employed. For Laplace's equation on a rectangular domain with Dirichlet boundary conditions, theories pioneered by Young (1954) and Frankel (1950) lead to an expression for the optimum w (hereafter denoted by wop,). First, defining σ as:
and the optimum w is given by:
To illustrate the procedure, consider the solution to Laplace’s equation on a square domain with Dirichlet boundary conditions using the five-point scheme. If we agree to start at the bottom of the square and sweep up by rows, we could write, for the general point:
Alternatively, the over-relaxation parameter w can be introduced prior to solution of the simultaneous algebraic equations. This is done by substituting the over-relaxation scheme introduced earlier into the above equation. The result is:
The above equation is then solved for each row by Thomas algorithm. The graphic depiction of the procedure can be seen in the figure below: 9
For the specified conditions of our problem, β is equal to 1 and wopt is calculated to be 1.989. The MATLAB code can be seen in the following pages: %Heat equation solution by SLOR method %by Ata Ghasemi Esfahani 92129071 clear all clc %Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% omega=1.989; a=-omega/4; m=10; iteration=1000;
Solution domain increments
%Initial conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for t=1:iteration for j=1:m for i=1:m te(t,j,i)=0; end end end %Boundary conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for t=1:iteration+1 for i=1:m+1 te(t,1,i)=300; te(t,m+1,i)=0; end end 10
for t=1:iteration+1 for j=2:m+1 te(t,j,1)=100; te(t,j,m+1)=100; end end te2(:,:)=te(iteration,:,:); te2(:,:) %formation of matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=1:m for i=1:m if (i==j) ma(j,i)=1; end if (i==j+1) ma(j,i)=a; end if (i==j-1) ma(j,i)=a; end end end
Formation of coefficients matrix
for t=1:iteration for j=2:m for i=2:m k(i-1)=(1-omega)*te(t,j,i-1)+te(t,j+1,i-1)+te(t+1,j-1,i1); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Thomas Algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% md(1,1)=1; for j1=1:m-1 for i=1:m-1 if (j1~=1) if (i==j1) md(j1,i)=ma(j1,i)-(a^2)/md(j1-1,i); end end if (i==j1+1) md(j1,i)=a; end end end for i=2:9 k(i)=k(i)-k(i-1)*a/md(i-1,i-1); end te(t,j,m)=k(m-1)/md(m-1,m-1); for i=1:m-1 11
i=m-i; te(t,j,i)=(k(i)-a*te(t,j,i+1))/md(i,i); end end end te1(:,:)=te(iteration,:,:);
4. Results and discussion 4.1 ADI method The results of the numerical procedure are presented in this section. Temperature distributions in the domain are presented so that the reader is able to compare the consequences of introducing various changes to the problem. In what follows n denotes the time level.
Fig. 1 Temperature distribution, ADI method, Δx = 0.1, n = 5, Δt = 0.1
12
Fig. 2 Temperature distribution, ADI method, Δx = 0.1, n = 5, Δt = 0.5
Fig. 3 Temperature distribution, ADI method, Δx = 0.2, n = 5, Δt = 0.1
13
Fig. 3 Temperature distribution, ADI method, Δx = 0.2, n = 5, Δt = 0.5
Fig. 4 Temperature distribution, ADI method, Δx = 0.1, n = 30, Δt = 0.1 14
Fig. 4 Temperature distribution, ADI method, Δx = 0.1, n = 50, Δt = 0.1
Fig. 5 Temperature distribution, ADI method, Δx = 0.1, n = 100, Δt = 0.1
15
Fig. 6 Temperature distribution, ADI method, Δx = 0.1, n = 1000, Δt = 0.1
Fig. 7 Temperature distribution, ADI method, Δx = 0.2, n = 1000, Δt = 0.1
16
Fig. 8 Temperature distribution, ADI method, Δx = 0.2, n = 1000, Δt = 0.5
300 300 300 300 300 300 300 300 300 300 300 100
0
0
0
0
0
0
0
0
0
100
100
0
0
0
0
0
0
0
0
0
100
100
0
0
0
0
0
0
0
0
0
100
100
0
0
0
0
0
0
0
0
0
100
100
0
0
0
0
0
0
0
0
0
100
100
0
0
0
0
0
0
0
0
0
100
100
0
0
0
0
0
0
0
0
0
100
100
0
0
0
0
0
0
0
0
0
100
100
0
0
0
0
0
0
0
0
0
100
100
0
0
0
0
0
0
0
0
0
100
Fig. 9 Temperature distribution, ADI method, Δx = 0.1, n = 1, Δt = 0.1, solution matrix
17
Columns 1 through 7
300.0000 300.0000 300.0000 300.0000 300.0000 300.0000 300.0000 100.0000 274.8239 297.8872 299.8227 299.9850 299.9975 299.9850 100.0000 49.8107 49.6214 49.4322 49.2429 49.0536 48.8643 100.0000 49.8107 49.6214 49.4322 49.2429 49.0536 48.8643 100.0000 49.8107 49.6214 49.4322 49.2429 49.0536 48.8643 100.0000 49.8107 49.6214 49.4322 49.2429 49.0536 48.8643 100.0000 49.8107 49.6214 49.4322 49.2429 49.0536 48.8643 100.0000 49.8107 49.6214 49.4322 49.2429 49.0536 48.8643 100.0000 49.8107 49.6214 49.4322 49.2429 49.0538 48.8666 100.0000 49.8107 49.6214 49.4322 49.2429 49.0538 48.8666 100.0000
0
0
0
0
0
0
Columns 8 through 11
300.0000 300.0000 300.0000 300.0000 299.8227 198.0254 274.8239 100.0000 48.6750 102.3387 49.2429 100.0000 48.6750 58.3989 49.2429 100.0000 48.6750 49.9288 49.2429 100.0000 48.6750 48.8138 49.2429 100.0000 48.6750 48.6585 49.2429 100.0000 48.6750 48.3911 49.2429 100.0000 48.7019 46.3765 49.2695 100.0000 48.7019 34.8840 49.2695 100.0000 0
0
0 100.0000
Fig. 10 Temperature distribution, ADI method, Δx = 0.1, n = 1000, Δt = 0.1, solution matrix 18
Columns 1 through 7
300.0000 300.0000 300.0000 300.0000 300.0000 300.0000 300.0000 100.0000 292.8531 299.8297 299.9959 299.9999 300.0000 299.9999 100.0000 49.9931 49.9863 49.9794 49.9726 49.9657 49.9588 100.0000 49.9931 49.9863 49.9794 49.9726 49.9657 49.9588 100.0000 49.9931 49.9863 49.9794 49.9726 49.9657 49.9588 100.0000 49.9931 49.9863 49.9794 49.9726 49.9657 49.9588 100.0000 49.9931 49.9863 49.9794 49.9726 49.9657 49.9588 100.0000 49.9931 49.9863 49.9794 49.9726 49.9657 49.9588 100.0000 49.9931 49.9863 49.9794 49.9726 49.9657 49.9588 100.0000 49.9931 49.9863 49.9794 49.9726 49.9657 49.9588 100.0000
0
0
0
0
0
0
Columns 8 through 11
300.0000 300.0000 300.0000 300.0000 299.9959 264.7318 292.8531 100.0000 49.9520 70.9772 49.9726 100.0000 49.9520 50.9738 49.9726 100.0000 49.9520 49.9887 49.9726 100.0000 49.9520 49.9531 49.9726 100.0000 49.9520 49.9518 49.9726 100.0000 49.9520 49.9442 49.9726 100.0000 49.9520 49.7358 49.9726 100.0000 49.9520 45.4135 49.9726 100.0000 0
0
0 100.0000
Fig. 11 Temperature distribution, ADI method, Δx = 0.2, n = 1000, Δt = 0.1, solution matrix Columns 1 through 7
300.0000 300.0000 300.0000 300.0000 300.0000 300.0000 300.0000 19
100.0000 219.6133 278.4530 294.1989 298.3425 299.1713 298.3425 100.0000 48.7815 47.5629 46.3444 45.1258 43.9073 42.6887 100.0000 48.7815 47.5629 46.3444 45.1258 43.9073 42.6887 100.0000 48.7815 47.5629 46.3444 45.1258 43.9073 42.6887 100.0000 48.7815 47.5629 46.3444 45.1258 43.9073 42.6887 100.0000 48.7815 47.5629 46.3444 45.1258 43.9073 42.6887 100.0000 48.7815 47.5629 46.3444 45.1258 43.9073 42.6887 100.0000 48.7821 47.5655 46.3540 45.1619 44.0419 43.1912 100.0000 48.7821 47.5655 46.3540 45.1619 44.0419 43.1912 100.0000
0
0
0
0
0
0
Columns 8 through 11
300.0000 300.0000 300.0000 300.0000 294.1989 132.3841 219.6133 100.0000 41.4702 99.9029 45.1258 100.0000 41.4702 84.3623 45.1258 100.0000 41.4702 60.1294 45.1258 100.0000 41.4702 47.9626 45.1258 100.0000 41.4702 42.2473 45.1258 100.0000 41.4702 37.2052 45.1258 100.0000 43.3455 28.5179 46.8756 100.0000 43.3455 16.2208 46.8756 100.0000 0
0
0 100.0000
Fig. 12 Temperature distribution, ADI method, Δx = 0.1, n = 1000, Δt = 0.5, solution matrix
20
4.2 SLOR method
Fig. 13 Temperature distribution, SLOR method, Δx = 0.1, iteration level = 1000, w=1.989
Fig. 14 Temperature distribution, SLOR method, Δx = 0.05, iteration level = 1000, w=1.989
21
Fig. 14 Temperature distribution, SLOR method, Δx = 0.0025, iteration level = 1000, w=1.989
Fig. 15 Temperature distribution, SLOR method, Δx = 0.05, iteration level = 1000, w=1.5 22
Fig. 16 Temperature distribution, SLOR method, Δx = 0.05, iteration level = 1000, w=1.2
Fig. 17 Temperature distribution, SLOR method, Δx = 0.05, iteration level = 1000, w=1.15 23
4.3 Discussion Although both methods are stable, careful examination of the ADI method results indicates that: 1. Increasing the time steps would speed up the convergence to final state 2. Increasing the size of the grid prevents us from seeing the accurate temperature distributions that we would see otherwise. 3. With the same grid size and time and same time level, bigger time steps yield different results. 4. The code is biased from left to right, in other words the boundary conditions from the left side of the domain propagate through the domain faster than the boundary conditions on the right. 5. The computation time does not increase significantly with time steps that are in the order of convergence time step. It must be mentioned that the time step n was changed manually to observe temperature variations but the code converged at 976 steps.
The SLOR method, on the other hand, presented physically unacceptable results. It can be observed that: 1. At the same iteration level, using smaller mesh net results in the diminishment of the region that has a higher temperature than that of the boundaries. 2. Although the optimum relaxation factor calculated according to [1] is 1.989, reducing the relaxation factor results in more acceptable results in terms of the physical reality of the problem. 3. Increasing the iteration levels exponentially increases the computation time, although even at 1000 iterations, the code has not converged yet. Therefore, by comparing the results obtained from ADI and SLOR methods it is found that the ADI method, although a little bit harder to code, 1. Yields realistic results 2. Significantly reduces the computation time 3. Converges much more quickly
24
5. References [1] Computational Fluid Mechanics and Heat Transfer, John C. Tannehil, Dale A. Anderson, Richard H. Pletcher, Seond Edition, Taylor and Francis, 1997 [2] Computational Fluid Dynamics: The Basics with Applications, John D. Anderson, McGraw-Hill, 1995
25
