* A. Salih, Dept. of Mechanical Engg., NIT - Trichy, India. *
*************************************************************************
* PROGRAM FOR STEADY, ONE-DIMENSIONAL CONVECTION-DIFFUSION EQUATION *
* USING FINITE VOLUME METHOD (UNIFORM CV) *
*************************************************************************
*
program CONVEC_DIFFUSE_1D
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/t.dat',status='unknown')
c
pi = 4*datan(1d0)
c
CALL READ_IN ! read-in and write-out data
close (unit=11)
CALL GRID ! setting up the grid points in the domain
CALL BOUNDARY_COND ! setting up boundary values
CALL SOLVE ! setting up system of equations
CALL PRINTOUT ! printing out the 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,*) u ! value of advection velocity
read(11,*) alpha ! thermal diffusivity
read(11,*) conv_scheme ! flag for convection discretization
! 1-CDS, 2-Upwind scheme
c
c calculation of Peclet number grid Peclet number
dx = alx /(m-1)
PeL = u*alx/alpha
Pe = u*dx/alpha
if ((conv_scheme == 1).and.(Pe > 2d0)) then
print*, 'Warning: The grid Peclet number is > 2, wiggles may ap
$pear in the solution.'
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'
c
write(21,*) 'Linear, 1-dimensional, steady convection-diffusion:'
c
if (conv_scheme == 1) then
write(21,*) 'CDS is used for discretization of convection term.'
elseif (conv_scheme == 2) then
write(21,*) 'Upwind scheme is used for discretization of convec
$tion term.'
endif
c
write(21,*) 'Advection velocity, u =', u
write(21,*) 'thermal diffusivity =', alpha
write(21,*) 'Peclet number, Pe =', PeL
write(21,*) 'Grid Peclet number =', Pe
write(21,*) 'Number of control volumes:', mx
write(21,*) 'dx =', dx
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
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
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
TMP(1) = tw
TMP(m) = te
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'
c
c CDS - CDS
if (conv_scheme == 1) then
c construction of coefficient matrix
do i = 2, m-1
AW(i) = 1d0 + Pe/2
AE(i) = 1d0 - Pe/2
AP(i) = -(AW(i) + AE(i))
C(i) = 0
enddo
c incorporation of boundary conditions
c west
AW(2) = 8d0/3d0 + Pe
AE(2) = 4d0/3d0 - Pe/2
AP(2) = -(AW(2) + AE(2))
C(2) = C(2) - AW(2)*TMP(1)
AW(2) = 0
c east
AW(m-1) = 4d0/3d0 + Pe/2
AE(m-1) = 8d0/3d0 - Pe
AP(m-1) = -(AW(m-1) + AE(m-1))
C(m-1) = C(m-1) - AE(m-1)*TMP(m)
AE(m-1) = 0
elseif (conv_scheme == 2) then
c construction of coefficient matrix
do i = 2, m-1
AW(i) = 1d0 + max(Pe,0d0)
AE(i) = 1d0 - min(Pe,0d0)
AP(i) = -(AW(i) + AE(i))
C(i) = 0
enddo
c incorporation of boundary conditions
c west
AW(2) = 8d0/3d0 + max(Pe,0d0)
AE(2) = 4d0/3d0 - min(Pe,0d0)
AP(2) = -(AW(2) + AE(2))
C(2) = C(2) - AW(2)*TMP(1)
AW(2) = 0
c east
AW(m-1) = 4d0/3d0 + max(Pe,0d0)
AE(m-1) = 8d0/3d0 - min(Pe,0d0)
AP(m-1) = -(AW(m-1) + AE(m-1))
C(m-1) = C(m-1) - AE(m-1)*TMP(m)
AE(m-1) = 0
endif
c solution of system of equations
CALL TDMA(AW,AP,AE,C,TMP)
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),U(0:m1),D(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: 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,102) (TMP(i), i=1,m)
101 format(e13.6)
102 format(e13.6)
c
return
end
|
