#!/usr/bin/env python """ Use Maccormack's technique to solve the one-dimensional subsonic-supersonic isentropic nozzle flow. Provision of nozzle profile with external file is possible. Calculate analytic solution for ideal-gas and non-ideal-gas gamma. Use conservative form of flow equations. Author: Bob Wimmer Date: January 27 & 29, 2009 """ from numpy import * import Gnuplot, Gnuplot.funcutils import time #global definitions gamma = 1.4 #ratio of specific heats, gamma=5/3 for ideal mono-atomic gas C = 0.5 #Courant number (0.5 is a good choice) nx = 30 #number of grid intervals L = 3. #Length of nozzle x = linspace(0.,L,nx+1) #Array along the nozzle dx = L/float(nx) #x-interval along nozzle A = zeros(nx+1) #Profile of the nozzle's cross section R = 8.314 #Universal gas constant in J/(mol K) T_env = 293.15 #degrees Kelvin (room temperature) t_scale = L/sqrt(gamma*R*T_env) #time scale res_acc = 1.e-5 #Accuracy in residuals to which calculations are carried out mass_flux_const = 0.59 #estimate for the final value of the mass flux used for initial guess #Provide cross section profile via external file #infilename = 'A_prof_1.dat' #infile = open(infilename,'r') #i = 0 #for line in infile: # #print line # bla = line.split(' ') # A[i] = bla[1] # i = i+1 #infile.close() #Provide cross section profile via analytic expression for i in xrange(0,nx+1,1): A[i] = 1. + 2.2*(x[i] - 1.5)**2 #print x[i], A[i] #A[i] = 1. + x[i]**1.5+2.2*(x[i] - 1.5)**2 - 0.5*(x[i]-1.5)**4 A_min = A.min() #This is the minimum cross section area, the throat of the nozzle for i in xrange(0,nx+1,1): if (A[i] == A_min): j_min = i #print "j_mi = " + str(j_min) class A_of_Mach(): """Calculate cross section as a function of Mach number, M. Calculate derivative with respect to Mach number, M. This function is defined as a class because we need to pass cross sectional area as a parameter. """ def __init__(self,A=1.): self.A = A #This is how I can change the parameter def __call__(self,M=1.): M_sq = M**2 gexp = (gamma + 1.)/(gamma - 1.) paren = 1. + (gamma - 1.)/2.*M_sq f_of_M = (paren**gexp)/M_sq - ((gamma + 1.)/2.)**gexp*self.A**2 df = -2.*paren**gexp/M_sq/M + (gamma +1.)/M*paren**(2./(gamma -1.)) #print f_of_M, df return f_of_M, df def rtsafe(func,x1,x2,xacc): """Safe version of a Newton-Raphson algorithm that combines this with bisection. Call with function func to be solved (we want to solve f(x) = 0 for x)) x1 and x2 are values that bracket the solution and must be passed by user to initialize xacc gives the accuracy that should be attempted. Don't be too greedy! """ maxit = 100 fl, df = func(x1) fh, df = func(x2) #print 333, x1, fl, x2, fh, fl*fh if (fl*fh > 0): print "error with guess_lo = ", x1, ", f= ", fl, " guess_hi = ", x2, ", f= ",fh return "This will not work as initial guesses do not bracket solution" if (fl < 0): xl = x1 xh = x2 else: xh = x1 xl = x2 rtsafe = 0.5*(xl+xh) dxold = fabs(xh-xl) dx = dxold f, df = func(rtsafe) for i in xrange(0,maxit,1): if ((((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) > 0.0) or (fabs(2.0*f) > fabs(dxold*df))): """use bisection""" dxold = dx dx = 0.5*(xh - xl) rtsafe = xl+dx if (xl == rtsafe): return rtsafe else: """use Newton-Raphson step""" dxold=dx dx=f/df temp=rtsafe rtsafe -= dx if (temp == rtsafe): return rtsafe #no appreciable change ===> convergence if (fabs(dx) < xacc): return rtsafe #convergence criterion met f, df = func(rtsafe) if (f < 0.): xl = rtsafe else: xh = rtsafe print "maximum iterations reached in rtsafe! value for x will not be good" return rtsafe #Calculate analytic solution #declare arrays M_an = zeros(nx+1) M_an_sq = zeros(nx+1) T_an = M_an.copy() v_an = M_an.copy() p_an = M_an.copy() rho_an = M_an.copy() j = 0 #determine analytic solutions for gamma = 5/3 if (gamma == 5./3.): p = [1., 12., 54., 0, 81.] #polynomial of order 4 has four (complex) solutions for i in xrange(0,nx+1,1): p[3]=(108. - 256.*(A[i]/A_min)**2) sol= roots(p) sol= real_if_close(sol)#this ensures that we dont have numerically insignificant imginary parts in the two real solutions #Make sure we use the right solution. Need the larger one after the throat #if A[i] == A_min: j = 1 M_an_sq[i] = real(sol[3-j]) M_an[i] = sqrt(M_an_sq[i]) check = 1./M_an_sq[i]*(0.75*(1.+1./3.*M_an_sq[i]))**4 #print i + 1, j, A[i], sqrt(check), M_an[i] rho_an[i] = (1.+1./3.*M_an_sq[i])**(-3./2.) T_an[i] = 1./(1.+1./3.*M_an_sq[i]) p_an[i] = (1.+1./3.*M_an_sq[i])**(-5./2.) else: AM = A_of_Mach() guess_lo = 0.01#0.01*rtsafe(AM,1.e-3,1.,1.e-4) guess_hi = 0.99#200.*guess_lo #print guess_lo, guess_hi for i in xrange(0,nx+1,1): AM.A = A[i] agl = AM(guess_lo)[0] agh = AM(guess_hi)[0] k = 0 while ((agl*agh > 0) and (k < 10)): guess_hi = 1.5*guess_hi guess_lo = 0.95*guess_lo agl = AM(guess_lo)[0] agh = AM(guess_hi)[0] k = k+1 #print 11111, guess_lo, agl, guess_hi, agh #print i, A[i], guess_lo, agl, guess_hi, agh #print i, A[i], rtsafe(AM,guess_lo,guess_hi,1.e-4) if (A[i] == A_min): M_an[i] = 1. else: M_an[i] = rtsafe(AM,guess_lo,guess_hi,1.e-4) guess_lo = M_an[i] guess_hi = 1.5*guess_lo #print M_an[i] if (i> 2): guess_hi = M_an[i] + 1.5*(M_an[i]-M_an[i-1]) agl = AM(guess_lo)[0] agh = AM(guess_hi)[0] #print guess_lo, agl, guess_hi, agh M_an_sq[i] = M_an[i]**2 p_an[i] = (1.+(gamma-1.)/2.*M_an_sq[i])**(-gamma/(gamma - 1.)) rho_an[i] = (1.+(gamma-1.)/2.*M_an_sq[i])**(-1./(gamma - 1.)) T_an[i] = (1.+(gamma-1.)/2.*M_an_sq[i])**-1. a = Gnuplot.Gnuplot() data_rho_an = Gnuplot.Data(x,rho_an,using=(1,2)) data_p_an = Gnuplot.Data(x,p_an,using=(1,2)) data_T_an = Gnuplot.Data(x,T_an,using=(1,2)) data_M_an = Gnuplot.Data(x,M_an,using=(1,2)) a.plot(data_M_an,data_rho_an,data_p_an,data_T_an) def solver(user_action=None): t0 = time.clock() #start timer to find out how much CPU time is used in solver """Define arrays etc.""" U1 = zeros(nx+1) U2 = zeros(nx+1) U3 = zeros(nx+1) F1 = zeros(nx+1) F2 = zeros(nx+1) F3 = zeros(nx+1) J2 = zeros(nx+1) du1_dt = U1.copy() du2_dt = U1.copy() du3_dt = U1.copy() U1_p = U1.copy() U2_p = U1.copy() U3_p = U1.copy() F1_p = U1.copy() F2_p = U1.copy() F3_p = U1.copy() J2_p = U1.copy() rho_p = U1.copy() v_p = U1.copy() T_p = U1.copy() du1_p_dt = U1.copy() du2_p_dt = U1.copy() du3_p_dt = U1.copy() U1_pp = U1.copy() U2_pp = U1.copy() U3_pp = U1.copy() F1_pp = U1.copy() F2_pp = U1.copy() F3_pp = U1.copy() J2_pp = U1.copy() rho_pp = U1.copy() v_pp = U1.copy() T_pp = U1.copy() #above only needed for conservative form rho = zeros(nx+1) v = rho.copy() T = v.copy() dta = v.copy() cs = v.copy() """Set up initial conditions""" for i in xrange(0,nx+1,1): if (x[i] <= 0.5): rho[i] = 1. T[i] = 1. elif (x[i] > 0.5 and x[i] < 1.5): rho[i] = 1. - 0.366*(x[i] - 0.5) T[i] = 1. - 0.167*(x[i] - 0.5) else: rho[i] = 0.634 - 0.3879*(x[i] - 1.5) T[i] = 0.833 - 0.3507*(x[i] - 1.5) v[0:nx+1] = mass_flux_const/rho[0:nx+1]/A[0:nx+1] #constant mass flow cs[0:nx+1] = sqrt(T[0:nx+1]) U1[0:nx+1] = rho[0:nx+1]*A[0:nx+1] U2[0:nx+1] = U1[0:nx+1]*v[0:nx+1] U3[0:nx+1] = U1[0:nx+1]*(T[0:nx+1]/(gamma - 1.) + gamma/2.*v[0:nx+1]**2) F1[0:nx+1] = U2[0:nx+1] F2[0:nx+1] = U2[0:nx+1]**2/U1[0:nx+1] + (gamma - 1.)/gamma*(U3[0:nx+1] - gamma/2.*U2[0:nx+1]**2/U1[0:nx+1]) F3[0:nx+1] = gamma*U2[0:nx+1]*U3[0:nx+1]/U1[0:nx+1] - gamma*(gamma-1.)/2.*U2[0:nx+1]**3/U1[0:nx+1]**2 J2[0:nx] = (gamma - 1.)/gamma* (U3[0:nx]- gamma/2.*U2[0:nx]**2/U1[0:nx])* (log(A[1:nx+1])- log(A[0:nx]))/dx J2[nx] = (gamma - 1.)/gamma *(U3[nx] - gamma/2.*U2[nx]**2/U1[nx])*(log(A[nx]) - log(A[nx-1]))/dx # print '-------------------------------------------------------------------' # print 'initial values' # print '-------------------------------------------------------------------' # for i in xrange(0,nx+1,1): # print i, x[i], rho[i], v[i], T[i], U1[i], U2[i], U3[i], v[i] + sqrt(T[i]) # raw_input('Please press the return key to continue...\n') tmp = Gnuplot.Gnuplot() t = 0. dres = 1.e6 oldres = 1. # while dres > 5e5: # dres = 5.e-6 while dres > res_acc: """compute time step""" for i in xrange(0,nx+1,1): dta[i] = fabs(C*dx/(cs[i]+v[i])) # print 'time step calculation: ', i, t, dta[i], dx, cs[i], v[i], dx/(cs[i]+v[i]) dt = dta.min() #dt = 0.0267 #print 'final dt is: ', dta, '\n******', dt #print 'final dt is: ', dt t_old = t; t += dt * t_scale #raw_input('Please press the return key to continue...\n') """predictor step""" #prepare derivatives using forward differences for i in xrange(0,nx,1): du1_dt[i] = - (F1[i+1] - F1[i])/dx du2_dt[i] = - (F2[i+1] - F2[i])/dx + J2[i] du3_dt[i] = - (F3[i+1] - F3[i])/dx #Now calculate predictor values for i in xrange(0,nx,1): U1_p[i] = U1[i] + du1_dt[i]*dt U2_p[i] = U2[i] + du2_dt[i]*dt U3_p[i] = U3[i] + du3_dt[i]*dt F1_p[i] = U2_p[i] v_p[i] = U2_p[i]/U1_p[i] aux1 = U2_p[i]*v_p[i] F2_p[i] = aux1 + (gamma - 1.)/gamma *(U3_p[i] - gamma/2.*aux1) F3_p[i] = gamma*v_p[i]*U3_p[i] - gamma*(gamma - 1.)/2. * v_p[i]**2*U2_p[i] J2_p[i] = (gamma - 1.)/gamma *(U3_p[i] - gamma/2.*aux1)*(log(A[i+1]) - log(A[i]))/dx rho_p[i] = U1_p[i]/A[i] T_p[i] = (gamma - 1.)*(U3_p[i]/U1_p[i] - gamma/2.*v_p[i]**2) """Update boundary conditions""" rho_p[0] = 1. U1_p[0] = rho_p[0]*A[0] U2_p[0] = 2.*U2_p[1] - U2_p[2] v_p[0] = U2_p[0]/U1_p[0] T_p[0] = 1. U3_p[0] = U1_p[0]*(T[0]/(gamma -1.)+gamma/2*v_p[0]**2) F1_p[0] = U2_p[0] aux1 = U2_p[0]**2/U1_p[0] F2_p[0] = aux1 + (gamma - 1.)/gamma *(U3_p[0] - gamma/2.*aux1) F3_p[0] = gamma*v_p[0]*U3_p[0] - gamma*(gamma - 1.)/2. * v_p[0]**2*U2_p[0] J2_p[0] = (gamma - 1.)/gamma *(U3_p[0] - gamma/2.*aux1)*(log(A[1]) - log(A[0]))/dx U1_p[nx] = 2.*U1_p[nx-1] - U1_p[nx-2] U2_p[nx] = 2.*U2_p[nx-1] - U2_p[nx-2] U3_p[nx] = 2.*U3_p[nx-1] - U3_p[nx-2] rho_p[nx]= U1_p[nx]/A[nx] v_p[nx] = U2_p[nx]/U1_p[nx] T_p[nx] = (gamma - 1.)*(U3_p[nx]/U1_p[nx] - gamma/2.*v_p[nx]**2) F1_p[nx] = U2_p[nx] aux1 = U2_p[nx]**2/U1_p[nx] F2_p[nx] = aux1 + (gamma - 1.)/gamma *(U3_p[nx] - gamma/2.*aux1) F3_p[nx] = gamma*v_p[nx]*U3_p[nx] - gamma*(gamma - 1.)/2. * v_p[nx]**2*U2_p[nx] J2_p[nx] = (gamma - 1.)/gamma *(U3_p[nx] - gamma/2.*aux1)*(log(A[nx]) - log(A[nx-1]))/dx #F1_p[nx] = 2.*F1_p[nx-1] - F1_p[nx-2] #F2_p[nx] = 2.*F2_p[nx-1] - F2_p[nx-2] #F3_p[nx] = 2.*F3_p[nx-1] - F3_p[nx-2] #J2_p[nx] = 2.*J2_p[nx-1] - J2_p[nx-2] """corrector step""" #prepare derivatives using rearward differences for i in xrange(1,nx+1,1): du1_p_dt[i] = - (F1_p[i] - F1_p[i-1])/dx du2_p_dt[i] = - (F2_p[i] - F2_p[i-1])/dx + J2_p[i] du3_p_dt[i] = - (F3_p[i] - F3_p[i-1])/dx #print 'dx_p_dt values at i = 16: ', drho_p_dt[15], dv_p_dt[15], dT_p_dt[15] """time propagator""" for i in xrange(1,nx+1,1): U1_pp[i] = U1[i] + 0.5*(du1_dt[i] + du1_p_dt[i])*dt U2_pp[i] = U2[i] + 0.5*(du2_dt[i] + du2_p_dt[i])*dt U3_pp[i] = U3[i] + 0.5*(du3_dt[i] + du3_p_dt[i])*dt F1_pp[i] = U2_pp[i] v_pp[i] = U2_pp[i]/U1_pp[i] aux1 = U2_pp[i]**2*v_pp[i] F2_pp[i] = aux1 + (gamma - 1.)/gamma *(U3_pp[i] - gamma/2.*aux1) F3_pp[i] = gamma*v_pp[i]*U3_pp[i] - gamma*(gamma - 1.)/2. * v_pp[i]**2*U2_pp[i] rho_pp[i] = U1_pp[i]/A[i] T_pp[i] = (gamma - 1.)*(U3_pp[i]/U1_pp[i] - gamma/2.*v_pp[i]**2) for i in xrange(1,nx,1): J2_pp[i] = (gamma - 1.)/gamma *(U3_pp[i] - gamma/2.*U2_pp[i]**2/U1_pp[i])*(log(A[i+1]) - log(A[i]))/dx """insert boundary conditions""" rho_pp[0] = 1. U1_pp[0] = rho_pp[0]*A[0] U2_pp[0] = 2.*U2_pp[1] - U2_pp[2] v_pp[0] = U2_pp[0]/U1_pp[0] T_pp[0] = 1. U3_pp[0] = U1_pp[0]*(T[0]/(gamma -1.)+gamma/2*v_pp[0]**2) F1_pp[0] = U2_pp[0] aux1 = U2_pp[0]**2/U1_pp[0] F2_pp[0] = aux1 + (gamma - 1.)/gamma *(U3_pp[0] - gamma/2.*aux1) F3_pp[0] = gamma*v_pp[0]*U3_pp[0] - gamma*(gamma - 1.)/2. * v_pp[0]**2*U2_pp[0] J2_pp[0] = (gamma - 1.)/gamma *(U3_pp[0] - gamma/2.*aux1)*(log(A[1]) - log(A[0]))/dx U1_pp[nx] = 2.*U1_pp[nx-1] - U1_pp[nx-2] U2_pp[nx] = 2.*U2_pp[nx-1] - U2_pp[nx-2] U3_pp[nx] = 2.*U3_pp[nx-1] - U3_pp[nx-2] rho_pp[nx]= U1_pp[nx]/A[nx] v_pp[nx] = U2_pp[nx]/U1_pp[nx] T_pp[nx] = (gamma - 1.)*(U3_pp[nx]/U1_pp[nx] - gamma/2.*v_pp[nx]**2) F1_pp[nx] = U2_pp[nx] aux1 = U2_pp[nx]**2/U1_pp[nx] F2_pp[nx] = aux1 + (gamma - 1.)/gamma *(U3_pp[nx] - gamma/2.*aux1) F3_pp[nx] = gamma*v_pp[nx]*U3_pp[nx] - gamma*(gamma - 1.)/2. * v_pp[nx]**2*U2_pp[nx] J2_pp[nx] = (gamma - 1.)/gamma *(U3_pp[nx] - gamma/2.*aux1)*(log(A[nx]) - log(A[nx-1]))/dx #F1_pp[nx] = 2.*F1_pp[nx-1] - F1_pp[nx-2] #F2_pp[nx] = 2.*F2_pp[nx-1] - F2_pp[nx-2] #F3_pp[nx] = 2.*F3_pp[nx-1] - F3_pp[nx-2] #J2_pp[nx] = 2.*J2_pp[nx-1] - J2_pp[nx-2] print 'pressure at end of nozzle is: ', rho_pp[nx], T_pp[nx], rho_p[nx]*T_pp[nx] #raw_input('Please press the return key to continue...\n') """The above values are the new ones (at the next time step)""" rho, v, T = rho_pp, v_pp, T_pp U1, U2, U3 = U1_pp, U2_pp, U3_pp F1, F2, F3 = F1_pp, F2_pp, F3_pp J2 = J2_pp cs[0:nx+1] = sqrt(T[0:nx+1]) tmp.reset() tmp('set ylabel "y-axis [arb. units]"') #this is how you access gnupot commands tmp('set xlabel "x-axis [arb. units]"') tmp('set yrange [-4:6]') string = 'set title "time: ' + str(t) +'"' tmp(string) data_u1 = Gnuplot.Data(x,U1,using=(1,2), title = 'U1') data_u2 = Gnuplot.Data(x,U2,using=(1,2), title = 'U2') data_u3 = Gnuplot.Data(x,U3,using=(1,2), title = 'U3') data_f1 = Gnuplot.Data(x,F1,using=(1,2), title = 'F1') data_f2 = Gnuplot.Data(x,F2,using=(1,2), title = 'F2') data_f3 = Gnuplot.Data(x,F3,using=(1,2), title = 'F3') data_j2 = Gnuplot.Data(x,J2,using=(1,2), title = 'J2') tmp.plot(data_u1, data_u2, data_u3, data_f1, data_f2, data_f3, data_j2) # print '-------------------------------------------------------------------' # print 'end of time step values' # print '-------------------------------------------------------------------' # for i in xrange(0,nx+1,1): # print i, rho[i], v[i], T[i], U1[i], U2[i], U3[i], F1[i], F2[i], F3[i], J2[i] # raw_input('Please press the return key to continue...\n') """Determine residuals""" U1_res = 0. U2_res = 0. U3_res = 0. for i in xrange(0,nx+1,1): U1_res += (0.5*(du1_dt[i] + du1_p_dt[i]))**2 U2_res += (0.5*(du2_dt[i] + du2_p_dt[i]))**2 U3_res += (0.5*(du3_dt[i] + du3_p_dt[i]))**2 #print 'residuals: ', rho_res, v_res, T_res dres = fabs(oldres - U1_res - U2_res - U3_res) oldres = U1_res + U2_res + U3_res if user_action is not None: user_action(rho, v, T, U1_res, U2_res, U3_res, dres, x, t) #I can do whatever I'd like here t1 = time.clock() return t1-t0 , t #solver does NOT return solution u! #But it returns how long it worked (t1-t0) and the current model time # ---------------- end of solver ------------------------------------------- def test_solver(N): """ Very simple test case. Store the solution at every N time levels. Measure how long the solver actually works """ g = Gnuplot.Gnuplot() rho_solutions = [] v_solutions = [] T_solutions = [] timesteps = [] rho_residuals = [] v_residuals = [] T_residuals = [] # Need time_level_counter as global variable since # it is assigned in the action function (that makes # a variable local to that block otherwise). global time_level_counter time_level_counter = 0 def action(rho, v, T, rho_res, v_res, T_res, dres, x, t): #this is where I can access solution u global time_level_counter if time_level_counter % N == 0: #only do this every N times timesteps.append(t) rho_residuals.append(rho_res) v_residuals.append(v_res) T_residuals.append(T_res) rho_solutions.append(rho.copy()) v_solutions.append(v.copy()) T_solutions.append(T.copy()) p = zeros(nx+1) p[0:nx+1] = rho[0:nx+1]*T[0:nx+1] M = zeros(nx+1) M[0:nx+1] = v[0:nx+1]/sqrt(T[0:nx+1]) g.reset() string = 'set title "time: ' + str(t) +' s"' g(string) g('set ylabel "y-axis [arb. units]"') #this is how you access gnupot commands g('set xlabel "x-axis [arb. units]"') data_rho = Gnuplot.Data(x,rho,using=(1,2), title = 'rho') data_v = Gnuplot.Data(x,v,using=(1,2), title = 'v') data_T = Gnuplot.Data(x,T,using=(1,2), title = 'T') data_p = Gnuplot.Data(x,p,using=(1,2), title = 'p') data_M = Gnuplot.Data(x,M,using=(1,2), title = 'M') data_rho_an = Gnuplot.Data(x,rho_an,using=(1,2), with_ = 'line') data_p_an = Gnuplot.Data(x,p_an,using=(1,2), with_ = 'line') data_T_an = Gnuplot.Data(x,T_an,using=(1,2), with_ = 'line') data_M_an = Gnuplot.Data(x,M_an,using=(1,2), with_ = 'line') g.plot(data_rho,data_v,data_p,data_T,data_M,data_rho_an,data_p_an,data_T_an, data_M_an) #g.plot(data_rho,data_v,data_T) #z = time.sleep(0.1) #wait a little so you can look at it #print rho_res, dres time_level_counter += 1 #call the solver here, note that the solution, u, is not returned #but cpu time is cpu, t = solver(user_action=action) print 'CPU time:', cpu, ' sec' print 'time: ', t, 'seconds physical time' return timesteps, rho_solutions, v_solutions, T_solutions, rho_residuals, v_residuals, T_residuals #test_solver has access to solutions through 'action' #and returns solution array to main # --------------------------------- end test_solver ----------------------------- def main(): # open output file #output = open("1dwave-sol.dat","w") timesteps, rho_solutions, v_solutions, T_solutions, rho_residuals, v_residuals, T_residuals = test_solver(1) #this is where I call the solver #print solutions #output.write(str(solutions)) #note: the output is an array #output.close() raw_input('Please press the return key to continue...\n') main()