* A. Salih, Dept. of Mechanical Engg., NIT - Trichy, India. *
*************************************************************************
* THIS PROGRAM SOLVES TRANSIENT HEAT EQUATION IN A 2-D DOMAIN *
* USING FINITE VOLUME METHOD (UNIFORM CV) *
* ALTERNATING DIRECTION IMPLICIT METHOD (PEACEMAN and RACHFORD) *
*************************************************************************
*
program ADI
implicit doubleprecision (a-h,o-z)
include 'scalars/integers.inc'
include 'scalars/parameters.inc'
open(unit=11,file='input.dat',status='unknown')
open(unit=21,file='output.dat',status='unknown')
open(unit=31,file='outputs/x.dat',status='unknown')
open(unit=32,file='outputs/y.dat',status='unknown')
open(unit=33,file='outputs/tmp.dat',status='unknown')
open(unit=34,file='outputs/time.dat',status='unknown')
open(unit=35,file='outputs/time2.dat',status='unknown')
open(unit=36,file='outputs/t_trans.dat',status='unknown')
c
pi = 4*datan(1d0)
c
CALL READ_IN ! read-in data
close (unit=11)
CALL GRID ! setting up the grid points in the domain
CALL INITIAL_COND ! setting up initial condition
CALL BOUNDARY_COND ! setting up boundary values
CALL SOLVE ! setting up system of equations
CALL PRINTOUT ! printing out the computational results
CALL WRITE_OUT ! write-out the basic data used for computation
c
stop
end
*
***SUBROUTINE: READING IN THE BASIC DATA *******************************
*
subroutine READ_IN
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
c
read(11,*) alx ! actual length of domain
read(11,*) aly ! actual length of domain
read(11,*) alpha ! thermal diffusivity
read(11,*) Fo ! value of mean grid Fourier number for
! calculating the time step
read(11,*) tolstdy ! convergence criterion for steady-state
read(11,*) tfinal ! time at which the solution is desired
read(11,*) print_freq !
c
c calculation of time-step dt based on the stability condition
dx = alx /(m-2)
dy = aly /(n-2)
dt = 2*Fo *(dx**2 * dy**2)/(dx**2 + dy**2) /alpha
maxntimestp = tfinal /dt
dt = tfinal /maxntimestp ! revised time-step
Fox = alpha*dt/dx**2
Foy = alpha*dt/dy**2
Fo = (Fox + Foy)/2d0 ! revised mean Fourier number
if (tfinal < dt) then
print*, 'Warning: final time is less than time-step!'
stop
endif
c
return
end
*
***SUBROUTINE: WRITING OUT THE BASIC DATA USED FOR COMPUTATION**********
*
subroutine WRITE_OUT
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/boundaryvalues.inc'
c
write(21,*) 'Dirichlet bc is specified on all boundaries'
write(21,*) 'Thermal diffusivity', alpha
write(21,*) 'Mean grid Fourier number, Fo =', Fo
write(21,*) 'Length of the domain in x-direction', alx
write(21,*) 'Length of the domain in y-direction', aly
write(21,*) 'Number of control volumes:', nx, ' X', ny
write(21,*) 'dx =', dx
write(21,*) 'dy =', dy
write(21,*) 'dt =', dt
write(21,*) 'No. of Time step = ', ntimestp
write(21,*) 'Final time = ', time
write(21,*) 'convergence criterion for steady-state solution:',
$tolstdy
c
return
end
*
****SUBROUTINE: GRID****************************************************
*
SUBROUTINE GRID
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'arrays/grid.inc'
c
c m-2 = number of control volumes the in x-direction
c n-2 = number of control volumes the in y-direction
dx = alx /(m-2)
x(1) = 0
x(2) = x(1) + dx/2d0
do i = 3, m-1
x(i) = x(i-1) + dx
enddo
x(m) = x(m-1) + dx/2d0
c
dy = aly /(n-2)
y(1) = 0
y(2) = y(1) + dy/2d0
do i = 3, n-1
y(i) = y(i-1) + dy
enddo
y(n) = y(n-1) + dy/2d0
c
return
end
*
***SUBROUTINE: BOUNDARY CONDITIONS FOR TEMPERATURE**********************
*
SUBROUTINE BOUNDARY_COND
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/boundaryvalues.inc'
include 'scalars/parameters.inc'
include 'arrays/temperature.inc'
include 'arrays/grid.inc'
c
c defining the boundary conditions
include 'bc/bc.inc'
c
c boundary_type: - dirichlet: all
c west and east
do j = 2, n-1
TMP(1,j) = BTW(x(1),y(j))
TMP(m,j) = BTE(x(m),y(j))
enddo
c south and north
do i = 2, m-1
TMP(i,1) = BTS(x(i),y(1))
TMP(i,n) = BTN(x(i),y(n))
enddo
c corner boundary nodes
TMP(1,1) = (BTW(x(1),y(1)) + BTS(x(1),y(1)))/2d0
TMP(m,1) = (BTE(x(m),y(1)) + BTS(x(m),y(1)))/2d0
TMP(1,n) = (BTW(x(1),y(n)) + BTN(x(1),y(n)))/2d0
TMP(m,n) = (BTE(x(m),y(n)) + BTN(x(m),y(n)))/2d0
c
return
end
*
***SUBROUTINE: INITIAL CONDITIONS FOR TEMPERATURE**********************
*
SUBROUTINE INITIAL_COND
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/parameters.inc'
include 'arrays/temperature.inc'
include 'arrays/grid.inc'
c
c defining the initial conditions
TINIT(x,y,t) = 0 ! dsin(pi*x/alx)
c
do j = 2, n-1
do i = 2, m-1
TMP(i,j) = TINIT(x(i),y(j),0)
enddo
enddo
c
return
end
*
***SUBROUTINE: CONSTRUCTION OF SYSTEM OF EQUATION AND ITS SOLUTION******
*
subroutine SOLVE
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
parameter( mn = max(m,n)+1 )
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/boundaryvalues.inc'
include 'scalars/parameters.inc'
include 'arrays/grid.inc'
include 'arrays/temperature.inc'
include 'arrays/coeff.inc'
doubleprecision TN(0:m1,0:n1), TNN(0:m1,0:n1), maxdif
doubleprecision T(0:m1), TT(0:m1)
c
include 'arithmeticfunctions/source.inc' ! source.inc file
! contains the definition of source term
c
c storing the values of temperature at nth time level
TNN = TMP
TN = TMP
c
ntimestp = 0
maxdif = 1d9
time = 0d0
write(34,101) time ! continuous time
write(35,101) time ! periodic time
monitorx = nx/4
monitory = ny/4
write(36,*) TMP(monitorx,monitory)
c
Fox = Fox/2d0
Foy = Foy/2d0
10 if ((maxdif > tolstdy).and.(ntimestp < maxntimestp)) then
c row-sweep
c south control volumes
j = 2
c south-west
i = 2
AW(i) = -8d0/3d0 *Fox
AE(i) = -4d0/3d0 *Fox
AP(i) = 1d0 - (AW(i) + AE(i))
C(i) = TN(i,j) + Foy*(8*TN(i,j-1) - 12*TN(i,j) +
$ 4*TN(i,j+1)) + S(x(i),y(j))*dt/2d0
C(i) = C(i) - AW(i)*TMP(i-1,j)
AW(i) = 0
c south
do i = 3, m-2
AW(i) = -Fox
AE(i) = -Fox
AP(i) = 1d0 - (AW(i) + AE(i))
C(i) = TN(i,j) + Foy*(8*TN(i,j-1) - 12*TN(i,j) +
$ 4*TN(i,j+1)) + S(x(i),y(j))*dt/2d0
enddo
c south-east
i = m-1
AW(i) = -4d0/3d0 *Fox
AE(i) = -8d0/3d0 *Fox
AP(i) = 1d0 - (AW(i) + AE(i))
C(i) = TN(i,j) + Foy*(8*TN(i,j-1) - 12*TN(i,j) +
$ 4*TN(i,j+1)) + S(x(i),y(j))*dt/2d0
C(i) = C(i) - AE(i)*TMP(i+1,j)
AE(i) = 0
c
do i = 1, m
T(i) = TMP(i,j)
enddo
CALL TDMA(AW,AP,AE,C,T)
do i = 2, m-1
TMP(i,j) = T(i)
enddo
c
c interior control volumes
do j = 3, n-2
c west
i = 2
AW(i) = -8d0/3d0 *Fox
AE(i) = -4d0/3d0 *Fox
AP(i) = 1d0 - (AW(i) + AE(i))
C(i) = TN(i,j) + Foy*(TN(i,j-1) -2*TN(i,j) + TN(i,j+1)) +
$ S(x(i),y(j))*dt/2d0
C(i) = C(i) - AW(i)*TMP(i-1,j)
AW(i) = 0
c inerior
do i = 3, m-2
AW(i) = -Fox
AE(i) = -Fox
AP(i) = 1d0 - (AW(i) + AE(i))
C(i) = TN(i,j) + Foy*(TN(i,j-1) - 2*TN(i,j) +
$ TN(i,j+1)) + S(x(i),y(j))*dt/2d0
enddo
c east
i = m-1
AW(i) = -4d0/3d0 *Fox
AE(i) = -8d0/3d0 *Fox
AP(i) = 1d0 - (AW(i) + AE(i))
C(i) = TN(i,j) + Foy*(TN(i,j-1) -2*TN(i,j) + TN(i,j+1)) +
$ S(x(i),y(j))*dt/2d0
C(i) = C(i) - AE(i)*TMP(i+1,j)
AE(i) = 0
c
do i = 1, m
T(i) = TMP(i,j)
enddo
CALL TDMA(AW,AP,AE,C,T)
do i = 2, m-1
TMP(i,j) = T(i)
enddo
enddo
c north control volumes
j = n-1
c north-west
i = 2
AW(i) = -8d0/3d0 *Fox
AE(i) = -4d0/3d0 *Fox
AP(i) = 1d0 - (AW(i) + AE(i))
C(i) = TN(i,j) + Foy*(4*TN(i,j-1) - 12*TN(i,j) +
$ 8*TN(i,j+1)) + S(x(i),y(j))*dt/2d0
C(i) = C(i) - AW(i)*TMP(i-1,j)
AW(i) = 0
c north
do i = 3, m-2
AW(i) = -Fox
AE(i) = -Fox
AP(i) = 1d0 - (AW(i) + AE(i))
C(i) = TN(i,j) + Foy*(4*TN(i,j-1) - 12*TN(i,j) +
$ 8*TN(i,j+1)) + S(x(i),y(j))*dt/2d0
enddo
c north-east
i = m-1
AW(i) = -4d0/3d0 *Fox
AE(i) = -8d0/3d0 *Fox
AP(i) = 1d0 - (AW(i) + AE(i))
C(i) = TN(i,j) + Foy*(4*TN(i,j-1) - 12*TN(i,j) +
$ 8*TN(i,j+1)) + S(x(i),y(j))*dt/2d0
C(i) = C(i) - AE(i)*TMP(i+1,j)
AE(i) = 0
c
do i = 1, m
T(i) = TMP(i,j)
enddo
CALL TDMA(AW,AP,AE,C,T)
do i = 2, m-1
TMP(i,j) = T(i)
enddo
c updating the temperature for column sweep
TN = TMP
c column-sweep
c west control volumes
i = 2
c south-west
j = 2
AS(j) = -8d0/3d0 *Foy
AN(j) = -4d0/3d0 *Foy
APP(j) = 1d0 - (AS(j) + AN(j))
CC(j) = TN(i,j) + Fox*(8*TN(i-1,j) - 12*TN(i,j) +
$ 4*TN(i+1,j)) + S(x(i),y(j))*dt/2d0
CC(j) = CC(j) - AS(j)*TMP(i,j-1)
AS(j) = 0
c west
do j = 3, n-2
AS(j) = -Foy
AN(j) = -Foy
APP(j) = 1d0 - (AS(j) + AN(j))
CC(j) = TN(i,j) + Fox*(8*TN(i-1,j) - 12*TN(i,j) +
$ 4*TN(i+1,j)) + S(x(i),y(j))*dt/2d0
enddo
c north-west
j = n-1
AS(j) = -4d0/3d0 *Foy
AN(j) = -8d0/3d0 *Foy
APP(j) = 1d0 - (AS(j) + AN(j))
CC(j) = TN(i,j) + Fox*(8*TN(i-1,j) - 12*TN(i,j) +
$ 4*TN(i+1,j)) + S(x(i),y(j))*dt/2d0
CC(j) = CC(j) - AN(j)*TMP(i,j+1)
AN(j) = 0
c
do j = 1, n
TT(j) = TMP(i,j)
enddo
CALL TDMA(AS,APP,AN,CC,TT)
do j = 2, n-1
TMP(i,j) = TT(j)
enddo
c
c interior control volumes
do i = 3, m-2
c south
j = 2
AS(j) = -8d0/3d0 *Foy
AN(j) = -4d0/3d0 *Foy
APP(j) = 1d0 - (AS(j) + AN(j))
CC(j) = TN(i,j) + Fox*(TN(i-1,j) -2*TN(i,j) + TN(i+1,j)) +
$ S(x(i),y(j))*dt/2d0
CC(j) = CC(j) - AS(j)*TMP(i,j-1)
AS(j) = 0
c inerior
do j = 3, n-2
AS(j) = -Foy
AN(j) = -Foy
APP(j) = 1d0 - (AS(j) + AN(j))
CC(j) = TN(i,j) + Fox*(TN(i-1,j) - 2*TN(i,j) +
$ TN(i+1,j)) + S(x(i),y(j))*dt/2d0
enddo
c north
j = n-1
AS(j) = -4d0/3d0 *Foy
AN(j) = -8d0/3d0 *Foy
APP(j) = 1d0 - (AS(j) + AN(j))
CC(j) = TN(i,j) + Fox*(TN(i-1,j) -2*TN(i,j) + TN(i+1,j)) +
$ S(x(i),y(j))*dt/2d0
CC(j) = CC(j) - AN(j)*TMP(i,j+1)
AN(j) = 0
c
do j = 1, n
TT(j) = TMP(i,j)
enddo
CALL TDMA(AS,APP,AN,CC,TT)
do j = 2, n-1
TMP(i,j) = TT(j)
enddo
enddo
c north control volumes
i = m-1
c south-east
j = 2
AS(j) = -8d0/3d0 *Foy
AN(j) = -4d0/3d0 *Foy
APP(j) = 1d0 - (AS(j) + AN(j))
CC(j) = TN(i,j) + Fox*(4*TN(i-1,j) - 12*TN(i,j) +
$ 8*TN(i+1,j)) + S(x(i),y(j))*dt/2d0
CC(j) = CC(j) - AS(j)*TMP(i,j-1)
AS(j) = 0
c east
do j = 3, n-2
AS(j) = -Foy
AN(j) = -Foy
APP(j) = 1d0 - (AS(j) + AN(j))
CC(j) = TN(i,j) + Fox*(4*TN(i-1,j) - 12*TN(i,j) +
$ 8*TN(i+1,j)) + S(x(i),y(j))*dt/2d0
enddo
c north-east
j = n-1
AS(j) = -4d0/3d0 *Foy
AN(j) = -8d0/3d0 *Foy
APP(j) = 1d0 - (AS(j) + AN(j))
CC(j) = TN(i,j) + Fox*(4*TN(i-1,j) - 12*TN(i,j) +
$ 8*TN(i+1,j)) + S(x(i),y(j))*dt/2d0
CC(j) = CC(j) - AN(j)*TMP(i,j+1)
AN(j) = 0
c
do j = 1, n
TT(j) = TMP(i,j)
enddo
CALL TDMA(AS,APP,AN,CC,TT)
do j = 2, n-1
TMP(i,j) = TT(j)
enddo
c
time = time + dt
ntimestp = ntimestp + 1
c
maxdif = 0
do j = 2, n-1
do i = 2, m-1
dif = TMP(i,j) - TNN(i,j)
if (dabs(dif) > maxdif) then
maxdif = dabs(dif)
endif
enddo
enddo
c
c updating the temperatures
TNN = TMP
TN = TMP
c
write(34,101) time ! continuous time
write(6,102) ntimestp, maxdif
c
c printing out transient temperature at selected point of domain
j = jcounter + 1
jcounter = mod(j,print_freq)
if (jcounter == 0) then
write(35,101) (time + dt) ! periodic time
write(36,*) TMP(monitorx,monitory)
endif
go to 10
endif
c
101 format(e14.7)
102 format(i7,5x,e14.7)
c
return
end
*
***SUBROUTINE: TDMA*****************************************************
*
subroutine TDMA(L,D,U,C,X)
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
doubleprecision L(0:m1),D(0:m1),U(0:m1),C(0:m1),X(0:m1)
doubleprecision P(-1:m1),Q(-1:m1)
c
P(1) = 0d0
Q(1) = 0d0
c
c forward elimination
do i = 2, m-1
denom = D(i) + L(i)*P(i-1)
P(i) = -U(i) /denom
Q(i) = (C(i) - L(i)*Q(i-1)) /denom
enddo
c
c back substitution
do i = m-1, 2, -1
X(i) = P(i)*X(i+1) + Q(i)
enddo
c
return
end
*
***SUBROUTINE: TDMA*****************************************************
*
subroutine TDMA2(L,D,U,C,X)
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
doubleprecision L(0:n1),D(0:n1),U(0:n1),C(0:n1),X(0:n1)
doubleprecision P(-1:n1),Q(-1:n1)
c
P(1) = 0d0
Q(1) = 0d0
c
c forward elimination
do i = 2, n-1
denom = D(i) + L(i)*P(i-1)
P(i) = -U(i) /denom
Q(i) = (C(i) - L(i)*Q(i-1)) /denom
enddo
c
c back substitution
do i = n-1, 2, -1
X(i) = P(i)*X(i+1) + Q(i)
enddo
c
return
end
*
***SUBROUTINE: Printing-out the output data*****************************
*
SUBROUTINE PRINTOUT
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'arrays/temperature.inc'
include 'arrays/grid.inc'
c
write(31,101) (x(i), i=1,m)
write(32,101) (y(i), i=1,n)
write(33,102) ((TMP(i,j), i=1,m), j=1,n)
101 format(e13.6)
102 format(34(e13.6,1x))! format number = no. of CV in x-direction + 2
c
return
end
|
