<jupyter_start><jupyter_text>Analysis of Quasi-One-Dimensional Nozzle Flow with Analytical and Numerical MethodsThis notebook covers an in-depth analysis of an isotropic flow through a converging diverging nozzle for both the subsonic and the supersonic regime based on the problem addressed in chapter 7 of <NAME>'s book *Computational Fluid Dynamics: The Basics with Applications*. The inlet of the nozzle is fed from a reservoir and will have a constant pressure and temperature $(p_0, T_0)$. The flow will be subsonic in the converging portion of the nozzle, sonic at the throat (M=1) and supersonic at the exit of the nozzle.  Analytical SolutionThis problem is solved analytically which gives us a great reference point to try out numerical methods. Let's first look at the analytical solution and then move on to numerical methods to solve the problem using CFD. The analytical solutions for continuity, energy and momentum are as follows. Continuity: $\rho_1 V_1 A_1 = \rho_2 V_2 A_2 $Momentum: $p_1 A_1 + \rho_1 V_1^2 A_1 + \int_{A_1}^{A_2}pdA = p_2 A_2 + \rho_2 V_2^2 A_2$ Energy: $h_1 + \frac{v_1^2}{2} = h_2 + \frac{V_2^2}{2}$The different subscripts denote the location of the flow inside of the nozzle. In addition to the flow equations, we also have the ideal gas equation:$$p = \rho R T$$This flow consists of a calorically perfect gas:$$h = c_p T$$The Mach number in the nozzle is determined by the area ratio $A/A^*$ where $A^*$ is the area of the nozzle at the choke point (M=1).$$\Big (\frac{A}{A^*} \Big)^2 = \frac{1}{M^2}\Big[\frac{2}{\gamma +1}\Big(1+\frac{\gamma-1}{2}M^2\Big)\Big]^{(\gamma+1)/(\gamma-1)}$$Gamma ($\gamma$) is defined as the ratio of specific heats and for air, the standard value is 1.4. The rest of the parameters can be solved in terms of the mach number:$$\frac{p}{p_0}  = \Big(1+\frac{\gamma-1}{2}M^2\Big)^{(-\gamma)/(\gamma-1)}$$$$\frac{\rho}{\rho_0}  = \Big(1+\frac{\gamma-1}{2}M^2\Big)^{-1/(\gamma-1)}$$$$\frac{T}{T_0}  = \Big(1+\frac{\gamma-1}{2}M^2\Big)^{-1}$$To view these relationships, we can plot the analytical solution.<jupyter_code>import numpy as np                 #loading in numpy
from matplotlib import pyplot      #load matplotlib
import math                        # load math in 
#this makes matplotlib plots appear in the notebook (instead of a separate window)
%matplotlib inline 
# Mach number vs Area Ratio
M = np.linspace(0.4,2,100)
aRatio = (1/M**2)* ((2/(2.4))*(1+(.4/2)*M**2))**((2.4)/(0.4))
pyplot.plot(aRatio,M)
pyplot.xlabel('Area Ratio')
pyplot.ylabel('Mach Number')
M = np.linspace(0,3,100)
# Pressure Ratio
pRatio = (1+ (0.4/2)*M**2)**(-1.4/0.4)
pyplot.plot(M,pRatio)
pyplot.xlabel('Mach Number')
pyplot.ylabel('Pressure Ratio p/p0')
# Density Ratio
rhoRatio = (1+ (0.4/2)*M**2)**(-1/0.4)
pyplot.plot(M,rhoRatio)
pyplot.xlabel('Mach Number')
pyplot.ylabel('Density Ratio rho/rho0')
# Temperature Ratio
TRatio = (1+ (0.4/2)*M**2)**(-1)
pyplot.plot(M,TRatio)
pyplot.xlabel('Mach Number')
pyplot.ylabel('Temperature Ratio T/T0')<jupyter_output><empty_output><jupyter_text>Numerical Solution The first step of any numerical solution is to first write down the governing equations for the problem. In this case, we will be looking for a continuity, momentum and energy equation that will describe the flow. Continuity EquationFor the continuity equation, we will use the integral finite volume method.$$\frac{\partial}{\partial t} \iiint \limits_V \rho dV + \frac{\partial}{\partial t} \iint \limits_S \rho \textbf{V} \cdot \textbf{dS} = 0$$This equation means that the total rate of change of mass inside of the volume is equal to the amount of stuff crossing the surface of the volume. To relate the volume integral to out physical problem, we can write the following equation:$$\frac{\partial}{\partial t} \iiint \limits_V \rho dV = \frac{\partial}{\partial t} (\rho A dx)$$In this equation the value of dx*A represents an elemental volume inside of our nozzle. It is like a cross sectional area multiplied by some small length of dx. Now for the surface integral, an equation can be written as follows. $$\frac{\partial}{\partial t} \iint \limits_S \rho \textbf{V} \cdot \textbf{dS} = -\rho V A + (\rho + d \rho)(V+ dV)(A+dA)$$Expanding this equation and combining it with the previous one, we have a form of the contiuity equation in terms of differential elements relating to our system.$$\frac{\partial (\rho A)}{\partial t} + \frac{\partial (\rho A V)}{\partial x} = 0$$ Momentum EquationThe next step in describing our system is to utilize the momentum equation. Neglecting viscous  stress terms and body forces we can represent the momentum equation as shown below.$$\frac{\partial}{\partial t} \iiint \limits_V (\rho u) dV + \frac{\partial}{\partial t} \iint \limits_S (\rho u \textbf{V}) \cdot \textbf{dS} = \iint \limits_S (p dS)_x$$The first two terms can be transformed in a similar way as shown in the previous section but with an additional velocity component. $$\frac{\partial}{\partial t} \iiint \limits_V (\rho u) dV = \frac{\partial}{\partial t} (\rho V A dx)$$$$\frac{\partial}{\partial t} \iint \limits_S (\rho u \textbf{V}) \cdot \textbf{dS} = -\rho V^2 A + (\rho + d \rho)(V+ dV)^2(A+dA)$$The term on the right side of the equal sign represents the pressure around our element of fluid. This term can be represented with the following expression.$$\iint \limits_S (p dS)_x = -pA + (p+dp)(A+dA) - 2p\Big(\frac{dA}{2}\Big)$$Combining all of these terms we get the following equation.$$\frac{\partial (\rho V A)}{\partial t} + \frac{\partial (\rho V^2 A)}{\partial x} = -A \frac{\partial p}{\partial x}$$This is the conservation form of this equation, but we want a non-conservation form for our numerical analysis. WE can do this by multiplying the continuity equation by V resulting in the following equation. $$V\frac{\partial (\rho A)}{\partial t} + V\frac{\partial (\rho A V)}{\partial x} = 0$$Now we can subtract this equation from the momentum equation we derived just before. With some manipulation, we are left with the final momentum equation. $$\rho \frac{\partial V}{\partial t} + \rho u \frac{\partial u}{\partial x} = -\frac{\partial p}{\partial x}$$ Energy EquationThe integral form of the energy equation with no energy source terms is as follows.$$\frac{\partial}{\partial t} \iiint \limits_V \rho \Big(e+\frac{V^2}{2}\Big) dV + \iint \limits_S \rho\Big(e+\frac{V^2}{2}\Big)\textbf{V} \cdot \textbf{dS} = \iint \limits_S (p \textbf{V}) \cdot \textbf{dS}$$With some manipulation, we are able to transform the equation into an equation that represents our system and is in non-conservation form.  $$\rho \frac{\partial e}{\partial t} + \rho V \frac{\partial e}{\partial x} = -p\frac{\partial V}{\partial x} - pV\frac{\partial(ln A)}{\partial x}$$There is one thing left to consider, and that is the heat coefficient for the fluid. In this case, we are dealing with a calorically perfect gas. $$e = c_v T$$Finally we have an energy equation that we can use for our numerical methods.$$\rho  c_v\frac{\partial T}{\partial t} + \rho V c_v \frac{\partial T}{\partial x} = -p\frac{\partial V}{\partial x} - pV\frac{\partial(ln A)}{\partial x}$$ Summary of EquationsAs we take a moment to examine our equations, we can see that the variables are $\rho$, V, p and T but with the use of the ideal gas equation, we can reduce the number of unknowns by one. $$p = \rho R T$$$$\frac{\partial p}{\partial x} = R \Big( \rho \frac{\partial T}{\partial x} + T \frac{\partial \rho}{\partial x}\Big )$$Now we can implement this equation and reduce the continuity, momentum, and energy equations to only have three unknown variables. The continuity equation is simply expanded and the other two are rearranged with the new definition. Continuity: $\frac{\partial (\rho A)}{\partial t} + \rho A\frac{\partial V}{\partial x}+ \rho V \frac{\partial A}{\partial x}+ VA\frac{\partial \rho}{\partial x} = 0$Momentum: $\rho \frac{\partial V}{\partial t} + \rho V \frac{\partial V}{\partial x} = -R\Big(\rho \frac{\partial T}{\partial x}+ T \frac{\partial \rho}{\partial x}\Big)$Energy: $\rho  c_v\frac{\partial T}{\partial t} + \rho V c_v \frac{\partial T}{\partial x} = -\rho R T\Big[\frac{\partial V}{\partial x} + V\frac{\partial(ln A)}{\partial x}\Big]$At this point, we have our flow equations and we only have a few more things to take care of before we can start writing some programs to solve the problem.  Nondimensional NumbersIn order to simplify our problem, we will use non-dimensional number for our calculations. We can write our non-dimensional temperature and density as follows.$$T' = \frac{T}{T_0}$$  and $$\rho' = \frac{\rho}{\rho_0}$$Additionally, we can define out non-dimensional length with the following equation.$$x' = \frac{x}{L}$$Where L is the length of the nozzle. The speed of sound is defined as $$a_0 = \sqrt{\gamma R T_0}$$and a nondimensional velocity of $$V' = \frac{V}{a_0}$$Time can also be put in a nondimensional form as shown in the following equation.$$t' = \frac{t}{L/a_0}$$Lastly, the area can be nondimensionalized by dividing by the throat area where the flow is at $M=1$. $$A' = \frac{A}{A^*}$$With all of the relationships we need, we can nondimensionalize our continuity, momentum and state equations.Continuity: $\frac{\partial \rho'}{\partial t'} = -\rho ' \frac{\partial V'}{\partial x'} -\rho 'V'\frac{\partial (ln A')}{\partial x'} - V'\frac{\partial \rho'}{\partial x'}$Momentum: $\frac{\partial V'}{\partial t'} = -V' \frac{\partial V'}{\partial x'} - \frac{1}{\gamma}\Big(\frac{\partial T'}{\partial x'} + \frac{T'}{\rho'}\frac{\partial \rho '} {\partial x'}\Big)$Energy: $\frac{\partial T'}{\partial t'} =-V' \frac{\partial T'}{\partial x'} - (\gamma -1)T'\Big[\frac{\partial V'}{\partial x'} +V'\frac{\partial (ln A')}{\partial x'}\Big]$Yayy! We finally made it to the end of the math part (well kinda...). The next step is to take the equations that describe our flow and put it in a form that the computer can understand. In this case, we will be using a finite difference method with MacCormack's explicit technique which uses a predictor-corrector method to help us get second order accuracy without doing an insane amount of algebra and book keeping. Because this is a quasi-one-dimensional problem, we only need to make a grid that is one dimensional.<jupyter_code># Generate our grid of points 
numPts = 14 # number of grid points
length = 3  # nondimensional length of out nozzle
grid = np.zeros(numPts)
pyplot.scatter(np.linspace(0,length,numPts),grid) #plotting the nodes 
pyplot.xlabel('Nondimensional Length (L\')')<jupyter_output><empty_output><jupyter_text>The plot above shows the uniform distribution of our points. Now we are ready to start developing the equations that we will use for the finite difference method. The first step is to execute the predictor step. To do this, we calculate the first order derivative at our point of interest at time $t$.$$\Big(\frac{\partial \rho}{\partial t}\Big)_i^t= -\rho_i^t \frac{V^t_{i+1} - V^t_i}{\Delta x} -\rho_i^t V^t_i\frac{ln A^t_{i+1} - ln A^t_i}{\Delta x} -V_i^t \frac{\rho^t_{i+1} - \rho^t_i}{\Delta x}$$The same is done for the velocity and temperature. $$\Big(\frac{\partial V}{\partial t}\Big)_i^t= -V_i^t \frac{V^t_{i+1} - V^t_i}{\Delta x} -\frac{1}{\gamma}\Big(\frac{T^t_{i+1} - T^t_i}{\Delta x} +\frac{T^t_i}{\rho^t_i} \frac{\rho^t_{i+1} - \rho^t_i}{\Delta x}\Big)$$$$\Big(\frac{\partial T}{\partial t}\Big)_i^t= -V_i^t\frac{T^t_{i+1} - T^t_i}{\Delta x} - (\gamma -1)T^t_i\Big(\frac{V^t_{i+1} - V^t_i}{\Delta x} + V^t_i \frac{ln A^t_{i+1} - ln A^t_i}{\Delta x}\Big)$$These values are then used to calculate the predicted step with the following equations$$\bar{\rho}_i^{t+\Delta t} =\rho_i^t +\Big(\frac{\partial \rho}{\partial t}\Big)_i^t \Delta t$$$$\bar{V}_i^{t+\Delta t} =V_i^t +\Big(\frac{\partial V}{\partial t}\Big)_i^t \Delta t$$$$\bar{T}_i^{t+\Delta t} =T_i^t +\Big(\frac{\partial T}{\partial t}\Big)_i^t \Delta t$$The second part of this technique is the corrector step. The corrector step gets the slope of the predicted value which will be used in the final step to determine $(\partial /\partial t)_{av}$$$\Big(\frac{\bar{\partial \rho}}{\partial t}\Big)_i^{t+\Delta t t}= -\bar{\rho}_i^{t+\Delta t} \frac{\bar{V}^{t+\Delta t}_{i} - \bar{V}^{t+\Delta t}_{i-1}}{\Delta x} -\bar{\rho}_i^{t+\Delta t} \bar{V}^{t+\Delta t}_i\frac{ln \bar{A}^{t+\Delta t}_{i} - ln \bar{A}^{t+\Delta t}_{i-1}}{\Delta x} -\bar{V}_i^{t+\Delta t} \frac{\bar{\rho}^{t+\Delta t}_{i} - \bar{\rho}^{t+\Delta t}_{i-1}}{\Delta x}$$The same is done for the velocity and temperature. $$\Big(\frac{\bar{\partial V}}{\partial t}\Big)_i^{t+\Delta t}= -\bar{V}_i^{t+\Delta t} \frac{\bar{V}^{t+\Delta t}_{i} - \bar{V}^{t+\Delta t}_{i-1}}{\Delta x} -\frac{1}{\gamma}\Big(\frac{\bar{T}^{t+\Delta t}_{i} - \bar{T}^{t+\Delta t}_{i-1}}{\Delta x} +\frac{\bar{T}^{t+\Delta t}_i}{\bar{\rho}^{t+\Delta t}_i} \frac{\bar{\rho}^{t+\Delta t}_{i} - \bar{\rho}^{t+\Delta t}_{i-1}}{\Delta x}\Big)$$$$\Big(\frac{\bar{\partial T}}{\partial t}\Big)_i^{t+\Delta t}= -\bar{V}_i^{t+\Delta t}\frac{\bar{T}^{t+\Delta t}_{i} - \bar{T}^{t+\Delta t}_{i-1}}{\Delta x} - (\gamma -1)\bar{T}^{t+\Delta t}_i\Big(\frac{\bar{V}^{t+\Delta t}_{i} - \bar{V}^{t+\Delta t}_{i-1}}{\Delta x} + \bar{V}^{t+\Delta t}_i \frac{ln \bar{A}^{t+\Delta t}_{i} - ln \bar{A}^{t+\Delta t}_{i-1}}{\Delta x}\Big)$$The partial derivative average can be calculated using the following equation$$\Big(\frac{\partial}{\partial t}\Big)_{av} = \frac{1}{2}\Big[\Big(\frac{\partial}{\partial t}\Big)_i^t+\Big(\frac{\bar{\partial}}{\partial t}\Big)_i^{t+\Delta t}\Big]$$Finally, the value we want can be calculated (with second order accuracy) using this equation$$\rho_i^{t+\Delta t} = \rho_i^t + \Big(\frac{\partial \rho}{\partial t}\Big)_{av}\Delta t$$$$V_i^{t+\Delta t} = V_i^t + \Big(\frac{\partial V}{\partial t}\Big)_{av}\Delta t$$$$T_i^{t+\Delta t} = T_i^t + \Big(\frac{\partial T}{\partial t}\Big)_{av}\Delta t$$ Now that we know how to solve the equations of our flow, let's write a function that can take information from some given nodes and return the values after one time iteration. Doing this now will save us time and space when we write the main code to solve our problem.<jupyter_code>def calcStepSuper(rho,vel,temp,area,gamma,deltat,deltax):
    """
    Import three numpy arrays that represent the 
    density velocity and temperature of each element
    with a supersonic exit flow BC
    rho : density array 
    vel : velocity array
    temp: temperature array
    deltat : time step for the elements
    deltax : space between each element
    """ 
    #loop through each element starting at index 1 and ending at index numPts-1 (N-1)
    numElm = len(rho)
    
    #Density
    fwddRhodt  = np.zeros(numElm)
    bkwddRhodt  = np.zeros(numElm)
    rhoBar = np.zeros(numElm)
    dRhoAv = np.zeros(numElm)
    #Velocity
    
    fwddvdx  = np.zeros(numElm)
    bkwddvdx  = np.zeros(numElm)
    VBar = np.zeros(numElm)
    dVAv = np.zeros(numElm)
    #Temperature 
    
    fwddTdx = np.zeros(numElm)
    bkwddTdx = np.zeros(numElm)
    TBar = np.zeros(numElm)
    dTAv = np.zeros(numElm)
    
    ##################
    ####Predictor Step
    ##################
    for i in range(1,numElm-1):
            #Predictor Step
            # Calculate first derivative of density
            fwddRhodt[i] = -rho[i] * ( (vel[i+1]-vel[i])/deltax) -(rho[i]*vel[i]*(math.log(area[i+1])-math.log(area[i]))/deltax) - (vel[i]*(rho[i+1]-rho[i])/deltax)
            rhoBar[i]  = rho[i]+  fwddRhodt[i] * deltat
            
            # Calculate first derivative of velocity
            fwddvdx[i] = -vel[i]*((vel[i+1]-vel[i])/deltax)-(1/gamma)*(((temp[i+1]-temp[i])/deltax)+(temp[i]/rho[i])*(rho[i+1]-rho[i])/deltax)
            VBar[i]  = vel[i] + fwddvdx[i] * deltat
                 
            # Calculate first derivative of temperatures
            a = -vel[i]*(temp[i+1]-temp[i])/deltax
            b = (gamma-1)*temp[i]
            c = (vel[i+1]-vel[i])/deltax
            d = vel[i]*((math.log(area[i+1])-math.log(area[i]))/deltax)
            fwddTdx[i] = a-b*(c+d)
            TBar[i] = temp[i] + (fwddTdx[i]) * deltat
            
    # Calculate BC's for the preictor step                                
    rhoBar[0] = 1
    VBar[0] = 2*VBar[1]-VBar[2]
    TBar[0] = 1
                                    
    
    # We have the values for the predicted values and the slopes
    ###################
    ####Corrector Step
    ###################
    for i in range(1,numElm-1):
        bkwddRhodt[i] = -rhoBar[i] * ( (VBar[i]-VBar[i-1])/deltax) - (rhoBar[i]*VBar[i]*(math.log(area[i])-math.log(area[i-1]))/deltax) -(VBar[i]*(rhoBar[i]-rhoBar[i-1])/deltax)
        
        bkwddvdx[i] = -VBar[i]*((VBar[i]-VBar[i-1])/deltax)-(1/gamma)*(((TBar[i]-TBar[i-1])/deltax)+(TBar[i]/rhoBar[i])*(rhoBar[i]-rhoBar[i-1])/deltax)
            
        a = -VBar[i]*(TBar[i]-TBar[i-1])/deltax
        b = (gamma-1)*TBar[i]
        c = (VBar[i]-VBar[i-1])/deltax
        d = VBar[i]*((math.log(area[i])-math.log(area[i-1]))/deltax)
        bkwddTdx[i] = a-b*(c+d)
        
        dRhoAv[i] = 0.500*(fwddRhodt[i]+bkwddRhodt[i])
        dVAv[i]   = 0.500*(fwddvdx[i]+bkwddvdx[i])
        dTAv[i]   = 0.500*(fwddTdx[i]+bkwddTdx[i])
        
        rho[i] = rho[i]   + dRhoAv[i] * deltat
        vel[i] = vel[i]   + dVAv[i]   * deltat
        temp[i] = temp[i] + dTAv[i]   * deltat
        
    #Set BC's for the final values 
    #Far left BC
    rho[0] = 1
    vel[0] = 2*vel[1]-vel[2]
    temp[0] = 1
    #Far right BC                                
    rho[numElm-1]= (2*rho[numElm-2])-rho[numElm-3]
    vel[numElm-1] = (2*vel[numElm-2])-vel[numElm-3]
    temp[numElm-1] = (2*temp[numElm-2])-temp[numElm-3]
    return rho,vel,temp,dRhoAv,dVAv #output new values and d/dt for residuals<jupyter_output><empty_output><jupyter_text>The code above implements the finite difference method we developed earlier. There were several lines of code in that function dealing with some material that we have yet to cover, boundary conditions. We need to establish the boundary conditions to anchor our mathematics to our physical problem. The left side of our grid represents the subsonic inflow of from a reservoir this is our very first node. From there we have our regular nodes governed by our flow equations and at the other end we have a supersonic our flow (we will address a subsonic exit flow later). The Nth node will be our outflow condition. The first step in setting our boundary conditions is to use the method of characteristics for an unsteady, one-dimensional flow. This type of flow is characterized by a hyperbolic PDE which means that the solution has two real characteristic lines. For our subsonic inlet, the left-running line propagates behind the inlet and the right-running characteristic line runs downstream toward the exit of the nozzle. One way to think about these lines is to consider how information travels in a flow. For the subsonic flow, information can travel upstream, because the velocity of the flow is less than the speed of sounds (how fast waves propagate in a flow). The super sonic outlet is moving at speeds greater than the speed of sound, so information can not propagate upstream. Because of this, we must allow one of the variables at the first node to float while the other two flow variables are set. As for the outflow boundary condition, all values are allowed to float due to the fact that the left and right-running characteristic lines are carried downstream.For inlet node: $$V_1 = 2V_2 - V_3$$$$\rho_1 =1$$ (Constant)$$T_1 = 1$$ (Constant)For outlet node:$$V_N = 2V_{N-1}-V_{N-2}$$$$\rho_N = 2\rho_{N-1}-\rho_{N-2}$$$$T_N = 2T_{N-1}-T_{N-2}$$  Time StepFor our time step, we must consider the Courant number$\Delta t = C\frac{\Delta x}{a+V}$Where a is the speed of sound locally and V is the local velocity. For every element, there will be a different value for a given Courant number and although we can use a local timestep, we will instead just use the minimum of the whole system at a given time step.  Nozzle ShapeFor our nozzle shpae, we will use the following equation for area$A = 1+2.2(x-1.5)^2$<jupyter_code>x=np.linspace(0,3,100)
A=1+2.2*(x-1.5)**2
y=(A/3.1415)**0.5
pyplot.plot(x,y,'black')
pyplot.plot(x,-y,'black')
pyplot.xlim([0,3])
pyplot.ylim([-1.5,1.5])<jupyter_output><empty_output><jupyter_text>Initial ConditionsIt is important to propose good initial conditions for a couple of reasons. One reason is that the closer you are to the final steady state solution, then the less iterations you will need to converge on a solution. Another reason is that if the guess is way off, then the simulation can blow up. The initial conditions for this simulation will be$$\rho = 1-0.3146x$$$$T=1-0.2314x$$$$V=(0.1+1.09x)T^{1/2}$$  Running the Simulation (Main Code)We have finally arrived at the spot where we can put it all together. The following code will be a culmination of all of the information we have gathered before. Let's rock and roll!<jupyter_code>## Main Script for quasi-1D-nozzle flow
def main(totalIter,numPts): # number of iterations, number of grid points
    #Initialize our gird
    length = 3  # nondimensional length of out nozzle
    gridLoc = np.linspace(0,length,numPts)
    area = 1+2.2*(gridLoc-1.5)**2 #initialize area

    #Initial conditions of flow
    gamma = 1.4
    rho = 1-0.3146*gridLoc
    temp = 1-0.2314*gridLoc
    vel = (0.1+1.09*gridLoc)*temp**0.5
    a = temp**0.5 #speed of sound

    #Iterations and plotting set up
    rhoRes = np.ones(numPts)
    velRes = np.ones(numPts)
    k = int(numPts/2)
    
    #Calculate first time step 
    cNum = 0.3 # Courant number
    deltax = length/(numPts-1)
    dtLocal = cNum *(deltax/(a+vel))
    deltat = min(dtLocal)

    #Initialize Matrix to store our output data
    # set up to be n point number of rows X total iterations
    rhoResults = np.zeros(shape=(numPts,totalIter+1))
    tempResults = np.zeros(shape=(numPts,totalIter+1))
    velResults = np.zeros(shape=(numPts,totalIter+1))
    rhoResidualResults = np.zeros(shape=(numPts,totalIter+1))
    velResidualResults = np.zeros(shape=(numPts,totalIter+1))
    
    rhoResults[:,0] = rho
    tempResults[:,0] = temp
    velResults[:,0] = vel
    rhoResidualResults[:,0] = rhoRes
    velResidualResults[:,0] = velRes
    
    iter=1
    #Iterate
    while iter <= totalIter:
        rho,vel,temp,rhoRes,velRes = calcStepSuper(rho,vel,temp,area,gamma,deltat,deltax)
        # Save info in matrix
        rhoResults[:,iter] = rho
        tempResults[:,iter] = temp
        velResults[:,iter] = vel
        rhoResidualResults[:,iter] = rhoRes
        velResidualResults[:,iter] = velRes 
        
        #Recalculate timestep
        a = temp**0.5 #speed of sound
        dtLocal = cNum *(deltax/(a+vel))
        deltat = min(dtLocal)
        # Add 1 to itteration
        iter=iter+1    
        
    return rhoResults,tempResults,velResults,rhoResidualResults,velResidualResults,area
rhoResults,tempResults,velResults,rhoResidualResults,velResidualResults,area = main(1500,31)
print('Done.')<jupyter_output>Done.<jupyter_text>Post Processing Initial RunNow that we have the results of our first simulation, let's look at the results and the iteration steps.<jupyter_code>#Lets plots out results 
x = np.linspace(0,3,len(rhoResults[:,-1]))
fig, ax1 = pyplot.subplots()
ax1.plot(x,rhoResults[:,-1],label='Nondim density')
ax1.plot(x,tempResults[:,-1],label='Nondim Temp')
ax1.plot(x,tempResults[:,-1]*rhoResults[:,-1],label='Nondim Pressure')

ax1.set_xlabel('Nondimensional length (x)')
# Make the y-axis label, ticks and tick labels match the line color.
ax1.set_ylabel('den,temp,pres')
ax1.tick_params('y')
ax1.legend(loc=6)
ax2 = ax1.twinx()
ax2.plot(x,velResults[:,-1]/tempResults[:,-1]**0.5,'y',label='Mach')
ax2.set_ylabel('Mach Number (M)')
ax2.tick_params('y')
ax2.legend(loc = 8)

fig.tight_layout()
pyplot.show()

# Density vs iteration at throat of nozzle
pyplot.plot(rhoResults[int(len(rhoResidualResults)/2),:])
pyplot.xlabel('Iteration')
pyplot.ylabel('Nondimensional density')
pyplot.title('Nondimensional Density vs Iteration')
# Temp vs iteration at throat of nozzle
pyplot.plot(tempResults[int(len(rhoResidualResults)/2),:])
pyplot.xlabel('Iteration')
pyplot.ylabel('Nondimensional temperature')
pyplot.title('Nondimensional Temperature vs Iteration')
# Pressure vs iteration at throat of nozzle
pyplot.plot(rhoResults[int(len(rhoResidualResults)/2)]*tempResults[int(len(rhoResidualResults)/2),:])
pyplot.xlabel('Iteration')
pyplot.ylabel('Nondimensional pressure')
pyplot.title('Nondimensional Pressure vs Iteration')
# Mach vs iteration at throat of nozzle
pyplot.plot(velResults[int(len(rhoResidualResults)/2),:]/tempResults[int(len(rhoResidualResults)/2),:]**0.5)
pyplot.xlabel('Iteration')
pyplot.ylabel('Mach Number')
pyplot.title('Mach Number vs Iteration')<jupyter_output><empty_output><jupyter_text>As we can see above, the values begin to converge around 600 iterations. By looking at the mach number plot, we can see that the selected location is $A^*$ because the mach number converges at $M=1$. As shown in all of the parameters, the values begin with a great rate of change in value. As the iterations continue, the solution settles around a single value. These big oscillations in the beginning are due to unsteady compression and expansion waves that move within the nozzle. We can also look at the residuals for our simulation as a function of iterations. The residuals represent the change in solution from iteration to iteration. This value can be obtained by plotting the value of $(\partial \rho /\partial t)_{av}$ and $(\partial V /\partial t)_{av}$ against the iterations.<jupyter_code>pyplot.plot(abs(rhoResidualResults[int(len(rhoResidualResults)/2),:]),label="Density Residual")
pyplot.plot(abs(velResidualResults[int(len(rhoResidualResults)/2),:]),label="Velocity Residual")
pyplot.yscale('log')
pyplot.xlabel('Iteration')
pyplot.ylabel('Residual')
pyplot.legend()
pyplot.title('Residuals')<jupyter_output><empty_output><jupyter_text>At 1500 iterations, the residual is still going down at a steady rate. If we were to increase the iterations to 5000 we would see the residual plateau. Another metric to look at is the mass flow parameter. We can calculate mass flux with the following equation.$$\dot{m} = \rho V A$$If our steady state solution represents physics well, then that value should be constant at any point in the nozzle. $$\rho V A = const$$<jupyter_code># Rerun the simulation with more iterations
rhoResults,tempResults,velResults,rhoResidualResults,velResidualResults,area = main(2000,61)
mdot = np.transpose(np.transpose(rhoResults)*np.transpose(velResults)*area)
pyplot.figure(figsize=(5, 8))
pyplot.plot(mdot[:,0],label='iter=0')
pyplot.plot(mdot[:,50],label='iter=50')
pyplot.plot(mdot[:,100],label='iter=100')
pyplot.plot(mdot[:,500],label='iter=500')
pyplot.plot(mdot[:,1000],label='iter=1000')
pyplot.plot(mdot[:,2000],label='iter=2000')
pts = 10
pyplot.scatter(np.linspace(0,60,pts),np.ones(pts)*0.579,label = 'analytic solution')
pyplot.xlim([0,60])
pyplot.legend()<jupyter_output><empty_output><jupyter_text>This plot shows us that as the iterations continue, the numerical solution gets closer to a flat line and approaches the analytical solution which turns out to be 0.579. To further compare our results to the analytical solution, we can plot our steady state results against analytical data points.<jupyter_code>x = np.linspace(0,3,len(rhoResults[:,-1]))
pyplot.plot(x,rhoResults[:,-1],label='Nondim density')
pyplot.scatter([0.0,.5,1,1.5,2,2.5,3],[0.995,0.983,0.920,0.634,0.258,0.105,0.052],color='r')
pyplot.xlim([0,3])
pyplot.xlabel('x/L')
pyplot.ylabel('Nondimensional density')
pyplot.title('Numerical vs Analytical Density')
pyplot.plot(x,velResults[:,-1]/tempResults[:,-1]**0.5,label='Mach')
pyplot.scatter([0.0,.5,1,1.5,2,2.5,3],[0.098,0.185,0.413,1,1.896,2.706,3.359],color='r')
pyplot.xlim([0,3])
pyplot.xlabel('x/L')
pyplot.ylabel('Mach Number')
pyplot.title('Numerical vs Analytical Mach Number')<jupyter_output><empty_output><jupyter_text>Solution of Purely Subsonic Isentropic Nozzle FlowThe difference between this version of the problem and the one we just solved is that for purely subsonic flow, there are an infinite number of solutions depending on the pressure ration between the inlet and outlet. For a subsonic flow, te fluid will speed up at it approaches the throat but because the flow has not reached the speed of sound, then the flow will slow down as the nozzle expands. The following equations govern the purely subsonic flow.$$\frac{p_e}{p_0} = \Big(1+ \frac{\gamma-1}{2}M_e^2\Big)^{\gamma/(\gamma-1)}$$Here, $M_e$ and $p_e$ represent the mach number and pressure at the exit respectively. $p_0$ is the pressure at the inlet. $$\frac{A_e}{A^*} = \frac{1}{M_e^2}\Big[\frac{2}{\gamma-1}\Big(1+\frac{\gamma-1}{2}M_e^2\Big)\Big]^{(\gamma+1)/(\gamma-1)}$$$A^*$ is a reference value for the above equation. There is no physical $A^*$, because that is where the flow would theoretically go to $M=1$. This theoretical values will provide the value for the local mach number using the equation below.$$\Big (\frac{A}{A^*} \Big)^2 = \frac{1}{M^2}\Big[\frac{2}{\gamma +1}\Big(1+\frac{\gamma-1}{2}M^2\Big)\Big]^{(\gamma+1)/(\gamma-1)}$$The rest of the local parameters can be solved just as before:$$\frac{p}{p_0}  = \Big(1+\frac{\gamma-1}{2}M^2\Big)^{(-\gamma)/(\gamma-1)}$$$$\frac{\rho}{\rho_0}  = \Big(1+\frac{\gamma-1}{2}M^2\Big)^{-1/(\gamma-1)}$$$$\frac{T}{T_0}  = \Big(1+\frac{\gamma-1}{2}M^2\Big)^{-1}$$For this simulation, we will model the nozzle with the following equations.$$\frac{A}{A_t} = 1+2.2\Big(\frac{x}{L} - 1.5\Big)^2$ for $0\leq\frac{x}{L}\leq1.5$$$$\frac{A}{A_t} = 1+2.2223\Big(\frac{x}{L} - 1.5\Big)^2$ for $1.5\leq\frac{x}{L}\leq3.0$$<jupyter_code>x = np.linspace(0,3,100)
area = np.zeros(len(x))
index = 0
for i in x:
    if i<1.5:
        area[index] = 1+2.2*(i-1.5)**2
    else:
        area[index] = 1+0.2223*(i-1.5)**2
    index =index+1
    
y=(area/3.1415)**0.5
pyplot.figure(figsize=(8, 4))
pyplot.plot(x,y,'black')
pyplot.plot(x,-y,'black')
pyplot.xlim([0,3])
pyplot.ylim([-1.5,1.5])<jupyter_output><empty_output><jupyter_text>The governing equations for the simulation are the same but the main difference will be the boundary conditions at the inlet and outlet of the nozzle. The inlet BC is handled the same way as before, but this time we must address the outflow boundary conditions. Because the flow is subsonic, information can travel upstream and therefor we need to let two variables float. For this reason, we will only set one variable, and that will be the exit pressure. As you may have noticed, we have been dealing with density, velocity and temperature, so to set the value of p we use the following equation. $$p = \rho RT$$In non-dimensional terms, it reduces to $$p'_e = \rho'_eT'_e = specified$$To apply this boundary condition to our flow, we first extrapolate temperature or density, and then use that extrapolated value with the relationship for pressure and get the other value. you can start with temperature and get density or visa versa, but the end result should ultimately be the same. Here is what it would look like if we first extrapolated temperature. First extrapolate to get the value of $T'$ at the outflow boundary.$$T'_N = 2T'_{N-1}-T'_{N-2}$$With the new value for T, calculate density with the set pressure boundary condition.$$\rho'_N = \frac{p'_N}{T'_N} = \frac{specified value}{T'_N}$$Lastly, we can extraploate the value for velocity as such:$$V'_N = V'_{N-1} - V'_{N-2}$$Now let us modify the funciton we made before with the new boundary condition. The input to the function will be the same except for a set value for the exit pressure.<jupyter_code>def calcStepSub(exitPres,rho,vel,temp,area,gamma,deltat,deltax):
    """
    Import three numpy arrays that represent the 
    density velocity and temperature of each element
    with a supersonic exit flow BC
    rho : density array 
    vel : velocity array
    temp: temperature array
    deltat : time step for the elements
    deltax : space between each element
    """ 
    #loop through each element starting at index 1 and ending at index numPts-1 (N-1)
    numElm = len(rho)
    
    #Density
    fwddRhodt  = np.zeros(numElm)
    bkwddRhodt  = np.zeros(numElm)
    rhoBar = np.zeros(numElm)
    dRhoAv = np.zeros(numElm)
    #Velocity
    
    fwddvdx  = np.zeros(numElm)
    bkwddvdx  = np.zeros(numElm)
    VBar = np.zeros(numElm)
    dVAv = np.zeros(numElm)
    #Temperature 
    
    fwddTdx = np.zeros(numElm)
    bkwddTdx = np.zeros(numElm)
    TBar = np.zeros(numElm)
    dTAv = np.zeros(numElm)
    
    ##################
    ####Predictor Step
    ##################
    for i in range(1,numElm-1):
            #Predictor Step
            # Calculate first derivative of density
            fwddRhodt[i] = -rho[i] * ( (vel[i+1]-vel[i])/deltax) -(rho[i]*vel[i]*(math.log(area[i+1])-math.log(area[i]))/deltax) - (vel[i]*(rho[i+1]-rho[i])/deltax)
            rhoBar[i]  = rho[i]+  fwddRhodt[i] * deltat
            
            # Calculate first derivative of velocity
            fwddvdx[i] = -vel[i]*((vel[i+1]-vel[i])/deltax)-(1/gamma)*(((temp[i+1]-temp[i])/deltax)+(temp[i]/rho[i])*(rho[i+1]-rho[i])/deltax)
            VBar[i]  = vel[i] + fwddvdx[i] * deltat
                 
            # Calculate first derivative of temperatures
            a = -vel[i]*(temp[i+1]-temp[i])/deltax
            b = (gamma-1)*temp[i]
            c = (vel[i+1]-vel[i])/deltax
            d = vel[i]*((math.log(area[i+1])-math.log(area[i]))/deltax)
            fwddTdx[i] = a-b*(c+d)
            TBar[i] = temp[i] + (fwddTdx[i]) * deltat
    
                           
    rhoBar[0] = 1
    TBar[0] = 1
    VBar[0] = 2*vel[1]-vel[2]
    
    # We have the values for the predicted values and the slopes
    ###################
    ####Corrector Step
    ###################
    for i in range(1,numElm-1):
        
        bkwddRhodt[i] = -rhoBar[i] * ( (VBar[i]-VBar[i-1])/deltax) - (rhoBar[i]*VBar[i]*(math.log(area[i])-math.log(area[i-1]))/deltax) -(VBar[i]*(rhoBar[i]-rhoBar[i-1])/deltax)
        
        bkwddvdx[i] = -VBar[i]*((VBar[i]-VBar[i-1])/deltax)-(1/gamma)*(((TBar[i]-TBar[i-1])/deltax)+(TBar[i]/rhoBar[i])*(rhoBar[i]-rhoBar[i-1])/deltax)
            
        a = -VBar[i]*(TBar[i]-TBar[i-1])/deltax
        b = (gamma-1)*TBar[i]
        c = (VBar[i]-VBar[i-1])/deltax
        d = VBar[i]*((math.log(area[i])-math.log(area[i-1]))/deltax)
        bkwddTdx[i] = a-b*(c+d)
    
  
        dRhoAv[i] = 0.500*(fwddRhodt[i]+bkwddRhodt[i])
        dVAv[i]   = 0.500*(fwddvdx[i]+bkwddvdx[i])
        dTAv[i]   = 0.500*(fwddTdx[i]+bkwddTdx[i])
        
        rho[i] = rho[i]   + dRhoAv[i] * deltat
        vel[i] = vel[i]   + dVAv[i]   * deltat
        temp[i] = temp[i] + dTAv[i]   * deltat
        
    #Set BC's for the final values 
    #Far left BC
    rho[0] = 1
    vel[0] = 2*vel[1]-vel[2]
    temp[0] = 1
    
    #Far right BC (DIFFERENT FOR SUBSONIC OUTFLOW)                               
    rho[numElm-1] = 2*rho[numElm-2]-rho[numElm-3]
    temp[numElm-1] = exitPres/rho[numElm-1]
    vel[numElm-1] = 2*vel[numElm-2]-vel[numElm-3]

    return rho,vel,temp,dRhoAv,dVAv #output new values and d/dt for residuals<jupyter_output><empty_output><jupyter_text>Main Code for Purely Subsonic Flow<jupyter_code>## Main Script for quasi-1D-nozzle flow
def mainSub(totalIter,numPts,exitPres): # number of iterations, number of grid points
    #Initialize our gird
    length = 3  # nondimensional length of out nozzle
    gridLoc = np.linspace(0,length,numPts)
    
    #initialize area
    area=np.zeros(numPts)
    index = 0
    for i in gridLoc:
        if i<1.5:
            area[index] = 1+2.2*(i-1.5)**2
        else:
            area[index] = 1+0.2223*(i-1.5)**2
        index =index+1
    
    #Initial conditions of flow
    gamma = 1.4
    rho = 1-0.023*gridLoc
    temp = 1-0.009333*gridLoc
    vel = 0.05+0.11*gridLoc
    a = temp**0.5 #speed of sound

    #Residual Array allowcation
    rhoRes = np.ones(numPts)
    velRes = np.ones(numPts)
    
    #Calculate first time step 
    cNum = 0.55 # Courant number
    deltax = length/(numPts-1)
    dtLocal = cNum *(deltax/(a+vel))
    deltat = min(dtLocal)

    #Initialize Matrix to store our output data
    # set up to be n point number of rows X total iterations
    rhoResults = np.zeros(shape=(numPts,totalIter+1))
    tempResults = np.zeros(shape=(numPts,totalIter+1))
    velResults = np.zeros(shape=(numPts,totalIter+1))
    rhoResidualResults = np.zeros(shape=(numPts,totalIter+1))
    velResidualResults = np.zeros(shape=(numPts,totalIter+1))
    
    rhoResults[:,0] = rho
    tempResults[:,0] = temp
    velResults[:,0] = vel
    rhoResidualResults[:,0] = rhoRes
    velResidualResults[:,0] = velRes
    
    iter=1
    #Iterate
    while iter <= totalIter:
        rho,vel,temp,rhoRes,velRes = calcStepSub(exitPres,rho,vel,temp,area,gamma,deltat,deltax)
        # Save info in matrix
        rhoResults[:,iter] = rho
        tempResults[:,iter] = temp
        velResults[:,iter] = vel
        rhoResidualResults[:,iter] = rhoRes
        velResidualResults[:,iter] = velRes 
        
        #Recalculate timestep
        a = temp**0.5 #speed of sound
        dtLocal = cNum *(deltax/(a+vel))
        deltat = min(dtLocal)
        # Add 1 to itteration
        iter=iter+1    
        
    return rhoResults,tempResults,velResults,rhoResidualResults,velResidualResults,area
# Run simulation
rhoResults,tempResults,velResults,rhoResidualResults,velResidualResults,area = mainSub(5000,61,0.93)

print('Done.')<jupyter_output>Done.<jupyter_text>Now with the new simulation complete, let's look at our residuals<jupyter_code>pyplot.plot(abs(rhoResidualResults[int(len(rhoResidualResults)/2),:]),label="Density Residual")
pyplot.plot(abs(velResidualResults[int(len(rhoResidualResults)/2),:]),label="Velocity Residual")
pyplot.yscale('log')
pyplot.xlabel('Iteration')
pyplot.ylabel('Residual')
pyplot.legend()
pyplot.title('Residuals')<jupyter_output><empty_output><jupyter_text>For this simulation, the residuals showed a constant downward trend. This tells us that our simulation did not blow up, but is only converged to around E-3 accuracy. To understand more about how the solution converged lets take a look at mass flow rate as the iterations progressed.<jupyter_code>mdot = np.transpose(np.transpose(rhoResults)*np.transpose(velResults)*area)
pyplot.figure(figsize=(6, 6))
x = np.linspace(0,3,61)
pyplot.plot(x,mdot[:,0],label='iter=0')
pyplot.plot(x,mdot[:,500],label='iter=500')
pyplot.plot(x,mdot[:,5000],label='iter=5000')
pyplot.legend(loc=1)
pyplot.xlabel('x/L')
pyplot.ylabel('mass flow rate')<jupyter_output><empty_output><jupyter_text>As one can see from this plot, the values of mass flow rate converge on a constant value as expected. This value seems to match up well with the CFD results from Anderson's book. Additionally, we can examine how the pressure changed over the time steps.<jupyter_code>pressure = rhoResults*tempResults
x = np.linspace(0,3,61)
pyplot.plot(x,rhoResults[:,500],label='500 iter')
pyplot.plot(x,rhoResults[:,1000],label='1000 iter')
pyplot.plot(x,rhoResults[:,5000],label='5000 iter')
pyplot.xlim([0,3])
pyplot.legend()
pyplot.xlabel('x/L')
pyplot.ylabel('Pressure (nondimensionalized)')
#Lets plots out results 
x = np.linspace(0,3,len(rhoResults[:,-1]))

fig, ax1 = pyplot.subplots(figsize=(10, 6))
ax1.plot(x,rhoResults[:,-1],label='Nondim density')
ax1.plot(x,tempResults[:,-1],label='Nondim Temp')
ax1.plot(x,tempResults[:,-1]*rhoResults[:,-1],label='Nondim Pressure')

ax1.set_xlabel('Nondimensional length (x)')
# Make the y-axis label, ticks and tick labels match the line color.
ax1.set_ylabel('den,temp,pres')
ax1.tick_params('y')
ax1.legend(loc=6)
ax2 = ax1.twinx()
ax2.plot(x,velResults[:,-1]/tempResults[:,-1]**0.5,'y',label='Mach')
ax2.set_ylabel('Mach Number (M)')
ax2.tick_params('y')
ax2.legend(loc = 1)
fig.tight_layout()

pyplot.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15357053993530728, 'fraction_numerical': 0.022070166708136354, 'mean_word_length': 3.963077303037787, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 60, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '497978', 'license': 'MIT', 'n_tokens_mistral': 14346, 'n_tokens_neox': 13621, 'n_words': 4076, 'path': 'nozzleFlow/nozzleNonconservationForm.ipynb', 'repo_name': 'bjreddish/CFD', 'size': '40200'}