PROGRAM normic ! Computes a uniformly-spaced, 3d grid for a 2.0 ft long straight ! duct with a diameter of 1.0 ft and initializes a flow with a normal ! shock midway through the duct. ! !----------------------------------------------------------------------- Implicit none Integer, parameter :: nim = 202 Integer, parameter :: njm = 11 Integer, parameter :: nkm = 5 Integer i, j, k, m Integer nc1 Integer ni, nj, nk Integer nblock Real q (nim,njm,nkm,5) Real x (nim,njm,nkm) Real y (nim,njm,nkm) Real z (nim,njm,nkm) Real a1, et1, m1, p1, r1, t1, u1, v1, w1 Real a2, et2, m2, p2, r2, t2, u2, v2, w2 Real gam, gm1, gp1 Real gascon Real pi Real alpha, reyn, time Real xl, xr, xd, rad Real thet Real p2p1, r2r1, u2u1 !----------------------------------------------------------------------- pi = 4.0 * atan(1.0) !..Set gas variables. gam = 1.4 gp1 = gam + 1.0 gm1 = gam - 1.0 gascon = 1716.0 !..Set conditions on supersonic side of shock. m1 = 1.3 p1 = 10.0 * 144.0 t1 = 520.0 r1 = p1 / ( t1 * gascon ) a1 = sqrt( gam * gascon * t1 ) u1 = m1 * a1 v1 = 0.0 w1 = 0.0 et1 = p1 / gm1 + 0.5 * r1 * u1 * u1 !..Create conditions on subsonic side of shock. r2r1 = gp1 * m1*m1 / ( 2.0 + gm1 * m1*m1 ) u2u1 = 1.0 / r2r1 p2p1 = 1.0 + 2.0*gam*(m1*m1-1.0)/gp1 m2 = sqrt( ( 1.0 + 0.5*gm1*m1*m1 ) / ( gam*m1*m1 - 0.5*gm1 ) ) p2 = p2p1 * p1 r2 = r2r1 * r1 u2 = u2u1 * u1 v2 = 0.0 w2 = 0.0 et2 = p2 / gm1 + 0.5 * r2 * u2 * u2 !..Set duct size and shock location. xl = 0.0 xr = 2.0 xd = 1.0 rad = 0.5 !..Generate an evenly spaced, 3D grid. ni = 202 nj = 11 nk = 5 do i = 1, ni do j = 1, nj do k = 1, nk thet = 10.0 * ( k - 1.0 ) / ( nk - 1.0 ) * pi / 180.0 x(i,j,k) = xl + ( xr - xl ) * ( i - 1.0 ) / ( ni - 1.0 ) y(i,j,k) = rad * ( j - 1.0 ) / ( nj - 1.0 ) * cos(thet) z(i,j,k) = rad * ( j - 1.0 ) / ( nj - 1.0 ) * sin(thet) enddo enddo enddo open ( unit=7, file='norm.x', form='unformatted', status='unknown' ) nblock = 1 write(7) nblock write(7) ni, nj, nk write(7) ((( x(i,j,k), i=1,ni), j=1,nj), k=1,nk), & ((( y(i,j,k), i=1,ni), j=1,nj), k=1,nk), & ((( z(i,j,k), i=1,ni), j=1,nj), k=1,nk) !..Load the initial solution using state 1 as reference. do i = 1, ni do j = 1, nj do k = 1, nk if ( x(i,j,k) .lt. xd ) then q(i,j,k,1) = r1 / r1 q(i,j,k,2) = r1 * u1 / ( r1 * a1 ) q(i,j,k,3) = r1 * v1 / ( r1 * a1 ) q(i,j,k,4) = r1 * w1 / ( r1 * a1 ) q(i,j,k,5) = et1 / ( r1 * a1 * a1 ) else q(i,j,k,1) = r2 / r1 q(i,j,k,2) = r2 * u2 / ( r1 * a1 ) q(i,j,k,3) = r2 * v2 / ( r1 * a1 ) q(i,j,k,4) = r2 * w2 / ( r1 * a1 ) q(i,j,k,5) = et2 / ( r1 * a1 * a1 ) endif enddo enddo enddo !..Output the solution to file norm.q (3D plot3d unformated) open ( unit=8, file='norm.q', form='unformatted', status='unknown' ) alpha = 0.0 reyn = 1.0 time = 0.0 write (8) nblock write (8) ni, nj, nk write (8) m1, alpha, reyn, time write (8) (((( q(i,j,k,m), i=1,ni), j=1,nj), k=1,nk), m=1,5 ) !..Output the NPARC restart file. nc1 = 0 write (2) nc1 write (2) ni, nj, nk write (2) ((( x(i,j,k), i=1,ni), j=1,nj), k=1,nk), & ((( y(i,j,k), i=1,ni), j=1,nj), k=1,nk), & ((( z(i,j,k), i=1,ni), j=1,nj), k=1,nk) write (2) (((( q(i,j,k,m), i=1,ni), j=1,nj), k=1,nk), m=1,5 ) !----------------------------------------------------------------------- stop end