* 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) *
* FORWARD EULER (FTCS) METHOD *
*************************************************************************
*
program FTCS2
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,*) boundary_type ! flag for type of boundary conditions
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
if (boundary_type == 1) then
write(21,*) 'Dirichlet bc is specified on all boundaries'
elseif (boundary_type == 2) then
write(21,*) 'Dirichlet on W,E,N Neumann on S'
elseif (boundary_type == 3) then
write(21,*) 'Neumann on W,E,N Dirichlet on S'
endif
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
*
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: 1 - dirichlet: all
c boundary_type: 2 - dirichlet: W,E,N; neumann: S;
c boundary_type: 3 - neumann: W,E,N; dirichlet: S;
c
if (boundary_type == 1) then
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))
TMP(m,1) = BTE(x(m),y(1))
TMP(1,n) = BTW(x(1),y(n))
TMP(m,n) = BTE(x(m),y(n))
elseif (boundary_type == 2) then
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
qs = BQS(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))
TMP(m,1) = BTE(x(m),y(1))
TMP(1,n) = BTW(x(1),y(n))
TMP(m,n) = BTE(x(m),y(n))
elseif (boundary_type == 3) then
c west and east
do j = 2, n-1
qw = BQW(x(1),y(j))
qe = BQE(x(m),y(j))
enddo
c south and north
do i = 2, m-1
TMP(i,1) = BTS(x(i),y(1))
qn = BQN(x(i),y(n))
enddo
c corner boundary nodes
TMP(1,1) = BTS(x(1),y(1))
TMP(m,1) = BTS(x(m),y(1))
endif
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'
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), maxdif
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
TN = TMP
c inclusion of source term
do j = 2, n-1
do i = 2, m-1
C(i,j) = S(x(i),y(j))
enddo
enddo
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
10 if ((maxdif > tolstdy).and.(ntimestp < maxntimestp)) then
do j = 3, n-2
do i = 3, m-2
dTdx_west = (TN(i,j) - TN(i-1,j))/dx
dTdx_east = (TN(i+1,j) - TN(i,j))/dx
dTdy_south = (TN(i,j) - TN(i,j-1))/dy
dTdy_north = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_east - dTdx_west)/dx +
$ (dTdy_north - dTdy_south)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
enddo
enddo
c
c incorporation of boundary conditions:
c
c boundary_type: 1 - dirichlet: all
c boundary_type: 2 - dirichlet: W,E,N; neumann: S;
c boundary_type: 3 - neumann: W,E,N; dirichlet: S;
c
if ((boundary_type == 1).or.(boundary_type == 2)) then
c west
i = 2
do j = 3, n-2
dTdx_w = (-TN(i+1,j) + 9*TN(i,j) - 8*TN(i-1,j))/(3*dx)
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
enddo
c east
i = m-1
do j = 3, n-2
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = (8*TN(i+1,j) - 9*TN(i,j) + TN(i-1,j))/(3*dx)
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
enddo
c south
j = 2
do i = 3, m-2
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (-TN(i,j+1) + 9*TN(i,j) - 8*TN(i,j-1))/(3*dy)
if (boundary_type == 2) then
dTdy_s = qs
endif
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
enddo
c north
j = n-1
do i = 3, m-2
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = (8*TN(i,j+1) - 9*TN(i,j) + TN(i,j-1))/(3*dy)
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
enddo
c south-west
i = 2
j = 2
dTdx_w = (-TN(i+1,j) + 9*TN(i,j) - 8*TN(i-1,j))/(3*dx)
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (-TN(i,j+1) + 9*TN(i,j) - 8*TN(i,j-1))/(3*dy)
if (boundary_type == 2) then
dTdy_s = qs
endif
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
c south-east
i = m-1
j = 2
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = (8*TN(i+1,j) - 9*TN(i,j) + TN(i-1,j))/(3*dx)
dTdy_s = (-TN(i,j+1) + 9*TN(i,j) - 8*TN(i,j-1))/(3*dy)
if (boundary_type == 2) then
dTdy_s = qs
endif
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
c north-west
i = 2
j = n-1
dTdx_w = (-TN(i+1,j) + 9*TN(i,j) - 8*TN(i-1,j))/(3*dx)
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = (8*TN(i,j+1) - 9*TN(i,j) + TN(i,j-1))/(3*dy)
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
c north-east
i = m-1
j = n-1
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = (8*TN(i+1,j) - 9*TN(i,j) + TN(i-1,j))/(3*dx)
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = (8*TN(i,j+1) - 9*TN(i,j) + TN(i,j-1))/(3*dy)
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
elseif (boundary_type == 3) then
c west
i = 2
do j = 3, n-2
dTdx_w = qw
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
enddo
c east
i = m-1
do j = 3, n-2
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = qe
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
enddo
c south
j = 2
do i = 3, m-2
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (-TN(i,j+1) + 9*TN(i,j) - 8*TN(i,j-1))/(3*dy)
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
enddo
c north
j = n-1
do i = 3, m-2
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = qn
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
enddo
c south-west
i = 2
j = 2
dTdx_w = qw
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (-TN(i,j+1) + 9*TN(i,j) - 8*TN(i,j-1))/(3*dy)
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
c south-east
i = m-1
j = 2
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = qe
dTdy_s = (-TN(i,j+1) + 9*TN(i,j) - 8*TN(i,j-1))/(3*dy)
dTdy_n = (TN(i,j+1) - TN(i,j))/dy
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
c north-west
i = 2
j = n-1
dTdx_w = qw
dTdx_e = (TN(i+1,j) - TN(i,j))/dx
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = qn
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
c north-east
i = m-1
j = n-1
dTdx_w = (TN(i,j) - TN(i-1,j))/dx
dTdx_e = qe
dTdy_s = (TN(i,j) - TN(i,j-1))/dy
dTdy_n = qn
dTdt = alpha*((dTdx_e - dTdx_w)/dx +
$ (dTdy_n - dTdy_s)/dy) + C(i,j)
TMP(i,j) = TN(i,j) + dTdt*dt
endif
c
time = time + dt
ntimestp = ntimestp + 1
c
maxdif = 0
do j = 2, n-1
do i = 2, m-1
dif = TMP(i,j) - TN(i,j)
if (dabs(dif) > maxdif) then
maxdif = dabs(dif)
endif
enddo
enddo
c
c updating the temperature
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
c calculation of boundary temperatures where Neumann BC is specified
if (boundary_type == 1) then
!
elseif (boundary_type == 2) then
c south
j = 1
do i = 1, m
TMP(i,j) = (9*TMP(i,j+1) - TMP(i,j+2) - 3*qs*dy) /8d0
enddo
elseif (boundary_type == 3) then
c west
i = 1
do j = 1, n
TMP(i,j) = (9*TMP(i+1,j) - TMP(i+2,j) - 3*qw*dx) /8d0
enddo
c east
i = m
do j = 1, n
TMP(i,j) = (9*TMP(i-1,j) - TMP(i-2,j) + 3*qe*dx) /8d0
enddo
c north
j = n
do i = 1, m
TMP(i,j) = (9*TMP(i,j-1) - TMP(i,j-2) + 3*qn*dy) /8d0
enddo
endif
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
|
