program quasi_mac c *** This program solves the quasi-1D Euler equations using MacCormack's c *** method. The equations are solved in non-dimensional form with the c *** inlet density, temperature, and velocity as well as the duct length c *** used as scaling parameters. Thus the code solves for c *** rho*=rho/rho1, u*=u/u1, p*=p/(rho1*u1^2), Et*=Et/(rho1*u1^2) c *** as a function of x*=x/L. The duct area is scaled with the throat area, c *** i.e. A*=A/At. parameter( Mx=1000 ) parameter( gamma_gas=1.4 ) real M_inlet, Nx_inv real cons(Mx,3), R(Mx,3), R_old(Mx,3), A(Mx), A_inv(Mx), dAdx(Mx), & rho(Mx), u(Mx), p(Mx), gamma(4), zeta(4), & rho_exact(Mx), u_exact(Mx), p_exact(Mx) c *** Open files and write headers. open(unit=2,file='rms_error.dat') open(unit=3,file='final_error.dat') open(unit=4,file='solution.dat') open(unit=7,file='p.dat') write(2,2) 2 format('# time error rho error u error p') write(3,3) 3 format('# x/L error rho error u error p') write(4,4) 4 format('# x/L rho/rho1 u/u1 p/p1 ') write(7,7) 7 format('# x/L p/p1 (p/p1)exact ') c *** Prompt user for inputs. print *, 'Enter M_inlet, p_ratio, Nx, CFL, and t_max' read(5,*) M_inlet, p_ratio, Nx, CFL, t_max c *** Set parameters. xL = 1.0 dx = xL/float(Nx) dx_inv = 1.0/dx Nx_inv = 1.0/float(Nx) dt = dx*CFL/( 1.0 + 1.0/M_inlet**2 ) nt_max = nint(t_max/dt) dt_print = .05 n_skip = nint(dt_print/dt) nrk_max = 2 gamma(1) = 1.0 zeta( 1) = 0.0 gamma(2) = 0.5 zeta( 2) = -0.5 c0 = gamma_gas - 1.0 c0_inv = 1.0/c0 p_fac = gamma_gas*M_inlet**2 p_fac_inv = 1.0/p_fac p_exit = p_ratio*p_fac_inv c *** Compute and store the exact solution (complete with the shock). call quasi_exact( M_inlet, p_ratio, Nx, rho_exact, u_exact, & p_exact, x_s ) c *** Compute the isentropic solution which will serve as an initial conditon. call quasi_exact( M_inlet, 1.0, Nx, rho, u, p, x_s ) c *** Initialize the conservative variables and the duct area arrays. do i=1, Nx+1 cons(i,1) = rho(i) cons(i,2) = rho(i)*u(i) cons(i,3) = p(i)*c0_inv + 0.5*rho(i)*u(i)**2 x = float(i-1)*dx A(i) = Area(x) A_inv(i) = 1.0/A(i) dAdx(i) = dAreadx(x) end do c *** Enter the time stepping loop, which consists of an outer loop over c *** the number of time steps and an inner loop over the number of c *** Runge-Kutta substeps. i_bias is an index shift that will c *** control the one-sided differencing for MacCormack's method. R_old = 0.0 i_bias = 0 do nt=1, nt_max do nrk=1, nrk_max if( i_bias .eq. 0 ) then i_bias = 1 else i_bias = 0 end if c *** Compute the primitive variables. u_max = 0.0 do i=1, Nx+1 rho(i) = cons(i,1) u( i) = cons(i,2)/cons(i,1) p( i) = c0*( cons(i,3) - 0.5*rho(i)*u(i)**2 ) u_max = max( u_max, abs(u(i)) ) end do if( u_max .gt. 2.0 ) then print *, 'ERROR: solution blew up' goto 90 end if c *** Compute the right hand side. do i=2, Nx ie = i + i_bias iw = ie - 1 F1e = cons(ie,2) F1w = cons(iw,2) F2e = rho(ie)*u(ie)**2 + p(ie) F2w = rho(iw)*u(iw)**2 + p(iw) F3e = ( cons(ie,3) + p(ie) )*u(ie) F3w = ( cons(iw,3) + p(iw) )*u(iw) R(i,1) = -(A(ie)*F1e - A(iw)*F1w)*dx_inv R(i,2) = -(A(ie)*F2e - A(iw)*F2w)*dx_inv + p(i)*dAdx(i) R(i,3) = -(A(ie)*F3e - A(iw)*F3w)*dx_inv end do c *** Perform the solution update. do n=1, 3 do i=2, Nx cons(i,n) = cons(i,n) + dt*( gamma(nrk)*R( i,n) + & zeta( nrk)*R_old(i,n) & )*A_inv(i) end do end do c *** Apply the exit boundary conditions. i=Nx+1 c_exit = sqrt( gamma_gas*p(i)/rho(i) ) drhodx = ( rho(i) - rho(i-1) )*dx_inv dudx = ( u( i) - u( i-1) )*dx_inv dpdx = ( p( i) - p( i-1) )*dx_inv rho_exit = rho(i) - u(i)*dt*( drhodx - dpdx/c_exit**2 ) u_exit = u(i) - (u(i)+c_exit)*dt* & ( dudx + dpdx/(rho(i)*c_exit) ) - & dt*u(i)*c_exit*dAdx(i)*A_inv(i) cons(i,1) = rho_exit cons(i,2) = rho_exit*u_exit cons(i,3) = p_exit*c0_inv + 0.5*rho_exit*u_exit**2 c *** Copy the current right hand side to the R_old array do n=1, 3 do i=2, Nx R_old(i,n) = R(i,n) end do end do end do time = time + dt c *** Done with time step. Compute and write the error. if( mod(nt,n_skip) .eq. 0 .or. nt .eq. 1 ) then sum1 = 0.0 sum2 = 0.0 sum3 = 0.0 do i=1, Nx+1 sum1 = sum1 + ( rho(i) - rho_exact(i) )**2 sum2 = sum2 + ( u( i) - u_exact( i) )**2 sum3 = sum3 + ( p( i) - p_exact( i) )**2 end do error_rho = sqrt( sum1*Nx_inv ) error_u = sqrt( sum2*Nx_inv ) error_p = sqrt( sum3*Nx_inv ) write(2,20) time, error_rho, error_u, error_p 20 format(1p6e12.4) end if end do c *** Done with time stepping loop. Write the final solution. 90 continue do i=1, Nx+1 x = (i-1)*dx write(3,20) x, rho(i)-rho_exact(i), u(i)-u_exact(i), & p(i)-p_exact(i) write(4,20) x, rho(i), u(i), p(i)*p_fac write(7,20) x, p(i)*p_fac, p_exact(i)*p_fac end do stop end subroutine quasi_exact( M_inlet, p_ratio, Nx, rho, u, p, x_s ) c *** Solve the quasi-1D equations exactly. c *** Arguments: c *** M_inlet - inlet Mach number (input) c *** p_ratio - exit to inlet pressure ratio (input) c *** Nx - number of points in x (input) c *** u(Nx+1) - velocity scaled by inlet velocity (output) c *** u(Nx+1) - pressure scaled by inlet density times inlet velocity^2 c *** (output) c *** rho(Nx+1) - density scaled by inlet density (output) c *** xs - shock location scaled by duct length (output) parameter( gamma=1.4 ) real M_inlet, M_save, M, M1, M1_sq, M2_sq real rho(Nx+1), u(Nx+1), p(Nx+1) c *** Define functions. f(M) = (A/A*)^2, df(m) = df/dM f( M) = ( ( c2*( 1.0 + c1*M**2 ) )**c3 )/M**2 df(M) = 2.0/M*( ( c2*( 1.0 + c1*M**2 ) )**c4 - f(M) ) c *** Define frequently used constants. c1 = 0.5*( gamma - 1.0 ) c2 = 2.0/( gamma + 1.0 ) c3 = ( gamma + 1.0 ) / ( gamma - 1.0 ) c4 = c3 - 1.0 c5 = gamma / ( gamma - 1.0 ) c6 = 1.0 / ( gamma - 1.0 ) c7 = 1.0/(gamma*M_inlet**2) p_exit = c7*p_ratio dx = 1.0/float(Nx) c *** Enter the iteration loop for the shock position. c *** The shock position, x_s, is set to an initial guess. if( p_ratio .eq. 0.0 ) then x_s = 10000.0 else x_s = 0.7654321 end if do k=1, 10 c *** Set reference quantities to values at the inlet. x_ref = 0.0 A_ratio_sq_ref = f(M_inlet) T_factor = 1.0 + c1*M_inlet**2 p_factor = T_factor**c5 M_save = M_inlet c *** Enter the loop to march down the inlet. do i=1, Nx+1 x = (i-1)*dx c *** Compute the local area ratio. A_ratio_sq = A_ratio_sq_ref*(Area(x)/Area(x_ref))**2 c *** Use Newton iteration to determine the local Mach number. c *** As an initial guess, use the Mach number from the prior x station. M = M_save do j=1, 10 error = f(M) - A_ratio_sq M = M - error/df(M) end do c *** Compute the remainder of the flow variables. T = T_factor / ( 1.0 + c1*M**2 ) pp = p_factor / ( 1.0 + c1*M**2 )**c5 rho(i) = pp/T p(i) = c7*pp u(i) = M/M_inlet*sqrt( T ) M_save = M c *** Compute the jumps across the shock wave. c *** Note that x_s need not lie on a mesh point. if( x .le. x_s .and. x_s .lt. x+dx .and. M .gt. 1.0 ) then x = x_s A_ratio_sq = A_ratio_sq_ref*(Area(x)/Area(x_ref))**2 M1 = M do j=1, 10 error = f(M1) - A_ratio_sq M1 = M1 - error/df(M1) end do M1_sq = M1**2 M2_sq= ( 1.0 + c1*M1_sq )/( gamma*M1_sq - c1 ) p_jump = 1.0 + c2*gamma*( M1_sq - 1.0 ) rho_jump = ( M1_sq/c2 ) / ( 1.0 + c1*M1_sq ) p0_jump= ( rho_jump**gamma/p_jump )**c6 c *** Adjust the pressure and area ratio reference values c *** to account for the loss in total pressure. p_factor= p_factor*p0_jump A_ratio_sq_ref = A_ratio_sq_ref*p0_jump**2 M_save = sqrt( M2_sq ) end if end do c *** Use secant bisection to update the shock position. if( p_ratio .eq. 0.0 ) p_exit = p(Nx+1) error = p(Nx+1) - p_exit if( k .eq. 1 ) then if( abs(error) .lt. 1.e-7 ) return x_s_old = x_s x_s = x_s + 0.1 else derror = ( error - error_old )/( x_s - x_s_old ) x_s_old = x_s x_s = x_s - error/derror if( abs(error) .lt. 1.e-7 ) return end if error_old = error end do c *** If the program gets here the shock position iteration failed print *, 'WARNING shock position iteration did not converge' print *, 'final error = ', error return end function Area( x ) c *** Compute duct area. c *** Arguments: c *** x - x/L (input) c *** Area - A/At (output) Area = 1 + 4*( x - 0.5 )**2 return end function dAreadx( x ) c *** Compute derivative of duct area. c *** Arguments: c *** x - x/L (input) c *** dAreadx - dAdx/(At/L) (output) dAreadx = 8*( x - 0.5 ) return end