7 Mar 2013 1:51
#11670
Hi all can you please guide me on following I am using in CFD course I am getting error stating Error 94, Unit has neither been OPENed nor preconnected
!!!!!!!******* MAE 561 Homework 3 *******!!!!!!!!!
!!!!!!!******* Fortran95 Code for Problem 2 *******!!!!!!!!!
!!!!!!!******* ANANDRAO BIRADAR, 1206263652 *******!!!!!!!!!
program problem_2_FTCS_BTCS_explicit_method
implicit none
real*8, dimension (0:40,0:540) :: u, a, b, c
real*8, dimension (0:40) :: y
real*8, dimension (39,0:540) :: d
real*8, dimension (3) :: p
real*8 :: h, v, den, z, q
integer :: i, k, j, m
den = 800.0d0 !!!!! Density = d
z = 0.002 !!!!! z = delta(t)
h = 0.001 !!!!! grid spacing
!!!! Applying initial condition !!!!
y(0) = 0 !! boudary conditions for y
y(40)= 0.04
do i = 1,39
y(i) = y(i-1)+h
end do
open (10, file = 'Values of y.txt')
write (10,1) y
1 format(1f40.8)
close (10)
do i = 1,40
u(i,0) = 0
end do
v = 0.000217 !!!! v = kinematic viscosity
do i = 0,540
u(0,i) = 40.0d0 !!!! As lower plate start moving with u = 40 m/s
u(40,i)= 0 !!!! Upper plate is stationary
end do
!!!!!!!!!!!!!!!!!*************** Using FTCS Method with given conditions ****************!!!!!!!!!!!!!!!!
do j = 1,3
p(1) = 0.0d0 !!!!! p = dp/dx
p(2) = 20000.0d0
p(3) = -30000.0d0
do k = 0,539
do i = 1,39
u(i,k+1)= u(i,k)+((v*z)/h**2)*(u(i+1,k)-2.0d0*u(i,k)+u(i-1,k))-((z/den)*p(j)) !!!!!!!! The FTCS Formula
end do
!!!!!!!!! Saving the results in different files !!!!!!
if (j == 1) then !!!!!!! for dp/dx = 0
if ( k == 0 ) then
print*, 'At t = 0.00'
open (20, file = 'Values of u at dp-dx = 0.txt')
do i = 0,40
write (20,2) u(i,k)
end do
2 format(1f40.20)
else if ( k == 89 ) then
print*, 'At t = 0.18'
do i = 0,40
write (20,2) u(i,k)
end do
else if ( k == 179 ) then
print*, 'At t = 0.36'
do i = 0,40
write (20,2) u(i,k)
end do
else if ( k == 269 ) then
print*, 'At t = 0.54'
do i = 0,40
write (20,2) u(i,k)
end do
else if ( k == 359 ) then
print*, 'At t = 0.72'
do i = 0,40
write (20,2) u(i,k)
end do
else if ( k == 449 ) then
print*, 'At t = 0.9'
do i = 0,40
write (20,2) u(i,k)
end do
else if ( k == 539 ) then
print*, 'At t = 1.08'
do i = 0,40
write (20,2) u(i,k)
end do
close(20)
end if
else if (j == 2) then !!!!!!! for dp/dx = 2.0*10^4
if ( k == 0 ) then
print*, 'At t = 0.00'
open (30, file = 'Values of u at dp-dx = 20000.txt')
do i = 0,40
write (30,2) u(i,k)
end do
else if ( k == 89 ) then
print*, 'At t = 0.18'
do i = 0,40
write (30,2) u(i,k)
end do
else if ( k == 179 ) then
print*, 'At t = 0.36'
do i = 0,40
write (30,2) u(i,k)
end do
else if ( k == 269 ) then
print*, 'At t = 0.54'
do i = 0,40
write (30,2) u(i,k)
end do
else if ( k == 359 ) then
print*, 'At t = 0.72'
do i = 0,40
write (30,2) u(i,k)
end do
else if ( k == 449 ) then
print*, 'At t = 0.9'
do i = 0,40
write (30,2) u(i,k)
end do
else if ( k == 539 ) then
print*, 'At t = 1.08'
do i = 0,40
write (30,2) u(i,k)
end do
close(30)
end if
else if (j == 3) then !!!!!!! for dp/dx = - 3.0*10^4
if ( k == 0 ) then
print*, 'At t = 0.00'
open (40, file = 'Values of u at dp-dx = - 30000.txt')
do i = 0,40
write (40,2) u(i,k)
end do
else if ( k == 89 ) then
print*, 'At t = 0.18'
do i = 0,40
write (40,2) u(i,k)
end do
else if ( k == 179 ) then
print*, 'At t = 0.36'
do i = 0,40
write (40,2) u(i,k)
end do
else if ( k == 269 ) then
print*, 'At t = 0.54'
do i = 0,40
write (40,2) u(i,k)
end do
else if ( k == 359 ) then
print*, 'At t = 0.72'
do i = 0,40
write (40,2) u(i,k)
end do
else if ( k == 449 ) then
print*, 'At t = 0.9'
do i = 0,40
write (40,2) u(i,k)
end do
else if ( k == 539 ) then
print*, 'At t = 1.08'
do i = 0,40
write (40,2) u(i,k)
end do
close(40)
end if
end if
end do
end do
!$$$$$$ print*, 'Printing the whole matrix'
!$$$$$$ open (541, file = 'Printing the whole matrix.txt')
!$$$$$$ do i = 0,40
!$$$$$$ do k = 0,540
!$$$$$$ read (541,9) u(i,k)
!$$$$$$ end do
!$$$$$$ end do
!$$$$$$ write (541,9) u
!$$$$$$ 9 format(541f40.20)
!$$$$$$ close (541)
!!!!!!!!!!!!!!!!!*************** Using BTCS Method with given conditions ****************!!!!!!!!!!!!!!!!
! z(1) = 0.01 !! z = delta(t) z(2) = 0.002 has already been defined initially
!! Other initial conditions
do i = 1,40
u(i,0) = 0
end do
v = 0.000217 !!!! v = kinematic viscosity
do i = 0,540
u(0,i) = 40.0d0 !!!! As lower plate start moving with u = 40 m/s
u(40,i)= 0 !!!! Upper plate is stationary
end do
do m = 1,3
if ( m == 1 ) then
z = 0.01
q = 0.0 !!!! q = dp/dx
else if ( m == 2) then
z = 0.01
q = 2000.0d0
else if ( m == 3) then
z = 0.002
q = 2000.0d0
end if
!! Defining Elements
do k = 1,539
do i = 0,40
if ( i == 0 ) then
b(i,k) = -(1+(2.0d0*v*z/h**2))
c(i,k) = (v*z)/h**2
else if (i == 40) then
a(i,k) = (v*z)/h**2
b(i,k) = -(1+(2.0d0*v*z/h**2))
else
a(i,k) = (v*z)/h**2
b(i,k) = -(1+(2.0d0*v*z/h**2))
c(i,k) = (v*z)/h**2
end if
end do
!!!! Initiating & calculating d !!!!!
do i = 1,39
if ( i == 1 ) then
d(i,k) = - u(i,k) + (0.002)*(q/den) - a(i,k)*u(i-1,k+1)
else if ( i == 39 ) then
d(i,k) = - u(i,k) + (0.002)*(q/den) - c(i,k)*u(i+1,k+1)
else
d(i,k) = - u(i,k) + (0.002)*(q/den)
end if
/ end do
!!!! Gauss Elimination !!!!
do i = 2,39
b(i,k) = b(i,k) - c(i-1,k)*a(i,k)/b(i-1,k)
d(i,k) = d(i,k) - d(i-1,k)*a(i,k)/b(i-1,k)
end do
!end do
!!!!! Back Substition !!!!
do i = 39 ,1,-1
if (i == 39) then
d(39,k) = d(39,k)/b(39,k)
else
d(i,k) = (d(i,k)-c(i,k)*d(i+1,k))/b(i,k)
end if
end do
do i = 1,39
u(i,k) = d(i,k)
end do
!!!!! Saving the values to different files !!!!!!!!!
!!!!!!! for (I) delta(t) = 0.01s, dp/dx = 0 !!!!!!!
if (m == 1) then
if ( k == 0 ) then
print*, 'At t = 0.00'
open (50, file = 'Values of u at delta(t)=0_dp-dx = 0.txt')
do i = 0,40
write (50,2) u(i,k)
end do
print*,'u(2,89)',u(2,89)
else if ( k == 89 ) then
print*, 'At t = 0.18'
do i = 0,40
write (50,2) u(i,k)
end do
else if ( k == 179 ) then
print*, 'At t = 0.36'
do i = 0,40
write (50,2) u(i,k)
end do
else if ( k == 269 ) then
print*, 'At t = 0.54'
do i = 0,40
write (50,2) u(i,k)
end do
else if ( k == 359 ) then
print*, 'At t = 0.72'
do i = 0,40
write (50,2) u(i,k)
end do
else if ( k == 449 ) then
print*, 'At t = 0.9'
do i = 0,40
write (50,2) u(i,k)
end do
else if ( k == 539 ) then
print*, 'At t = 1.08'
do i = 0,40
write (50,2) u(i,k)
end do
close(50)
end if
!!!!!!! for (II) delta(t) = 0.01s, dp/dx = 2.0*10^4 N/m^2/m !!!!!!!
else if (m == 2) then
if ( k == 0 ) then
print*, 'At t = 0.00'
open (60, file = 'Values of u at delta(t)=0_dp-dx = 2000.txt')
do i = 0,40
write (60,2) u(i,k)
end do
else if ( k == 89 ) then
print*, 'At t = 0.18'
do i = 0,40
write (60,2) u(i,k)
end do
else if ( k == 179 ) then
print*, 'At t = 0.36'
do i = 0,40
write (60,2) u(i,k)
end do
else if ( k == 269 ) then
print*, 'At t = 0.54'
do i = 0,40
write (60,2) u(i,k)
end do
else if ( k == 359 ) then
print*, 'At t = 0.72'
do i = 0,40
write (60,2) u(i,k)
end do
else if ( k == 449 ) then
print*, 'At t = 0.9'
do i = 0,40
write (60,2) u(i,k)
end do
else if ( k == 539 ) then
print*, 'At t = 1.08'
do i = 0,40
write (60,2) u(i,k)
end do
close(60)
end if
!!!!!!! for (III) delta(t) = 0.002s, dp/dx = 2.0*10^4 N/m^2/m !!!!!!!
else if (m == 3) then
if ( k == 0 ) then
print*, 'At t = 0.00'
open (70, file = 'Values of u at delta(t)=0 dp-dx = 2000.txt')
do i = 0,40
write (70,2) u(i,k)
end do
else if ( k == 89 ) then
print*, 'At t = 0.18'
do i = 0,40
write (70,2) u(i,k)
end do
else if ( k == 179 ) then
print*, 'At t = 0.36'
do i = 0,40
write (70,2) u(i,k)
end do
else if ( k == 269 ) then
print*, 'At t = 0.54'
do i = 0,40
write (70,2) u(i,k)
end do
else if ( k == 359 ) then
print*, 'At t = 0.72'
do i = 0,40
write (70,2) u(i,k)
end do
else if ( k == 449 ) then
print*, 'At t = 0.9'
do i = 0,40
write (70,2) u(i,k)
end do
else if ( k == 539 ) then
print*, 'At t = 1.08'
do i = 0,40
write (70,2) u(i,k)
end do
close(70)
end if
end if
end do
end do
end program problem_2_FTCS_BTCS_explicit_method
[/quote]